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,981if (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;elsereturn 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;
}