POJ 2429 GCD & LCM Inverse（素数判定Miller-Rabin+素因子分解Pollard-rho）

Description

给出gcd(a,b)和lcm(a,b)，求a和b，如果存在多组方案则输出a+b最小的那一组

Input

两个整数gcd(a,b)和lcm(a,b)，数值均不超过2^63，保证有解

Output

输出满足条件的a和b(a<=b)，如果有多组方案则输出a+b最小的那一组

Sample Input

3 60

Sample Output

12 15

Solution

令g=gcd(a,b),l=lcm(a,b)，则有gcd(a/g,b/g)=1,(a/g)*(b/g)=l/g，所以问题转化为求t=l/g所有的素因子和每个素因子对应的幂指数，但由于t可能会非常大，所以正常的素因子分解会超时，此处需要用Pollard-rho来分解大数，而在素数判定时，同样的由于数会很大所以要用Miller-Rabin判定，对于比较小的数可以首先打一个素数表来判断，当对t素因子分解后，2^res枚举a/g应含的素因子（res为t的素因子数）更新最优解即可

Code

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<vector>
#include<map>
#include<ctime>
using namespace std;
#define maxn 200000
typedef long long ll;
bool is_prime[maxn];
vector<int>prime;
map<ll,int>factor;
void get_prime()
{
for(int i=0;i<maxn;i++)is_prime[i]=1;
is_prime[0]=is_prime[1]=0;
for(int i=2;i<maxn;i++)
if(is_prime[i])
{
prime.push_back(i);
for(int j=i;j<maxn;j+=i)is_prime[j]=0;
}
}
ll gcd(ll a,ll b)
{
if(b==0)return a;
return gcd(b,a%b);
}
ll mod_mul(ll a,ll b,ll p)
{
ll ans=0ll;
a%=p,b%=p;
if(a>b) swap(a,b);
while(b)
{
if(b&1) ans=(ans+a)%p;
a=(a+a)%p;
b>>=1;
}
return ans;
}
ll mod_pow(ll a,ll b,ll p)
{
ll ans=1ll;
a%=p;
while(b)
{
if(b&1) ans=mod_mul(ans,a,p);
a=mod_mul(a,a,p);
b>>=1;
}
return ans;
}
int Times=20;
bool witness(ll a,ll n)
{
ll m=n-1;
int j=0;
while((m&1)==0)
{
j++;
m>>=1;
}
ll x=mod_pow(a,m,n);
if(x==1||x==n-1) return false;
while(j--)
{
x=mod_mul(x,x,n);
if(x==n-1) return false;
}
return true;
}
bool Miller_Rabin(ll n)
{
srand(time(0));
if(n<2) return 0;
if(n==2) return 1;
if((n&1)==0) return 0;
for(int i=0;i<Times;i++)
{
ll a=rand()%(n-1)+1;
if(witness(a,n)) return false;
}
return true;
}
ll Pollard_Rho(ll n,int c)
{
ll x=2,y=2,d=1;
while(d==1)
{
x=mod_mul(x,x,n)+c;
y=mod_mul(y,y,n)+c;
y=mod_mul(y,y,n)+c;
d=gcd((x-y>=0?x-y:y-x),n);
}
if(d==n) return Pollard_Rho(n,c+1);
return d;
}
bool Is_Prime(ll n)
{
if(n<maxn&&is_prime
||n>=maxn&&Miller_Rabin(n))return 1;
return 0;
}
void Find_Factor(ll n)
{
if(Is_Prime(n))
{
factor
++;
return;
}
for(int i=0;i<prime.size()&&prime[i]<=n;i++)
if(n%prime[i]==0)
{
while(n%prime[i]==0)
{
factor[prime[i]]++;
n/=prime[i];
}
}
if(n!=1)
{
if(Is_Prime(n))
factor
++;
else
{
ll p=Pollard_Rho(n,1);
Find_Factor(p);
Find_Factor(n/p);
}
}
}
int main()
{
ll l,g;
get_prime();
while(~scanf("%lld%lld",&g,&l))
{
ll t=l/g;
factor.clear();
Find_Factor(t);
ll ans=1e18,x=1,y=t;
vector<ll>fact;
for(map<ll,int>::iterator it=factor.begin();it!=factor.end();it++)
{
ll temp=1ll;
while(it->second)
{
temp*=it->first;
it->second--;
}
fact.push_back(temp);
}
for(int i=0;i<(1<<fact.size());i++)
{
ll tx=1ll,ty;
for(int j=0;j<fact.size();j++)
if(i&(1<<j))tx*=fact[j];
ty=t/tx;
if(tx<ty&&ans>tx+ty)
ans=tx+ty,x=tx,y=ty;
}
printf("%lld %lld\n",x*g,y*g);
}
return 0;
}