Miller Rabin 算法簡介

0.1 一些閒話

最近一次更新是在2019年11月12日。以前的文章有不少問題：當我把個人代碼交到LOJ上，發現只有60多分。我調了一個晚上，嘗試用{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 61, 24251, 2147483647, 998244353}這麼一大串數做爲基底，而後左改右改，總算過去了。特別感謝 @騙分過樣例 的提醒，如今張貼的代碼應該是值得信賴的了。

以前個人同窗好像就指出過個人文章的不少問題。好比說我以前寫到，Miller Rabin在面對合數$46856248255981$時，若是用{2,3,7,61,24251}作基底，是會判斷錯誤的。但實際上，他說他對着個人代碼寫了一遍，發現這個數是能夠判掉的。

OI中的數學須要細心。相比其餘算法方面，數學真的很差調試——一個公式算錯了，一個下標寫反了，程序就錯了。並且複雜的數學代碼很難用gdb查錯，只能反覆本身的公式是否寫對，而且在轉換成代碼的時候是否有差錯。編程須要細緻和求實精神。當你在寫代碼，亦或是在寫題解時，多問本身一個：有沒有問題？是否是哪裏寫錯了？這裏爲何要這麼寫？可不能夠造數據Hack？尤爲是在寫博客的時刻，每個OIer都須要作到足夠細緻——由於這些文章不是寫出來好看，讓同窗膜拜的，而是真的要幫到網絡另外一端，須要幫助的人。

這是我尤其欠缺的。最近個人博文被指出了不少錯誤，我也意識到本身的不細緻會給本身帶來這麼多麻煩。也但願我能做爲反面教材，警醒後人吧。

1.1 問題引入

判斷一個數$n$是否是質數.

當n不太大時,咱們直接用$O(\sqrt{n})$的試除法就能夠了.據說還能夠利用素數定理,只用$O(ln(x))$就能夠了?這種說法我還暫時沒試過.

不過,對於$n \geqslant 10^{18}$時,這個算法是確定不可取的.咱們須要一個更加快速的算法.誰不但願全部的問題被$O(1)$解決呢?

一想到$O(1)$,咱們不妨試一下隨機算法.在以後介紹Pollard Rho算法後,也會用到這種思想.

1.2 隨機的依據

按理來說,咱們應該是在$[1,\sqrt{n}]$裏面隨機,再經過某些方式來提升正確率的.可是這樣作太慢!咱們要確保每個在$[1,\sqrt{n}]$中的數都不是n的因子,顯然用隨機試的算法是行不通的.那麼,咱們該怎麼辦?

數論中,咱們有Fermat(也就是費馬)小定理:

$ p \in \text{Prime},\quad x^p \equiv p \pmod{p}$ 或 $ x^{p-1} \equiv 1 \pmod{p}$.

這個定理對於任何質數成立;可是反過來不必定.知足Fermat小定理的數不必定是質數.事實上,若是對任意$ b \in \mathbb{N}, b^{p-1} \equiv 1\pmod{p}$,咱們就稱p爲Carmichael數.運用Fermat小定理，咱們能夠斷定一個數不是質數，可是不能斷定一個數是質數.

不過,這給咱們的隨機枚舉提供了一個不錯的啓示:不是隨機枚舉$n$的因子,而是一個基底$b$.若對於不一樣的基底,都有$b^{p-1} \equiv 1 \pmod{p}$,p是質數的可能性就更大了.事實上,若是p知足$b^{p-1} \equiv 1 \pmod{p}$,那麼p就叫作"基於b的僞素數".

若是可能的話,這裏會張貼費馬小定理的證實.

1.3 算法的改進: 二次探測定理

這個定理叫作"二次"是有緣由的.其一,它是費馬小定理探測質數的輔助工具,做爲第二道屏障;其二,它和一個"二次同餘方程"有關.定理以下:

若$x^2 \equiv 1 \pmod{p}$,則 $x \equiv 1 \pmod{p}$ 或 $x \equiv p-1 \pmod{p}$

這個定理比較好證實.這裏是證實過程.

對於一個質數p,只要p和二次探測定理不符,那麼咱們就能夠確定p不是質數.事實上,咱們經常直接用這個定理去斷定質數,而不是用費馬小定理.

2.1 基本實現

不管是費馬小定理仍是二次探測定理,咱們都須要選取若干個合適的基底來進行測試.通常而言,咱們會選取這五個基底:$2,3,7,61,24251$。當須要測試的數據規模在$10^{16}$左右時，它的表現效果仍是能夠的。但當數據規模更大些時，你必須考慮將基底擴大一些。在這個基底的基礎上，能夠考慮選取前$10\sim 15$大的質數，從而解決幾乎$100%$的數據。

對於每個基底$b$和一個待判斷的數$x$,咱們進行以下測試: 1.用費馬小定理斷定一下這個數是否是合數. 2.若是$x-1$是偶數,咱們能夠對費馬小定理作以下變形:

$b^{x-1}\equiv 1 \pmod{x} \implies b^{\frac{x-1}{2}\cdot 2} \equiv 1 \pmod{x}$ 這樣就能夠用二次探測定理看一下有沒有$b^{\frac{x-1}{2}} \equiv 1 \text{或} x-1 \pmod{x}$ 了.若是都不知足,x就必定不是質數.

若是$\frac{x-1}{2}$一樣是一個偶數,咱們能夠繼續將它拆成$\frac{x-1}{4} \cdot 2$,而後再用一次二次探測定理來作.

3.對於二次探測定理$X^k \equiv 1 \pmod{x}$,若是$k$是奇數或者$X \equiv x-1 \pmod{x}$,此時咱們沒有辦法再拿這個數作文章了.咱們只能暫時認定$x$是質數.

注意到這就是Miller Rabin算法具備隨機性的證據之二:對質數的斷定不充分.但事實上,用上面這種方法已經能夠成功地斷定許多質數了,在實際應用中是值得信賴的.

上面這個算法的時間複雜度是$O((\log n)^2)$的：瓶頸在於快速冪，和分解待驗證的數；但實際運行時，速度是至關快的。

(以前這裏是有代碼的，可是好像有點bug。這裏就不張貼了。)

2.2 算法的優化

顯然,咱們不難看出上述代碼有大量須要優化的地方.就好比說,開始的那次費馬小定理的斷定是徹底沒有必要的.咱們能夠直接把這一過程留到二次探測定理去執行. 其次,咱們使用二次探測定理時每次都要探測$\log n$次,這個次數是能夠稍微有所優化的.

對於前者,咱們直接刪去開頭的費馬小定理判斷就能夠了.對於後者,網上廣爲流傳這樣一種優化技巧:

設$k=x-1=2^p \cdot d$,$d$是一奇數,那麼以前二次探測定理$X^2 \equiv 1 \pmod{x}$檢測的數$X$分別是以下幾個數:$b^{2^p*d},b^{2^{p-1}*d},b^{2^{p-2}*d},\dots,b^{d}$. 若是咱們從後往前檢測,若是其中一個數$X$經過了二次探測定理,就直接斷定$x$是質數.

我會嘗試着證實這個的.往後可能會貼在這裏.

這個新算法的流程以下: 1.按上面的方法計算出d和p.記$cur=b^d$.若是$cur\equiv 1 \text{或} x-1 \pmod{x}$就直接斷定x是質數. 2.每次把$cur$賦值爲$cur^2 \pmod{x}$直到$cur=b^k$.這個過程等價與將$cur$從$b^d$往$b^k$的方向掃描. 3.若是$cur = x-1$則直接斷定x爲質數.不然轉2. 4.若是上述測試都沒有斷定x爲質數,則直接斷定x爲合數.

bool mr(ll x,ll b)
{
    ll d,p=0;
    d=x-1;
    while((d&1)==0)
        d>>=1,++p;
    int i;
    ll cur=qp(b,d,x);//快速冪
    if(cur==1 || cur==x-1)
        return true;
    for(i=1;i<=p;++i)
    {
        cur=cur*cur%x; //爲了防止溢出,這裏最好用快速乘或者直接用__int128轉化一下
        if(cur==x-1)
            return true;
    }
    if(i>p)
        return false;
}
bool mr(ll x)
{
    if(x==46856248255981 || x<2)
        return false;
    if(x==2 || x==3 || x==7 || x==61 || x==24251)
        return true;
    return mr(x,2)&&mr(x,3)&&mr(x,7)&&mr(x,61)&&mr(x,24251);
}

算法的具體實現我也是參考了洛谷P4718大佬@Great_Influence 的題解,不太清楚背後的緣由.能夠參考一下他的代碼.