1. 程式人生 > 其它 >Miller Rabin演算法詳解

Miller Rabin演算法詳解

原文

演算法作用

判斷一個數是否是素數

演算法依據

費馬小定理

如果P是素數,且整數a不是p的倍數有:

\[a^{p-1}\equiv 1 \pmod P \]

費馬定理只是n是素數的必要條件。即費馬定理不成立,n一定是合數;費馬定理成立,n可能是素數。

點選檢視證明

性質 1:$p-1$個整數 $a,2a,3a,...(p-1)a$中沒有一個是 $p$的倍數 

性質 2:$a,2a,3a,...(p-1)a$中沒有任何兩個同餘與模 $p$的

所以 $a,2a,3a,...(p-1)a$對模 $p$的同餘既不為零,也沒有兩個同餘相同

因此,這 $p-1$個數模 $p$的同餘一定是 $a,2a,3a,...(p-1)a$的某一種排列

即 $a*2a*3a*...*(p-1)a \equiv {1*2*3*...*(p-1)} \pmod p$

化簡為

$a^{p-1}*(p-1)! \equiv {p-1}! \pmod p$

根據威爾遜定理可知 $(p-1)!$與 $p$互質,所以同時約去 $(p-1)!$

即得到 $a^{p-1}\equiv 1 \pmod P$

二次探測定理

\(p\)為素數,a為任意整數有:

\[a^{2}\equiv 1\pmod P \]

易得

\[a\equiv \pm 1\pmod P \]
點選檢視證明

$a^{2}\equiv 1\pmod P$

$a^{2}-1\equiv 0\pmod P$

$(a+1)*(a-1)\equiv 0\pmod P$

那麼

$(a+1)\equiv 0\pmod P$

或者

$(a-1)\equiv 0\pmod P$

(此處可根據唯一分解定理證明)

即

$a\equiv \pm 1\pmod P$
$a^{2}\equiv 1\pmod P$

$a^{2}-1\equiv 0\pmod P$

$(a+1)*(a-1)\equiv 0\pmod P$

那麼

$(a+1)\equiv 0\pmod P$

或者

$(a-1)\equiv 0\pmod P$

(此處可根據唯一分解定理證明)

即

$a\equiv \pm 1\pmod P$

推理定理

首先,根據 Miller Rabin 演算法的過程

假設需要判斷的數是 \(p\)

我們把 \(p-1\)分解為 \(2^k*t\)的形式

\(p\)是素數,有 \(a ^ {2^k * t} \equiv 1 \pmod p\)

然後隨機選擇一個數 \(a\),計算出 \(a^t \pmod p\)

讓其不斷的自乘,同時結合二次探測定理進行判斷

如果我們自乘後的數 \(\pmod p = 1\),但是之前的數 \(\pmod p \not = \pm 1\)

那麼這個數就是合數 (違背了二次探測定理)

這樣乘 \(k\)次,最後得到的數就是 \(a^{p-1}\)

那麼如果最後計算出的數不為 \(1\)

,這個數也是合數 (費馬小定理)

正確性判斷

老祖宗告訴我們,若 \(p\)通過一次測試,則 \(p\)不是素數的概率為 \(25\)%

那麼經過 \(t\)輪測試,\(p\)不是素數的概率為 \(\dfrac {1}{4^{t}}\)

我習慣用 \(2,3,5,7,11,13,17,19\)這幾個數進行判斷

在資訊學範圍內出錯率為 \(0\)%(不帶高精)

程式碼

注意在進行素數判斷的時候需要用到快速冪。。

#include<cstdio>
#define LL long long 
inline int read() {
    char c = getchar(); int x = 0, f = 1;
    while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();}
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x * f;
}
int N, M, Test[10] = {2, 3, 5, 7, 11, 13, 17};
int pow(int a, int p, int mod) {
    int base = 1;
    for(; p; p >>= 1, a = (1ll * a * a) % mod) 
        if(p & 1) base = (1ll * base * a) % mod;
    return base % mod;
}
bool Query(int P) {
    if(P == 1) return 0;
    int t = P - 1, k = 0;
    while(!(t & 1)) k++, t >>= 1;
    for(int i = 0; i < 4; i++) {
        if(P == Test[i]) return 1;
        LL a = pow(Test[i], t, P), nxt = a;
        for(int j = 1; j <= k; j++) {
            nxt = (a * a) % P;
            if(nxt == 1 && a != 1 && a != P - 1) return 0;
            a = nxt;
        }
        if(a != 1) return 0;
    }
    return 1;
}
main() { 
    N = read(); M = read();    
    while(M--) puts(Query(read()) ? "Yes" : "No");
}