米勒-拉宾检验

米勒-拉宾质数判定法（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;
```