POJ 1811 （Miller_rabin+pollard_rho）

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#include<cstdlib>
#include<cmath>
#define maxn 0x7fffffff

using namespace std;

typedef long long ll;

int cnt;
ll fac[100];
ll aa,bb;
ll Max;
ll f[100];

ll random(ll k){
    return (ll)((double)rand()/RAND_MAX*k);
}

ll mulit(ll a,ll b,ll m){
    ll ans=0;
    while(b){
        if(b&1) ans=(ans+a)%m;
        b>>=1;
        a=(a+a)%m;
    }
    return ans;
}

ll quick_mod(ll a,ll b,ll m){
    ll ans=1;
//    a%=m;
    while(b){
        if(b&1) ans=mulit(ans,a,m);
        b>>=1;
        a=mulit(a,a,m);
    }
    return ans;
}

bool witness(ll k,ll n){
    ll m=n-1;
    int j=0;
    while(!(m&1)){
        m>>=1;
        j++;
    }
    ll x=quick_mod(k,m,n);
    if(x==1||x==n-1)
        return false;
    while(j--){
        x=mulit(x,x,n);
        if(x==n-1)
            return false;
    }
    return true;
}

bool miller(ll n){
    if(n==2) return true;
    if(n<2||!(n&1)) return false;
    for(int i=0;i<12;i++){
        ll t=random(n-2)+1;
        if(witness(t,n))
            return false;
    }
    return true;
}

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

ll pollard_rho(ll n,int c){
    ll x,y,i=1,j=2;
    x=random(n-2)+1;
    y=x;
    while(1){
        i++;
        x=(mulit(x,x,n)+c)%n;
        ll d=gcd((x-y+n)%n,n);
        if(d>1&&d<n) return d;
        if(y==x)        //失败，成为一个环
            return n;
        //Brent判环
        if(i==j){
            j<<=1;
            y=x;
        }
    }
}

void find(ll n,int k){
    if(n==1)
        return ;
    if(miller(n)){
        fac[cnt++]=n;
        return ;
    }
    ll p=n;
    while(p>=n){
        p=pollard_rho(p,k++);
    }
    find(p,k);
    find(n/p,k);
}

void dfs(ll a,ll b,int c,int len){
    if(a+b>=Max) return;
    if(c==len+1){
        if(a+b<Max){
            Max=a+b;
            aa=a;
            bb=b;
        }
        return ;
    }
    dfs(a*f[c],b,c+1,len);
    dfs(a,b*f[c],c+1,len);
}

int main()
{
    ll a,b;
    ll gcd,lcm;
    while(scanf("%lld%lld",&gcd,&lcm)!=EOF){
        memset(fac,0,sizeof(fac));
        memset(f,0,sizeof(0));
        cnt=0;
        ll n=lcm/gcd;
        if(miller(n)){
            a=1;
            b=n;
            printf("%lld %lld\n",a*gcd,b*gcd);
        }
        else if(lcm==gcd)
            printf("%lld %lld\n",gcd,gcd);
        else{
            find(n,1);
            sort(fac,fac+cnt);
            f[0]=fac[0];
            int j=0;
            for(int i=1;i<cnt;i++){
                if(fac[i]==fac[i-1])
                    f[j]*=fac[i];
                else{
                    f[j+1]=fac[i];
                    j++;
                }
            }
            aa=f[0];
            bb=n/f[0];
            Max=aa+bb;
            dfs(f[0],1,1,j);
            if(aa>bb)
                swap(aa,bb);
            printf("%lld %lld\n",aa*gcd,bb*gcd);
        }

    }
    return 0;
}