MillerRabin 素性测试

拉宾-米勒素性测试（Miller-Rabin 质数测试）（参考Silverman《数论概论》和 Wikipedia）

转载请注明出处：http://www.cnblogs.com/luna-lovegood/archive/2012/07/13/2590925.html

　　Miller-Rabin素性测试是一个随机算法，它能以很高的概率指出某个数可能是素数，或者判定某个数一定不是素数。Miller-Rabin判定为合数是确定的，但是判决出是素数则是以大概率的。

　　Miller-Rabin测试的正确性或者叫有效性依赖于以下关于素数的一个性质：

　　设p是一个奇素数，a是与p互素的整数，p-1 可以写为 (2k )*q的形式，其中q是奇数。则下列两种情况之一发生：

　　　　1. aq = 1 (mod p)

　　　　2. a(2^i)*q = p-1 (mod p)

　　证明：根据费马小定理有 ap-1 = 1 (mod p) , 即 a2^k *q = 1 (mod p) 。

　　考虑序列：aq , a2*q, a2^2*q, ... , a2^k-1*q, a2^k *q,每一项都是前一项的平方。

　　如果第一项aq mod p为1，则能满足费马小定理。因为此时aq = n*p + 1, 那么a2*q = (n*p + 1)^2 = (n*n*p + 2*n )*p + 1 ,可以推出a2*q = 1 (mod p)，这个关系可以一直推至a2^k*q。

　　如果第一项不为1 (mod p)，在最后一项或最后一项之前某一项要为1 (mod p)，才能满足费马小定理。假设第一个a2^i*q = 1(mod p) ，是第一个(mod p) 为1的项，且i != 0，即不为第一项，那么a2^(i-1)*q 必然与 -1 模p同余。假设a2^(i-1)*q = n*p + m，a2^i*q = (n*p + m)^2 = (n*n*p +2*n)*p + m*m。根据a2^i*q = 1(mod p) ，有m*m = 1(mod p)，即p整除m*m - 1，或者写成p 整除(m+1)*(m-1)，由此可得 p整除 m+1，或p整除m-1，其中1<m<p，只有m == p-1才能满足要求。所以a2^(i-1)*q = p-1 (mod p) 。

　　现在描述Miller-Rabin素性测试：

　　P是我们要测试的数，p可以写为(2k )*q的形式，其中q是奇数。对于与p互素的整数a，如果同时满足下面两点：

　　　　1. aq != 1 (mod p)

　　　　2. a2^i*q != p-1 (mod p) , 对所有0<=i<k

　　则p一定不是素数，否则p是可能的素数（概率大于等于3/4）。如果我们用多个a测试p，都发现p是可能的素数，则p为素数的可能性会很大，一句小概率事件不发生原理，我们就可以认定p是素数。

代码(mul_mod函数用二进制实现)：

 1 //zzy2012.7.13
 2  #include<cstdio>
 3  #include<iostream>
 4  #define ll long long
 5 
 6  using namespace std;
 7 
 8 
 9  ll mul_mod(const ll &a, ll b, const ll &n)
10 {
11     ll back(0), temp(a % n);
12     b %= n;
13     while ( b > 0 )
14     {
15         if ( b & 0x1 )
16         {
17             if ( (back = back + temp) > n )
18             back -= n;
19         }
20         if ( (temp <<= 1) > n )
21         temp -= n;
22         b >>= 1;
23     }
24     return back;
25 }
26 
27 ll pow_mod(const ll &a, ll b, const ll &n)
28 {
29     ll d(1), dTemp(a % n);//当前二进制位表示的是进制数值
30     while ( b > 0 )
31     {
32         if ( b & 0x1 ){
33             d = mul_mod(d, dTemp, n);
34             b ^= 1;
35         }
36         else{
37             dTemp = mul_mod(dTemp, dTemp, n);
38             b >>= 1;
39         }
40     }
41     return d;
42 }
43 
44  bool MillerRabin(long long p,long long a){
45      long long k=0,q=p-1,m;
46      while(q%2 == 0){
47          k++;
48          q/=2;
49      }
50      m = pow_mod(a,q,p);
51      if(m==1 || m==p-1)
52          return true;
53      for(long long i=1; i<k; i++){
54          m = pow_mod(m,2,p);
55          if(m == p-1)
56              return true;
57      }
58      return false;
59  }
60 
61  int main()
62  {
63      long long n,i;
64      bool sign;
65      cin>>n;
66      sign = true;
67      for(long long i=n-1; i >= n-100LL && i >= 2LL; i--){
68          if(MillerRabin(n,i)==false){
69              sign = false;
70              break;
71          }
72      }
73      if(sign == false)
74          printf("%lld is not a prime.\n",n);
75      else
76          printf("%lld is a prime.\n",n);
77      return 0;
78  }