백준 20226번 Luggage

문제 요약

숫자 $p$가 주어진다. $p=whd$로 나타낼 수 있는 $(w, h, d)$중에서 $w+h+d$의 값이 최소가 되는 $(w,h,d)$의 $w+h+d$를 출력하시오.

주어지는 테스트케이스의 최대 숫자는 300개이며 $p \le 10^{15}$이다.

풀이

가장 간편한 소인수분해인 $O(\sqrt{N})$짜리 소인수분해 방법으론 대략 3000만 곱하기 300으로 시간초과를 받을거 같다. 그래서 소인수분해 하려면 더 빠른 방법을 써야만 한다. 나는 폴라드 로를 사용했다.

이제 $p$를 소인수 분해해서 $p$의 약수를 전부 만들어준다. 약수가 제일 많은 수라고 해봤자 3만개를 넘지 않는다.

이제 이 약수들 중에서 하나를 고정한 다음에 나머지 숫자들을 찾아보도록 하자. 고정하는 수는 $d$로 하자.

$d$를 고정했으니 $wh = \frac{p}{d}$인 $w, h$ 중에서 $w+h+d$가 최소인 친구를 찾아야 한다. 일단 $w, h$가 $p$의 약수이면서 $p/d$의 약수여야 되는 조건을 떼고 생각을 해보자. 그리고 $k = p/d$로 표현하자.

$w, h$는 양수라는 조건을 넣으면 $\frac{w+h}{2} \ge \sqrt{wh}$가 성립하고, 이 부등식은 $w=h$일 때 등호가 성립하는 것이 알려져 있다. 그리고 그 값은 $w=h=\sqrt{k}$이다.

이제 다시 원래 문제로 돌아오자. $wh=k$이면서 $w+h$가 최소가 되는 자연수 $w, h$를 찾고 싶다. $\sqrt{k}$에서 최솟값을 가지기 때문에 $\sqrt{k}$와 가장 가까운 $w, h$에서 그 최솟값을 가지게 될 것이다.

따라서, $d$를 고정한 뒤에 $\sqrt{\frac{p}{d}}$에 가장 가까운 $w, h$를 찾고 확인하면 $d$에 대해서 나올 수 있는 최솟값을 확인한 것이다. 따라서 모든 $d$에 대해서 이 과정을 수행하면 된다.

시간복잡도는 폴라드 로가 $O(p^{1/4}\log p)$정도 일 것이고, 약수의 수를 $p$에 대한 식으로 어떻게 나타날지 모르겠지만 이걸 $n$이라고 하자. 하나의 $d$에 대해서 $O(\log n)$ 정도가 걸리니까 $O(n \log n)$이다. 따라서, 하나의 케이스당 $O(p^{1/4} \log p + n \log n)$정도로 통과한다.

#include<bits/stdc++.h>
using namespace std;
#define all(x) begin(x), end(x)
typedef long long ll;
typedef pair<int, int> pii;
typedef vector<int> vi;
typedef unsigned long long ull;

// from https://lemire.me/blog/2013/12/26/fastest-way-to-compute-the-greatest-common-divisor/
ull bgcd(ull u, ull v) {
  int shift;
    if (u == 0) return v;
    if (v == 0) return u;
    shift = __builtin_ctz(u | v);
    u >>= __builtin_ctz(u);
    do {
        v >>= __builtin_ctz(v);
        if (u > v) {
            ull t = v;
            v = u;
            u = t;
        }  
        v = v - u;
    } while (v != 0);
    return u << shift;
}
// from https://github.com/kth-competitive-programming/kactl/blob/master/content/number-theory/ModMulLL.h
ull modmul(ull a, ull b, ull M) {
    ll ret = a * b - M * ull(1.L / M * a * b);
    return ret + M * (ret < 0) - M * (ret >= (ll)M);
}
ull modpow(ull b, ull e, ull mod) {
    ull ans = 1;
    for (; e; b = modmul(b, b, mod), e /= 2)
        if (e & 1) ans = modmul(ans, b, mod);
    return ans;
}
// from https://github.com/kth-competitive-programming/kactl/blob/master/content/number-theory/MillerRabin.h
bool isPrime(ull n) {
    if (n < 2 || n % 6 % 4 != 1) return (n | 1) == 3;
    static ull A[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
    ull s = __builtin_ctzll(n-1), d = n >> s;
    for (ull a : A) {   // ^ count trailing zeroes
        ull p = modpow(a%n, d, n), i = s;
        while (p != 1 && p != n - 1 && a % n && i--)
            p = modmul(p, p, n);
        if (p != n-1 && i != s) return 0;
    }
    return 1;
}
// from https://github.com/kth-competitive-programming/kactl/blob/master/content/number-theory/Factor.h
ull pollard(ull n) {
    auto f = [&n](ull x) { return modmul(x, x, n) + 1; };
    ull x = 0, y = 0, t = 30, prd = 2, i = 1, q;
    while (t++ % 40 || bgcd(prd, n) == 1) {
        if (x == y) x = ++i, y = f(x);
        if ((q = modmul(prd, max(x,y) - min(x,y), n))) prd = q;
        x = f(x), y = f(f(y));
    }
    return bgcd(prd, n);
}
vector<ull> factor(ull n) {
    if (n == 1) return {};
    if (isPrime(n)) return {n};
    ull x = pollard(n);
    auto l = factor(x), r = factor(n / x);
    l.insert(l.end(), all(r));
    return l;
}

ull n;
vector<pair<ull, int>> factors;
vector<ull> divisors;

void dfs(int idx, ull x) {
    if(idx == factors.size()) {
        divisors.push_back(x);
        return ;
    }
    ull t = 1;
    for(int i=0;i<=factors[idx].second;++i) {
        dfs(idx+1, x*t);
        t *= factors[idx].first;
    }
}

int main() {
    while(true) {
        cin >> n;
        if(!n) break;
        auto v = factor(n);
        map<ull, int> mp;
        for(ull i : v) {
            if(mp.find(i) == mp.end()) mp[i] = 0;
            mp[i]++;
        }
        factors.clear();
        for(auto it : mp) factors.emplace_back(it.first, it.second);
        divisors.clear();
        dfs(0, 1);
        sort(divisors.begin(), divisors.end());
        assert(divisors.size() <= 30000);
        ull ans = 1e18;
        for(int i=0;i<divisors.size();++i) {
            ull d = n/divisors[i];
            ull x = divisors[i];
            ull opt = ceil(2.0*sqrt((long double)x));
            if(opt + d > ans) break;
            int idx = lower_bound(divisors.begin(), divisors.end(), (ull)sqrt((long double)x)) - divisors.begin();
            int diff = 0;
            while(idx - diff >= 0 || idx + diff < divisors.size()) {
                bool flag = false;
                if(idx-diff>=0 && x % divisors[idx-diff] == 0) {
                    ull w = divisors[idx-diff];
                    ull h = x/w;
                    ans = min(ans, w+h+d);
                    flag = true;
                }
                if(idx+diff < divisors.size() && x % divisors[idx+diff] == 0) {
                    ull w = divisors[idx+diff];
                    ull h = x/w;
                    ans = min(ans, w+h+d);
                    flag = true;
                }
                if(flag) break;
                ++diff;
            }
        }
        cout << ans << '\n';
    }
    return 0;
}