Miller-Rabin 素性检验和 Pollard-Rho 算法

Miller-Rabin 素性检验 ​

常规确定性素性检验的时间复杂度是 $O(\sqrt{n})$，但是对于 $2^{64}$ 规模的数据，该素性检验显得有些捉襟见肘了，于是我们引入 Miller-Rabin 素性检验法。该算法是一个随机化算法，能够在可接受的时间复杂度内，进行错误率在可接受范围内的素性检验。除此之外，可以证明，在一定规模内，选取合理的一组数据，可以实现确定性素性检验^[1]，在下面会进一步提到。

费马素性检验 ​

由费马小定理，可以得到：若 $p$ 是一个素数，则对于所有 $a\perp p$，有

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

该定理的逆否命题：若存在 $1\le a<n$，满足 $a^{n-1}\not\equiv 1\pmod{n}$，则 $p$ 一定不是素数。这是一个真命题（逆否命题的性质）。

借此启发我们随机找一些底数 $a$，计算 $a^{n-1}\bmod{n}$ 的结果，如果该数字不为 $1$，则可判定其为合数。

然而遗憾的是，有一些合数能够通过某些费马素性验证。例如当 $a=2,n=341$ 时，$2^{340}\equiv 1\pmod{341}$，然而 $341=11\times 31$ 是一个合数。此时 $n$ 称为 以 $a$ 为底的伪素数。

更为惊人的是，有一些满足对于所有 $a\perp n$，都有 $a^{n-1}\equiv 1\pmod{n}$ 成立，这样的数被称为 Carmichael 数。

所以我们需要对费马素性检验进行改进。

二次探测定理 ​

内容 ​

如果 $p$ 是奇素数，那么 $x^2\equiv 1\pmod{p}$ 的解为 $x\equiv 1\pmod{p}$ 和 $x\equiv -1\pmod{p}$。

证明 ​

由原式可以导出 $x^2-1\equiv0\pmod{p}$，运用平方差公式得到 $(x+1)(x-1)\equiv0\pmod{p}$，也就是说 $p\mid (x+1)(x-1)$。

由于 $p$ 是一个素数，所以 $p$ 的倍数必然含有素因子 $p$，所以 $x$ 必然满足 $x+1\equiv p\pmod{p}$ 和 $x-1\equiv p\pmod{p}$ 其中一个。容易验证两者都成立，证毕。

Miller-Rabin 素性测试 ​

结合费马素性检验和二次探测定理可以得到最终的素性检验。

我们称 $x^2\equiv 1\pmod{p}$ 的解 $x\equiv 1\pmod{p}$ 和 $x\equiv -1\pmod{p}$ 为其平凡解。

如果我们找到 $n$ 的一个非平凡解，则可以判定 $n$ 不是素数。

一个偶数的素性测试是容易的。

对于一个奇数 $n$，将 $n - 1$ 分解成 $n=2^k\times u$，其中 $u$ 是一个奇数。

对于一个随机的 $a$，先求 $a^{u}\bmod n$ 的值，然后对其平方 $k$ 次，若发现非平凡平方根时，该数字不是素数。如果最终结果不为 $1$，则该数字不是素数。

时间复杂度分析 ​

由于对 $n$ 的规模没有限制，我们并不能默认单次乘法的时间复杂度为 $O(1)$。常规情况下对两个小于 $n$ 的数单次乘法的时间复杂度为 $O(\log^{2} n)$。单次校验需要使用一次复杂度为 $O(\log^{3} n)$ 的快速幂算法，再进行 $O(\log n)$ 次乘法。故单次校验总时间复杂度为 $O(\log^{3} n)$。

如果对乘法使用快速傅里叶变换算法进行优化，单次乘法操作复杂度降为 $O(\log n\log\log n)$，于是总复杂度变为 $O(\log^{2}n\log\log n)=\widetilde{O}(\log^{2}n)$。

假设使用朴素算法进行 $k$ 次判断，则时间复杂度为 $O(k\log^{3} n)$。

错误率分析 ​

该算法的错误率取决于有多少个底数，使得一个合数被错误地判定为素数。

而这样的底数占总数的比率不超过 $\frac{1}{4}$ ^[2]，所以单次错误率是 $\frac{1}{4}$，故进行 $k$ 轮校验得到的错误率是 $O(4^{-k})$ 的，当 $k=8$ 时错误率约为十万分之，这是非常可以接受的。

实现 ​

这个程序实现了对 long long 范围内的整数进行 miller-rabin 测试，测试 TEST_TIME 中设定的次数次。

std::mt19937_64 Rad(std::chrono::steady_clock::now().time_since_epoch().count());

using ll = long long;

ll qpow(ll a, ll b, ll p) {
  ll res = 1;
  for (; b; b >>= 1, a = (__int128)a * a % p)
    if (b & 1) res = (__int128)res * a % p;
  return res;
}

ll rad(ll l, ll r) {
  return std::uniform_int_distribution<ll>(l, r)(Rad);
}

const int TEST_TIME = 8;

bool miller_rabin(ll n) {
  if (n < 3) return n == 2;
  if (n == 3) return true; // 判掉一些边边角角的情况
  ll u = n - 1, t = 0;
  while (!(u & 1)) ++t, u >>= 1;
  for (int i = 1; i <= TEST_TIME; ++i) {
    ll a = rad(2, n - 2);
    ll v = qpow(a, u, n);
    if (v == 1) continue; // 如果 v 为 1，直接通过测试
    int s;
    for (s = 0; s < t; ++s) {
      if (v == n - 1) break; // 在某一次 v 为 n - 1 时，以后 v 都是 1，通过测试
      v = v * v % n;
    }
    if (s == t) return false; // 当 s == t 时，未通过费马素性检验
  }
  return true;
}

拓展-确定性判素 ​

假设广义 Riemann 猜想成立，那么进行确定性判素只需要选择 $[2,\lfloor2(\log{n})^{2}\rfloor]$ 为底，进行 Miller-Rabin 测试^[1:1]。

但是即使广义 Riemann 猜想成立，这样判定时间复杂度变为 $O(\log^{5}n)$，此时运行时间已然不可接受。

但是在算法竞赛中，通常我们只需要对 long long 范围内的素数进行判定。在一定范围内，我们可以选用若干个固定底数进行确定性判素。例如^[3]：

• 对于 $2^{32}$ 范围内的判素，选取 $2,7,61$ 三个数为底。

• 对于 $2^{64}$ 范围内的判素，选取 $2,325,9375,28178,450775,9780504,1795265022$ 七个数为底。

• 如果记不住上面七个数，可以使用 $2,3,5,7,11,13,17,19,23,29,31,37$ 十二个数（即前 12 个素数）为底，可以适用于 $2^{78}$ 之内的判素。

Pollard-Rho 算法 ​

分解质因数有 $O(\sqrt{n})$ 的朴素算法。

同样地，对于更大规模的数据，朴素算法显得力不从心了，我们引入 Pollard-Rho 算法，该算法可以在期望 $O(n^{\frac{1}{4}})$ 的时间复杂度内求得 $n$ 的其中一个 非平凡因子。

生日悖论 ​

假设每年都是 365 天，求一个房间中至少多少人，才能使其中两个人生日相同的概率达到 $50\%$ ?

设房间里有 $n(n\le 365)$ 个人，有两个人同一天生日的概率为 $P$，那么没有两个人同一天生日的概率可以这样计算：

$$1-P=\frac{365}{365}\times \frac{364}{365}\times\cdots\times\frac{365-n+1}{365}=\frac{365!}{365^{n}(365-n)!}$$

可以代入几个数字得到：

可以求得若有 23 人，即可使得两个人生日相同概率达到 $50\%$。对于一个 41 人的房间，存在两个人生日相同的概率高达 $90\%$。对于 57 人的房间，存在两个人生日相同的概率到了惊人的 $99\%$。如果有 100 人，那么 没有 两个人同一天生日的概率达到了千万分之三。

这样的数学事实与我们的观念差异很大，所以称之为悖论。

这启示我们，在一个区间内随机选取一些数字，其中含有相同数字的概率很高。

更具体地，当在 $n$ 个数字中随机选取大约 $\sqrt{2n\ln 2}$ 个数字^[4]时，有 $50\%$ 的概率选择到两个相同的数字。也就是选择大约 $O(\sqrt{n})$ 个数字，有 $50\%$ 的概率选择到两个相同的数字。

如果在 $N$ 个数字中选出两个不同的数字，期望需要进行 $O(\sqrt{N})$ 次选取。^[5]

最大公约数求约数 ​

如果我们直接在 $(1,n)$ 中随机选出一个数，这个数字是 $n$ 的约数的概率是较小的。

我们巧妙地利用最大公约数的性质：假设有一个数字 $x$ 满足 $1<\gcd(x,n)$，那么 $\gcd(x,n)$ 一定是 $n$ 的一个非平凡因子。

满足 $1<\gcd(x,n)$ 的数字是多的：每一个质因子的倍数都满足条件。这可以大大提高我们的成功概率。

然而直接找满足条件的 $x$ 概率也并不太高，设 $n$ 的最小质因子为 $p$，假如能找到 $i\not=j$ 且 $i\equiv j\pmod{p}$，那么 $\vert i-j\vert$ 是一个满足条件的 $x$。

注意到模 $p$ 剩余系的大小不超过 $p$，并且 $p$ 不超过 $\sqrt{n}$。根据生日悖论，在 $(1,n)$ 中期望随机 $O(n^{\frac{1}{4}})$ 个数字，就能找到 $i\not=j$ 且 $i\equiv j\pmod{p}$，利用 $\gcd(\vert i-j\vert,n)$ 就能找出一个非平凡约数。

然而如果我们真的随机 $O(n^{\frac{1}{4}})$ 个数字，两两进行求差，求 gcd 操作，时间复杂度为 $O(\sqrt{n}\log n)$，比朴素做法更慢。所以我们使用两种方法进行实现：

Floyd 判环 ​

下面将使用多项式函数 $f(x)=(x^{2}+c)\bmod n$ 生成随机序列，其中 $c$ 是一个随机参数，这个函数能够保证，对于 $x\equiv y\pmod{p}$，有 $f(x)\equiv f(y)\pmod{p}$，其中 $p$ 是 $n$ 的约数。

可以感受一下，这个性质保证了若 $x\equiv y\pmod{p}$，此时能够找到一个因子，则 $f(x)\equiv f(y)\pmod{p}$ 仍然能够找到一个因子。

将任意一个数代入数列的第一项，就可以生成一个随机数列。但是由于模 $n$ 的结果是有限的，所以数列必然出现循环节，也就是环。

考虑两个在操场上跑圈的人，当跑得快的人追上跑得慢的人时，至少遍历了整个环一圈。

设 $a=f(1),b=f(1)$，每一次 $a=f(a),b=f(f(b))$，也就是模拟 $b$ 是跑得快的人，而 $a$ 是跑得慢的人。每一次计算 $\gcd(\vert a-b\vert,n)$ 的值，如果不为 $1$，那么找到一个非平凡因子。如果 $a=b$，此时已经进入了一个环，说明这个环中找不到答案，退出重新随机 $c$。

至于到底因为是环上任意两个位置都不能找到答案，还是因为我们使用的方法只检测了部分的位置，导致我们找不到答案，我就不知道了。

下面给出 模板题 的示例代码，由于该做法常数较大，无法通过。

#include <chrono>
#include <cstdio>
#include <random>

std::mt19937_64 Rad(std::chrono::steady_clock::now().time_since_epoch().count());

using ll = long long;

ll rad(ll l, ll r) {
  return std::uniform_int_distribution<ll>(l, r)(Rad);
}

const int TEST_TIME = 8;

ll qpow(ll a, ll b, ll p) {
  ll res = 1;
  for (; b; b >>= 1, a = a * a % p)
    if (b & 1) res = res * a % p;
  return res;
}
bool miller_rabin(ll n) {
  if (n < 3) return n == 2;
  if (n == 3) return true;
  ll u = n - 1, t = 0;
  while (!(u & 1)) ++t, u >>= 1;
  for (int i = 1; i <= TEST_TIME; ++i) {
    ll a = rad(2, n - 2);
    ll v = qpow(a, u, n);
    if (v == 1) continue;
    int s;
    for (s = 0; s < t; ++s) {
      if (v == n - 1) break;
      v = v * v % n;
    }
    if (s == t) return false;
  }
  return true;
}
ll f(ll x, ll a, ll n) {
  return ((__int128)x * x + x + a) % n;
}
ll gcd(ll a, ll b) {
  return b ? gcd(b, a % b) : a;
}
ll pollard_rho(ll n) {
  ll a = rad(1, n);
  ll x1, x2;
  x1 = x2 = rad(0, n - 1);
  while (true) {
    x1 = f(x1, a, n);
    x2 = f(f(x2, a, n), a, n);
    if (x1 == x2) {
      x1 = x2 = rad(0, n - 1);
      a = rad(1, n);
      continue;
    }
    ll d = gcd(abs(x1 - x2), n);
    if (d > 1 && d < n) return d;
  }
}
ll ans;
void find(ll n) {
  if (miller_rabin(n)) {
    ans = std::max(ans, n);
    return;
  }
  ll p = pollard_rho(n);
  find(p);
  find(n / p);
}

int main() {
  int tt;
  scanf("%d", &tt);
  while (tt--) {
    ll n;
    scanf("%lld", &n);
    if (miller_rabin(n)) {
      puts("Prime");
      continue;
    }
    ans = 0;
    find(n);
    printf("%lld\n", ans);
  }
}

倍增优化 ​

单次调用 $\gcd$ 函数的复杂度为 $O(\log n)$，经常调用会导致算法低效。注意到将若干个 $\vert s-t\vert$ 乘起来，对 $n$ 取模。如果存在一个 $\gcd(\vert s-t\vert,n)\not=1$，那么乘积与 $n$ 的 $\gcd$ 也不为 $1$。

所以我们可以累积一些乘积之后再求 $\gcd$。例如选择 128。

接着使用倍增进行优化。枚举一个步数 $2^{i}$，固定起点 $s$，从起点一步一步走 $2^{i}$ 步，将差值累乘起来，每隔 $128$ 步计算一次 $\gcd$。再将起点更新为终点，倍增步数。

下面给出模板的代码：

#include <chrono>
#include <cstdio>
#include <random>

std::mt19937_64 Rad(std::chrono::steady_clock::now().time_since_epoch().count());

using ll = long long;

ll rad(ll l, ll r) {
  return std::uniform_int_distribution<ll>(l, r)(Rad);
}

const int TEST_TIME = 8;

ll qpow(ll a, ll b, ll p) {
  ll res = 1;
  for (; b; b >>= 1, a = (__int128)a * a % p)
    if (b & 1) res = (__int128)res * a % p;
  return res;
}
bool miller_rabin(ll n) {
  if (n < 3) return n == 2;
  if (n == 3) return true;
  ll u = n - 1, t = 0;
  while (!(u & 1)) ++t, u >>= 1;
  for (int i = 1; i <= TEST_TIME; ++i) {
    ll a = rad(2, n - 2);
    ll v = qpow(a, u, n);
    if (v == 1) continue;
    int s;
    for (s = 0; s < t; ++s) {
      if (v == n - 1) break;
      v = (__int128)v * v % n;
    }
    if (s == t) return false;
  }
  return true;
}
ll f(ll x, ll a, ll n) {
  return ((__int128)x * x + x + a) % n;
}
ll gcd(ll a, ll b) {
  return b ? gcd(b, a % b) : a;
}
ll pollard_rho(ll n) {
  ll a = rad(1, n);
  ll s = 0, t = 0, val = 1;
  for (ll goal = 1; ; goal <<= 1, s = t, val = 1) {
    for (ll step = 1; step <= goal; ++step) {
      t = f(t, a, n);
      val = (__int128)val * abs(t - s) % n;
      if (!val) return n;
      if (step % 128 == 0) {
        ll d = gcd(val, n);
        if (d > 1) return d;
      }
    }
    ll d = gcd(val, n);
    if (d > 1) return d;
  }
}
ll ans;
void find(ll n) {
  if (n == 1) return;
  if (miller_rabin(n)) {
    ans = std::max(ans, n);
    return;
  }
  ll p = pollard_rho(n);
  find(p);
  find(n / p);
}

int main() {
  int tt;
  scanf("%d", &tt);
  while (tt--) {
    ll n;
    scanf("%lld", &n);
    if (miller_rabin(n)) {
      puts("Prime");
      continue;
    }
    ans = 1;
    find(n);
    printf("%lld\n", ans);
  }
}

参考资料与链接 ​