poj 2429 GCD & LCM Inverse miller_rabin素数判定和pollard_rho因数分解

题意：

给gcd(a,b)和lcm(a,b),求a+b最小的a和b。

分析：

miller_rabin素数判定要用费马小定理和二次探测定理。pollard_rho因数分解算法导论上讲的又全又好，网上的资料大多讲不清楚。

代码：

//poj 2429
//sep9
#include <iostream>
#include <map>
#include <vector>
#define gcc 10007
#define max_prime 200000
using namespace std;
typedef long long ll;
vector<int> primes;
vector<bool> is_prime;
ll gcd(ll a,ll b)
{
	ll m=1;
	if(!b)	return a;
	while(m){
		m=a%b;
		a=b;
		b=m;
	}
	return a;
}

ll multi_mod(ll a,ll b,ll mod)
{
	ll sum=0;
	while(b){
		if(b&1)	sum=(sum+a)%mod;
		a=(a+a)%mod;
		b>>=1;
	}
	return sum;
}

ll pow_mod(ll a,ll b,ll mod)
{
	ll sum=1;
	while(b)
	{
		if(b&1) sum=multi_mod(sum,a,mod);
		a=multi_mod(a,a,mod);
		b>>=1;
	}
	return sum;
} 

bool miller_rabin(ll n,int times)
{
	if(n<2) return false;
	if(n==2) return true;
	if(!(n&1)) return false;
	ll q=n-1;
	int k=0;
	while(!(q&1)){
		++k;
		q>>=1;
	}
	for(int i=0;i<times;++i){
		ll a=rand()%(n-1)+1;
		ll x=pow_mod(a,q,n);
		if(x==1) continue;
		for(int j=0;j<k;++j){
			ll buf=multi_mod(x,x,n);
			if(buf==1&&x!=1&&x!=n-1)
				return false;
			x=buf;
		}
		if(x!=1)
			return false;
	}
	return true;
}

ll pollard_rho(ll n)
{
	ll d,i,k,x,y;
	x=rand()%(n-1)+1;
	y=x;
	i=1;
	k=2;
	do
	{
		++i;
		d=gcd(n+y-x,n);
		if(d>1&&d<n)
			return d;
		if(i==k)
			y=x,k*=2;
		x=((multi_mod(x,x,n)-gcc)%n+n)%n;
	}while(y!=x);
	return n;
}

bool isPrime(ll n)
{
	if(n<=max_prime) return is_prime[n];
	return miller_rabin(n,20);
}

void factorize(ll n,map<ll,int>& factors)
{
	if(isPrime(n))
		++factors[n];
	else{
		for(int i=0;i<primes.size();++i){
			int p=primes[i];
			while(n%p==0){
				++factors[p];
				n/=p;
			}
		}
		if(n!=1){
			if(isPrime(n))
				++factors[n];
			else{
				ll d=pollard_rho(n);
				factorize(d,factors);
				factorize(n/d,factors);
			}
		}
	}
}

pair<ll,ll> solve(ll a,ll b)
{
	ll c=b/a;
	map<ll,int> factor;
	factorize(c,factor);
	vector<ll> mul_factors;
	for(map<ll,int>::iterator it=factor.begin();it!=factor.end();++it){
		ll mul=1;
		while(it->second){
			mul*=it->first;
			it->second--;
		}
		mul_factors.push_back(mul);
	}

	ll best_sum=1e18,best_x=1,best_y=c;
	for(int mask=0;mask<(1<<mul_factors.size());++mask){
		ll x=1;
		for(int i=0;i<mul_factors.size();++i)
			if((mask>>i)&1) x*=mul_factors[i];
		ll y=c/x;
		if(x<y&&x+y<best_sum){
			best_sum=x+y;
			best_x=x;
			best_y=y;
		}
	}
	return make_pair(best_x*a,best_y*a);
}

void ini_primes()
{
	is_prime=vector<bool>(max_prime+1,true);
	is_prime[0]=is_prime[1]=false;
	for(int i=2;i<=max_prime;++i){
		if(is_prime[i]){
			primes.push_back(i);
			for(int j=i*2;j<=max_prime;j+=i)
				is_prime[j]=false;
		}
	}
}

int main()
{
	ll a,b;
	ini_primes();
	while(scanf("%I64d%I64d",&a,&b)==2){
		pair<ll,ll> ans=solve(a,b);
		printf("%I64d %I64d\n",ans.first,ans.second);
	}
	return 0;
}