Miller-Rabin 与 Pollard-Rho

历史记录

清除记录

猜你想搜

AcWing热点
App
登录/注册

Miller-Rabin 与 Pollard-Rho

作者：

天元之弈 , 2025-06-11 17:28:23 · 福建 , 所有人可见 , 阅读 16

2

1

更好阅

前言

本文记录个人对于 $\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$ 是困难的，所以要通过合理的随机化算法。

$\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)$，只需要改一改即可。

喵喵喵喵喵喵，全文完，喵喵喵喵喵喵。

0 评论

App 内打开