质数检验有不少算法，一般使用的质数检验复杂度是\(O(\sqrt{n})\)；

又如线性筛可以在\(O(n)\)的时间内求出所有1~n的质数

但是，当n非常大，连\(O(\sqrt{n})\)的复杂度也难以接受时，上述算法便不能满足要求

这篇blog记录了一些关于Miller_Rabin算法的内容

\[a^{p-1}\equiv1\pmod p\]

其中\(a,p\)互质

我们猜想，任意选取\(a\)，如果一个数\(p\)满足以上式子，那么它就很有可能是一个质数

但是这个猜想很容易找到反例：\(a=2\)，\(p=341,561, 645, 1105\)时，费马小定理逆命题不成立，人们把\(a^{p-1}\equiv1\pmod p\)的合数称为以\(a\)为底的“伪素数”

在\([1,10^9]\)中，质数一共有\(50847534\)个，而满足\(2^{p-1}\equiv1\pmod p\)的合数\(p\)有\(5597\)个，算法出错的概率太高了

一个想法是，同时使用多个a来进行判断，例如同时检验\(a=2,3\)

在\([1,10^9]\)中，同时以\(2,3\)为底的伪素数只有\(1272\)个，不到以2为底的\(\frac{1}{4}\)

选取的\(a\)越多时，算法越准确，这便是费马素性检验

（以下内容引自Matrix67的博客）

人们自然会想，如果考虑了所有小于n的底数a，出错的概率是否就可以降到0呢？没想到的是，居然就有这样的合数，它可以通过所有a的测试（这个说法不准确，详见我在地核楼层的回复）。Carmichael第一个发现这样极端的伪素数，他把它们称作Carmichael数。你一定会以为这样的数一定很大。错。第一个Carmichael数小得惊人，仅仅是一个三位数，561。前10亿个自然数中Carmichael数也有600个之多。Carmichael数的存在说明，我们还需要继续加强素性判断的算法。

由此观之，费马素性检验仍不能满足我们的要求，我们需要改进素性检验算法

对于\(\forall a<p\)，\(p\)是质数

若\(a^2\equiv1\pmod p\)

那么有\(a=1\)或\(a=p-1\)

证明：

\[a^2\equiv1\pmod p\]

\[\Rightarrow a^2-1=(a+1)*(a-1)\equiv0\pmod p\]

Miller_Rabin素数测试的方法是，对于待检验数\(x\)，不断地提取\(x-1\)中的因子\(2\)，把\(x-1\)表示成\(2^r*d\)的形式，那么我们需要计算的东西变成了\(a^{2^r*d}\mod x\)

于是，如果\(x\)是质数，\(a^{2^r*d}\mod x\)要么等于\(1\)，要么等于\(x-1\)

如果\(a^{2^{r-1}*d}\mod x =1\)，定理同样适用于\(a^{2^{r-2}*d}\mod x\)

不断地开方，直到存在一个\(i\in[0,r)\)，使得\(a^{2^i*d}\mod x=x-1\)

至此，费马素性检验被强化为如下形式：

尽可能提取因子\(2\)，将待检验数\(x\)表示为\(2^r*d\)的形式，其中\(d\)是一个奇数

如果\(a^d\mod n=1\)，或者找到一个\(i\in[0,r)\)，使得\(a^{2^i*d}\mod x=x-1\)，那么\(x\)极有可能是一个质数，否则\(x\)就不是一个质数

上文提到极有可能，是因为Miller_Rabin同样是不确定算法，我们称能够通过以\(a\)为底的Miller_Rabin检验的合数称为以\(a\)为底的“强伪素数”，第一个以\(2\)为底的强伪素数为\(2047\)，而第一个以\(2\)和\(3\)为底的强伪素数则为1373653

如果选用\(2,3,7,61\)和\(24251\)作为\(a\)，那么\(10^{16}\)内唯一的强伪素数为\(46\space856\space248\space255\space981\)，正确率可以接受

Miller_Rabin核心代码

这里，我采用先将n-1的所有因子2提取出来，再不断地平方回到原来的n-1的枚举方法

typedef long long ll;
const ll a[10]={0,2,3,7,61,24251};

ll fastpow(ll a,ll b,ll p)
{
    ll re=1,base=a;
    while(b)
    {
        if(b&1)
            re=re*base%p;
        base=base*base%p;
        b>>=1;
    }
    return re;
}

bool Miller_Rabin(ll x,ll a)
{
    if(x==a)
        return true;
    if(x%a==0 || x<2)
        return false;
    ll d=x-1,r=0;
    while(!(d&1))
    {
        d>>=1;
        ++r;
    }
    ll k=fastpow(a,d,x);
    if(k==1 || k==x-1)
        return true;
    for(ll i=1;i<r;++i)
    {
        k=k*k%x;
        if(k==x-1)
            return true;
    }
    return false;
}

bool test(ll x)
{
    for(int i=1;i<=5;++i)
        if(!Miller_Rabin(x,a[i]))
            return false;
    return true;
}

题外话

由于水平及时间有限，博客中可能存在较多疏漏，希望读者能够指出，非常感谢

撰写过程中，很多内容参考了matrix67的博客，在此表示感谢

参考出处：http://www.matrix67.com/blog/archives/234

原文地址：https://www.cnblogs.com/lizbaka/p/Miller_Rabin.html