MILLER_RABIN素数判定法

 1 //MILLER_RABINE素数判别法
 2 #include<iostream>
 3 using namespace std;
 4 
 5 #include<vector>
 6 #include<cstdlib>
 7 #include<cmath>
 8 
 9 int MODULAR_EXPONENTIATION(int a,int b,int n)//a^b mod n
10 {
11     int c = 0,d = 1;
12     vector<int> vb;
13     while(b)
14     {
15         vb.push_back(b%2);
16         b /= 2;
17     }
18     int k = vb.size()-1;
19     for(int i = k;i >= 0;--i)
20     {
21         c *= 2;
22         d = (d*d)%n;
23         if(vb[i] == 1)
24         {
25             ++c;
26             d = (d*a)%n;
27         }
28     }
29     vb.clear();
30     return d;
31 }
32 
33 bool WITNESS(int a,int n)
34 {
35     int t = 0,u,k = n;
36     --k;
37     while(k%2 == 0)
38     {
39         k /= 2;
40         ++t;
41     }
42     u = k;
43     int *x = new int[t+1];
44     x[0] = MODULAR_EXPONENTIATION(a,u,n);
45     for(int i = 1;i <= t;++i)
46     {
47         x[i] = x[i-1]*x[i-1]%n;
48         if(x[i] == 1 && x[i-1] != 1 && x[i-1] != n-1)
49             return true;
50     }
51     if(x[t] != 1) return true;
52     return false;
53 }
54 
55 bool MILLER_RABIN(int n,int s)
56 {
57     int a;
58     srand(0);
59     for(int j = 1;j <= s;++j)
60     {
61         a = rand();
62         a %= n;
63         if(a < n-1) ++a;
64         if(WITNESS(a,n)) return false;
65     }
66     return true;
67 }
68 
69 int main()
70 {
71     vector<int> vi;
72     vi.push_back(2);
73     for(int i = 3;i <= 32768;i += 2)
74     {
75         if(MILLER_RABIN(i,10)) vi.push_back(i);
76     }
77 
78     return 0;
79 }