首页 > 代码库 > POJ 2429

POJ 2429

使用Pollard_rho算法分解lcm/gcd的质因数,原因不说也明白了。然后排序,把相同的质因子合并,因为如果相同的质因子分落在两个因数,会使ab的GCD值改变。

然后,枚举各种组合,呃。。。。这个实在想不到好方法,只好枚举了,真想不明白,那么两位数时间的是怎么样做到的。

 

#include <iostream>#include <cstdio>#include <algorithm>#include <cstring>#include <stdlib.h>#include <time.h>#define LL __int64#define C 201using namespace std;const int Max=100000;int atop;LL ans[100];bool se[Max];bool cmp(LL a,LL b){	if(a<b) return true;	return false;}LL gcd(LL a,LL b){	if(b==0) return a;	return gcd(b,a%b);}LL random(LL n){	return (LL)((double)rand()/RAND_MAX*n+0.5);}LL multi(LL a,LL b,LL m){	LL ret=0;	while(b>0){		if(b&1) ret=(ret+a)%m;		b>>=1;		a=(a<<1)%m;	}	return ret;}LL quick(LL a,LL b,LL m){	LL ans=1;	a%=m;	while(b){		if(b&1){			ans=multi(ans,a,m);		}		b>>=1;		a=multi(a,a,m);	}	return ans;}LL Pollard_rho(LL n,int c){	LL x,y,d,i=1,k=2;	x=random(n-1)+1;	y=x;	while(true){		i++;		x=(multi(x,x,n)+c)%n;		d=gcd(y-x,n);		if(d>1&&d<n) return d;		if(y==x) return n;		if(i==k){			y=x;			k=k<<1;		}	}}bool Witness(LL a,LL n){	LL m=n-1;	int j=0;	while(!(m&1)){		j++;		m=m>>1;	}	LL x=quick(a,m,n);	if(x==1||x==n-1)	return false;	while(j--){		x=multi(x,x,n);		if(x==n-1) return false;	}	return true;}bool Miller_Rabin(LL n){	if(n<2) return false;	if(n==2) return true;	if(!(n&1)) return false;	for(int i=1;i<=12;i++){		LL a=random(n-2)+1;		if(Witness(a,n)) return false;	}	return true;}void find(LL n, int k){	if(n==1) return ;	if(Miller_Rabin(n)){		ans[atop++]=n;		return ;	}	LL p=n;	while(p>=n)	p=Pollard_rho(p,k--);	find(p,k);	find(n/p,k);}LL getrt(int g){	LL ret=1;	for(int i=0;i<atop;i++){		if(g&(1<<i)){			ret*=ans[i];		}	}	return ret;}int main(){	LL gcd,lcm;	srand(time(0));	while(scanf("%I64d%I64d",&gcd,&lcm)!=EOF){		lcm/=gcd; atop=0;		memset(se,false,sizeof(se));		find(lcm,C);		sort(ans,ans+atop,cmp);		int m=atop;		atop=1;		ans[0]=ans[0];		for(int i=1;i<m;i++){			if(ans[i]==ans[i-1]){				ans[atop-1]*=ans[i];			}			else if(ans[i]!=ans[i-1]){				ans[atop++]=ans[i];			}		}	//	cout<<atop<<endl;		int all=(1<<atop)-1;		LL at=-1;		LL sum=(1LL<<63)-1;	//	cout<<sum<<endl;		LL ret;		for(int i=0;i<=all;i++){			if(se[i]) continue;			se[i]=se[all^i]=true;			ret=getrt(i);		//	cout<<ret<<endl;			if(sum>ret+lcm/ret){			//	cout<<ret<<endl;				at=ret;				sum=ret+lcm/ret;			}		}		LL a=gcd*at,b=gcd*(lcm/at);		printf("%I64d %I64d\n",a<b?a:b,a>b?a:b);	}	return 0;}

  

POJ 2429