POJ 2429 GCD & LCM Inverse ★(pollard-ρ && DFS枚举)

题目链接：http://poj.org/problem?id=2429 题目大意：给定gcd(a,b)和lcm(a,b)（<2^63），求a和b，如果有多种情况，输出和最小的情况.   首先gcd(a,b) * lcm(a,b) = a*b，但是如果我们直接从a*b中分解因子的话，a*b是可能超过long long的，这样就不好处理了. 我们可以先把gcd(a,b)都分给a,b，因为他们的因子中都要有gcd(a,b).于是现在还剩下lcm(a,b)/gcd(a,b)了，于是我们先用pollard-rho给他分解因子. 那么还有一个问题，能随便分么？分出来的a,b虽然能保证a*b，但是能保证他们的gcd和lcm都是给定的么？不一定. 所以我们还需要注意分因数的过程中还要保证gcd和lcm的正确性. 但这个问题其实不难，因为a*b一定是不变的，所以我们只要保证gcd和lcm其中一个不变，另一个也就自然不变了.显然保证gcd比较简单. 我们知道lcm/gcd = p1^q1 * p2 ^q2 *……* pn^qn，其中p1,p2,……,pn是因子.我们只要保证某个数把某个pi全部取走，这样他们除了先前取走的gcd外再无公因数，则可保gcd正确.这样的话，2^63内的数所有不同的因子个数最多也就十几个，枚举无压力.  

#include 
#include 
#include 
using namespace std;

//return a * b % m
unsigned long long mul_mod(unsigned long long a, unsigned long long b, unsigned long long m){
    //为了防止long long型a * b溢出，有时需要把乘法变加法
    //且因为暴力加法会超时要使用二分快速乘法模（模仿二分快速幂模……）
    unsigned long long res = 0, tmp = a % m;
    while(b){
        if (b & 1)
        {
            res = res + tmp;
            res = (res >= m ? res - m : res);
        }
        b >>= 1;
        tmp <<= 1;
        tmp = (tmp >= m ? tmp - m : tmp);
    }
    return res;
}

//return a ^ b % m
long long exp_mod(long long a, long long b, long long m){
    long long res = 1 % m, tmp = a % m;
    while(b){
        if (b & 1){
            //如果m在int范围内直接用下一式乘就可以，否则需要用下二式把乘法化加法，用快速乘法模
            //res = (res * t) % m;
            res = mul_mod(res, tmp, m);
        }
        //同上
        //t = t * t % m;
        tmp = mul_mod(tmp, tmp, m);

        b >>= 1;
    }
    return res;
}

/*-------------Miller-Rabin 素数测试 部分(用到上面mul_mod和exp_mod   素数return true)--------------*/
bool Miller_Rabin(long long n){
    int a[5] = {2, 3, 7, 61, 24251};
    //一般Miller_Rabin素数测试是随机选择100个a,这样的错误率为0.25^100
    //但在OI&&ACM中，可以使用上面一组a，在这组底数下，10^16内唯一的强伪素数为46,856,248,255,981

    if (n == 2)
        return true;
    if (n == 1 || (n & 1) == 0)
        return false;

    long long b = n - 1;
    for (int i = 0; i < 5; i ++){
        if (a[i] >= n)
            break;
        while((b & 1) == 0)    b >>= 1;
        long long t = exp_mod(a[i], b, n);
        while(b != n - 1 && t != 1 && t != n - 1){
            t = mul_mod(t, t, n);
            b <<= 1;
        }
        if (t == n - 1 || (b & 1))
            continue;
        else
            return false;
    }
    return true;
}
/*-------------Miller-Rabin 素数测试 部分--------------*/

/*-------------pollard-rho 大整数n因子分解 部分(用到mul_mod()和Miller-Rabin测试)--------------*/
long long factor[100];      //存n的素因子
long long nfactor, minfactor;

long long gcd(long long a, long long b){
    return b ? gcd(b, a%b) : a;
}
void Factor(long long n);
void pollard_rho(long long n){
    if (n <= 1)
        return ;
    if (Miller_Rabin(n)){
        factor[nfactor ++] = n;
        if (n < minfactor)
            minfactor = n;
        return ;
    }
    long long x = 2 % n, y = x, k = 2, i = 1;
    long long d = 1;
    while(true){
        i ++;
        x = (mul_mod(x, x, n) + 1) % n;
        d = gcd((y - x + n) % n, n);
        if (d > 1 && d < n){
            pollard_rho(d);
            pollard_rho(n/d);
            return ;
        }
        if (y == x){
            Factor(n);
            return ;
        }
        if (i == k){
            y = x;
            k <<= 1;
        }
    }
}
void Factor(long long n){
    //有时候RP不好 or n太小(比如n==4就试不出来……)用下面的pollard_rho没弄出来，则暴力枚举特殊处理一下
    long long d = 2;
    while(n % d != 0 && d * d <= n)
        d ++;
    pollard_rho(d);
    pollard_rho(n/d);
}
/*-------------pollard-rho 大整数n因子分解 部分--------------*/

vector  > vfac;
void find(long long n, long long &a, long long &b){
    long long sum = (1LL << 62);
    long long suma;
    for (int i = 0; i < (1 << vfac.size()); i ++){
        long long res = 1;
        for (size_t k = 0; k < vfac.size(); k ++){
            if (i & (1 << k)){
                for (int p = 0; p < vfac[k].second; p ++){
                    res *= vfac[k].first;
                }
            }
        }
        long long remain = n / res;
        if (res + remain < sum){
            sum = res + remain;
            a = suma = res;
            b = remain;
        }
        if (res + remain == sum){
            if (res < suma){
                a = suma = res;
                b = remain;
            }
        }
    }
    return ;
}
int main(){

    long long g, l, n;
    while(cin >> g >> l){
        n = l / g;
        vfac.clear();
        nfactor = 0;
        pollard_rho(n);
        long long tmp = n;
        for(int i = 0; i < nfactor; i ++){
            int facnum = 0;
            while(tmp % factor[i] == 0){
                tmp /= factor[i];
                facnum ++;
            }
            vfac.push_back(make_pair(factor[i], facnum));
        }
        long long a, b;
        find(n, a, b);
        cout << a * g << " " << b * g << endl;
    }
    return 0;
}