正文

Miller_Rabbin算法判断大素数

IP云
IP云V管理员/6 阅读
09/10
提示：扫一扫查出行
【扫一扫了解最新限行尾号】
复制提示

普通的素数测试我们有O(√ n)的试除算法。事实上，我们有O(s*log³n)的算法。

下面就介绍一下Miller_Rabbin算法思想：

定理一：假如p是质数，且(a,p)=1，那么a^(p-1)≡1(mod p)。即假如p是质数，且a,p互质，那么a的(p-1)次方除以p的余数恒等于1。（费马小定理）

定理二：如果p是一个素数，那么对于x(0<x<p)，若x^2 mod p 等于1，则x=1或p-1。

它利用了费马小定理，即：如果p是质数，且a，p互质，那么a^(p-1) mod p恒等于1。也就是对于所有小于p的正整数a来说都应该复合a^(p-1) mod p恒等于1。那么根据逆否命题，对于一个p，我们只要举出一个a（a<p）不符合这个恒等式（就是a^(p-1) mod p恒等于1式子），则可判定p不是素数。Miller-rabin算法就是多次用不同的a来尝试p是否为素数。

但是每次尝试过程中还做了一个优化操作，以提高用少量的a检测出p不是素数的概率。这个优化叫做二次探测。它是根据一个定理：如果p是一个素数，那么对于x(0<x<p)，若x^2 mod p 等于1，则x=1或p-1。逆否命题：如果对于x(0<x<p)，若x^2 mod p 不等于1，则p不是素数。根据这个定理，我们要计算a^(p-1) mod p是否等于1时，可以这样计算，设p-1=(2^t) * k。我们从a^k开始，不断将其平方直到得到a^(p-1)，一旦发现某次平方后mod p等于1了，那么说明符合了二次探测定理的逆否命题使用条件，立即检查x是否等于1或p-1，如果不是则可直接判定p为合数。

一些问题：

1、对于p-1=(2^t) * k这个t的找法，可以记录一下x（x表示(p-1)）的二进制形式下末尾有多少0，比如二进制(1000)B代表的是十进制(8)D,8的二进制数形式下末尾0有3个，那么也就是8=(2^3)*1;再例如二进制(11000)B代表的是十进制(24)D,24的二进制数形式下末尾0有3个，那么也就是24=(2^3)*3

2、对于随机数a的生成，因为我们这个算法判断的是p这个数是不是素数，那么素数和任意一个数都互质，题目又有要求a<p所以就随机生成一个[1,p-1]范围内的数就行

3、这个算法的复杂度是O(s*log³n)，这个s是我们提前定好的，因为上面说了我们用的是费马小定理的逆否命题，但是该定理的逆命题是不一定成立的，但是令人可喜的是大多数情况是成立的。所以我们可以多找几个数a来验证一下p是否是素数，找几个数，这个多少数就是s

可以证明，使用以上两个定理以后，检验s次出错的概率至多为2^(-s)，所以这个算法是很可靠的。

1 #include<stdio.h>
2 #include<string.h>
3 #include<iostream>
4 #include<algorithm>
5 #include<queue>
6 #include<map>
7 #include<vector>
8 #include<math.h>
9 #define mem(a,x) memset(a,x,sizeof(a))
10 using namespace std;
11 typedef long long LL;
12 const int maxn=50005;
13 const int mod=26;
14 const int INF=0x3f3f3f3f;
15 const int Times = 10;
16 const int N = 5500;
17 LL ct, cnt;
18 LL fac[N], num[N];
19 LL gcd(LL a, LL b) //求两数最大公因子
20 {
21 return b? gcd(b, a % b) : a;
22 }
23 LL multi(LL a, LL b, LL m) //快速乘
24 {
25 LL ans = 0;
26 a %= m;
27 while(b)
28 {
29 if(b & 1)
30 {
31 ans = (ans + a) % m;
32 b--;
33 }
34 b >>= 1;
35 a = (a + a) % m;
36 }
37 return ans;
38 }
39 LL pow(LL a, LL b, LL m) //快速幂
40 {
41 LL ans = 1;
42 a %= m;
43 while(b)
44 {
45 if(b & 1)
46 {
47 ans = multi(ans, a, m);
48 b--;
49 }
50 b >>= 1;
51 a = multi(a, a, m);
52 }
53 return ans;
54 }
55 bool Miller_Rabin(LL n) //判断是不是素数
56 {
57 if(n == 2) return true;
58 if(n < 2 || !(n & 1)) return false;
59 LL m = n - 1;
60 int k = 0;
61 while((m & 1) == 0)
62 {
63 k++; //这个k就是我们讲的时候的t
64 m >>= 1; //这个m就是k
65 }
66 for(int i=0; i<Times; i++) //Times就是我们事先定义的s（要找a的个数）
67 {
68 LL a = rand() % (n - 1) + 1; //找一个[1,n-1]内的任意数
69 LL x = pow(a, m, n);
70 LL y = 0;
71 for(int j=0; j<k; j++)
72 {
73 y = multi(x, x, n);
74 if(y == 1 && x != 1 && x != n - 1) return false;
75 x = y;
76 }
77 if(y != 1) return false;
78 }
79 return true;
80 }
81 int main()
82 {
83 LL n;
84 while(cin>>n)
85 {
86 if(Miller_Rabin(n))
87 printf("是素数\n");
88 else printf("不是素数\n");
89 }
90 return 0;
91 }