與數論的愛恨情仇--01：判斷大素數的Miller-Rabin

你以為你以為的就是你以為的嗎發表於2019-04-15

　　在我們需要判斷一個數是否是素數的時候，最容易想到的就是那個熟悉的O(√n)的演算法。那個演算法非常的簡單易懂，但如果我們仔細想想，當n這個數字很大的時候，這個演算法其實是不夠用的，時間複雜度會相對比較高。

　　怎麼解決呢？我們先來了解一下“費馬小定理”。假設我們有一個素數p，且另一個數a和p互素，就可以得到ap-1≡1(mod p)。這個定理很巧妙啊，有人就想了，能不能通過費馬小定理來判斷一個數是否是素數呢？也就是說，當我們判斷一個數p是否是素數時，只需要判斷ap-1≡1(mod p)是否成立即可。這裡的a因為是任意數，乾脆就讓它等於2，那麼判斷一個數p是否是素數就轉化成了判斷2p-1≡1(mod p)是否成立。乍一看，這好像沒什麼問題。當這個式子不成立時，p一定是一個合數，這沒毛病；但是當式子成立的時候p就一定是素數嗎？我們舉個反例。當p = 341時，2340≡1(mod 341)成立，然而很不巧，341 = 11 * 31是一個純正的合數。這就是數學中所說的，對於所有的a都存在對應的偽素數。(ps:這個問題並不能完全通過改變基數解決)

　　那我們該怎麼辦呢？其實很簡單，只需要進行“二次探測”把偽素數揪出來就ok了。這可不是亂改，是有依據的：當p為素數時，方程x2≡1(mod p)有兩個根 x = 1 和 x = p - 1。這兩個根被我們賦予了一個奇怪的名字：平凡平方根。那麼，判斷一個數p是否是素數這個問題就又被我們轉化，變成了判斷在模p意義下是否存在1的非平凡平方根，若存在則p為合數，反之則為素數。這一測試被我們“親切地”稱為“Miller - Rabin測試”。

　　具體操作步驟如下：

　　①選取多個基數a進行測試；

　　②尋找模p為1的非平凡平方根；令p - 1 = 2t*u(t >= 1, u為奇數)，ap-1 = a2t*u = a2t  ，先算出x=au (mod p)，再把 x 平方 t 次，每次模上 p，這樣我們就得到了一個長度為 t + 1 的序列。我們希望這個序列以 1 結尾。若中間某一項為 1，則這一項的前一項必須為 1 或 n - 1，否則p就是合數。

　　在Miller - Rabin測試中進行s次測試，這並不代表這項測試是簡單地驗證費馬小定理，它大大降低了出錯的概率，研究表明，Miller - Rabin測試的出錯概率至多為 2-s 這可以說是非常小了。所以不用擔心它的準確度和嚴謹性。

　　程式碼如下：(你沒看錯就是rand())

#include<bits/stdc++.h>
using namespace std;

int pow(long long a, long long b, long long n) {
    long long ans = 1;
    while (b) {
        if (b & 1) {
            ans = ans * a % n;
        }
        a = a * a % n;
        b >>= 1;
    }
    return ans;
}

bool judge(long long n, long long a) {
    long long u = 0, t = n - 1;
    while (t % 2 == 0) {u++; t /= 2;}
    long long x = pow(a, t, n);
    
    for (int i = 1; i <= u; i++) {
        long long nx = x * x % n;
        if ((nx == 1) && (x != 1) && (x != n - 1))
            return true;
        x = nx;
    }
    if (x != 1) return true;
    return false;
}

bool miller(long long n, int s) {
    if (n == 2) return true;
    if (n < 2 || n % 2 == 0) return false;
    for (int i = 1; i <= s; i++) {
        long long a = rand() % (n - 2) + 2;
        if (judge(n, a)) return false;
    }
    return true;
}

int main() {
    int t;
    cin >> t;
    long long a;
    while (t--) {
        cin >> a;
        if (miller(a, 10)) {
            cout << "yes" << endl;
        }
        else cout << "no" << endl;
    }
    return 0;
}

相關文章