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 200000using 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;}