Miller-Rabin随机性素数测试法

 

Miller–Rabin 随机性素数测试：

 

前言：

我们普通的判素数的方法一般就是for循环找因子、打素数表判断因子，这样的复杂度下限差不多也就O(sqrt(n))了，而对于比较大的n就难以处理了。这里介绍一下Miller-Robin测试法，虽然该算法是一种随机算法，即无法保证判断结果百分之百正确，但是该算法在绝大多数时候都表现地很好。

 

费马小定理：

首先要说明一下费马小定理：如果p是素数，那么对于小于p的任意正整数a都有：

                                                  a^(p-1)%p=1                              --①

那么也就是说如果我们能找到任意一个小于p的正整数a不满足①式，那么我们就可以说p是一个合数。而要判断素数，我们更关心的是它的逆命题是否成立，也就是当所有的小于p的正整数都满足①时，是否可以判定p是一个素数。结论是否定的，并且当①式成立且p为合数时，称p是一个基为a的伪素数。

但是乐观的是我们可以认为它“几乎成立”，例如假定a=2，然后用①式去测试一个数n是否为素数，出错的概率为多少呢？在小于10000的值中，只有22个值的判定会产生错误，最靠前的四个值分别为341，561，645，1105（来自黑书）。然而这样的错误概率在实际的IO中已经无法接受，所以我们需要更有效的算法。

上面提到只用a=2去对p进行测试虽然会出错，但是得到出错的概率已经“较小”了，那么我们是否能通过取不同的基数a去对p进行测试，从而避免出错的概率呢？回答是不能，虽然取不同的基底可以减小出错的概率，但是有这么一类数，就算你将a取遍(0，p)的所有值，①式均成立，但是p却是一个合数，这样的一类数被称为Carmichael数。前三个Carmichael数分别是561、1105、1729。Carmichael数极少，在小于100000000的数中，只有255个，尽管如此，这样的出错的概率还是不能接受。而Miller-Rabin测试法则可以解决这一问题。

 

二次探测：

在Miller-Rabin测试中不仅仅用了费马小定理进行判定，更重要的是对p的判定进行了“二次探测”，具体来说是用到了另外一条定理：如果p是一个奇素数且e>=1，则方程：x^2≡1(mod p^e)，仅有两个解，即x=1和x=-1(证明略)；且称x=1和x=-1为以n为模的1的两个“平凡”平方根。

由此定理可以得出一条推论：如果对模n存在1的非平凡平方根(x!=1&&x!=-1)，则n是合数。

其实上述的结论是不难得出的，若x^2=1(mod p)，那么有p|(x+1)或p|(x-1)，由于p为素数，故x=±1，这只是一个不规范直观理解，严谨的证明请阅读更权威的资料（黑书！）。

然后我们怎么进行二次探测呢？原来，我们只进行a^(p-1)%p=1进行验证，现在，我们先将p-1分解成2^t*u的形式，其中u为奇数(p为奇素数，故p-1一定可以得到这种形式)，我们计算X[0]=a^u%p，然后对于i∈[1，t]，我们计算X[i]=(a^u)^(2*i)%p，即有X[i]=X[i-1]^2%p，这样一来，更具上述定理及其推论，如果X[i]==1，而X[i-1]不等于±1(实际上是判断X[i]!=1&&X[i]!=-1)，我们就可以得出p是合数的结论。

 

举例：

再举一个具体的例子(来自黑书)来说明该做法排除Carmichael数561的过程：首先p-1=560=2^4*35，故t=4，u=35。假设选定a=7，则有X[0]=7^35%561=241，接下来依次可得到X[]的序列：X=<241，298，166，67，1>，我们注意到最后一次运算中得到的结果为X[4]=7^560%561=1，但是它的前一个值X[3]=7^280%561=67≠±1，故a=7可以说明561为合数。

 

正确性：

然而前面早已经说明了Miller-Rabin测试是一种随机的算法，它还是有一定的错误概率，但是这个错误的概率是由你选择的基数a的值及其个数确定的，并且其错误的概率可以证明不超过2^(-s)（实际上表现得比这要好）。我们总是不喜欢跟“不确定”的事件打交道，那么这里就给出一个关于选取基数的值和个数与其错误概率关系的结论(来自Wiki)：

• if n < 2,047, it is enough to test a = 2;

• if n < 1,373,653, it is enough to test a = 2 and 3;

• if n < 9,080,191, it is enough to test a = 31 and 73;

• if n < 25,326,001, it is enough to test a = 2, 3, and 5;

• if n < 3,215,031,751, it is enough to test a = 2, 3, 5, and 7;

• if n < 4,759,123,141, it is enough to test a = 2, 7, and 61;

• if n < 1,122,004,669,633, it is enough to test a = 2, 13, 23, and 1662803;

• if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11;

• if n < 3,474,749,660,383, it is enough to test a = 2, 3, 5, 7, 11, and 13;

• if n < 341,550,071,728,321, it is enough to test a = 2, 3, 5, 7, 11, 13, and 17.

• if n < 3,825,123,056,546,413,051, it is enough to test a = 2, 3, 5, 7, 11, 13, 17, 19, and 23.

• if n < 18,446,744,073,709,551,616 = 264, it is enough to test a = 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, and 37.

• if n < 318,665,857,834,031,151,167,461, it is enough to test a = 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, and 37.

• if n < 3,317,044,064,679,887,385,961,981, it is enough to test a = 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, and 41.

根据上述数据，我们在平时的OI处理中应用Miller-Rabin算法就已经够用了。

至此，Miller-Rabin随机性素数测试法介绍完毕，最后给出该算法的具体实现：

 

代码：

#include<iostream>
#include<cstring>
#include<string>
#include<cstdio>
#include<algorithm>
#define LL long long
using namespace std;
int prime[15]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47};
LL minfactor;
LL Mult(LL x,LL y,LL mod){        //避免两个long long相乘溢出
    LL ans=0;ans%=mod;
    if(!y)return 0;
    while(y){
        if(y&1){
            ans+=x;
            if(ans>=mod)ans-=mod;
        }
        y>>=1;x<<=1;
        if(x>=mod)x-=mod;
    }
    return ans;
}
LL Kpow(LL x,LL y,LL mod){        //快速幂
    LL ans=1;
    if(!y)return 1;
    while(y){
        if(y&1){
            ans=Mult(ans,x,mod);
        }
        y>>=1;x=Mult(x,x,mod);
    }
    return ans;
}
bool is_prime(LL x,LL mod){       //取一次基数
    LL mid=Kpow(x,mod-1,mod);
    if(mid!=1)return 0;
    LL n=mod-1;
    while(!(n%2)&&mid==1){
        n>>=1;
        mid=Kpow(x,n,mod);
    }
    if(mid==1||mid==mod-1)return 1;
    return 0;
}
bool Miller_Rabin(LL x){          //Miller-Rabin测试
    if(x<=47){
        for(int i=0;i<15;i++)if(prime[i]==x)return 1;
        return 0;
    }
    for(int i=0;i<15;i++){
        if(!is_prime(prime[i],x))return 0;
    }
    return 1;
}
int main(){
    LL T,n;
    cin>>T;
    while(T--){
        cin>>n;
        if(Miller_Rabin(n))cout<<"Prime"<<endl;
        else cout<<"Composite"<<endl;
    }
    return 0;
}

 

总结：

本人数学渣，这次将该算法总算是好好地学了一遍（以前直接套的模板＞︿＜），收获挺大。

学习过程中难免出现各种错误，欢迎路过的各位朋友批评斧正。^_^