Miller-Rabin素数测试

 1 /*
 2 long long 范围的Miller-Rabin素数测试
 3 */
 4 
 5 #include <iostream>
 6 #include <stdio.h>
 7 #include <string.h>
 8 #include <stdlib.h>
 9 #include <time.h>
10 #define ll long long
11 using namespace std;
12 
13 //计算a*b%mod
14 ll mod_mul(ll a, ll b, ll n){
15     ll cnt = 0LL;
16     while(b){
17         if(b&1LL) cnt = (cnt+a)%n;
18         a=(a+a)%n;
19         b >>= 1LL;
20     }
21     return cnt;
22 }
23 //计算a^b%mod
24 ll mod_exp(ll a, ll b, ll n){
25     ll res = 1LL;
26     while(b){
27         if(b&1LL) res = mod_mul(res,a,n);
28         a = mod_mul(a,a,n);
29         b >>= 1LL;
30     }
31     return res;
32 }
33 //Miller-Rabin测试，测试n是否为素数
34 bool Miller_Rabin(ll n, int respat){
35     if(n==2LL || n==3LL || n==5LL || n==7LL || n==11LL) return true;
36     if(n==1 || !(n%2) || !(n%3) || !(n%5) || !(n%7) || !(n%11)) return false;
37     
38     int k = 0;
39     ll d = n-1; //要求x^u%n 不为1一定不是素数
40     
41     while(!(d&1LL)){
42         k++; d >>= 1LL;
43     }
44     srand((ll)time(0));
45     for(int i = 0; i < respat; i ++) {
46         ll a = rand()%(n-2)+2;     //在[2,n)中取整数
47         ll x = mod_exp(a,d,n);
48         ll y = 0LL;
49         for(int j = 0; j < k; j ++){
50             y = mod_mul(x,x,n);
51             if(1LL==y && 1LL!=x && n-1LL!=x)return false; //二次探测定理，这里如果y = 1则x必须等于 1
52                                                           //或n-1,否则可以判断不是素数
53             x = y;
54         }
55         if(1LL != y) return false;     //费马小定理
56     }
57     return true;
58 }
59 int main(){
60     int n;
61     cin>>n;
62     if(Miller_Rabin(n,5))cout << "Yes\n";
63     else cout << "No\n";
64     return 0;
65 }

  1 #include <cstdio>
  2 #include <cstring>
  3 #include <cstdlib>
  4 
  5 #define MAXL 4
  6 #define M10 1000000000
  7 #define Z10 9
  8 
  9 const int zero[MAXL - 1] = {0};
 10 
 11 struct bnum{
 12     int data[MAXL]; //  断成每截9个长度
 13     
 14     //  读取字符串并转存
 15     void read(){
 16         memset(data, 0, sizeof(data));
 17         char buf[32];
 18         scanf("%s", buf);
 19         int len = (int)strlen(buf);
 20         int i = 0, k;
 21         while (len >= Z10){
 22             for (k = len - Z10; k < len; ++k){
 23                 data[i] = data[i] * 10 + buf[k] - ‘0‘;
 24             }
 25             ++i;
 26             len -= Z10;
 27         }
 28         if (len > 0){
 29             for (k = 0; k < len; ++k){
 30                 data[i] = data[i] * 10 + buf[k] - ‘0‘;
 31             }
 32         }
 33     }
 34     
 35     bool operator == (const bnum &x) {
 36         return memcmp(data, x.data, sizeof(data)) == 0;
 37     }
 38     
 39     bnum & operator = (const int x){
 40         memset(data, 0, sizeof(data));
 41         data[0] = x;
 42         return *this;
 43     }
 44     
 45     bnum operator + (const bnum &x){
 46         int i, carry = 0;
 47         bnum ans;
 48         for (i = 0; i < MAXL; ++i){
 49             ans.data[i] = data[i] + x.data[i] + carry;
 50             carry = ans.data[i] / M10;
 51             ans.data[i] %= M10;
 52         }
 53         return  ans;
 54     }
 55     
 56     bnum operator - (const bnum &x){
 57         int i, carry = 0;
 58         bnum ans;
 59         for (i = 0; i < MAXL; ++i){
 60             ans.data[i] = data[i] - x.data[i] - carry;
 61             if (ans.data[i] < 0){
 62                 ans.data[i] += M10;
 63                 carry = 1;
 64             }
 65             else{
 66                 carry = 0;
 67             }
 68         }
 69         return ans;
 70     }
 71     
 72     //  assume *this < x * 2
 73     bnum operator % (const bnum &x){
 74         int i;
 75         for (i = MAXL - 1; i >= 0; --i){
 76             if (data[i] < x.data[i]){
 77                 return *this;
 78             }
 79             else if (data[i] > x.data[i]){
 80                 break;
 81             }
 82         }
 83         return ((*this) - x);
 84     }
 85     
 86     bnum & div2(){
 87         int  i, carry = 0, tmp;
 88         for (i = MAXL - 1; i >= 0; --i){
 89             tmp = data[i] & 1;
 90             data[i] = (data[i] + carry) >> 1;
 91             carry = tmp * M10;
 92         }
 93         return *this;
 94     }
 95     
 96     bool is_odd(){
 97         return (data[0] & 1) == 1;
 98     }
 99     
100     bool is_zero(){
101         for (int i = 0; i < MAXL; ++i){
102             if (data[i]){
103                 return false;
104             }
105         }
106         return true;
107     }
108 };
109 
110 void mulmod(bnum &a0, bnum &b0, bnum &p, bnum &ans){
111     bnum tmp = a0, b = b0;
112     ans = 0;
113     while (!b.is_zero()){
114         if (b.is_odd()){
115             ans = (ans + tmp) % p;
116         }
117         tmp = (tmp + tmp) % p;
118         b.div2();
119     }
120 }
121 
122 void powmod(bnum &a0, bnum &b0, bnum &p, bnum &ans){
123     bnum tmp = a0, b = b0;
124     ans = 1;
125     while (!b.is_zero()){
126         if (b.is_odd()){
127             mulmod(ans, tmp, p, ans);
128         }
129         mulmod(tmp, tmp, p, tmp);
130         b.div2();
131     }
132 }
133 
134 bool MillerRabinTest(bnum &p, int iter){
135     int  i, small = 0, j, d = 0;
136     for (i = 1; i < MAXL; ++i){
137         if (p.data[i]){
138             break;
139         }
140     }
141     if (i == MAXL){
142         // small integer test
143         if (p.data[0] < 2){
144             return  false;
145         }
146         if (p.data[0] == 2){
147             return  true;
148         }
149         small = 1;
150     }
151     if (!p.is_odd()){
152         return false;   //  even number
153     }
154     bnum a, s, m, one, pd1;
155     one = 1;
156     s = pd1 = p - one;
157     while (!s.is_odd()){
158         s.div2();
159         ++d;
160     }
161     
162     for (i = 0; i < iter; ++i){
163         a = rand();
164         if (small){
165             a.data[0] = a.data[0] % (p.data[0] - 1) + 1;
166         }
167         else{
168             a.data[1] = a.data[0] / M10;
169             a.data[0] %= M10;
170         }
171         if (a == one){
172             continue;
173         }
174         
175         powmod(a, s, p, m);
176         
177         for (j = 0; j < d && !(m == one) && !(m == pd1); ++j){
178             mulmod(m, m, p, m);
179         }
180         if (!(m == pd1) && j > 0){
181             return false;
182         }
183     }
184     return true;
185 }
186 
187 int main(){
188     bnum x;
189     
190     x.read();
191     puts(MillerRabinTest(x, 5) ? "Yes" : "No");
192     
193     return 0;
194 }