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

Tangninghaha · 2025-01-10 08:09:32 · 算法·理论

Miller-Rabin 素性检验

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

费马素性检验

由费马小定理，可以得到：若 p 是一个素数，则对于所有 a\perp 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} [^3]，所以单次错误率是 \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 测试[^2]。

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

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

• 对于 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，那么没有两个人同一天生日的概率可以这样计算：

可以代入几个数字得到：

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

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

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

更具体地，当在 n 个数字中随机选取大约 \sqrt{2n\ln 2} 个数字[^1]时，有 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=x'p+r,y=y'p+r，那么 f(x)=(x'^{2}p^{2}+r^{2}+2x'pr+c)\bmod{n}=x'^{2}p^{2}+r^{2}+2x'pr+c-kn，同理得到 f(y)=y'^{2}p^{2}+r^{2}+2y'pr+c-kn。

实际上任何多项式函数都满足这个性质。但是我们希望项数比较少以更快地得到结果。此外，经过本人测试，如果直接使用 f(x)=(x^{2}+c)\bmod n 在对 4 进行分解时会陷入死循环。所以下文实际使用的是 f(x)=(x^{2}+x+c)\bmod 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);
  }
}

参考资料与链接

[^1]: https://oi-wiki.org/math/number-theory/pollard-rho/#%E7%94%9F%E6%97%A5%E6%82%96%E8%AE%BA 分解质因数 - OI wiki

[^2]: https://www.zhihu.com/question/308322307/answer/574767625 如何编程判断一个数是否是质数？ - Wiener ES的回答 - 知乎

[^3]: https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Accuracy Miller–Rabin primality test - Wikipedia

[^4]: https://www.luogu.com.cn/article/jr27vtto 题解 P4718/论 Miller-Rabin 算法的确定性化 - 洛谷专栏

[^5]: https://en.wikipedia.org/wiki/Birthday_problem#Probability_of_a_shared_birthday_\(collision) Birthday problem - wikipedia