GCD & LCM Inverse（poj 2429）

题意：给出$x$和$y$的$gcd$和$lcm$，让你求出$x+y$最小的$x，y$

思路：首先$\frac{(xy)}{gcd}=lcm$,可以构造$\frac{x}{gcd}\times \frac{y}{gcd}=\frac{lcm}{gcd}$，其中$\frac{x}{gcd},\frac{y}{gcd}$互质，假设$n=\frac{lcm}{gcd}$，那么只需分解n的因子，通过n的因子构造两个互质的数即可，因为两个数互质所以两个数不可能有公共的质因子，用dfs枚举取x+y最小的情况即可，这题的范围1e18，分解质因数需要使用miller_robin和pollard算法。

代码：

#include <iostream>
#include <algorithm>
#include <math.h>
#include <cstring>
#include <stdio.h>

using namespace std;
typedef unsigned long long  ll;
const int maxx=1e6+7,times=11;
ll a[maxx],cnt,minn=1e18,k=201;
ll random(ll n)
{
    return (ll)((double)rand()/RAND_MAX*n+0.5);
}
ll multi(ll a,ll b,ll mod)
{
    ll ans=0;
    while(b)
    {
        if(b&1)
            ans=(ans+a)%mod;
        b>>=1;
        a=(a<<1)%mod;
    }
    return ans;
}
ll qpow(ll a,ll b,ll mod)
{
    ll ans=1;
    while(b)
    {
        if(b&1)
            ans=multi(ans,a,mod);
        a=multi(a,a,mod);
        b>>=1;
    }
    return ans;
}
bool witness(ll a,ll mod)
{
    ll m=mod-1;
    int j=0;
    while(!(m&1))
    {
        j++;
        m>>=1;
    }
    ll x=qpow(a,m,mod);
    if(x==1||x==mod-1)
        return 0;
    while(j--)
    {
        x=x*x%mod;
        if(x==mod-1)
            return 0;
    }
    return 1;
}
bool miller_rabin(ll n)
{
    if(n<2) return 0;
    if(n==2) return 1;
    if(!(n&1)) return 0;
    for(int i=1;i<=times;i++)
    {
        ll a=random(n-2)+1;
        if(witness(a,n)) return 0;
    }
    return 1;
}
ll pollard_rho(ll n,ll c)
{
    ll x,y,d,i=1,k=2;
    x=random(n-1)+1;
    y=x;
    while(1)
    {
        i++;
        x=(multi(x,x,n)+c)%n;
        d=__gcd(y-x,n);
        if(1<d&&d<n)
            return d;
        if(x==y)
            return n;
        if(i==k)
            y=x,k<<=1;
    }
}
void finds(ll n,ll k)
{
    if(n==1) return;
    if(miller_rabin(n))
    {
        a[++cnt]=n;
        if(n<minn)
            minn=n;
        return;
    }
    ll p=n;
    while(p>=n)
        p=pollard_rho(p,k--);
    finds(p,k);
    finds(n/p,k);
}
ll p[maxx],len,x=2*1e18,y=2*1e18;
void dfs(int now,ll a,ll b)
{
    if(a+b>x+y) return;
    if(now==len+1)
    {
        if(a+b<x+y)
            x=a,y=b;
        return;
    }
    dfs(now+1,a*p[now],b);
    dfs(now+1,a,b*p[now]);
}
int main()
{
    ll n,gcd,lcm;
    while(cin>>gcd>>lcm)
    {
        len=0,cnt=0,x=2*1e18,y=2*1e18;
        n=lcm/gcd;
        finds(n,k);
        sort(a+1,a+1+cnt);
        for(int i=1;i<=cnt;i++)
        {
            if(a[i]==a[i-1])
                p[len]*=a[i];
            else p[++len]=a[i];
        }
        //for(int i=1;i<=len;i++)
        //    cout<<p[i]<<endl;
        dfs(1,1,1);
        if(x>y) swap(x,y);
        cout<<x*gcd<<" "<<y*gcd<<endl;
    }
}