SCU 2016 GCD & LCM Inverse(素性测试+DFS）

Description
We are all familiar with the Euclid's algorithm. Given two positive integers a, b. we can easily find the greatest common divisor (GCD) and least common multiple (LCM) of a, b through this method. But what about
the inverse? That is: given GCD and LCM of a and b, can you find a, b ?
Input
The input will contain multiple test cases, each of which contains two positive integers: GCD and LCM. You can safely assume these integers will be less
than 2^63-1.
Output
For each test case, output a, b in ascending order. If there are multiple solutions, output the pair with smallest a+b.
Sample Input

3 60

Sample Output

12 15

Source
分析：设a=gcd*k1,b=gcd*k2,那么k1,k2互质我们可以有：k1*k2=lcm/gcd
2^63分解后最多的素数个数也在20以内，搜一下就行了

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
#include<string>
#include<iostream>
#include<queue>
#include<cmath>
#include<map>
#include<set>
#include<time.h>
using namespace std;
#define REPF( i , a , b ) for ( int i = a ; i <= b ; ++ i )
#define REP( i , n ) for ( int i = 0 ; i < n ; ++ i )
#define CLEAR( a , x ) memset ( a , x , sizeof a )
typedef long long LL;
typedef pair<int,int>pil;
const int maxn=1e6+100;
const int mod=1e9+7;
//****************************************************************
// Miller_Rabin 算法进行素数测试
//速度快，而且可以判断 <2^63的数
//****************************************************************
const int S=20;
LL mult_mod(LL a,LL b,LL c)
{
    a%=c;
    b%=c;
    LL ret=0;
    while(b)
    {
        if(b&1){ret+=a;ret%=c;}
        a<<=1;
        if(a>=c)a%=c;
        b>>=1;
    }
    return ret;
}
LL pow_mod(LL x,LL n,LL mod)//x^n%c
{
    if(n==1)return x%mod;
    x%=mod;
    LL tmp=x;
    LL ret=1;
    while(n)
    {
        if(n&1) ret=mult_mod(ret,tmp,mod);
        tmp=mult_mod(tmp,tmp,mod);
        n>>=1;
    }
    return ret;
}
bool check(LL a,LL n,LL x,LL t)
{
    LL ret=pow_mod(a,x,n);
    LL last=ret;
    for(int i=1;i<=t;i++)
    {
        ret=mult_mod(ret,ret,n);
        if(ret==1&&last!=1&&last!=n-1) return true;
        last=ret;
    }
    if(ret!=1) return true;
    return false;
}
bool Miller_Rabin(LL n)
{
    if(n<2)return false;
    if(n==2)return true;
    if((n&1)==0) return false;
    LL x=n-1;
    LL t=0;
    while((x&1)==0){x>>=1;t++;}
    for(int i=0;i<S;i++)
    {
        LL a=rand()%(n-1)+1;
        if(check(a,n,x,t))
            return false;
    }
    return true;
}
LL fac[100];
int tot;
LL gcd(LL a,LL b)
{
    if(a==0)return 1;
    if(a<0) return gcd(-a,b);
    while(b)
    {
        LL t=a%b;
        a=b;
        b=t;
    }
    return a;
}

LL Pollard_rho(LL x,LL c)
{
    LL i=1,k=2;
    LL x0=rand()%x;
    LL y=x0;
    while(1)
    {
        i++;
        x0=(mult_mod(x0,x0,x)+c)%x;
        LL d=gcd(y-x0,x);
        if(d!=1&&d!=x) return d;
        if(y==x0) return x;
        if(i==k){y=x0;k+=k;}
    }
}
void findfac(LL n)
{
    if(Miller_Rabin(n))//素数
    {
        fac[tot++]=n;
        return;
    }
    LL p=n;
    while(p>=n)p=Pollard_rho(p,rand()%(n-1)+1);
    findfac(p);
    findfac(n/p);
}
int vis[100],ans[100],num[100];;
LL sum,ans1,ans2,d;
void dfs(int pos,int s)
{
    if(s)
    {
        LL x=1;
        for(int i=0;i<s;i++)
        {
            for(int j=0;j<num[ans[i]];j++)
                x=x*fac[ans[i]];
        }
        LL y=d/x;
        if(x+y<sum)
        {
            sum=x+y;
            ans1=x;
            ans2=y;
        }
    }
    for(int i=pos;i<tot;i++)
    {
        if(!vis[i])
        {
            vis[i]=1;
            ans[s]=i;
            dfs(i+1,s+1);
            vis[i]=0;
        }
    }
}
int main()
{
    LL a,b;
    while(~scanf("%lld%lld",&a,&b))
    {
        if(a==b)
        {
            printf("%lld %lld\n",a,b);
            continue;
        }
        tot=0;
        findfac(b/a);
        CLEAR(num,0);
        LL s=b/a;d=b/a;
        for(int i=0;i<tot;i++)
        {
            while(s%fac[i]==0)
            {
                num[i]++;
                s/=fac[i];
            }
        }
        sum=1e18;
        dfs(0,0);
        printf("%lld %lld\n",min(ans1,ans2)*a,max(ans1,ans2)*a);
    }
    return 0;
}