当前位置: 代码迷 >> 综合 >> poj 2429 GCD LCM Inverse miller_rabin素数判定和pollard_rho因数分解
  详细解决方案

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

热度:97   发布时间:2024-01-19 05:46:34.0

题意:

给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;	
}