POJ 2429 GCD & LCM Inverse(Pollard_rho法质因数分解）

http://poj.org/problem?id=2429

题意：已知两个数的最大公约数和最小公倍数，求原来这两个数。如果有多种可能，输出两数之和最小的组合。

把lcm/gcd分解成两个和最小的互质的数，这两个数再各自乘以gcd就是答案。

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<ctime>
#include<algorithm>
#define LL unsigned long long // 这道题说给出的gcd & lcm 都小于 2^63，为了避免运算过程中溢出，开到最大
#define INF ((LL)1<<61) // ((LL)(1<<61)) 这样是错的
#define N 1000
using namespace std;
LL ans[N], cnt, C[N], Ans[N], G, L, X, Y, Min;
LL gcd(LL a, LL b)
{
    if ( b == 0 )
        return a;
    return gcd ( b, a % b );
}
LL mod_mult (LL a, LL b, LL n)
{
    LL ret = 0;
    a = a % n;
    while ( b >= 1 )
    {
        if ( b & 1 )
		{
            ret += a;
            if ( ret >= n ) ret -= n;
        }
        a = a << 1;
        if ( a >= n ) a -= n;
        b = b >> 1;
    }
    return ret;
}
LL mod_exp (LL a, LL b, LL n)
{
    LL ret = 1;
    a = a % n;
    while ( b >= 1 )
    {
        if ( b & 1 )
            ret =  mod_mult(ret,a,n);
        a = mod_mult(a,a,n);
        b = b >> 1;
    }
    return ret;
}
bool witness (LL a, LL n)
{
	int i, t = 0;
    LL m = n - 1, x, y;
    while ( m % 2 == 0 ) { m >>= 1; t++; }
    x = mod_exp (a, m, n);
    for ( i = 1; i <= t; i++ )
    {
        y = mod_exp ( x, 2, n );
        if( y==1 && x!=1 && x!=n-1 )
            return true;
        x = y;
    }
    if ( y != 1 ) return true;
    return false;
}
bool miller_rabin (LL n, int times = 10 ) // 判断一个较大的n是不是素数
{
	if ( n == 2 ) return true;
    if ( n == 1 || n % 2 == 0 ) return false;
   // srand ( time(NULL) ); POJ上提交的时候必须把这句话删掉
    for ( int i = 1; i <= times; i++ )
    {
        LL a = rand() % (n-1) + 1;
        if ( witness(a,n) ) return false;
    }
    return true;
}
LL rho (LL n, int c)
{
    LL i, k, x, y, d;
 // srand ( time(NULL) );
    i = 1;  k = 2;
    y = x = rand() % n;
    while (true)
    {
        i++;
        x = (mod_mult(x,x,n)+c) % n;
        d = gcd (y - x, n);
        if (d > 1 && d < n) return d;
        if (y == x) break;
        if (i == k) {y = x; k *= 2;}
	}
	return n;
}
void pollard (LL n, int c) // 将n进行质因数分解，分解出来的数总共有cnt个，保存在ans（0~cnt-1)中，不按照升序排列
{
    if (n == 1)  return;
    if (miller_rabin(n)) {ans[cnt++] = n; return;}
    LL m = n;
    while (m >= n)
        m = rho (m, c--);
    pollard (m, c);
    pollard (n / m, c);
}
void dfs(LL a, LL b, LL p)
{
    if (a+b >= Min) return; //不断地把因子加到两堆中的某一堆，如果当前两堆的和已经超过最小值，那么后面的只可能更大
    if (p == cnt)
    {
        if (a+b < Min)
        {
            Min = a+b;
            X = a*G;
            Y = b*G;
            if (X > Y) swap(X, Y);
        }
        return;
    }
    dfs(a*Ans[p], b, p+1);
    dfs(a, b*Ans[p], p+1);
}
int main()
{
    while (scanf("%llu%llu",&G,&L) != EOF)
    {
        if (L == G)
        {
            printf("%llu %llu\n",L,G);
            continue;
        }
        L /= G;
        cnt = 0;
        pollard (L, 107);
        sort (ans, ans + cnt);
        LL j = 0;
        Ans[j] = ans[j];
        for (int i = 1; i < cnt; i++)
            if (ans[i] == ans[i-1])
        {
            Ans[j] *= ans[i];
        }
        else
        {
            j++;
            Ans[j] = ans[i];
        }
        cnt = j+1;//上一段的目的是，为了保证分解出来的两个数是互质的，所以把相同的质因子乘起来，让它们不可能在dfs过程中被分到两堆
        C[0] = -1;
        Min = INF;
        dfs(1, 1, 0);
        printf("%llu %llu\n", X, Y);
    }
    return 0;
}