Miller-Rabin素性测试

今天分享一个特别牛的判断一个大数是否为素数的方法，该方法基本可以通吃所有的关于判断素数的问题，它不像是传统的素数判定方法一样只适用于较小素数的判断，反之，数越大，判断正确率越高。

但美中不足的是仍然存在少量的Carmichael数无法准确判断，比如561、1105等，但这种数很少，1~一亿之间只有255个，关键是准确性高且效率高。

咖啡你冲不冲？冲~冲~冲~

那废话不多说，进入今天的重头戏。

一、二次探测定理：

我们使用Miller-Rabin素性测试方法，首先要明白什么是二次探测定理，所谓二次探测定理是指：如果p是一个奇素数，且e≥1，则方程仅有两个解：x=1和x=-1。当e=1时该方程仅有两个解：x=1和x=p-1。这两个解称为x对模p来说1的平凡平方根。

根据这个定理我们可以这样说：如果对模p存在1的非平凡平方根，那么p为合数！！（你好好悟这句话）

二、步骤

（1）输入一个数n（n>2），且n为奇数，测试它是否为素数。

不成立，那么a为合数－－－－费马定理

（2）由于我们想要解决的是测试大数，所以n－1是很大的，这个地方我们可以把n－1表示成幂的形式， 借助快速幂方法求解。

令（其中u为奇数，t为正整数）

n－1的二进制表示是u的二进制表示后边加上t个0

所以（

接下里我们就可以先计算a^u（mod n),然后连续平方t次取模就可以了

（3）判断

如果模运算结果不为1，说明不满足费马定理，n为合数。

如果模运算结果等于1，但是出现了1的非平凡平方根，说明不满足二次探测定理，n为合数。

（4）测试

选取多个（一般大于50次）就可以保证出错率忽略不计。

三、代码：

本代码依据hdu题库：How many prime numbers

#include<iostream>
using namespace std;
long long int fast_pow(long long int a,long long int u,long long int n)
{
    if(u==0)
        return 1;//任何数的0次幂都为1
    long long int x=fast_pow(a,u/2,n);//折半计算，提高效率
    long long int sum=x*x%n;//因为进行了折半，所以要左右两部分相乘取模
    if(u%2==1)//如果n为奇数，那么二分之后中间的a会少乘一次
        sum=sum*a%n;
    return sum;
}
bool witness(long long int a,long long int n)
{
    long long int u=n-1;
    int t=0;
    while(u&1==0)//u&1=0说明u的末尾数字为0，此时u需要右移一位，同时t+1
    {
        u=u>>1;
        t++;
    }
    long long int x1,x2;//用x1记录每次模计算结果，用x2记录平凡平方根
    x1=fast_pow(a,u,n);//先算a^u（mod n)
    while(t--)
    {
        x2=fast_pow(x1,2,n);
        if(x2==1&&x1!=1&&x1!=n-1)//不满足上述第二句判断条件，即不满足二次探测定理，n为合数
            return true;
        x1=x2;
    }
    if(x1!=1)//不满足上述第一句判断条件，即不满足费马定理，n为合数
        return true;
    return false;
}
int miller_rabin(long long int n,int s)
{
    if(n<2)return 0;//2是最小的素数
    if(n==2)return 1;
    if(n%2==0)
        return 0;//大于2的偶数一定为合数
    for(int i=1;i<n&&i<s;i++)
    {
        long long int a=rand()%(n-1)+1;//a随机取一个小于n的数
        if(witness(a,n))//判断是否为合数
            return 0;
    }
    return 1;
}
int main()
{
    int k;//k次测试
    while(scanf("%d",&k)!=EOF)
    {
        int ans=0;//记录素数的个数
        for(int i=1;i<=k;i++)
        {
            long long int n;
            cin>>n;
            int s=50;//a试值次数设置为50次就可以保证几乎不出错
            ans+=miller_rabin(n,s);
        }
        cout<<ans<<endl;
    }
    return 0;
}