转自 http://www.cnblogs.com/luna-lovegood/archive/2012/07/15/2591982.html
参考:
1.CLRS《算法导论》
2.http://www.csh.rit.edu/~pat/math/quickies/rho/#algorithm
Pollard rho方法是随机算法,《算法导论》上说是启发式方法。期望的时间为n的四次方根。关于rho的实现细节《算法导论》和参考2的文章有些不同。《算法导论》计算出x1,x2,x3,...,xk ... 序列,每个数xi都要参与验证,即计算gcd( xj - xi , n) ,for all i,其中xj是下标j满足小于等于i,且j为2的幂次的最大j对应的xj。参考2文章里也是要一次计算出x1,x2,x3,...,xk ... 序列,不过只对下标为偶数的i验证(下标从1开始),即计算gcd( x[i/2] - xi, n )。
Pollard rho方法实现需要考究的地方有两个:
一是判断出现循环,这需要记录算出的xi序列,并要对新计算的xk+1和前面k个值比对,这比较耗空间和时间,所以我采取了另一种方法,设定循环次数限制,但这会导致时间效率和有效性下降,好处是实现简单。
二是初始时种子的选取,我固定选取2为种子。种子可以选取多个,即在用某个种子无法分解时,用不同种子再试,实现会复杂一些,不过程序有效性可能有所提高。另外我在实现mul_mod函数时采用高精度,可以对17到18位左右的大数分解。
代码:
//zzy2012.7.14
#include<cstdio>
#include<iostream>
#include<cstdlib>
#define ll long long
using namespace std;
long long factor[1000],fnum;
ll mul_mod(const ll &a, ll b, const ll &n)
{
ll back(0), temp(a % n);
b %= n;
while ( b > 0 )
{
if ( b & 0x1 )
{
if ( (back = back + temp) > n )
back -= n;
}
if ( (temp <<= 1) > n )
temp -= n;
b >>= 1;
}
return back;
}
ll pow_mod(const ll &a, ll b, const ll &n)
{
ll d(1), dTemp(a % n);//当前二进制位表示的是进制数值
while ( b > 0 )
{
if ( b & 0x1 ){
d = mul_mod(d, dTemp, n);
b ^= 1;
}
else{
dTemp = mul_mod(dTemp, dTemp, n);
b >>= 1;
}
}
return d;
}
bool MillerRabin(long long p,long long a){
long long k=0,q=p-1,m;
while(q%2 == 0){
k++;
q/=2;
}
m = pow_mod(a,q,p);
if(m==1 || m==p-1)
return true;
for(long long i=1; i<k; i++){
m = pow_mod(m,2,p);
if(m == p-1)
return true;
}
return false;
}
bool isPrime(long long n){
for(long long i=1; i<=20; i++){
if(MillerRabin(n, rand()%(n-1)+1)==false)
return false;
}
return true;
}
long long gcd(long long a,long long b){
if(b==0LL)
return a;
return gcd(b,a%b);
}
ll Pollard_rho(const ll &c,const ll &num)
{
int i=1, k=2;
ll x = rand() % num;
ll y = x, comDiv;
do
{
++i;
if ( (x = mul_mod(x, x, num) - c) < 0 )
x += num;
if ( x == y )
break;
comDiv = gcd((y-x+num)%num, num);
if ( comDiv > 1 && comDiv < num )
return comDiv;
if ( i == k )
{
y = x;
k <<= 1;
}
}while ( true );
return num;
}
void FindFac(const ll &num)
{
if (isPrime(num)==true)
{
factor[fnum++] = num;
return;
}
ll fac;
do
{
fac=Pollard_rho(rand()%(num-1)+1,num);
}while (fac>=num);
FindFac(fac);
FindFac(num/fac);
}
int main()
{
long long n;
cin>>n;
if(isPrime(n) == true)
cout<<n<<endl;
else{
fnum = 0;
FindFac(n);
for(long long i=0 ;i<fnum; i++)
cout<<factor[i]<<' ';
cout<<endl;
}
return 0;
}