Miller–Rabin 素性测试

整理了 Miller–Rabin 素性测试的原理及实现方法.

为了探究一奇数 $n>1$ 是否为素数, 将其用 $e$ 和 $k$ 来表示, 使得 $n-1=2^ek$, 其中 $k$ 为奇数.

根据 Keith Conrad^[1] 的推导, 多项式 $x^{n-1}-1=x^{2^ek}-1$ 可以被分解, 直到指数不是 $2$ 的倍数：

$$\begin{aligned}x^{n-1}-1&=x^{2^ek}-1\\&=(x^{2^{e-1}k}-1)(x^{2^{e-1}k}+1)\\&=(x^{2^{e-2}k}-1)(x^{2^{e-2}k}+1)(x^{2^{e-1}k}+1)\\&\;\;\vdots\\&=(x^k-1)(x^k+1)(x^{2k}+1)(x^{4k}+1)\cdots(x^{2^{e-1}k}+1).\end{aligned}$$

若 $n$ 为素数, 且 $1\leq a\leq n-1$, 那么根据费马小定理, 可知 $a^{n-1}-1\equiv0\;({\rm mod}\;n)$, 再根据上式的分解过程, 可以得出

$$(a^k-1)(a^k+1)(a^{2k}+1)(a^{4k}+1)\cdots(a^{2^{e-1}k}+1)\equiv0\;({\rm mod}\;n).$$

所以上式中的某一项必须等于 $0\;({\rm mod}\;n)$, 也就是

$$a^k\equiv1\;({\rm mod}\;n)\;\text{or}\;a^{2^ik}\equiv-1\;({\rm mod}\;n)\;\text{for some}\;i\in\{0,\ldots,n-1\}.$$

像 Fermat test、Miller–Rabin test 等基于概率的素性测试算法, 目的是找出能够证明 $n$ 是合数的证据. 若找不到这样的证据, 那么 $n$ 很可能是素数. 用这类方法找出的素数称为伪素数(pseudo prime).

对于一个奇数 $n>1$, 在 $\{1,\ldots,n-1\}$ 中取一整数 $a$, 如果某个 $a$ 使上式不成立, 也就是

$$a^k\not\equiv1\;({\rm mod}\;n)\;\text{and}\;a^{2^ik}\not\equiv-1\;({\rm mod}\;n)\;\text{for all}\;i\in\{0,\ldots,n-1\},\quad(*)$$

那么称这个 $a$ 为一个 Miller–Rabin witness, 在素性测试中, 术语 “witness” 意为某个能够证明 $n$ 为合数的数.

可以证明^[1], 若一奇数是合数, 那么 $\{2,\ldots,n-2\}$ 中超过 $75\%$ 的数都是 Miller–Rabin witness. 所以, 若该算法对一个数进行 $k$ 轮检测, 并判定其为素数, 那么该数确实是素数的可能性至少有 $1-4^{-k}$.

我在 SICP 的 练习 1.28 中用 Scheme 实现了该素性测试, 其中用非平凡平方根(nontrivial square root) 的概念巧妙地找出了 Miller–Rabin witness, 关键点是修改后的 expmod 过程, 在本来执行 square 的地方加入了检测非平凡平方根的过程.

我也用 C++ 实现了一遍, 当然不再是递归版的, 思考方式变化极大, 不过更加直接地体现了 $(*)$ 式.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

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

ll expmod(ll base, ll ex, ll mod)
{
    ll ans = 1, x = base;
    while (ex) {
        if (ex & 1)
            ans = (ans * x) % mod;
        x = (x * x) % mod;
        ex >>= 1;
    }
    return ans;
}

bool miller_rabin(ll n, int rounds)
{
    ll k = n - 1, e = 0;
    while (!(k & 1))
        k >>= 1, ++e;
    while (rounds--) {
        ll a = (rand() % (n - 1)) + 1, b = expmod(a, k, n);
        if (b == 1 || b == n - 1)
            continue;
        bool composite = true;
        for (int i = 1; composite && i < e; ++i)
            if ((b = (b * b) % n) == n - 1)
                composite = false;
        if (composite)
            return false;
    }
    return true;
}

参考资料

[1] Keith Conrad, “The Miller–Rabin Test,” https://kconrad.math.uconn.edu/blurbs/ugradnumthy/millerrabin.pdf.