POJ2429 GCD & LCM Inverse pollard_rho大整数分解

Given two positive integers a and b, we can easily calculate the greatest common divisor (GCD) and the least common multiple (LCM) of a and b. But what about the inverse? That is: given GCD and LCM, finding a and b.

Input

The input contains multiple test cases, each of which contains two positive integers, the GCD and the LCM. You can assume that these two numbers are both less than 2^63.

Output

For each test case, output a and b in ascending order. If there are multiple solutions, output the pair with smallest a + b.

Sample Input

3 60

Sample Output

12 15

题意：给出两个数的最大公约数和最小公倍数，让你找出满足条件的两个数，使他们的和最小。

题解：

对于两个数a,b和他们的最大公约数gcd以及最小公倍数lcm，有lcm=a*b/gcd，进一步变形可以得到：（a/gcd * b/gcd）*gcd=lcm，即(a/gcd * b/gcd)=lcm/gcd，我们令m=a/gcd,n=b/gcd，那么问题就变为了：找出两个互素的整数，使他们的乘积为key=lcm/gcd。我们可以用Pollard_rho算法进行质因素分解的方法找出lcm/gcd的所有素因子，然后dfs找出任意几个素数因子之积(设为x)，设y=key/x，找出最小的x+y就行了

代码：

#include <cstdio>

#include <iostream>

#include <cstdlib>

#include <cmath>

#include <algorithm>

 

#define times 10

#define N 501

using namespace std;

typedef unsigned long long LL;

const LL INF=(LL)1<<61;

LL key,ct,cnt,gd,lm,resa,resb,mini;

LL fac[N],num[N];

 

LL gcd(LL a,LL b)

{

    return b?gcd(b,a%b):a;

}

 

LL multi(LL a,LL b,LL m)

{

    LL ans=0;

    a%=m;

    while(b)

    {

        if(b&1)

        {

            ans=(ans+a)%m;

            b--;

        }

        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=multi(ans,a,m);

            b--;

        }

        b>>=1;

        a=multi(a,a,m);

    }

    return ans;

}

 

bool Miller_Rabin(LL n)

{

    if(n==2) return true;

    if(n<2||!(n&1)) return false;

    LL m=n-1;

    int k=0;

    while(!(m&1))

    {

        k++;

        m>>=1;

    }

    for(int i=0;i<times;i++)

    {

        LL a=rand()%(n-1)+1;

        LL x=quick_mod(a,m,n);

        LL y=0;

        for(int j=0;j<k;j++)

        {

            y=multi(x,x,n);

            if(y==1&&x!=1&&x!=n-1) return false;

            x=y;

        }

        if(y!=1) return false;

    }

    return true;

}

 

LL Pollard_rho(LL n,LL c)

{

    LL i=1,k=2;

    LL x=rand()%(n-1)+1;

    LL y=x;

    while(true)

    {

        i++;

        x=(multi(x,x,n)+c)%n;

        LL d=gcd((y-x+n)%n,n);

        if(1<d&&d<n) return d;

        if(y==x) return n;

        if(i==k)

        {

            y=x;

            k<<=1;

        }

    }

}

 

void Find(LL n,LL c)

{

    if(n==1) return ;

    if(Miller_Rabin(n))

    {

        fac[ct++]=n;

        return ;

    }

    LL p=n;

    LL k=c;

    while(p>=n) p=Pollard_rho(p,c--);

    Find(p,k);

    Find(n/p,k);

}

 

void dfs(LL dept,LL product)

{//dept为递归深度，product为其中一个因子

    if(dept==cnt)

    {

        LL a=product;

        LL b=key/a;

        if(gcd(a,b)==1)

        {

            a*=gd;

            b*=gd;

            if(a+b<mini)

            {

                mini=a+b;

                resa=a;

                resb=b;

            }

        }

        return ;

    }

    for(int i=0;i<=num[dept];i++)

    {

        if(product>mini) return ;

        dfs(dept+1,product);

        product*=fac[dept];

    }

}

 

 

void Solve(LL n)

{

    ct=0;

    Find(n,120);

    sort(fac,fac+ct);

    num[0]=1;

    int k=1;

    for(int i=1;i<ct;i++)

    {

        if(fac[i]==fac[i-1])

            num[k-1]++;

        else

        {

            num[k]=1;

            fac[k++]=fac[i];

        }

    }

    cnt=k;

    dfs(0,1);

    if(resa>resb) swap(resa,resb);

}

 

int main()

{

    while(cin>>gd>>lm)

    {

        if(gd==lm)

        {

            printf("%llu %llu
",gd,lm);

            continue;

        }

        mini=INF;

        key=lm/gd;

        Solve(key);

        printf("%llu %llu
",resa,resb);

    }

    return 0;

}