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函数用二进制实现)：

 

//zzy2012.7.13
 #include<cstdio>
 #include<iostream>
 #define ll long long

 using namespace std;


 ll mul_mod(const ll &a, ll b, const ll &n)
{
    ll back(0), temp(a % n);
    b %= n;
    while ( b > 0 )
    {
        if ( b & 0x1 )
        {
            if ( (back = back + temp) > n )
            back -= n;
        }
        if ( (temp <<= 1) > n )
        temp -= n;
        b >>= 1;
    }
    return back;
}

ll pow_mod(const ll &a, ll b, const ll &n)
{
    ll d(1), dTemp(a % n);//当前二进制位表示的是进制数值
    while ( b > 0 )
    {
        if ( b & 0x1 ){
            d = mul_mod(d, dTemp, n);
            b ^= 1;
        }
        else{
            dTemp = mul_mod(dTemp, dTemp, n);
            b >>= 1;
        }
    }
    return d;
}

 bool MillerRabin(long long p,long long a){
     long long k=0,q=p-1,m;
     while(q%2 == 0){
         k++;
         q/=2;
     }
     m = pow_mod(a,q,p);
     if(m==1 || m==p-1)
         return true;
     for(long long i=1; i<k; i++){
         m = pow_mod(m,2,p);
         if(m == p-1)
             return true;
     }
     return false;
 }

 int main()
 {
     long long n,i;
     bool sign;
     cin>>n;
     sign = true;
     for(long long i=n-1; i >= n-100LL && i >= 2LL; i--){
         if(MillerRabin(n,i)==false){
             sign = false;
             break;
         }
     }
     if(sign == false)
         printf("%lld is not a prime.\n",n);
     else
         printf("%lld is a prime.\n",n);
     return 0;
 }