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重复上述过程
通过这个方式,我们只需要知道x1和c就可以生成一系列数值,c可以取任意给定值,一般取c=1。
若序列出现循环,则退出。
计算p=gcd(xi-1-xi, 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 }