C++:Miller-Rabin素数(质数)检测算法

2.1 Miller-Rabin理论基础

Fermat定理 若n是奇素数，a是任意正整数(1≤ a≤ n?1)，则 a^n?1≡1pmod n。[2]

Miller-Rabin算法的理论基础 如果n是一个奇素数，将n?1表示成2^s*r的形式，r是奇数，a与n是互素的任何整数，那么a^r≡1pmod n或者对某个j(0 ≤ j≤ s?1, j∈Z)等式a^2jr≡?1 pmod n成立。[2]

这个理论是由Fermat定理推导而来：n是一个奇素数，则方程 x²≡1pmod n只有±1两个解。

定理 3 设x，y和n是整数，如果x² = y²pmod n, 但x≠± ypmod n，则(x?y)和n的公约数中有n的非平凡因子。

2.2 Miller-Rabin算法描述

输入：一个大于3的奇整数n和一个大于等于1的安全参数t（用于确定测试轮数）。

输出：返回n是否是素数（概率意义上的，一般误判概率小于((1/2))⁸⁰)即可)。

1. 将n-1表示成2^sr

2. 对i从1到t做循环做以下操作：

   1. 选择一个随机整数a（2 ≤ a ≤ n?2）

   2. 计算y ← a^r bmod n

   3. 如果y≠1并且y≠ n?1循环做下面的操作，否则转3：

      1. j ← 1

      2. 当j ≤ s?1并且y ≠ n?1循环做下面操作，否则跳到（iv.）

      3. 计算y ← y² bmod n，如果y = 1返回“合数”，否则j ← j + 1

      4. 如果y ≠ n?1则返回“合数”

3. 返回“素数”

其中第(vi.)步是基于定理3得来的。 [2]

2.3 Miller-Rabin算法的误判概率

经过独立的t轮Miller-Rabin算法将一个合数误判为素数的可能性不大于4^?t。这个概率是给予Fermat定理的算法中是最好的。这个误判概率是利用下面的两个定理证明的。

定理 1：设d = gcd(k,m)那么在有限群{g1,g + 2,···,gm = 1}中（g是有限群的生成元，m是有限群的阶）有d个元素满足方程x^k = 1。

定理 2： 设p是一个素奇数，p?1 = 2^sh（h是奇数），那么在乘法群(Z/pZ)^*中满足方程 x^2rt = ?1pmod p（t是奇数）的元素个数是：0，如果r≥ s；2^rgcd(h,t)如果r<s。

利用这两个定理，对算法输入n分3种情况讨论：

1. n是可以被某个素数的平方整除的时候；

2. n是两个不同的素数的乘积的时候；

3. n是两个以上不同素数乘积的时候。

这样就可以证明Miiler-Rabin算法的误判概率上界。

2.4 Miller-Rabin算法的时间复杂度

Miller-Rabin算法是一个概率算法，算法的计算集中于(b)步和(c)步的循环中，最坏情况是(iv.)的循环没有中途推出，则一轮Miller-Rabin算法的最坏情况复杂度为(1 + O(1))log2(n)(以模n乘法为基本操作)。如果以单精度乘法操作作为时间复杂度的衡量，则一轮优化的Miller-Rabin算法的最坏情况时间复杂度是O(log2³(n))。从时间复杂度来看Miller-Rabin算法的性能是很好的。在实际应用中，Miller-Rabin算法的实际执行速度也很快。

代码实现：

点击(此处)折叠或打开

1. #include <iostream>
   
2. 
   
3. using namespace std;
   
4. 
   
5. typedef unsigned __int64 llong;
   
6. 
   
7. llong mod_pro(llong x,llong y,llong n)
   
8. {
   
9. llong ret=0,tmp=x%n;
   
10. while(y)
    
11. {
    
12. if(y&0x1)if((ret+=tmp)>n)ret-=n;
    
13. if((tmp<<=1)>n)tmp-=n;
    
14. y>>=1;
    
15. }
    
16. return ret;
    
17. }
    
18. llong mod(llong a,llong b,llong c)
    
19. {
    
20. llong ret=1;
    
21. while(b)
    
22. {
    
23. if(b&0x1)ret=mod_pro(ret,a,c);
    
24. a=mod_pro(a,a,c);
    
25. b>>=1;
    
26. }
    
27. return ret;
    
28. }
    
29. llong ran()
    
30. {
    
31. llong ret=rand();
    
32. return ret*rand();
    
33. }
    
34. bool is_prime(llong n,int t)
    
35. {
    
36. if(n<2)return false;
    
37. if(n==2)return true;
    
38. if(!(n&0x1))return false;
    
39. llong k=0,m,a,i;
    
40. for(m=n-1;!(m&1);m>>=1,k++);
    
41. while(t--)
    
42. {
    
43. a=mod(ran()%(n-2)+2,m,n);
    
44. if(a!=1)
    
45. {
    
46. for(i=0;i<k&&a!=n-1;i++)
    
47. a=mod_pro(a,a,n);
    
48. 
    
49. if(i>=k)return false;
    
50. }
    
51. }
    
52. return true;
    
53. }
    
54. int main()
    
55. {
    
56. llong n;
    
57. while(scanf("%I64u",&n)!=EOF)
    
58. if(is_prime(n,3))
    
59. cout<<"YES\n";
    
60. else
    
61. cout<<"NO\n";
    
62. return 0;
    
63. }