快速检验素数：Miller-Rabin 素性检验

本文最后更新于 2025年9月2日星期二 01:57

背景

Miller–Rabin 素性检验(Miller–Rabin primality test)是一种素数判定法则，利用随机化算法判断一个数是合数还是可能是素数。1976 年，卡内基梅隆大学的计算机系教授 Gary Lee Miller 首先提出了基于广义黎曼猜想的确定性算法，由于广义黎曼猜想并没有被证明，1980 年，以色列耶路撒冷希伯来大学的 Michael O. Rabin 教授作出修改，提出了不依赖于该假设的随机化算法。

前置知识

费马小定理

如果 \(p\) 是一个质数，而整数 \(a\) 不是 \(p\) 的倍数，则有 \(a^{p-1}\equiv 1 \pmod p\) 或 \(a^p\equiv a \pmod p\)。

二次探测定理

如果 \(p\) 是一个素数，那么 \(x^2\equiv 1 \pmod p\)的解只有两种可能，\(x = 1\) 或 \(x = p-1\)。

证明

\(x^2\equiv 1 \pmod p\) 可以转换为 \(x^2-1 \equiv 0\pmod p\)，即 \((x-1)(x+1) = 0 \pmod p\)，因此 \(p|(x-1)(x+1)\)，所以 \(x = 1\) 或 \(x = p-1\)。

Miller-Rabin 算法的证明

假设 \(n(n>2)\) 为素数，那么满足方程 \(x^2\equiv 1\pmod n\) 的解为 \(x = 1\) 或 \(x = n-1\)。

已知 \(n\) 为奇数，那么 \(n-1\) 一定为偶数，我们可以设 \(n-1 = m\cdot2^t\)，其中 \(m\) 为奇数。其实只要将 \(n-1\) 不断地除以 2，直到变为 \(m\)。

接下来我们取一个随机数 \(a\in[1, n)\)，且 \(a^{n-1}= a^{m\cdot2^t}\)，又由二次探测定理可知

\[ a^{m\cdot2^{t-1}}\equiv 1\pmod n \]

或

\[ a^{m\cdot2^{t-1}}\equiv n-1\pmod n \]

不断地用二次探测定理递归下去，直到最后 \(a^m\equiv 1\pmod n\) 或 \(a^m\equiv n-1\pmod n\)，中间只要有一次不满足二次探测定理，那么这个数就一定是合数，否则就有 3/4 的概率为素数。

Miller-Rabin 算法的正确率为 75%，所以要进行约 4~8 次来提高正确率。

相较于 Lucas–Lehmer 检验和费马素性检验，Miller-Rabin 素性检验的速度更快，误判率也更低。

算法流程

记要判断的数为 \(n\)。

\(n=2\) 返回 true
\(n<2\) 或者为偶数，返回 false
令 \(n-1 = m\cdot2^t\)，求出 \(m\) 和 \(t\)
生成随机数 \(a\in[1, n)\)
令 \(x = a^m\mod n\)
进行 \(t\) 次循环，每次循环中都进行 \(y =x^2\mod n\)，\(x = y\) 的操作，最终会变为 \(a^{n-1}\)，即 \(a^{m\cdot2^t}\)。运算过程中如果 \(x^2\mod n = 1\)，也就是 \(y = 1\)，如果 \(n\) 是素数，那么 \(x == 1\) 或者 \(x == n-1\)，否则 \(n\) 就一定不是素数，返回 false
整个循环结束后，程序运行到最后 \(x = a^{n-1}\), 根据费马小定理，如果 \(x\neq 1\)，那么肯定不是素数，返回 false

C++ 代码

ll Mul(ll a, ll b, ll p) {
  ll t = 0;
  for (; b; b >>= 1, a = (a << 1) % p) {
    if (b & 1) t = (t + a) % p;
  }
  return t;
}

ll Pow(ll a, ll b, ll p) {
  ll t = 1;
  for (; b; b >>= 1, a = a * a % p) {
    if (b & 1) t = t * a % p;
  }
  return t;
}

bool Miller_Rabin(ll n, int s) {
  if (n == 2) return true;
  if (n < 2 || !(n & 1)) return false;
  int t = 0;
  ll x, y, u = n - 1;
  while ((u & 1) == 0) ++t, u >>= 1;
  for (int i = 0; i < s; ++i) {
    ll a = rand() % (n - 1) + 1;
    ll x = Pow(a, u, n);
    for (int j = 0; j < t; ++j) {
      ll y = Mul(x, x, n);
      if (y == 1 && x != 1 && x != n - 1) return false;
      x = y;
    }
    if (x != 1) return false;
  }
  return true;
}