Miller-Rabin素数测试
2016-04-05 20:53
267 查看
#include <iostream> #include <cstdio> #include <cstring> using namespace std; typedef long long LL; const LL prime[12]={2,3,5,7,11,13,17,19,23,29,31,37};//结论:对于Long Long范围的素数只需要这些 LL mul(LL a,LL b,LL n){ //********************蒙哥马利算法,解决a*b%n的问题,将b二进制拆开“快速加”,避免了a*b溢出 LL res=0; while (b){ if (b&1) res=(res+a)%n; a=(a+a)%n; b>>=1; } return res; } LL pow(LL a,LL m,LL n){ LL res=1; while (m){ if (m&1) res=mul(res,a,n); a=mul(a,a,n); m>>=1; } return res; } bool Miller_Rabin(LL p){ if (p==2) return 1; if (p==1 || p%2==0 ) return 0; LL k=0,e=p-1,s,i,j; while (e%2==0) k++,e/=2; for (i=0;i<12;i++){ if (prime[i]==p) return 1; s=pow(prime[i],e,p); for (j=1;j<=k;j++){ if (mul(s,s,p)==1 && !(s==1 || s==p-1)) return 0;//***************只有当X^2=1(MOD P)时才要二次剩余检验 s=mul(s,s,p); } if (s!=1) return 0;//fermat检验 } return 1; } int main() { int t; LL x; scanf("%d",&t); while (t--) { scanf("%lld",&x);//****************************注意读入类型 Miller_Rabin(x) ? printf("Yes"):printf("No"); printf("\n"); } return 0; }
误判概率
经过独立的t轮Miller-Rabin算法将一个合数误判为素数的可能性不大于 (1/4)t,(正确实现的算法不会误判素数为合 数),这个误判概率在基于Fermat定理的算法中是最好的。这个误判概率是利用下面两个定理证明的:
定 理 1 :设d=gcd(k,m) ,那么在有限群{g,g2,g3,?,gm=1}中(g是有限群的生成元,m是有限群的阶) 有d个元素x满足方程式xk=1。
定 理 2 : 设p是一个奇素数,p-1=2sh(h是奇数),那么在乘法群(Z/pZ)*中满足方程式 x2rt=-1mod p(t是奇数 )的元素的个数是:0,如果r ≥s;2rgcd(h,t),如果r<s。 利用这两个定理,对算法输入n分3种情况讨论:
(1)n可以被某个素数的平方整除的时候;
(2)n是两个不同素数的乘积的时候;
(3)n是两个以上不同素数乘积的时候。
这样可以证明Miller-Rabin算法的误判概率上界。
优化实现
Miller-Rabin算法最为耗时的步骤在2.2模幂操作和2.3.2 循环。对算法的优化实现主要集中在对这两部分运算的优 化。对模幂操作的优化有两种途径:减少模幂算法中的模乘 操作和优化模乘操作。在求模幂的过程中不能先求幂最后一次求模,这样会产生一个十分巨大的中间结果,造成实际的 不可操作,所以在求模幂的算法中用模乘代替乘法,使得中 间结果的长度不超过模的长度。对模幂算法的优化,我们使 用改进的滑动窗口算法结合Montgomery模乘和模平方算法。表1给出模幂算法的比较。
模幂算法 | 预先计算 | 模平方 | 模乘法 | ||
模平方 | 模乘法 | 最坏情况 | 平均情况 | ||
平方乘算法 滑动窗口类算法 改进的滑动窗口算法 | 0 1 1 | 0 2k -3 2k-1-1 | t t-(k-1)≤次数≤t t-(k-1)≤次数≤t | t (t/k)-1 (t/k)-1 | t/2 t/k(2k-1)/ 2k k ≤t/k(2 -1)/ |
输入: x,e=(e tet-1?e1e0)2,其中et=1,k≥1( 窗口大小)
输出: xe mod n
1、预计算
1.1、x1← MontMul(x, R2,n),
x2←MontSqu(x 1, n)
1.2、对i 从1 到2k-1-1计算x2i+1←MontMul(x2i-1, x2,n)
2、A←R,i ←t
3、 当i≥ 0时作下面的操作:
3.1、如果ei=0,A←MontSqu(A ,n),i← i-1
3.2、否则找到最长的位串eiei-1?es使得i-s+1≤k并且es=1,计算
3.2.1、A <-A2i-s+1 , (利 用MontSqu函数计算)
3.2.2、A <-A*X(ee ...e )2 ,(利 用MontMul函数计算)
3.2.3、i ←s-1
4、A←MontMul(A ,1 ,n)
5、返回A
其中MontMul(x,y,n) 是Montgomery模乘函数,函数输出 结果为x*y*R-1 mod n,MontSqu(x,n) 是Montgomery模平方函 数,输出结果为x2R-1 mod n。模乘算法如果采用大整数乘法 和除法求模乘,因为涉及耗时的除法操作,所以要相对较 慢,结合大整数乘法和Barrett求模算法可以用2(n2+3n+1) 步 单精度乘法完成。使用Montgomery求模算法结合大整数乘法 算法,可以 在 2n(n+1) 步单精度乘法内完成算法。 Montgomery模平方的操作可以在3n(n+1)
/2步单精度乘法内 完成,而Barrett模平方需要(3n(n+3)/2+1) 步单精度乘法。结 合改进的滑动窗口算法和Montgomery类算法,可以得到目前 非特殊情况下的最优的模幂算法。
在Miller-Rabin算法的2.3.2循环中的模平方操作我们没有 使用Montgomery模平方算法,因为该算法给出的结果带有R-1这个参数,在2.3.2循环中处理掉这个参数将占整个循环运 行时间中的很大部分,尤其是在循环的控制参数s 相对较小的时候。我们在这里使用大整数平方算法结合Barrett求模算 法,2.3.2的循环最坏情况需要(s-1)(3n(n+3)/2+1)步单精度乘法。
相关文章推荐
- 【BZOJ-3626】LCA 树链剖分
- 2.猜数字
- QT学习笔记——“const"
- #leetcode#283. Move Zeroes
- 分类-3-生成学习-1-问题引入
- P值已经被废,学的东西有多少是真的?
- 【原】Spark中Stage的提交源码解读
- 【转】深入剖析Java中的装箱和拆箱
- Codeforces Round #345 (Div. 1) & CodeForces 650A Watchmen
- android.view.GLES20RecordingCanvas.drawBitmap异常
- 第5周项目4-长方柱类
- HDU 4123(RMQ(o1)模板 , 树直径)
- 如何阻止移动设备(手机,pad)浏览器双击放大网页?
- Fresco-Facebook的图片加载框架的使用
- BZOJ 2209: [Jsoi2011]括号序列&&BZOJ 2329: [HNOI2011]括号修复
- codevs2370小机房的树
- 使用Classie.js让添加、删除和检查类更容易
- iOS安全–在非越狱平台进行越狱开发(附分析流程)
- 《LeetBook》leetcode题解(3):Longest Substring Without Repeating Characters[M]——哈希判断重复
- 《LeetBook》leetcode题解(3):Longest Substring Without Repeating Characters[M]——哈希判断重复