整数的素因子分解：Pollard rho method

参考：

　　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位左右的大数分解。

　　文章代码已经改过，对之前代码错误的问题表示歉意！！！

代码：

  1 //zzy2012.7.14
  2 #include<cstdio>
  3 #include<iostream>
  4 #include<cstdlib>
  5 #define ll long long
  6 using namespace std;
  7 
  8 long long factor[1000],fnum;
  9 
 10 ll mul_mod(const ll &a, ll b, const ll &n)
 11 {
 12     ll back(0), temp(a % n);
 13     b %= n;
 14     while ( b > 0 )
 15     {
 16         if ( b & 0x1 )
 17         {
 18             if ( (back = back + temp) > n )
 19             back -= n;
 20         }
 21         if ( (temp <<= 1) > n )
 22         temp -= n;
 23         b >>= 1;
 24     }
 25     return back;
 26 }
 27 
 28 ll pow_mod(const ll &a, ll b, const ll &n)
 29 {
 30     ll d(1), dTemp(a % n);//当前二进制位表示的是进制数值
 31     while ( b > 0 )
 32     {
 33         if ( b & 0x1 ){
 34             d = mul_mod(d, dTemp, n);
 35             b ^= 1;
 36         }
 37         else{
 38             dTemp = mul_mod(dTemp, dTemp, n);
 39             b >>= 1;
 40         }
 41     }
 42     return d;
 43 }
 44 
 45 bool MillerRabin(long long p,long long a){
 46     long long k=0,q=p-1,m;
 47     while(q%2 == 0){
 48         k++;
 49         q/=2;
 50     }
 51     m = pow_mod(a,q,p);
 52 
 53     if(m==1 || m==p-1)
 54         return true;
 55     for(long long i=1; i<k; i++){
 56         m = pow_mod(m,2,p);
 57         if(m == p-1)
 58             return true;
 59     }
 60     return false;
 61 }
 62 
 63 bool isPrime(long long n){
 64     for(long long i=1; i<=20; i++){
 65         if(MillerRabin(n, rand()%(n-1)+1)==false)
 66             return false;
 67     }
 68     return true;
 69 }
 70 
 71 long long gcd(long long a,long long b){
 72     if(b==0LL)
 73         return a;
 74     return gcd(b,a%b);
 75 }
 76 
 77 
 78 ll Pollard_rho(const ll &c,const ll &num)
 79 {
 80     int i=1, k=2;
 81     ll x = rand() % num;
 82     ll y = x, comDiv;
 83     do
 84     {
 85         ++i;
 86         if ( (x = mul_mod(x, x, num) - c) < 0 )
 87         x += num;
 88         if ( x == y )
 89         break;
 90         comDiv = gcd((y-x+num)%num, num);
 91         if ( comDiv > 1 && comDiv < num )
 92         return comDiv;
 93         if ( i == k )
 94         {
 95             y = x;
 96             k <<= 1;
 97         }
 98     }while ( true );
 99     return num;
100 }
101 void FindFac(const ll &num)
102 {
103     if (isPrime(num)==true)
104     {
105         factor[fnum++] = num;
106         return;
107     }
108      ll fac;
109      do
110      {
111             fac=Pollard_rho(rand()%(num-1)+1,num);
112      }while (fac>=num);
113      FindFac(fac);
114      FindFac(num/fac);
115 }
116 
117 int main()
118 {
119     long long n;
120     cin>>n;
121     if(isPrime(n) == true)
122         cout<<n<<endl;
123     else{
124         fnum = 0;
125         FindFac(n);
126         for(long long i=0 ;i<fnum; i++)
127             cout<<factor[i]<<' ';
128         cout<<endl;
129     }
130     return 0;
131 }