Pollard_Rho 整数分解法

如果要对比较大的整数分解，显然之前所学的筛选法和是试除法都将不再适用。所以我们需要学习速度更快的Pollard_Rho算法

pollard_rho 算法流程

Pollard_rho算法的大致流程是先判断当前数是否是素数（Miller_rabin）了，如果是则直接返回。如果不是素数的话，试图找到当前数的一个因子（可以不是质因子）。然后递归对该因子和约去这个因子的另一个因子进行分解。

那么自然的疑问就是，怎么找到当前数n的一个因子？当然不是一个一个慢慢试验，而是一种神奇的想法。其实这个找因子的过程我理解的不是非常透彻，感觉还是有一点儿试的意味，但不是盲目的枚举，而是一种随机化算法。我们假设要找的因子为p，他是随机取一个x1，由x1构造x2，使得｛p可以整除x1-x2 && x1-x2不能整除n｝则p=gcd（x1-x2，n），结果可能是1也可能不是1。如果不是1就找寻成功了一个因子，返回因子；如果是1就寻找失败，那么我们就要不断调整x2，具体的办法通常是x2=x2*x2+c（c是自己定的）直到出现x2出现了循环==x1了表示x1选取失败重新选取x1重复上述过程

通过这个方式，我们只需要知道x₁和c就可以生成一系列数值，c可以取任意给定值，一般取c=1。

若序列出现循环，则退出。

计算p=gcd(x_i-1-x_i, n)，若p=1，返回上一步，直到p>1为止；若p=n，则n为素数，否则p为一个约数并递归分解p和n/p。

  1 #include <iostream>
  2 #include <time.h>
  3 #include <algorithm>
  4 #include <stdio.h>
  5 
  6 typedef long long LL;
  7 
  8 using namespace std;
  9 
 10 const int times = 20;
 11 LL fac[1001];
 12 int cnt;
 13 
 14 LL mul(LL a,LL b,LL mod){
 15     LL ans = 0;
 16     while (b){
 17         if (b & 1){
 18             ans = (ans + a) % mod;
 19         }
 20         a = (a<<1) % mod;
 21         b >>= 1;
 22     }
 23     return ans;
 24 }
 25 
 26 
 27 LL pow(LL a,LL b,LL mod){
 28     LL ans = 1;
 29     while (b){
 30         if (b & 1){
 31             ans = mul(ans,a,mod);
 32         }
 33         b >>= 1;
 34         a = mul(a,a,mod);
 35     }
 36     return ans;
 37 }
 38 
 39 
 40 bool witness(LL a,LL n){
 41     LL temp = n - 1;
 42     int j = 0;
 43     while (temp % 2 == 0){  // 其实就是得到 m
 44         j++;
 45         temp /= 2;
 46     }
 47     LL x = pow(a,temp,n);
 48     if (x == 1 || x == n-1){   // 判断a^m
 49         return true;
 50     }
 51     while (j--){
 52         x = mul(x,x,n);  // 进一步判断 a^(2m)  a^(4m) ...
 53         if (x == n-1)
 54             return true;
 55     }
 56     return false;
 57 }
 58 
 59 bool miller_rabin(LL n){
 60     if (n == 2){ //  如果是2肯定是素数
 61         return true;
 62     }
 63     if (n<2 || n % 2 == 0){ //如果小于2或者是大于2的偶数肯定不是素数
 64         return false;
 65     }
 66     for (int i=0;i<times;i++){ //随机化检验
 67         LL a = rand() % (n-1) + 1;
 68         if (!witness(a,n))
 69             return false;
 70     }
 71     return true;
 72 }
 73 
 74 LL gcd(LL a,LL b){ // 这里的gcd和一般的gcd不一样
 75     if (a == 0){  // pollard_rho的需要
 76         return 1;
 77     }
 78     if (a < 0){ // 可能有负数
 79         return gcd(-a,b);
 80     }
 81     while (b){
 82         LL t = a % b;
 83         a = b;
 84         b = t;
 85     }
 86     return a;
 87 }
 88 
 89 LL pollard_rho(LL n,LL c){ // 找因子
 90     LL i = 1,k = 2; // 用来判断是否成环
 91     LL xx = rand() % n,y = xx;
 92     while (1){
 93         i++;
 94         xx = (mul(xx,xx,n) + c) % n;
 95         LL d = gcd(y-xx,n);
 96         if (1<d && d<n){ // 找到一个因数
 97             return d;
 98         }
 99         if (y == xx){ // 出现循环，那么查找失败
100             return n;
101         }
102         if (i == k){ // 相当一个优化？
103             y = xx;
104             k <<= 1;
105         }
106     }
107 }
108 
109 void find(LL n){ // 通过因数来找质因子
110     if (miller_rabin(n)){
111         fac[cnt++] = n; // 记录质因子
112         return ;
113     }
114     LL p = n;
115     while (p >= n)
116         p = pollard_rho(p,rand() % (n-1) + 1); //如果转了一圈还是p那么继续
117     find(p);
118     find(n/p);
119 }
120 
121 int main(){
122     srand(time(NULL));
123     int t;
124     scanf("%d",&t);
125     while (t--){
126         LL x;
127         scanf("%lld",&x);
128         if (miller_rabin(x)){
129             printf("Prime\n");
130             continue;
131         }
132         cnt = 0;
133         find(x);
134         LL ans = fac[0];
135         for (int i=1;i<cnt;i++){
136             if (fac[i]<ans)
137                 ans = fac[i];
138         }
139         printf("%lld\n",ans);
140     }
141     return 0;
142 }