【演算法】米勒-拉賓素性檢驗

問題描述

問題來源於 Sphere Online Judge (SPOJ) 網站的的“Prime or Not”題目：

Prime or Not: Given the number, you are to answer the question: "Is it prime?"

Input:
t – the number of test cases, then t test cases follows. [t ≤ 500]
Each line contains one integer: N [2 ≤ N ≤ 263-1]

Output: For each test case output string "YES" if given number is prime and "NO" otherwise.

Time limit: 21s

Source limit: 5,000 Bytes

這道題目要求我們判斷大約五百個給定的數是不是素數，其中每個數不超過 263-1 。時間限制是 21 秒。如果通過試除法進行因子分解，肯定會逾時。如果使用篩法，肯定會造成記憶體超限。

演算法描述

首先，讓我們來看看跟素數有關的費馬定理：

如果 p 是素數，則對於所有整數 a，ap ≡ a (mod p)。

根據上述的費馬定理，每當 p 是素數以及 a 不是 p 的倍數時，我們有 a p - 1 ≡ 1 (mod p) 。而且有有效方法計算 a n - 1 mod n ，且只需要 O(log n) 個模 n 乘法運算。因此我們可以確定，當這個關係不成立時，n 不是素數。對於一個給定的數的非素性來說，費馬定理是一個強有力的檢驗。當 n 不是素數時，總是有可能來求 a < n 的一個值，使得 a n - 1 ≠ 1 (mod n) 。事實上，經驗證明，這樣一個值幾乎總能非常快地求出。有某些稀少的 n 值，它們經常使得 a n - 1 ≡ 1 (mod n) ，但在此情況下，n 有小於 n 1/3 的因子。

我沒有找到確定一個很大的數是否素數的有效演算法。但是 Miller-Rabin primatlity test 演算法能夠以很高的機率來檢驗一個很大的數是否素數。該演算法描述如下：

Input: n > 3, an odd integer to be tested for primality;
Input: k, a parameter that determines the accuracy of the test
Output: composite if n is composite, otherwise probably prime
01: write n − 1 as 2 s· d with d odd by factoring powers of 2 from n − 1
02: LOOP: repeat k times:
03:   pick a randomly in the range [2, n − 2]
04:    x ← a d mod n
05:   if x = 1 or x = n − 1 then do next LOOP
06:   for r = 1 .. s − 1
07:      x ← x2 mod n
08:     if x = 1 then return composite
09:     if x = n − 1 then do next LOOP
10:   return composite
11: return probably prime

構成該演算法的思想是，如果 a d ≠ 1 (mod n) 以及 n = 1 + 2s · d 是素數，則值序列

a d mod n， a 2 d mod n， a 4 d mod n，…， a 2 s d mod n

將以 1 結束，而且在頭一個 1 的前邊的值將是 n – 1 (當 p 是素數時，對於 y 2 ≡ 1 (mod p) ，僅有的解是 y ≡ ±1 (mod p)，因為 (y + 1)(y - 1)必須是 p 的倍數)。注意，如果在該序列中出現了 n – 1，則該序列中的下一個值一定是 1。因為：(n – 1)2  ≡  n2 – 2n + 1  ≡  1  (mod n)。在該演算法中：

• 該演算法用於判斷一個大於 3 的奇數 n 是否素數。參數 k 用於決定 n 是素數的機率。
• 該演算法能夠肯定地判斷 n 是合數，但是只能說 n 可能是素數。
• 第 01 行，將 n – 1 分解為 2s·d  的形式，這裡 d 是奇數。
• 第 02 行，將以下步驟(第 03 到 10 行)迴圈 k 次。
• 第 03 行，◇在 [2, n - 2] 的範圍中獨立和隨機地選擇一個正整數 a 。
• 第 04 行，◇計算該序列的第一個值：x ← ad mod n 。
• 第 05 行，◇如果該序列的第一個數是 1 或者 n - 1，符合上述條件，n 可能是素數，轉到第 03 行進行一下次迴圈。
• 第 06 行，◇迴圈執行第 07 到 09 行，順序遍曆該序列剩下的 s – 1 個值。
• 第 07 行，◇◇計算該序列的下一個值：x ← x2 mod n 。
• 第 08 行，◇◇如果這個值是 1 ，但是前邊的值不是 n - 1，不符合上述條件，因此 n 肯定是合數，演算法結束。
• 第 09 行，◇◇如果這個值是 n - 1，因此下一個值一定是 1，符合上述條件，n 可能是素數，轉到第 03 行進行下一次迴圈。
• 第 10 行，◇發現該序列不是以 1 結束，不符合上述條件，因此 n 肯定是合數，演算法結束。
• 第 11 行，已經對 k 個獨立和隨機地選擇的 a 值進行了檢驗，因此判斷 n 非常有可能是素數，演算法結束。

在一次檢驗中，該演算法出錯的可能頂多是四分之一。如果我們獨立地和隨機地選擇 a 進行重複檢驗，一旦此演算法報告 n 是合數，我們就可以確信 n 肯定不是素數。但如果此演算法重複檢驗 25 次報告都報告說 n 可能是素數，則我們可以說 n “幾乎肯定是素數”。因為這樣一個 25 次的檢驗過程給出關於它的輸入的錯誤資訊的機率小於 (1/4)25。這種機會小於 1015 分之一。即使我們以這樣一個過程驗證了十億個不同的素數，預料出錯的機率仍將小於百萬分之一。因此如果真出了錯，與其說此演算法重複地猜測錯，倒不如說由於硬體的失靈或宇宙射線的原因，我們的電腦在它的計算中丟了一位。這樣的機率性演算法使我們對傳統的可靠性標準提出一個問號：我們是否真正需要有素性的嚴格證明。(以上文字引用自 Donald E. Knuth 所著的《電腦程式設計藝術 第2卷 半數值演算法(第3版)》第 359 頁“4.5.4 分解素因子”中的“演算法P(機率素性檢驗)”後面的說明) 

Ruby 程式

根據上述演算法，編寫出相應的 Ruby 程式如下：

def modPow a, b, m  v = 1  p = a % m  while b > 0 do    v = (v * p) % m if (b & 1) != 0    p = (p * p) % m    b >>= 1  end  venddef witness a, n  n1 = n - 1  s2 = n1 & -n1  x = modPow a, n1 / s2, n  return false if x == 1 || x == n1  while s2 > 1 do    x = (x * x) % n    return true if x == 1    return false if x == n1    s2 >>= 1  end  trueenddef probably_prime? n, k  # http://en.wikipedia.org/wiki/Miller-Rabin_primality_test  # n, an integer to be tested for primality  # k, a parameter that determines the accuracy of the test  return true if n == 2 || n == 3  return false if n < 2 || n % 2 == 0  k.downto(1) do    return false if witness rand(n - 3) + 2, n  end  trueend# http://www.spoj.pl/problems/PON/gets.to_i.downto(1) do  puts probably_prime?(gets.to_i, 1) ? 'YES' : 'NO'end

在該網站提交，結果是“accepted”，已耗用時間 0.23 秒，記憶體佔用 4.7 MB ，目前在 Ruby 語言中排名第三位。

C 程式

相應的 C 語言程式如下所示：

01:  #include <stdio.h>02:  #include <stdlib.h>03:  #include <time.h>04:  05:  typedef unsigned long long U8;06:  typedef int bool;07:  08:  const bool true = 1;09:  const bool false = 0;10:  11:  U8 modMultiply(U8 a, U8 b, U8 m)12:  {13:    return a * b % m;14:  }15:  16:  U8 modPow(U8 a, U8 b, U8 m)17:  {18:    U8 v = 1;19:    for (U8 p = a % m; b > 0; b >>= 1, p = modMultiply(p, p, m))20:      if (b & 1) v = modMultiply(v, p, m);21:    return v;22:  }23:  24:  bool witness(U8 a, U8 n)25:  {26:    U8 n1 = n - 1, s2 = n1 & -n1, x = modPow(a, n1 / s2, n);27:    if (x == 1 || x == n1) return false;28:    for (; s2 > 1; s2 >>= 1)29:    {30:      x = modMultiply(x, x, n);31:      if (x == 1) return true;32:      if (x == n1) return false;33:    }34:    return true;35:  }36:  37:  U8 random(U8 high)38:  {39:    // http://www.cppreference.com/wiki/c/other/rand40:    return (U8)(high * (rand() / (double)RAND_MAX));41:  }42:  43:  // http://en.wikipedia.org/wiki/Miller-Rabin_primality_test44:  // n, an integer to be tested for primality45:  // k, a parameter that determines the accuracy of the test46:  bool probablyPrime(U8 n, int k)47:  {48:    if (n == 2 || n == 3) return 1;49:    if (n < 2 || n % 2 == 0) return 0;50:    while (k-- > 0) if (witness(random(n - 3) + 2, n)) return false;51:    return true;52:  }53:  54:  // http://www.spoj.pl/problems/PON/55:  int main()56:  {57:    srand(time(NULL));58:    int t;59:    scanf("%d", &t);60:    while (t-- > 0)61:    {62:      U8 n;63:      scanf("%lu", &n);64:      puts(probablyPrime(n, 3) ? "YES" : "NO");65:    }66:    return 0;67:  }

在該網站提交，運行結果是“wrong answer”。這是由於來源程式中第 11 到 14 行的 modMultiply 函數運算溢出造成的，雖然該函數的參數和傳回值都能夠表示題目要求的數值範圍(不超過 64 bits 整數的範圍)，但是其中的乘法的中間的結果需要 128 bits 的整數才能表達，造成了溢出。而 C 語言中並沒有內建的大整數類型。

參考資料

1. Wikipedia: Fermat’s little theorem
2. Wikipedia: Miller-Rabin primality test
3. Wikipedia: Sieve of Eratosthenes
4. The Art of Computer Programming

演算法和資料結構目錄