题目链接

GCD & LCM Inverse POJ - 2429

思路

$$ \begin{aligned} 已知gcd和lcm，求a，b 满足 a + b最小 \\ 求一个w = lcm / gcd \\ 然后把w分解因子 \\ 这个题要用Pollard’s Rho 板子调了一天。。。 \end{aligned} $$

时间复杂度

$$ O(N^{1/4})N为要分解的数 $$

代码

#include <cstdio>
#include <vector>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include <map>

using namespace std;

class Factor {
public:
    vector<bool> is_prime;
    //pi : min prime factor
    vector<int> prime, pi;

    Factor () {
        sieve((int)1e6 + 10);
    }

    int sieve(int sz) {
        int cnt = 0;

        is_prime.resize(sz + 1, false);
        pi.resize(sz + 1, 0);
        for (int i = 1; i <= sz; i++) {
            pi[i] = i;
        }
        for (int i = 2; i <= sz; i++) {
            if (pi[i] == i) {
                prime.push_back(i);
                is_prime[i] = true;
                ++cnt;
            }
            for (int j = 0; j < cnt && 1LL * prime[j] * i <= sz; j++) {
                pi[prime[j] * i] = prime[j];
                if (i % prime[j] == 0) {
                    break;
                }
            }
        }
        return cnt;
    }

    inline long long mul_mod(long long a, long long b, long long mod) const {
        if (mod <= 1000000000) {
            return a * b % mod;
        } else if (mod <= 1000000000000LL) {
            return (((a * (b >> 20) % mod) << 20) + (a * (b & ((1 << 20) - 1)))) % mod;
        } else {
            long long d = (long long) floor(a * (long double)b / mod + 0.5);
            long long res = (a * b - d * mod) % mod;
            if (res < 0) {
                res += mod;
            }
            return res;
        }
    }

    long long pow_mod(long long a, long long b, long long mod) {
        long long res = 1;
        while (b) {
            if (b & 1) {
                res = mul_mod(res, a, mod);
            }
            a = mul_mod(a, a, mod);
            b >>= 1;
        }
        return res;
    }

    bool witness(long long a, long long n) {
        int t = 0;
        long long u = n - 1;
        while ((~u) & 1) {
            t++;
            u >>= 1;
        }
        long long x = pow_mod(a, u, n);
        long long xx = 0;
        while (t--) {
            xx = mul_mod(x, x, n);
            if (xx == 1 && x != 1 && x != n - 1) {
                return true;
            }
            x = xx;
        }
        return xx != 1;
    }

    inline long long randll() {
        return (long long)rand() << 30 | rand();
    }



    bool miller(long long n) {
        if (n < 2) {
            return false;
        }
        if (n == 2) {
            return true;
        }
        if (n < (int)is_prime.size()) {
            return is_prime[n];
        }
        if ((~n) & 1) {
            return false;
        }
        for (int i = 0; i <= 7; i++) {
            if (witness(randll() % (n - 1) + 1, n)) {
                return false;
            }
        }
        return true;
    }

    // fast gcd
    long long gcd(long long a, long long b) {
        long long res = 1;
        while (a != 0) {
            if (((~a) & 1) && ((~b) & 1)) {
                res <<= 1;
                a >>= 1;
                b >>= 1;
            } else if ((~a) & 1) {
                a >>= 1;
            } else if ((~b) & 1) {
                b >>= 1;
            } else {
                if (a < b) {
                    swap(a, b);
                }
                a -= b;
            }
        }
        return res * b;
    }

    long long rho(long long n) {
        static long long a[1001000];
        while (true) {

            long long X = rand() % n, Z, T = 1;

            long long *lY = a, *lX = lY;
            int tmp = 20;

            int C = rand() % 10 + 3;
            X = mul_mod(X, X, n) + C;
            *(lY++) = X;
            lX++;

            long long Y = mul_mod(X, X, n) + C;
            *(lY++) = Y;

            while (X != Y) {
                long long t = X - Y + n;
                Z = mul_mod(T, t, n);
                if (Z == 0) {
                    return gcd(T, n);
                }
                tmp--;
                if (tmp == 0) {
                    tmp = 20;
                    Z = gcd(Z, n);
                    if (Z != 1 && Z != n) {
                        return Z;
                    }
                }
                T = Z;
                Y = *(lY++) = mul_mod(Y, Y, n) + C;
                Y = *(lY++) = mul_mod(Y, Y, n) + C;
                X = *(lX++);
            }
        }
    }

    vector<long long> fac;
    void print_fac() {
        for (int i = 0; i < (int)fac.size(); i++) {
            printf("%lld ", fac[i]);
        }
        puts("");
    }

    void _factor(long long n) {
        int sz = (int)fac.size();
        for (int i = 0; i < sz; i++) {
            if (n % fac[i] == 0) {
                n /= fac[i];
                fac.push_back(fac[i]);
            }
        }
        if (n <= (int)pi.size()) {
            while (n != 1) {
                fac.push_back(pi[n]);
                n /= pi[n];
            }
            return;
        }
        if (miller(n)) {
            fac.push_back(n);
            return;
        } else {
            long long x = rho(n);
            _factor(x);
            _factor(n / x);
        }

    }

    // gen all prime factors
    vector<long long> fact_big(long long n) {
        fac.clear();
        _factor(n);
        return fac;
    }
};
Factor fa;
map<long long, int> cnt;
map<int, long long> id;
int main() {
    long long a, b;

    while (~scanf("%lld %lld", &a, &b)) {

        if (a == b) {
            printf("%lld %lld\n", a, a);
            continue;
        }
        b /= a;
        vector<long long> w = fa.fact_big(b);
        cnt.clear();
        id.clear();
        int index = 0;
        for (int i = 0; i < (int)w.size(); i++) {
            if (cnt.count(w[i]) == 0) {
                id[index] = w[i];
                index++;
            }
            cnt[w[i]]++;
        }
        long long x = 1, y = b;
        for (int i = 1; i <= (1 << index) - 2; i++) {
            long long cur = 1;
            for (int j = 0; j < index; j++) {
                if ((i >> j) & 1) {
                    int k = cnt[id[j]];
                    while (k--) {
                        cur *= id[j];
                    }
                }
            }
            if (cur + b / cur < x + y) {
                x = cur;
                y = b / cur;
            }
        }
        if (x > y) {
            swap(x, y);
        }
        printf("%lld %lld\n", x * a, y * a);
    }

    return 0;
}