OI-wiki Pollard-Rho 算法完全指南从生日悖论到 O(N^{1/4}) 的大整数质因数分解【免费下载链接】OI-wiki:star2: Wiki of OI / ICPC for everyone. 某大型游戏线上攻略内含炫酷算术魔法项目地址: https://gitcode.com/GitHub_Trending/oi/OI-wiki给定一个正整数 $N$试快速找到它的一个非平凡因数即不等于 $1$ 和 $N$ 自身的真因数。朴素试除需要 $O(\sqrt N)$ 时间当 $N \ge 10^{18}$ 时完全不可行。Pollard-Rho 算法以随机化为武器用生日悖论作理论支撑能在 $O(\sqrt p) O(N^{1/4})$ 的期望时间内找到一个非平凡因子是 OI / ICPC 中处理 $10^{18}$ 量级大整数质因数分解的标准算法。本文以 OI-wiki 的 Pollard-Rho 文档为主体结合仓库内完整实现代码与配套测试数据从朴素算法出发逐步推导生日悖论、Floyd 判环、Brent 判环与倍增优化并给出可运行的完整解法帮助你从原理到实现彻底掌握这一算法。引入为什么需要比 $O(\sqrt N)$ 更快设 $N$ 的最小素因子为 $p$则朴素找因子的复杂度为 $O(p) O(\sqrt N)$。当 $N \ge 10^{18}$ 时即使每秒执行 $10^9$ 次除法$\sqrt N \approx 10^9$ 次遍历也需要约 1 秒甚至更久且随着数据规模增大迅速不可接受。此时一个自然想法是随机猜测直接随机猜一个数是不是 $N$ 的因数运气好可以 $O(1)$ 出解。但对 $N \ge 10^{18}$ 的数据成功概率仅 $\frac{1}{10^{18}}$期望猜测次数高达 $10^{18}$ 次显然不现实即便把猜测范围限制在 $[2, \sqrt N]$ 内成功率依然过低。我们需要一种方法来优化猜测Pollard-Rho 算法正是为此而生。朴素算法试除分解最简单的算法即从 $[2, \sqrt N]$ 遍历把能整除 $N$ 的素数收集起来 Ccpp vectorint breakdown(int N) { vectorint result; for (int i 2; i * i N; i) { if (N % i 0) { // 如果 i 能够整除 N说明 i 为 N 的一个质因子 while (N % i 0) N / i; result.push_back(i); } } if (N ! 1) { // 说明再经过操作之后 N 留下了一个素数 result.push_back(N); } return result; } Pythonpython def breakdown(N): result [] for i in range(2, int(sqrt(N)) 1): if N % i 0: # 如果 i 能够整除 N说明 i 为 N 的一个质因子 while N % i 0: N // i result.append(i) if N ! 1: # 说明再经过操作之后 N 留下了一个素数 result.append(N) return result 可以证明result中的所有元素即为 $N$ 的全体素因数。证明分四步N的变化过程当循环进行到i结束时刚执行完while (N % i 0) N / ii不再整除N且每次除去一个因子都保证N仍整除初始的 $N$。因此当循环进行到i时N是 $N$ 的一个因子且不被任何小于i的整数整除。result中的元素均为 $N$ 的因子能存入i的条件是N % i 0即i整除N而N又是 $N$ 的因子故i是 $N$ 的因子循环结束时若N不为 1 也会存入此时它同样是 $N$ 的因子。result中均为素数假设存在合数 $K \in \text{result}$则必有不超过 $\sqrt K$ 的 $i$ 整除 $K$。这样的 $K$ 不可能作为循环中的某个i存入——因为由第 1 点循环到 $K$ 时N不被任何小于 $K$ 的i整除也不可能在循环结束后加入——因为循环退出的条件是i * i N此时所有不超过 $\sqrt K$ 的i都已遍历完且这些i均不整除当前的N即 $K$。所有 $N$ 的素因子必出现在result中假设 $p$ 是 $N$ 的素因子但未出现在result中则 $p$ 不可能是循环中出现过的i。设退出循环前最后的i$ 为i则i p而退出后的N$ 不被之前的i整除故 $p$ 整除N于是最后的N大于 1 且必为素数即N p会在最后加入result与假设矛盾。值得指出如果事先打好素数表时间复杂度将从 $O(\sqrt N)$ 下降到 $O\left(\frac{\sqrt N}{\ln N}\right)$可前往 筛法 查阅打表信息。Pollard-Rho 算法算法思想概述暴力算法获得一个非平凡因子的复杂度为 $O(p) O(\sqrt N)$$p$ 是 $N$ 的最小素因子。Pollard-Rho 算法是一种随机化算法可以在 $O(\sqrt p) O(N^{1/4})$ 的期望复杂度内获得一个非平凡因子。注意非平凡因子不一定是素因子。其核心想法是对于随机自映射 $f: \mathbb Z_p \to \mathbb Z_p$从任何一点 $x_1$ 出发迭代计算 $x_n f(x_{n-1})$将在 $O(\sqrt p)$ 期望时间内进入循环。如果能找到 $x_i \equiv x_j \pmod p$则 $p$ 整除 $\gcd(|x_i - x_j|, N)$这一最大公约数就是 $N$ 的一个非平凡因子。理解进入循环的期望时间为 $O(\sqrt p)$可以从生日悖论中获得启发。生日悖论为什么 $O(\sqrt N)$ 次抽样就够不考虑出生年份假设每年 365 天一个房间中至少多少人才能使其中两个人生日相同的概率达到 $50%$答案是 23 人而这个反直觉的数学事实就是生日悖论。推导如下假设一年有 $n$ 天房间中有 $k$ 人每个人的生日均匀分布于 $n$ 天之中且相互独立。设 $k$ 个人生日互不相同为事件 $A$则$$ P(A)\prod_{i0}^{k-1}\frac{n-i}{n} $$至少两个人生日相同的概率为 $P(\overline A)1-P(A)$。令 $P(\overline A)\ge\frac{1}{2}$即$$ P(A)\prod_{i0}^{k-1}\frac{n-i}{n} \le \frac{1}{2} $$由不等式 $1x\le \mathrm{e}^x$ 可得$$ P(A) \le \prod_{i1}^{k-1}\exp\left({-\frac{i}{n}}\right)\exp \left({-\frac{k(k-1)}{2n}}\right) $$于是$$ \exp\left({-\dfrac{k(k-1)}{2n}}\right) \le \frac{1}{2}\implies P(A) \le \frac{1}{2} $$将 $n365$ 代入解得 $k \ge 23$。当 $k 56$、$n365$ 时出现两个人同一天生日的概率将大于 $99%$。一般地在一年有 $n$ 天的情况下房间中约有 $\frac{1}{2}(\sqrt{8n\ln 21}1) \approx \sqrt{2n\ln 2}$ 人时至少两人同日生日的概率约为 $50%$。进一步地随机均匀地选取一列生日首次获得重复生日所需人数的期望也是 $O(\sqrt n)$。设人数为 $X$则$$ E(X) \sum_{x1}^{n1}P(X\ge x1) \sum_{x0}^n\frac{n!}{(n-x)!n^x} \sqrt{\frac{\pi n}{2}}-\frac13o(1). $$这启发我们随机选取一列数字出现重复数字所需抽样规模的期望是 $O(\sqrt n)$ 的。这正是 Pollard-Rho 算法能以 $O(N^{1/4})$ 期望复杂度工作的根本原因——它只需要在模 $p$ 意义下随机抽样 $O(\sqrt p)$ 次就大概率能撞上两个模 $p$ 同余的数。利用最大公约数求出一个约数实际构建一列模 $p$ 的随机数列并不现实因为 $p$ 正是需要求的未知量。所以算法通过 $f(x)(x^2c)\bmod N$ 生成伪随机数序列 ${x_i}$随机取 $x_1$令 $x_2f(x_1),\ x_3f(x_2),\ \dots,\ x_if(x_{i-1})$其中 $c \in [1,N)$ 是随机选取的常数。该函数容易计算且往往能生成相当随机的序列但它并非完全随机。例如设 $n50,\ c6,\ x_11$$f(x)$ 生成的数据为$$ 1, 7, 5, 31, 17, 45, 31, 17, 45, 31,\dots $$数据在 $x_4$ 以后都在 $31,17,45$ 之间循环。将这些数按下图方式排列图像酷似一个希腊字母 $\rho$算法也因此得名 rho更重要的是这样的函数确实提供了 $\mathbb Z_p$ 上的一个自映射它满足关键性质若 $x \equiv y \pmod p$则 $f(x) \equiv f(y) \pmod p$。证明若 $x \equiv y \pmod p$则 $x^2c \equiv y^2c \pmod p$。注意到 $f(x)x^2c-k_xN$$k_x$ 是依赖于 $x$ 的整数且 $p \mid N$所以 $f(x) \equiv x^2c \pmod p$因而 $f(x) \equiv f(y) \pmod p$。作为 $\mathbb Z_p$ 上的伪随机自映射反复迭代${x_n \bmod p}$ 在 $O(\sqrt p)$ 的期望时间内就会出现重复。只要观察到这样的重复 $x_i \equiv x_j \pmod p$就可以根据 $\gcd(|x_i-x_j|, N)$ 求出 $N$ 的一个非平凡因子。由于 $p$ 未知无法直接判断重复的发生一个简单的判断方法正是检验 $\gcd(|x_i-x_j|, N)$ 是否严格大于 1。算法并非总能成功$\gcd(|x_i-x_j|, N)$ 可能等于 $N$即 $x_i \equiv x_j \pmod N$——此时 ${x_n \bmod p}$ 首次发生重复时恰好 ${x_n}$ 本身也发生重复了没有得到非平凡因子。而且 ${x_n}$ 一旦开始循环继续迭代只会重复这一循环没有意义。此时算法应输出分解失败更换 $f(x)$ 中选取的 $c$ 重新分解。理论上任何满足 $\forall x \equiv y \pmod p,\ f(x) \equiv f(y) \pmod p$ 且能保证一定伪随机性的函数 $f(x)$例如某些多项式函数都可以用在此处。实践中主要使用 $f(x)x^2c$其中 $c \neq 0, -2$。实现把判等转化为判环需要在迭代过程中快速判断 ${x_n \bmod p}$ 是否已经出现重复。将 $f$ 看成以 $\mathbb Z_p$ 为顶点的有向图上的边实际要实现的其实是一个判环算法只是把判等改成了判断 $\gcd(|x_i-x_j|, N)$ 是否大于 1。Floyd 判环设想两个人在赛跑A 速度快B 速度慢经过一定时间后 A 一定会和 B 相遇且相遇时 A 跑过的总距离减去 B 跑过的总距离一定是圈长的倍数。设 $af(0),\ bf(f(0))$每一次更新 $af(a),\ bf(f(b))$只要检查更新过程中 $a$ 和 $b$ 是否相等即可——相等说明出现了环。每次令 $d\gcd(|x_i-x_j|, N)$判断 $d$ 是否满足 $1dN$满足则直接返回 $d$若 $dN$说明 ${x_i}$ 已形成环此时不能继续操作直接返回 $N$ 本身并在后续操作里调整随机常数 $c$ 重新分解??? note 基于 Floyd 判环的 Pollard-Rho 算法 C cpp ll Pollard_Rho(ll N) { if (N 4) return 2; // 因为一开始跳了两步所以需要特判一下 4 ll c rand() % (N - 1) 1; ll t f(0, c, N); ll r f(f(0, c, N), c, N); while (t ! r) { ll d gcd(abs(t - r), N); if (d 1) return d; t f(t, c, N); r f(f(r, c, N), c, N); } return N; } Python python import random def Pollard_Rho(N): if N 4: return 2 # 因为一开始跳了两步所以需要特判一下 4 c random.randint(1, N - 1) t f(0, c, N) r f(f(0, c, N), c, N) while t ! r: d gcd(abs(t - r), N) if d 1: return d t f(t, c, N) r f(f(r, c, N), c, N) return N 注意代码中N 4的特判由于一开始跳了两步t和r分别从 $f(0)$、$f(f(0))$ 出发对 $N4$ 会直接进入死循环需要单独处理。Brent 判环Floyd 判环算法在常数上还可以改进。Brent 判环从 $k1$ 开始递增 $k$在第 $k$ 轮让 A 等在原地B 向前移动 $2^k$ 步如果过程中 B 遇到了 A则说明已经得到环否则让 A 瞬移到 B 的位置继续下一轮。可以证明这样得到环之前需要调用 $f$ 的次数永远不大于 Floyd 判环算法原论文中的测试表明Brent 判环需要的平均时间相较于 Floyd 判环减少了 $24%$。倍增优化减少 gcd 调用次数无论是 Floyd 判环还是 Brent 判环迭代次数都是 $O(\sqrt p)$但每次迭代都用 $\gcd$ 判断是否成环会拖慢运行速度。可以通过乘法累积来减少求 $\gcd$ 的次数。核心性质如果 $\gcd(a, N) 1$那么 $\gcd(ab \bmod N, N) \gcd(ab, N) 1$ 对任意 $b \in \mathbb N_$ 都成立。也就是说如果计算得到 $\gcd\left(\prod |x_i-x_j| \bmod N, N\right) 1$那么必然存在其中一对 $(x_i, x_j)$ 满足 $\gcd(|x_i-x_j|, N) 1$。如果该乘积在某一时刻得到 0则分解失败退出并返回 $N$ 本身。如果每 $k$ 对计算一次 $\gcd$算法复杂度降低到 $O\left(\sqrt p k^{-1}\sqrt p \log N\right)$其中 $\log N$ 是单次计算 $\gcd$ 的开销。当 $k$ 和 $\log N$ 大致同阶时可以得到 $O(\sqrt p)$ 的期望复杂度。具体实现中大多选取 $k128$下述实现取 $127$ 步一判。以下给出 Brent 判环 倍增优化的 Pollard-Rho 算法实现??? note 实现 C cpp ll Pollard_Rho(ll x) { ll t 0; ll c rand() % (x - 1) 1; ll s t; int step 0, goal 1; ll val 1; for (goal 1;; goal 1, s t, val 1) { for (step 1; step goal; step) { t f(t, c, x); val val * abs(t - s) % x; // 如果 val 为 0退出重新分解 if (!val) return x; if (step % 127 0) { ll d gcd(val, x); if (d 1) return d; } } ll d gcd(val, x); if (d 1) return d; } } Python python from random import randint from math import gcd def Pollard_Rho(x): c randint(1, x - 1) s t f(0, c, x) goal val 1 while True: for step in range(1, goal 1): t f(t, c, x) val val * abs(t - s) % x if val 0: return x # 如果 val 为 0退出重新分解 if step % 127 0: d gcd(val, x) if d 1: return d d gcd(val, x) if d 1: return d s t goal 1 val 1 实现要点逐行解析s t是当前轮的锚点t不断前进val累积 $\prod |t-s| \bmod x$每goal步Brent 判环的 $2^k$ 段结束后s瞬移到t的位置goal翻倍val重置为 1每 127 步取一次 $\gcd(val, x)$避免每一步都做昂贵的 $\gcd$val一旦为 0模 $x$ 意义下乘积为零说明出现退化情形直接返回x让外层重试。复杂度分析Pollard-Rho 算法中的期望迭代次数为 $O(\sqrt p)$$p$ 是 $N$ 的最小素因子。无论采用 Floyd 判环还是 Brent 判环如果不使用倍增优化期望复杂度都是 $O(\sqrt p \log N)$加上倍增优化后可以近似得到 $O(\sqrt p)$ 的期望复杂度。值得一提的是上述分析基于完全随机的自映射函数而 Pollard-Rho 算法实际使用的是伪随机函数所以该算法并没有严格的复杂度分析实践中通常跑得较快。完整实战求一个数的最大素因子以 P4718【模板】Pollard-Rho 算法 为例该题数据规模庞大Floyd 判环方法不够用必须采用倍增优化。整体思路对 $n$ 先用 Miller-Rabin 素性测试判断是否为素数是则直接返回否则用 Pollard-Rho 找一个因子 $p$将 $n$ 除去因子 $p$再递归分解 $n$ 和 $p$用 Miller-Rabin 判断是否出现质因子并用max_factor更新求出最大质因子。仓库中提供了完整可运行的 C 实现 docs/math/code/pollard-rho/pollard-rho_1.cpp配套测试数据在 pollard-rho_1.in 与 pollard-rho_1.ans。#include algorithm #include cstdlib #include ctime #include iostream using namespace std; using ll long long; using ull unsigned long long; int t; ll max_factor, n; ll gcd(ll a, ll b) { if (b 0) return a; return gcd(b, a % b); } ll bmul(ll a, ll b, ll m) { // 快速乘 ull c (ull)a * (ull)b - (ull)((long double)a / m * b 0.5L) * (ull)m; if (c (ull)m) return c; return c m; } ll qpow(ll x, ll p, ll mod) { // 快速幂 ll ans 1; while (p) { if (p 1) ans bmul(ans, x, mod); x bmul(x, x, mod); p 1; } return ans; } bool Miller_Rabin(ll p) { // 判断素数 if (p 2) return false; if (p 2) return true; if (p 3) return true; ll d p - 1, r 0; while (!(d 1)) r, d 1; // 将d处理为奇数 for (ll k 0; k 10; k) { ll a rand() % (p - 2) 2; ll x qpow(a, d, p); if (x 1 || x p - 1) continue; for (int i 0; i r - 1; i) { x bmul(x, x, p); if (x p - 1) break; } if (x ! p - 1) return false; } return true; } ll Pollard_Rho(ll x) { ll s 0, t 0; ll c (ll)rand() % (x - 1) 1; int step 0, goal 1; ll val 1; for (goal 1;; goal * 2, s t, val 1) { // 倍增优化 for (step 1; step goal; step) { t (bmul(t, t, x) c) % x; val bmul(val, abs(t - s), x); if ((step % 127) 0) { ll d gcd(val, x); if (d 1) return d; } } ll d gcd(val, x); if (d 1) return d; } } void fac(ll x) { if (x max_factor || x 2) return; if (Miller_Rabin(x)) { // 如果x为质数 max_factor max(max_factor, x); // 更新答案 return; } ll p x; while (p x) p Pollard_Rho(x); // 使用该算法 while ((x % p) 0) x / p; fac(x), fac(p); // 继续向下分解x和p } int main() { cin t; while (t--) { srand((unsigned)time(NULL)); max_factor 0; cin n; fac(n); if (max_factor n) // 最大的质因数即自己 cout Prime\n; else cout max_factor \n; } return 0; }配套测试数据验证仓库 docs/math/examples/pollard-rho/ 提供了 6 组输入输出样例可用来直接验证上述实现的正确性输入 $n$期望输出说明2Prime素数直接判定13Prime素数直接判定13467$134 2 \times 67$889741合数最大素因子为 4112345676543214649超过 32 位的大合数10000000000005$10^{12} 2^{12} \times 5^{12}$最大素因子为 5其中1234567654321和1000000000000均已超出 32 位整数范围直观展示了算法处理 $10^{12}$ 量级数据的能力配合long long与快速乘bmul可进一步扩展至 $10^{18}$ 量级。关键实现细节快速乘bmul计算 $a \times b \bmod m$ 时直接乘法可能溢出 64 位。仓库实现使用long double预估商再校正的做法避免__int128依赖在评测环境受限时依然可用。Miller-Rabin 素性测试将 $p-1$ 分解为 $d \times 2^r$随机取 10 个底数 $a$ 检验通过二次探测与费马小定理判定素数。其完整原理可参考 Miller-Rabin 素性测试。fac的递归结构先剔除因子 $p$while (x % p 0) x / p再对 $x$ 和 $p$ 分别递归if (x max_factor)剪枝避免无效递归。随机常数cPollard_Rho每次被调用都会重新随机 $c$一旦返回x分解失败就在外层while (p x)中重新尝试直到获得真正的因子 $p x$。参考资料与链接生日悖论Reverse problemhttps://en.wikipedia.org/wiki/Birthday_problem#Reverse_problemMenezes, Alfred J.; van Oorschot, Paul C.; Vanstone, Scott A. (2001).Handbook of Applied Cryptography. Section 3.11 and 3.12.Brent, R. P. (1980),An improved Monte Carlo factorization algorithm, BIT Numerical Mathematics, 20(2): 176–184.【免费下载链接】OI-wiki:star2: Wiki of OI / ICPC for everyone. 某大型游戏线上攻略内含炫酷算术魔法项目地址: https://gitcode.com/GitHub_Trending/oi/OI-wiki创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考