算法思想
基本思路:
对于给定的一个整数 n,显然如果 n 为素数(Miller-Rabin 算法判断),那么算法结束,返回唯一素因子 n。
否则,pollard's rho 算 ***试着找到 n 的一个因子 a(a 并不一定是素数),然后递归 Pollard_rho( a ) 和 Pollard_rho( n/a ) ,即将 n 分为了两个部分,将原问题变成了规模更小的两个子问题。
如何求因子 a?
一般情况下可能会本能的想到枚举,但这样过于耗时。寻找因子是该算法中最重要的一步,该算法中采用随机化算法来查找因子 a。假设此时 n 仅有 2 个因子 p 和 q,那么如果我们用随机数来找 n 的因子,成功率为 1/n。接下来我们思路是如何提升成功率。
寻找一个因子 a,等价于寻找是否存在 k 个数,使得其中有 ,由于生日悖论 我们可以得知当 时,该概率是 50%,提高了成功率。但不幸的是对于 10 位的整数,k = 1e5 时,仍要做 次比较,我们仍然选取 k 个数: 但我们不再询问是否存在 可以整除n,转而询问是否有 的情况。因为如果要问有多少个数可以整除n,那么只有p,q。但是如果问有多少个数满足,就会有更多的数的组合满足了,从而也增加了成功率。而 就是一个n的因子。但是经过验证,还是需要选择 才能有相当高的成功率。而 数量仍然很大,以至于可能不方便存在内存中。为了解决无法储存太多数的问题,Pollard's rho Algorithm 只将两个数存放在内存中,具体思路是:并不随机生成 k 个数并两两比较,而是一个一个地生成并检查连续的两个数。反复执行这个步骤并希望能够得到我们想要的数。
我们使用一个函数来生成伪随机数。
使用的时候从x1=2开始,让x2 = f(X1), x3 = f(x2) ,xn+1=f(xn),于是顺着这个策略,不断判断 是否大于1即可。
利用这个函数因为x1和x2再调整时最终一定会出现循环,形成一个类似希腊字母rho的形状,故因此得名。于是问题变成了“如何判断环的出现”。一种方法是记录所有出现过的数,当然这会耗费大量内存,故舍弃;另一种方法是 Floyd 发明的算法,这里可以举一个有趣的例子说明“假设 A 和 B 在一个很长的圆形轨道上走,那么我们如何判断 B 是否走完一圈呢?我们可以让 B 的速度是 A 的二倍,他们同时出发,当 B 第一次追上 A,就知道 B 至少已经走了一圈”,同样的道理运用到该算法中,框架如下:
a = 2; b = 2; while(b != a){ a = f(a); //一倍速 b = f(f(b));//二倍速 p = GCD(|b-a| , n); if(p > 1) return "Found factor: p"; } return "Failed";
算法流程:
#include<bits/stdc++.h> using namespace std; typedef long long ll; //**************************************************************** // Miller_Rabin 算法进行素数测试 //速度快,而且可以判断 <2^63的数 //**************************************************************** ll prime[10] = {2,3,5,7,11,13,17,19,23,29}; ll Quick_Multiply(ll a, ll b, ll c) { ll ans = 0, res = a; while(b) { if(b & 1) ans = (ans + res) % c; res = (res + res) % c; b >>= 1; } return ans; } ll Quick_Power(ll a, ll b, ll m) { ll ans = 1; while(b) { if(b & 1) { ans = ans * a % m; } a = a * a % m; b >>= 1; } return ans; } bool Miller_Rabin(ll x) //判断素数 { ll i, j, k; ll s = 0, t = x-1; if(x == 2) return true; //2是素数 if(x < 2 || !(x & 1)) return false; //如果x是偶数或者是0,1,那它不是素数 while(!(t & 1)) //将x-1分解成(2^s)*t的样子 { s++; t >>= 1; } for(i = 0; i < 10 && prime[i] < x; i++) //随便选一个素数进行测试 { ll a = prime[i]; ll b = Quick_Power(a, t, x); //先算出a^t for(j = 1; j <= s; j++) //然后进行s次平方 { k = Quick_Multiply(b,b,x); //求b的平方 if(k == 1 && b != 1 && b != x-1) return false; b = k; } if(b != 1) return false; //用费马小定理判断 } return true; //如果进行多次测试都是对的,那么x就很有可能是素数 } //************************************************ //pollard_rho 算法进行质因数分解 //************************************************ ll factor[110]; //用来存放被分解的因子(这里也可以用map来存) int tot = 0; //因子个数 ll gcd(ll a, ll b) { /* 返回a和b的最大公约数 */ if(a == 0) return 1; if(a < 0) return gcd(-a, b); while(b) { ll t = a % b; a = b; b = t; } return a; } ll pollard_rho(ll x, ll c) { /* 返回 x 的一个因子或 x 本身 */ ll i = 1, k = 2; ll tx = rand()%x; ll y = tx; while(true) { i++; tx = (Quick_Multiply(tx, tx, x) + c) % x; ll d = gcd(y-tx, x); if(d != 1 && d != x) return d; //x是个合数,所以因子应该大于1,小于n if(y == tx) return x; //如果循环一圈了,仍未找到因子,返回到findFac函数重新取c if(i == k) y = tx, k += k; //保持y是1倍速,tx是二倍速 } } void findFac(ll n) { /* 对 n 进行质因数分解 */ if(Miller_Rabin(n)) { factor[++tot] = n; return; } ll p = n; while(p >= n) p = pollard_rho(p, rand()%(n-1)+1); //rand()%(n-1)+1 可以保证c值范围是在1~n-1之间 findFac(p); //递归分解 findFac(n/p); } int main() { ll x, T; scanf("%lld",&T); while(T--) { memset(factor, 0, sizeof(factor)); tot = 0; scanf("%lld",&x); tot = 0; if(Miller_Rabin(x)) { printf("isprime\n"); printf("%lld\n",x); } else { printf("noprime\n"); findFac(x); sort(factor+1, factor+tot+1); for(int i = 1; i < tot; i++) { if(factor[i] != factor[i+1]) printf("%lld ",factor[i]); } printf("%lld\n",factor[tot]); } } }
参考博客:
https://hacpai.com/article/1569220928268
https://www.cnblogs.com/fzl194/p/9047710.html