POJ 2429 GCD & LCM Inverse

题目链接：GCD & LCM Inverse

解题思路：这个给了a和b的GCD和LCM，让求出a和b，这里a和b有很多解，让输出a + b最小的一组。

gcd(a, b) = k

lcm(a, b) = m

a = pk, b = qk

a + b = (p + q)k

p * q * k = m

可以通过以上得出

p * q = m / k，gcd(p, q) = 1， 所以我们要求出q + p最小的。

我们可以知道当d是m / k的约数的时候并且gcd(d, m / k / d) = 1，那么我们这里可以看到，当p q的差距越来越小的时候p + q就越来越小，但是我们对于longlong没有sqrt，所以我们要另一种思路。

我们将m / k分解质因数，我们将相同的质因数相乘合成n个约数，这里为什么要合成，因为同一个约数既在p也在q的时候gcd(p, q) ！= 1，所以我们这样合成的都是两两互质的，所以搜索一遍就得出结果了。

这里涉及大整数分解，那就是Pollard rho算法，详细请见

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#include<ctime>
#include<algorithm>
#include<cstdlib>
#include<vector>
#define ll __int64
#define MAX 5500

const int S = 20;//素数测试误判率=1 / 2^S   S越大误差率越小

using namespace std;


ll gcd(ll a, ll b){
    return b == 0 ? a : gcd(b, a % b);
}


ll multi_mod(ll a, ll b, ll p){
    ll ret = 0, q = a % p;
    while(b){
        if(b & 1){
            ret = (ret + q) % p;
        }
        q = (q + q) % p;
        b >>= 1;
    }
    return ret % p;
}

ll pow_mod(ll a, ll b, ll p){
    ll ret = 1, q = a % p;
    while(b){
        if(b & 1){
            ret = multi_mod(ret, q, p);
        }
        q = multi_mod(q, q, p);
        b >>= 1;
    }
    return ret % p;
}

bool Miller_Rabin(ll p){
    int i, j, k;
    ll u, t, x, y;
    if(p == 2)  return true;
    if(p % 2 == 0 || p == 1) return false;
    u = p - 1, t = 0;
    while(u % 2 == 0){
        t++;
        u >>= 1;
    }
    for(int i = 0; i < S; i++){
        x = rand() % (p - 1) + 1;
        x = pow_mod(x, u, p);
        ll y = 0;
        for(j = 0; j < t; j++){
            y = multi_mod(x, x, p);
            if(y == 1 && x != 1 && x != p - 1){
                return false;
            }
            x = y;
        }
        if(y != 1)  return false;
    }
    return true;
}



ll fact[MAX];
ll tot = 0;

ll pollard_rho(ll n, ll c){
    ll i = 1, k = 2;
    ll x = rand() % (n - 1) + 1;
    ll y = x;
    while(true)
    {
        i++;
        x = (multi_mod(x, x, n) + c) % n;
        ll d = gcd((y - x + n) % n, n);
        if(1 < d && d < n) return d;
        if(y == x) return n;
        if(i == k)
        {
            y = x;
            k <<= 1;
        }
    }
}

void find(ll n, ll c){
    if(n == 1){
        return;
    }
    if(Miller_Rabin(n)){
        fact[tot++] = n;
        return;
    }
    ll d = n;
    ll k = c;
    while(d >= n){
        d = pollard_rho(n, c--);
    }
    find(d, k);
    find(n / d, k);
}


vector<ll> num;


int main(){
    int i, j, k;
    ll a, b, p;
    ll n;
    //srand(time(NULL));
    while(cin >> a >> b){
        ll gcd = a;
        n = b / a;
        tot = 0;
        num.clear();
        find(n, 120);
        sort(fact, fact + tot);
        for(i = 0; i < tot; i++){
            if(i == 0 || fact[i] != fact[i - 1]){
                num.push_back(fact[i]);
            }
            else{
                num[num.size() - 1] *= fact[i];
            }
        }
        p = n + 1;
        a = 1, b = n;
        ll tem, am;
        for(i = 0; i < (1 << num.size()); i++){
            int j = i;
            tem = 1, am = 0;
            while(j){
                if(j & 1){
                    tem *= num[am];
                }
                j >>= 1;
                am++;
            }
            ll tt = n / tem;
            //cout << tt << " " << tem << endl;
            if(tt + tem < p){
                p = tt + tem;
                a = min(tt, tem);
                b = max(tt, tem);
            }
        }
        cout << gcd * a << " " << gcd * b << endl;
    }
    return 0;
}