poj 2429 GCD & LCM Inverse

题目链接：http://poj.org/problem?id=2429

题意：给你两个数n和m,使得gcd(a,b)=n, lcm(a,b)=m. 让你求a,b.如果a,b有多组，输出a+b最小的一组。（gcd表示a,b的最大公约数，lcm表示a,b的最小公倍数）

思路：因为a,b的值都可能很大，所以要采用pollard_rho 算法进行质因素分解。

因为lcm(a,b)=a*b/gcd(a,b) ,  则a/gcd(a,b) * b/gcd(a,b)=lcm(a,b)/gcd(a,b);

可以设t=lcm(a,b)/gcd(a,b),需对t进行质因素分解，并且相同的质因子应该放在一起，如果相同的质因子不在一起，则a,b的最大公约数就不可能是n。然后对其进行任意组合，求出其中a+b最小的一组。

因为a*b是一定的为gcd(a,b)*lcm(a,b); 要求a+b最小，则a/gcd(a,b) 应该最接近sqrt(t)，并且小于sqrt(t)。

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<time.h>
#include<algorithm>
using namespace std;
const int S=8;//随机算法判定次数,一般8~10次就够了
long long factor[100];//质因素分解结果（刚返回时是无序的）
int tol;//质因素的个数,编号0~tol-1
long long mult_mod(long long a,long long b,long long c)
{
a%=c;
b%=c;
long long ret=0;
long long tmp=a;
while(b)
{
if(b&1)
{
ret+=tmp;
if(ret>c)
ret-=c;
}
tmp<<=1;
if(tmp>c)
tmp-=c;
b>>=1;
}
return ret;
}
//计算ret=(a^n)%mod
long long pow_mod(long long a,long long n,long long mod)
{
long long ret=1;
long long temp=a%mod;
while(n)
{
if(n&1)
ret=mult_mod(ret,temp,mod);
temp=mult_mod(temp,temp,mod);
n>>=1;
}
return ret;
}
//*********************************************
long long gcd(long long a,long long b)
{
long long t;
while(b)
{
t=a;
a=b;
b=t%b;
}
if(a>0)
return a;
else
return -a;
}
//************************************************

//找出一个因子
//**************************************************
long long pillard_rho(long long x,long long c)
{
long long i=1,k=2;
srand(time(NULL));
long long x0=rand()%(x-1)+1;
long long y=x0;
while(1)
{
i++;
x0=(mult_mod(x0,x0,x)+c)%x;
long long d=gcd(y-x0,x);
if(d!=1 && d!=x)
return d;
if(y==x0)
return x;
if(i==k)
{
y=x0;
k+=k;
}
}
}
//**************************************************
bool check(long long a,long long n,long long x,long long t)
{
long long ret = pow_mod(a,x,n);
long long 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;
else
return false;
}
//Miller Rabin 算法
//是素数返回true,(可能是伪素数)
//不是素数返回false
bool Miller_Rabin(long long n)
{
if(n<2)
return false;
if(n==2)
return true;
if((n&1)==0)
return false;//偶数
long long x=n-1;
long long t=0;
while((x&1)==0)
{
x>>=1;
t++;
}
srand(time(NULL));
for(int i =0;i<S;i++)
{
long long a=rand()%(n-1)+1;
if(check(a,n,x,t))
return false;
}
return true;
}
//*********************************************

//对n进行素因子分解，存入factor.k设置为107左右即可

//**************************************************
void findfac(long long n,int k)
{
if(n==1)
return ;
if(Miller_Rabin(n))
{
factor[tol++]=n;
return;
}
long long p=n;
int c=k;
while(p>=n)
p=pillard_rho(p,c--);//值变化，防止死循环k
findfac(p,k);
findfac(n/p,k);
}

long long pow(int a,int n)//快速取幂a^n
{
if(n == 0)  return 1;
long long  x = pow(a,n/2);
long long ans = x * x;
if(n%2 == 1)
ans = ans * a;
return ans;
}

//**************************************************
int main()
{
long long n,m,t;
long long a[17];
long long num,ans;
int k=107;
int flag[17];
int cont;
while(~scanf("%I64d %I64d",&n,&m))
{
memset(factor,0,sizeof(factor));
t=m/n;
if(t==1){
printf("%lld %lld\n",n,m); //m=n,则这两个数就为m,n;
continue;
}
tol=0;
findfac(t,k);
sort(factor,factor+tol);//排序
cont=0;
a[cont]=factor[0];
for(int i=1;i<tol;i++)
{
if(factor[i]==factor[i-1]) a[cont]*=factor[i];//将素因子相同的值的积统计并记录下来
else a[++cont]=factor[i];
}
int min=0;
int max=(int )(sqrt(t*1.0)+0.5);
num=pow(2,cont+1);//求2^cont+1
//num=1<<(cont+1);
for(int i=0;i<num;i++){
memset(flag,0,sizeof(flag));
int j=0;
int p=i;
while(p){
flag[j++]=p%2;
p/=2;
}
ans=1;
for(int x=0;x<j;x++){
if(flag[x]) ans*=a[x];
}
if(ans<=max && ans>=min) min=ans;
}
if(t/min<min) min=t/min;
printf("%lld %lld\n",min*n,t/min*n);
}
return 0;
}