POJ/PKU 2429 GCD & LCM

　　题目说给出最大公约数和最小公倍数，求a，b。如果存在多个解取最小a+b。

　　思路很简单，a/GCD和b/GCD是互质的，所以a/GCD和b/GCD分解以后产生的因子肯定不同。所以他们的乘积a/GCD*(b/GCD)=a*b/GCD/GCD=LCM/GCD因数分解以后产生的任何一种质因数只能属于a/GCD或b/GCD。很明显种数不会超过15，深搜或者位运算很容易求出所有组合，取其中a+b最小的a，b即可。

　　质因数分解直接用PollardRho的模板，但是在G++下仍然会re，可能是编译器的原因。

/*
 * File:   main.cpp
 * Author: acmer
 *
 * Created on June 11, 2011, 3:55 PM
 */

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <algorithm>

using namespace std;

#define Random(n) (rand()%(n+1))

typedef __int64 int64;

const int kMaxT(50);

int cnt;
int64 factor[207];

int64 Gcd(int64 a,int64 b){
    for(int64 t=a%b;t;a=b,b=t,t=a%b);
    return b>0?b:(b*(-1));
}

int64 MutiMod(int64 a,int64 b,int64 n){
    int64 exp(a%n),res(0);
    while(b){
        if(b&1){
            res+=exp;
            if(res>n)res-=n;
        }
        exp<<=1;
        if(exp>n)exp-=n;
        b>>=1;
    }
    return res;
}

int64 ExpMod(int64 a,int64 n,int64 b){
    int64 r(1),t(a%b);
    if(n==0) return 1%b;
    while(n>1){
        if(n&1)
            r=MutiMod(r,t,b);
        t=MutiMod(t,t,b);
        n>>=1;
    }
    return MutiMod(r,t,b);
}

bool MillerRabbin(int64 n){
    if(n==2)
        return true;
    if(n<2 || !(n&1))
        return false;
    int64 a,u(n-1),x,y;
    int t(0);
    while(u%2==0){
        t++;
        u>>=1;
    }
    srand(time(NULL));
    for(int i=1;i<=kMaxT;i++){
        a=Random(n-2)+1;
        x=ExpMod(a,u,n);
        for(int j=0;j<t;j++){
            y=MutiMod(x,x,n);
            if(y==1 && x!=1 && x!=n-1)
                return false;
            x=y;
        }
        if(y!=1)
            return false;
    }
    return true;
}

int64 PollardRho(int64 n,int c){
    int64 x(Random(n-2)+1),y(x),d,i(1),k(2);
    while (true){
        i++;
        x=(MutiMod(x,x,n)+c)%n;
        d=Gcd(y-x,n);
        if(d > 1 && d < n)
            return d;
        if(x==y)
            return n;
        if(i==k){
            y=x;
            k<<=1;
        }
    }
}

int64 temp;

//找出所有质因数
void FindFactor(int64 n, int k){
    if(n == 1)
        return;
    if(MillerRabbin(n)){
        //保存小于n的所有质因数
        if (n != temp){
            factor[cnt++] = n;
        }
        return;
    }
    int64 p(n);
    while (p >= n)
        p = PollardRho(p,k--);
    FindFactor(p, k);
    FindFactor(n/p,k);
}

int l;
int64 t[67], minf, a, n;

//求最小的a+b
void dfs(int i, int64 sum){
    if (i == l + 1){
        if (minf > (sum + n / sum)){
            minf = sum + n / sum;
            a = sum;
        }
        return;
    }
    dfs(i+1, sum * t[i]);
    dfs(i+1, sum * 1);
    return;
}

int main(){

    int64 LCM, GCD;

    while (cin>>GCD>>LCM){

        n = LCM / GCD;
        if (n == 1 || MillerRabbin(n)){
            cout<<GCD<<" "<<LCM<<endl;
            continue;
        }

        cnt = 0;
        temp = n;

        FindFactor(n, 107);
        sort(factor, factor + cnt);

        l = 0;
        t[l] = factor[0];
        int j = 0;
        for (int i = 1; i < cnt; i++)
            if (factor[i] == factor[j])
                t[l] *= factor[i];
            else {
                t[++l] = factor[i];
                j = i;
            }
        minf = 2e63;//科学计数法越界以后会自动换为此范围内最大数
        dfs(0, 1);

        cout<<min(a*GCD, LCM/a)<<" "<<max(a*GCD, LCM/a)<<endl;
    }
    return 0;
}