Miller-Rabin 与 Pollard-Rho
前言
本文记录个人对于 \(\mathcal{Miller-Rabin}\) 和 \(\mathcal{Pollard-Rho}\) 的一些理解，一些证明可能不是很完美，欢迎指出。
喵呜喵呜喵喵喵~
Miller-Rabin
二次探测定理
如果 \(p\) 是奇素数，则 \(x^2 \equiv 1 \pmod p\) 的解为 \(x \equiv 1 \pmod p\) 或 \(x \equiv p - 1 \pmod p\)。
证明是简单的，移项后用平方差公式展开即可。
算法
一种用 \(log(n)\) 的时间快速判断一个数是否为质数的算法。
假设我们已经排除了所有的偶数和 \(n < 3\) 的情况，那么对于要判定的 \(n\) 来说 \(n - 1\) 一定是偶数。一种直观的方法是根据费马小定理 \(a^{p - 1} \equiv 1 \pmod p\) 的逆定理来判定素数，但这么做会被 \(\mathcal{Carmichael}\) 数给卡掉（比如 \(561\)）。
接着介绍一个定理：
对于质数 \(p \geq 3\)，设 \(p - 1 = d \times 2^r\)（\(d\) 为奇数），对于任意 \(a \leq p\)，以下两个命题至少有一个成立：
- \(a^d \equiv 1 \pmod p\)
- \(\exists 0 \leq i < r\) 使得 \(a^{d \times 2^i} \equiv p - 1 \pmod p\)
证明：
设 \(s_i = a^{d \times 2 ^ i}\)，则有 \(s_{i} = s_{i - 1}^2\)，由于质数一定满足费马小定理，所以一定有 \(s_r \equiv 1 \pmod p\)。
那么假设最小的满足 \(s_i \equiv 1 \pmod p\) 的下标为 \(i\)，那么分类讨论：
- \(i = 0\)，则 \(a^d \equiv 1 \pmod p\) 成立
- $i \ge1 $，那么有 \(s_i \equiv 1 \pmod p\)，所以 \(s_{i - 1}^2 \equiv 1 \pmod p\)。根据二次探测定理，就有 \(s_{i - 1} \equiv p - 1 \pmod p\)。
综上，原定理成立。
知道这么个定理后 \(\mathcal{Miller-Rabin}\) 算法就很简单了。我们随机几个底数 \(a\)，对于每个底数 \(a\) 判断一下是否通过检验。
需要注意的是这个定理同样不能直接判定一个数是质数，只能保证概率正确。但对于 long long 范围内的整数，取底数为 \(\{2, 325, 9375, 28178, 450775, 9780504, 1795265022\}\) 可以保证 \(100\%\) 的正确率。或者选 \(\{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37\}\) 也同样可行。
模板题 代码：
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define rep(i, a, b) for (int i = (a); i <= (b); i ++)
#define fro(i, a, b) for (int i = (a); i >= b; i --)
#define INF 0x3f3f3f3f
#define eps 1e-6
#define lowbit(x) (x & (-x))
#define reg register
#define IL inline
typedef long long LL;
typedef std::pair<int, int> PII;
inline int read() {
int x = 0, f = 1;
char ch = getchar();
while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
while (ch >= '0' && ch <= '9') { x = (x << 1) + (x << 3) + (ch ^ 48); ch = getchar(); }
return x * f;
}
int primes[8] = {2, 3, 5, 7, 11, 13, 17, 37};
int T, n;
int mul(int a, int b, int p) {
int res = 0;
while (b) {
if (b & 1) res = (res + a) % p;
b >>= 1; a = (a + a) % p;
}
return res;
}
int qmi(int a, int b, int p) {
int res = 1;
while (b) {
if (b & 1) res = mul(res, a, p);
b >>= 1; a = mul(a, a, p);
}
return res;
}
bool Miller_Rabin(int a, int n) {
if (n < 3 || n % 2 == 0) return n == 2;
int d = n - 1, r = 0;
while (!(d & 1)) d >>= 1, r ++;
int t = qmi(a, d, n); if (t == 1) return true;
for (int i = 0; i < r; i ++) {
if (t == n - 1) return true;
t = mul(t, t, n);
}
return false;
}
bool isprimes(int n) {
if (n <= 1) return false;
for (int i = 0; i < 8; i ++) {
if (n == primes[i]) return true;
if (n % primes[i] == 0) return false;
if (!Miller_Rabin(primes[i], n)) return false;
}
return true;
}
signed main() {
int T = read();
while (T --) {
int n = read();
puts(isprimes(n) ? "YES" : "NO");
}
return 0;
}
Pollard-Rho
生日悖论
随机选取一列生日，首次获得重复生日需要的人数的期望是 \(\sqrt \frac{\pi n}{2}\)，即 \(O(\sqrt n)\)。
证明省略，可以自行查阅 oi-wiki。
算法
\(\mathcal{Pollard-Rho}\) 算法用来快速找出一个整数 \(n\) 的一个非 \(1, n\) 的因数。
首先先介绍一种找因数的方法，如果存在 \(d = \gcd(k, n) > 1\) 那么 \(d\) 显然是 \(n\) 的因数，但直接找 \(k\) 是困难的，所以要通过合理的随机化算法。
我们考虑找到 \(x_i \equiv x_j \pmod p\)，那么就有 \(p | \gcd(|x_i - x_j|, n)\)，如果 \(\gcd(|x_i - x_j|, n) > 1\) 那么 \(p\) 就是 \(n\) 的一个非平凡因子。
\(\mathcal{Pollard-Rho}\) 算法通过构造函数 \(F(x) = (x ^ 2 + c) \bmod n\) 来进行伪随机。通过一些观察会发现，这个函数会在某一时刻进入循环。原因是每个 \(0 \sim n - 1\) 的点会向其它点连边，最后就会形成一个形状酷似 \(\rho\) 的基环树。
同时我们还有一个很好的性质，如果 \(x \equiv y \pmod p\)，则 \(f(x) \equiv f(y) \pmod p\)，将式子展开即可证明。那么如果 \(x_i\) 和 \(x_j\) 在环上，对 \(|x_i - x_j|\) 进行一次判断往后的 \(x_{i + k}, x_{j + k}\) 就都不需要判断了。
所以我们考虑在找环的同时每次尽可能选 \(|x_i - x_j| \bmod p\) 不一样的点对进行判断。
那么有快慢指针的做法，即指针 \(x\) 每次前进一步，\(y\) 每次前进两步，如此既能够找环也可以保证每一步距离不一样。
注意如果出现了 \(|x_i - x_j| = 0\) 的情况我们需要重新选取 \(c\) 再做一遍。
IL int Pollard_Rho(int n) {
static mt19937_64 sj(chrono::steady_clock::now().time_since_epoch().count());
uniform_int_distribution<int> u0(1, n - 1);
int c = u0(sj);
auto f = [&](int x) { return (mul(x, x, n) + c) % n; };
int x = f(0), y = f(x);
while (x != y) {
int d = gcd(abs(x - y), n);
if (d != 1) return d;
x = f(x), y = f(f(y));
}
return n;
}
于是你得到了朴素的做法，这么做根据生日悖论找到重复的 \(\{x_n \bmod p\}\) 的期望复杂度为 \(O(\sqrt p)\)，其中 \(p\) 是最小的质因子，所以也就是 \(O(n^{\frac{1}{4}})\)。同时由于你做了很多 \(\gcd\) 操作还要再乘一个 \(log\)，并且常数很大。
考虑少做一些 \(\gcd\) 操作，这里使用倍增累积法，把 \(2^k\) 个需要判断的 \(|x - y|\) 都乘起来，然后同一计算。并且前人的经验告诉我们，如果累积的次数超过了 \(127\) 就单独计算一次。
IL int Pollard_Rho(int n) {
static mt19937_64 sj(chrono::steady_clock::now().time_since_epoch().count());
uniform_int_distribution<int> u0(1, n - 1);
int c = u0(sj);
auto f = [&](int x) { return (mul(x, x, n) + c) % n; };
int x = 0, y = 0, s = 1;
for (reg int k = 1; ; k <<= 1, y = x, s = 1) {
for (reg int i = 1; i <= k; i ++) {
x = f(x); s = mul(s, abs(x - y), n);
if (i % 127 == 0) {
int d = gcd(s, n);
if (d > 1) return d;
}
}
int d = gcd(s, n);
if (d > 1) return d;
}
return n;
}
应用
到这里你就已经学会求一个因子了，但如果要求分解质因数怎么办？
那么考虑每次找出 \(n\) 的因子 \(x\)，如果是质数就直接加入答案集合（使用 \(\mathcal{Miller-Rabin}\) 算法），否则递归处理 \(x\) 和 \(\frac{n}{x}\) 即可。
以 P4718 【模板】Pollard-Rho 为例，要求找出最大的质因子，给出代码。
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define rep(i, a, b) for (int i = (a); i <= (b); i ++)
#define fro(i, a, b) for (int i = (a); i >= b; i --)
#define INF 0x3f3f3f3f
#define eps 1e-6
#define lowbit(x) (x & (-x))
#define reg register
#define IL inline
typedef long long LL;
typedef std::pair<int, int> PII;
inline int read() {
int x = 0, f = 1;
char ch = getchar();
while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
while (ch >= '0' && ch <= '9') { x = (x << 1) + (x << 3) + (ch ^ 48); ch = getchar(); }
return x * f;
}
int primes[12] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
int ans;
IL int gcd(int a, int b) { return b ? gcd(b, a % b) : a; }
IL int mul(int a, int b, int p) { return (__int128)a * b % p; }
IL int qmi(int a, int b, int p) {
int res = 1;
while (b) {
if (b & 1) res = mul(res, a, p);
b >>= 1; a = mul(a, a, p);
}
return res;
}
IL bool Miller_Rabin(int a, int n) {
if (n < 3 || n % 2 == 0) return n == 2;
int d = n - 1, r = 0;
while (!(d & 1)) d >>= 1, r ++;
int t = qmi(a, d, n); if (t == 1) return true;
for (reg int i = 0; i < r; i ++) {
if (t == n - 1) return true;
t = mul(t, t, n);
}
return false;
}
IL bool isprimes(int n) {
if (n < 3 || n % 2 == 0) return n == 2;
for (reg int i = 0; i < 12; i ++) {
if (n == primes[i]) return true;
if (n % primes[i] == 0) return false;
if (!Miller_Rabin(primes[i], n)) return false;
}
return true;
}
IL int Pollard_Rho(int n) {
static mt19937_64 sj(chrono::steady_clock::now().time_since_epoch().count());
uniform_int_distribution<int> u0(1, n - 1);
int c = u0(sj);
auto f = [&](int x) { return (mul(x, x, n) + c) % n; };
int x = 0, y = 0, s = 1;
for (reg int k = 1; ; k <<= 1, y = x, s = 1) {
for (reg int i = 1; i <= k; i ++) {
x = f(x); s = mul(s, abs(x - y), n);
if (i % 127 == 0) {
int d = gcd(s, n);
if (d > 1) return d;
}
}
int d = gcd(s, n);
if (d > 1) return d;
}
return n;
}
IL void work(int n) {
if (n <= 1) return;
if (n == 4) { ans = max(ans, 2ll); return; }
if (isprimes(n)) { ans = max(ans, n); return; }
int x = n;
while (x == n) x = Pollard_Rho(n);
work(x); work(n / x);
}
void solve() {
ans = 0; int n = read();
if (isprimes(n)) puts("Prime");
else {
work(n);
printf("%lld\n", ans);
}
}
signed main() {
int T = read();
while (T --) solve();
return 0;
}
你甚至可以快速求出 \(\varphi(n)\)，只需要改一改即可。
浙公网安备 33010602011771号