質數檢驗有不少演算法，一般使用的質數檢驗複雜度是\(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\)

在\([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\)

若\(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