米勒-拉賓檢驗

米勒-拉賓質數判定法（Miller–Rabin primality test）是一種質數判定法則，利用隨機化算法判斷一個數是合數還是可能是素數。1976年，卡內基梅隆大學的電腦系教授蓋瑞·米勒 (電腦科學家)首先提出了基於廣義黎曼猜想的確定性算法，由於廣義黎曼猜想並沒有被證明，於1980年，由以色列耶路撒冷希伯來大學的麥可·拉賓教授作出修改，提出了不依賴於該假設的隨機化算法。

## 概念

首先介紹一個相關的引理。我們發現 $1^2 \bmod p$ 和 $(-1)^2\bmod p$ 總是得到 $1$，我們稱 $-1$ 和 $1$ 是 $1 \bmod p$ 的「平凡平方根」，當 $p$ 是素數且 $p > 2$ 時，不存在 $1 \bmod p$ 的「非平凡平方根」。為了證明該引理，首先假設 $x$ 是 $1 \bmod p$ 的平方根之一，於是有

$$x^2 \equiv 1 \pmod{p}$$

$$(x + 1)(x - 1) \equiv 0 \pmod{p}$$

也就是說，素數 $p$ 能夠整除 $(x - 1)(x + 1)$ 或者 $x = p - 1$，根據歐幾里得引理，$x - 1$ 或者 $x + 1$ 能夠被 $p$ 整除，即 $x \equiv 1 \pmod{p}$ 或 $x \equiv -1 \pmod{p}$,

即 $x$ 是 $1 \bmod p$ 的平凡平方根。

現在假設 $n$ 是一個素數，且 $n > 2$。於是 $n - 1$ 是一個偶數，可以被表示為 $2^s*d$ 的形式，$s$ 和 $d$ 都是正整數且$d$是奇數。對任意在 $(Z/nZ)^*$ 範圍內的 $a$，必須滿足以下兩種形式的一種：

$$a^d \equiv 1 \pmod{n} \text{ ① }$$ $$a^{2^rd} \equiv -1 \pmod{n}\text{ ② }$$

其中 $r$ 是某個滿足 $0 \leq r \leq s - 1$ 的整數。

因為由於 費馬小定理，對於一個素數 $n$，有

$$a^{n - 1} \equiv 1 \pmod{n}$$

又由於上面的引理，如果我們不斷對$a^{n - 1}$取平方根後，我們總會得到 $1$ 和 $-1$。如果我們得到了 $-1$，意味著②式成立，$n$ 是一個素數。如果我們從未得到 $-1$，那麼通過這個過程我們已經取遍了所有 $2$ 的冪次，即①式成立。

米勒-拉賓素性測試就是基於上述原理的逆否，也就是說，如果我們能找到這樣一個 $a$，使得對任意$0 \leq r \leq s - 1$以下兩個式子均滿足：

$$a^d \not\equiv 1 \pmod{n}$$ $$a^{2^rd} \not\equiv -1 \pmod{n}$$那麼 $n$ 就不是一個素數。這樣的 $a$ 稱為 $n$ 是合數的一個憑證（witness）。否則 $a$ 可能是一個證明 $n$ 是素數的「強偽證」（strong liar），即當$n$確實是一個合數，但是對當前選取的 $a$ 來說上述兩個式子均不滿足，這時我們認為$n$是基於$a$的大機率素數。

每個奇合數 $n$ 都有很多個對應的憑證 $a$，但是要生成這些 $a$ 並不容易。當前解決的辦法是使用機率性的測試。我們從集合 $(Z/nZ)^*$ 中隨機選擇非零數 $a$，然後檢測當前的 $a$ 是否是 $n$ 為合數的一個憑證。如果 $n$ 本身確實是一個合數，那麼大部分被選擇的 $a$ 都應該是 $n$ 的憑證，也即通過這個測試能檢測出 $n$ 是合數的可能性很大。然而，仍有極小機率的情況下我們找到的 $a$ 是一個「強偽證」（反而表明了 $n$ 可能是一個素數）。通過多次獨立測試不同的 $a$，我們能減少這種出錯的機率。

對於測試一個大數是否是素數，常見的做法隨機選取基數$a$，畢竟我們並不知道憑證和偽證在$[1, n - 1]$這個區間如何分布。典型的例子是 Arnault 曾經給出了一個397位的合數$n$，但是所有小於307的基數$a$都是$n$是素數的「強偽證」。不出所料，這個大合數通過了 Maple 程序的`isprime()` 函數（被判定為素數）。這個函數通過檢測 $a = 2,3,5,7,11$ 這幾種情況來進行素性檢驗。

## 例子

假設我們需要檢驗 $n = 221$ 是否是一個素數。首先$n - 1 = 220 = (2^2) * 55$，於是我們得到了$s = 2$和$d = 55$.我們隨機從$[1, n - 1]$中選取$a$，假設$a = 174$，於是有：

$$a^{2^0d} \bmod n = 174^{55} \bmod 221 = 47 \not= 1,-1$$ $$a^{2^1d} \bmod n = 174^{110} \bmod 221 = 220 = -1$$因為我們得到了一個 -1，所以要麼$n = 221$確實是一個素數，要麼$a = 174$是一個「強偽證」。我們再選取$a = 137$，於是有：

$$a^{2^0d} \bmod n = 137^{55} \bmod 221 = 188 \not= 1,-1$$ $$a^{2^1d} \bmod n = 137^{110} \bmod 221 = 205 \not= -1$$即$a = 137$是$n = 221$為合數的一個憑證，而$a = 174$是一個「強偽證」。

選取特定的整數可以在一定範圍內確定（而非單純基於機率猜測）某個整數是質數還是合數。對於小於$2^{32}$的情形，選取$2, 7, 61$共三個憑據可以做到這一點；對於小於$2^{64}$的情形，選取$2, 325, 9375, 28178, 450775, 9780504, 1795265022$共七個憑據可以做到這一點。

## 演算法複雜度

這一演算法可以被表示成如下偽代碼：

```
'''Input #1''': ''n'' > 3, an odd integer to be tested for primality;
'''Input #2''': ''k'', a parameter that determines the accuracy of the test
'''Output''': ''composite'' if ''n'' is composite, otherwise ''probably prime''

write ''n'' − 1 as 2<sup>''r''</sup>·''d'' with ''d'' odd by factoring powers of 2 from ''n'' − 1
WitnessLoop: '''repeat''' ''k'' times:
   pick a random integer ''a'' in the range [2, ''n'' − 2]
   ''x'' ← ''a''<sup>''d''</sup> mod ''n''
   '''if''' ''x'' = 1 or ''x'' = ''n'' − 1 '''then'''
      '''continue''' WitnessLoop
   '''repeat''' ''r'' − 1 times:
      ''x'' ← ''x''<sup>2</sup> mod ''n''
      '''if''' ''x'' = ''n'' − 1 '''then'''
         '''continue''' WitnessLoop
   '''return''' ''composite''
'''return''' ''probably prime''
```

使用模冪運算，這個演算法的時間複雜度是$O(k\log ^3 n)$，$k$是我們測試的不同的 $a$ 的值。也就是說，這是一個高效的多項式時間演算法。使用快速傅立葉變換能夠將這個時間推進到 O(k log2n log log n log log log n) = Õ(k log2n).

如果我們加入最大公因數演算法到上述演算法中，我們有時候可以得到 $n$ 的因數，而不僅僅是證明 $n$ 是一個合數。例如，若 $n$ 是一個基於$a$ 的可能的素數，但不是一個大機率素數，則$\gcd((a^d \bmod n ) - 1, n )$或$\gcd((a^{2^rd} \bmod n ) - 1, n )$將得到 $n$ 的因數。如果因式分解是必要的，這一$GCDs$演算法可以加入到上述的演算法中，代價是增加了一些額外的運算時間。

例如，假設 $n = 341$，則$n - 1 = 340 = 85*4$.於是$2^{85} \bmod 341 = 32$，這也告訴我們 $n$ 不是一個大機率素數，即 $n$ 是一個合數。如果這個時候我們求最大公因數，我們可以得到一個$n = 341$的因數：$\gcd((2^{85} \bmod 341) - 1, 341) = 31$.這時可行的，因為$n = 341$是一個基於2的偽素數，但不是一個「強偽素數」。

## 示例代碼

下面是根據以上定義而使用Magma語言編寫的「米勒-拉賓」檢驗程序。

```matlab
function ModPotenz(a,b,n)
erg:=1;
while (b ne 0) do
       b, rest:=Quotrem(b,2);
       if (rest ne 0) then erg:=((erg*a) mod n); end if;
       a:= (a^2 mod n);
     end while;
     return erg;
end function;
;
function Miller_rabin(n)
  if (n mod 2 ne 0) then
  d:=n-1; s:=0;t:=50;
  while (d mod 2) eq 0 do
    d:=d div 2;
    s:=s+1;
  end while;
  k:=0;
  while(k lt t) do
    a:=Random(1,n-1);
    k:=k+1;
    if GCD(a,n) ne 1 then
      continue;
    end if;
    x:=ModPotenz(a,d,n);
    if(x ne 1) then
      for r in [0,s-1] do
        x:=ModPotenz(a,2^r*d,n);
        if(x eq (n-1)) then
           return "probably prime";
        end if;
      end for;
    elif (x eq 1) then
      break;
    end if;
  end while;
  end if;
  return "composite";
end function;
```