Miller-Rabin素性测试算法

原理

由费马小定理可知,若p为质数,对于任意与p互质的整数a,有a^p-1≡1(modp),假设我们要测试的数是x,然后在1到p-1内随机生成一个数作为底数a,然后测试它是否符合费马小定理,如果不符合则一定不是素数,符合则有可能是素数.单纯用费马小定理检验素数出错率很高(指满足费马小定理又不是素数的数-这些数称为Carmichael数,也称弱可能素数),因此需要一个更强的检定方法.

 二次探测定理:p是质数,x²≡1(modp),x的解只能为1或p-1.利用这条定理就可以大大加强我们的素性检测.具体做法如下,对于被测数x,首先检验它是否能通过费马小定理,如果通过了,看p-1是否为偶数(只有p为2的情况不符合),如果是偶数这利用二次检测定理检测x^(p-1)/2是否为1或p-1,如果是p-1的话就不能检测下去了,如果是1的话,则有x^(p-1)/2≡1(modp),检查(p-1)/2是否为偶数,如果是偶数的话又能再次利用二次探测定理,检测x^(p-1)/4的值是否为1或p-1.这样重复检测下去直到xⁿ=p-1或n为奇数时才停止检测,在这期间任何一次检测不通过则可以断定x为合数.由于只用一个底数检测出错率还是不够低,因此一般是用多个底数进行判定,这样出错率才能降到令人接受的程度.由整个检测过程可知,我们刚开始生成的底数a的取值范围可以改进,想一下1ⁿ=1(modp),(p-1)ⁿ=1(modp)或(p-1)ⁿ=-1(modp)恒成立(p-1的用二项式拆开就证明).因此只要在2到p-2内随机生成一个数作为底数a即可.

不过经过一些大佬的验证,如果待测数是在1至1e15里面的话,用固定底数2,3,7,61,24251测速度会快很多,而且不会出错,这组底数能出错的第一个数是46856248255981.这是个RP算法,多次运行可以降低出错率,假设测的次数为n,出错率大约(1/4)ⁿ.

代码

 

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#include <utility>
#include <ctime>
#include <cmath>
#include <cstring>
#include <string>
#include <stack>
#include <queue>
#include <vector>
#include <set>
#include <map>

using namespace std;
typedef long long LL;
const int INF_INT=0x3f3f3f3f;
const LL INF_LL=0x3f3f3f3f3f3f3f3f;
//防爆乘法
LL q_muti(LL a,LL b,LL p)
{
    LL res=0;
    a%=p;
    while(b)
    {
        if(b&1) res=(res+a)%p;
        a=(a<<1)%p;
        b>>=1;
    }
    return res;
}

LL q_power(LL a,LL b,LL p)
{
    LL res=1;
    a%=p;
    while(b)
    {
        if(b&1) res=q_muti(res,a,p);
        a=q_muti(a,a,p);
        b>>=1;
    }
    return res;
}

bool m_l(LL x,LL a)
{
    if(q_power(a,x-1,x)!=1) return false;
    LL q=x-1;
    while(!(q&1))
    {
        q>>=1;
        LL t=q_power(a,q,x);
        if(t==1) continue;
        if(t==x-1) break;
        return false;
    }
    return true;
}

bool is_prime(LL x)
{
    if(x==0||x==1) return false;
    if(x==2||x==3) return true;
    else
    {

        int k=10; //设置测试底数个数
        for(int i=0;i<k;i++)
        {
            LL a=rand()%(x-3)+2;
            if(!m_l(x,a)) return false;
        }
        return true; //所有底数测试通过即可判为质数
//      if(x==2||x==3||x==7||x==61||x==24251) return true;
//      return m_l(x,2)&&m_l(x,3)&&m_l(x,7)&&m_l(x,61)&&m_l(x,24251);
    }
}

int main()
{
//    freopen("std.in","r",stdin);
//    freopen("std.out","w",stdout);
    srand(time(0));
    LL n;
    while(~scanf("%lld",&n))
    {
        if(is_prime(n)) printf("YES\n");
        else printf("NO\n");
    }
    return 0;
}

 

转载于:https://www.cnblogs.com/VBEL/p/11421434.html