说实话，第一眼看到 P4718 这个模板题，我是有点不屑的——Pollard-Rho 也好，Miller-Rabin 也好，板子不是一抓一大把？结果真上手写，WA 和 TLE 来回横跳了一下午，最后发现坑全藏在乘法溢出、gcd 调用频率、随机数取值这些“边角料”里。网上很多板子直接抄过来不一定能过，因为数据范围摆在那里：单个 n 能到 4×10^18 级别，任何二次方乘法都可能爆掉 long long。这篇就按我实际调试的顺序，把 Miller-Rabin 和 Pollard-Rho 的原理、优化、完整可提交模板一次讲清楚。
如果你是正在刷这题的选手，或者只是想把“会背板子”升级成“知道板子在干什么”，这篇应该都能帮到你。我会尽量把每个判断、每个优化背后的理由说透，而不是丢一段代码让你背。
1. 题目到底在考什么：Miller-Rabin与Pollard-Rho的分工
1.1 从“最大质因子”这句话读出什么
题意很直接：给你若干个可能高达 4×10^18 的数，对每个数输出它的最大质因子。注意不是“质因子个数”，也不是“是不是质数”，而是要把这个数真正拆开，找到所有质因子，然后取最大值。
暴力思路到 sqrt(n) 试除，理论上要跑到 2×10^9 次循环，绝对不可行。筛法只对 1e7 以内的数友好，在大数面前形同虚设。所以这条路从一开始就断了。
题目没有直接说“请实现 Pollard-Rho”，但所有标程都指向同一个组合：先判断这个数是不是质数，如果不是，就把它拆成两个因子，然后对每个因子继续判断、继续拆，直到所有因子都是质数。整个流程其实就是“素性判定 + 大数分解”的嵌套。
1.2 两个算法怎么配合
Miller-Rabin（下面简称 MR）负责回答“n 是不是素数”，答案只有 true/false。Pollard-Rho（下面简称 PR）负责回答“如果 n 是合数，给我一个非平凡因子”，这个因子可能是质数，也可能是合数。
这两个算法是上下级关系，不是并列关系。流程是这样的：
factor(n): 如果 n == 1，返回 如果 MR(n) 判定为质数：记录 n，返回 否则： d = PR(n) // 拿到一个因子 factor(d) // 递归分解 factor(n / d) // 递归分解另一半为什么 PR 返回的因子还要继续递归？因为 PR 只保证给你一个“非平凡因子”，不保证它是质数。比如 n=1009×1013，PR 第一次很可能直接命中 1009 或 1013 这种质数，但也可能返回 1009×1013 本身之外的某个合数因子。只有 MR 能最终裁决谁是质数。
第 6 章我会给出完整的 factor 代码，现在先别急，我们把 MR 和 PR 各自内部的门道拆开。
2. Miller-Rabin：二次探测是加锁，不是加戏
2.1 费马小定理会漏掉一类合数
MR 的起点是费马小定理：如果 n 是素数，那么对任意与 n 互质的 a，都有
a^(n-1) ≡ 1 (mod n)反过来用这个性质测试：随机取一些 a，算 a^(n-1) mod n，如果结果不是 1，那 n 肯定是合数；如果都是 1，那 n 大概率是素数。
但“大概率”三个字就是问题所在。存在一类合数叫卡迈克尔数，比如最小的 561 = 3×11×17，它对任意与 n 互质的 a 都满足 a^(n-1) ≡ 1 mod n。你要是只测费马小定理，561 会被当成素数放过去，而它显然不是。所以 MR 不能只用费马小定理这一道锁。
2.2 二次探测：非平凡平方根现场抓包
MR 真正加锁的地方是二次探测。核心引理很简单：如果 n 是奇素数，那么方程
x^2 ≡ 1 (mod n)的解只有 x ≡ 1 或 x ≡ -1（也就是 n-1）。这是模素数下多项式方程根数不超过次数的直接结论。
怎么利用它？先把 n-1 写成
n - 1 = d * 2^s其中 d 是奇数。然后从底数 a 开始，依次计算：
a^d, a^(2d), a^(4d), ..., a^(2^(s-1)d)每一项都是前一项的平方。对真正的素数来说，最终结果必须等于 1（费马小定理保证 a^(n-1) ≡ 1）。而在这个逐步平方的过程中，序列一旦出现“前一项不是 ±1，平方后变成 1”，就意味着我们找到了 1 的非平凡平方根，n 必为合数。
所以 MR 的判定逻辑是：
- 先算 x = a^d mod n；
- 如果 x == 1 或 x == n-1，这一轮直接通过；
- 否则连续平方最多 s-1 次，每次平方后检查是否等于 n-1，如果等于，通过；
- 如果一直没见到 n-1，判定合数。
为什么 x == 1 也直接通过？因为 x 已经是 1 了，后面再怎么平方都是 1，拿不到额外信息，只能视为“没发现破绽”。
2.3 写MR时的三个细节
细节一：拆分 n-1 时，d 必须保持为奇数，s 统计的是因子 2 的个数。写成while ((d & 1) == 0) d >>= 1, s++;就行。
细节二：基底选择。64 位整数范围内，用下面这组基底可以做到确定性判定：
{2, 325, 9375, 28178, 450775, 9780504, 1795265022}这组数是前人用程序暴力搜出来的“最小测试集”，覆盖 2^64 以内所有数，不会误判。另一个常见做法是直接用前 12 个小质数 {2,3,5,7,11,13,17,19,23,29,31,37}，同样有覆盖 64 位范围的结论。我习惯用魔数那组，因为只需要测 7 轮，常数小一些。
细节三：如果基底 a 恰好是 n 的倍数，即a % n == 0，这一轮要跳过。因为此时幂运算结果恒为 0，没有任何判别意义。
下面是最小可用的 MR 实现：
bool isPrime(ll n) { if (n < 2) return false; static const ll smallPrimes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}; for (ll p : smallPrimes) { if (n == p) return true; if (n % p == 0) return false; } ll d = n - 1, s = 0; while ((d & 1) == 0) d >>= 1, s++; static const ll bases[] = {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL}; for (ll a : bases) { if (a % n == 0) continue; ll x = qpow(a, d, n); if (x == 1 || x == n - 1) continue; bool ok = false; for (int r = 0; r < s - 1; r++) { x = mul_mod(x, x, n); if (x == n - 1) { ok = true; break; } } if (!ok) return false; } return true; }这里的qpow和mul_mod都还没写，第 5 章专门讲。先记着：所有乘法都过 mul_mod，不允许直接写a*b%mod。
3. Pollard-Rho原理：生日悖论怎样把因子撞出来
3.1 随机撒点：因子在碰撞里现身
Pollard-Rho 的数学直觉来自生日悖论。一个班里 365 个生日，只要大约 23 人，就有超过一半概率出现两个人生日相同——因为比较的是任意两人，而不是“有没有人和我生日相同”。
把同样的逻辑搬到模运算里。假设合数 n 有一个不太大的素因子 p，那么在区间 [0, p-1] 里随机撒点，大约撒 O(√p) 个点，就会出现重复。一旦两个随机点 x_i、x_j 在模 p 意义下相等，就说明它们的差 |x_i - x_j| 是 p 的倍数。这个差很可能不是 n 的倍数，于是
gcd(|x_i - x_j|, n)就会落在 1 和 n 之间——恭喜，你拿到了 n 的一个非平凡因子。
这就是 Pollard-Rho 的全部核心：不需要知道 p 是多少，只需要构造出“两个点在模 p 下撞车”。
3.2 f函数与Floyd判圈：不用存历史的随机源
问题是，如果完全用真随机数发生器，每生成一个新点都要和之前所有点比较，存储量会爆炸。Pollard-Rho 的巧妙之处是引入一个确定性的伪随机迭代函数：
f(x) = (x^2 + c) mod n给定一个起点 x0 和常数 c，序列 x0, f(x0), f(f(x0)), ... 完全确定。因为状态空间是有限的，这个序列最终一定会进入循环。而一旦两个点在模 p 下撞车，在循环结构上就表现为“轨迹形成了环”，这就是 Rho 这个名字的来源——轨迹长得像希腊字母 ρ。
为了不记录历史状态，用 Floyd 判圈：一个慢指针 x 每轮走一步，一个快指针 y 每轮走两步。如果轨迹里有环，两指针必然相遇。每轮算一个差值，然后求 gcd。
朴素写法如下：
ll pollard_rho_naive(ll n) { if (n % 2 == 0) return 2; ll c = rng() % (n - 1) + 1; ll x = rng() % (n - 2) + 2; ll y = x; ll d = 1; while (d == 1) { x = f(x, c, n); y = f(f(y, c, n), c, n); d = gcd(llabs(x - y), n); } if (d == n) return pollard_rho_naive(n); // 这次参数白给，重来 return d; }这段代码思路完全正确，但有两个隐患：
- 每轮循环都要做一次 gcd，而 gcd 在迭代次数很多时是明显的常数开销；
- 一旦出现 d == n，说明两指针在模 n 意义下重合了，本轮信息全部报废，只能换 c 重来。
第 4 章的批量优化就是同时解决这两个问题。
3.3 参数选择和失败重来
c 的取值范围一般是 [1, n-1]，初始点 x 在 [2, n-2]。为什么不能取 0 或 1？如果 c=0 且 x=0，序列恒为 0，整个算法直接卡死；即使只取 c=0，当 x 比较特殊时序列也可能迅速退化。随机数生成对 PR 影响很大，后面专门讲。
d==n 意味着什么？它说明 |x-y| 是 n 的倍数，两个指针在模 n 下已经重合。这本身不说明没有因子，只说明“当前这一轮没抓到”。朴素写法直接放弃重来，确实有概率把到手的信息扔掉。所以实际竞赛中很少有人直接用朴素版，大家更常用的是第 4 章的倍增批量版。
4. 批量GCD与倍增：把“偶然发现”升级成“必然效率”
4.1 问题：gcd调用太频繁
朴素版每轮迭代都算一次gcd(llabs(x-y), n)。gcd 本身复杂度是 O(log n)，单次看不算贵，但 PR 面对待分解的大合数时，迭代次数可能达到百万级，百万次 gcd 乘上几百个测试数据，差距就出来了。
优化的思路非常朴素：gcd 不用每轮做，攒一批差值，把这一批的差值全部乘起来，对乘积只做一次 gcd。
4.2 反直觉的事实：乘积一次gcd不会漏因子
这一步是全文最容易被误解的地方。有人会担心，万一因子被其他差值“掩盖”了怎么办？其实数学上不会漏：
如果某个差值 delta_i 与 n 的 gcd > 1，即 delta_i 含有 n 的某个因子 p，那么把这批差值全部乘起来得到 prod，prod 同样含有 p。因此
gcd(prod, n) >= p > 1也就是说，只要这批差值里存在任何与 n 不互质的一项，对乘积求 gcd 一定至少能发现这个因子。不会漏。
真正要小心的反而是另一种情况：gcd(prod, n) == n。这说明整批差值的乘积已经是 n 的倍数。可这并不能直接说明单个差值能给出因子，因为可能是好几个差值分别贡献了 n 的不同因子，乘在一起“拼凑”出完整的 n。这时候必须回退，逐个差值重新求 gcd，定位到具体是哪一个位置贡献了因子。
4.3 撞到n时的回退
批量版本的核心结构：用 step 控制每批检测的步数，step 从 1 开始指数增长，每批开始时记录起点 x，然后让 y 往前走 step 步，每走一步把差值乘进 prod。最后对 prod 求一次 gcd。
- 如果 gcd(prod, n) == 1：这批全部与 n 互质，没货，继续下一批；
- 如果 gcd(prod, n) == d，并且 1 < d < n：拿到因子 d，直接返回；
- 如果 gcd(prod, n) == n：这批里有货，但被乘积掩盖了，把 y 退回到这批起点，逐个重新求 gcd，找到第一个非平凡因子为止。
指数增长的 step 保证了单批批量不会太大，同时整体采样范围迅速扩大。相比朴素版每轮 gcd，batch 方法把 gcd 次数压缩到 O(log step) 级别，而乘法成本远低于 gcd。
4.4 完整pollard_rho代码与逐行解读
ll pollard_rho(ll n) { if (n % 2 == 0) return 2; ll c = rng() % (n - 1) + 1; ll y = rng() % (n - 2) + 2; ll d = 1; for (int step = 1; ; step <<= 1) { ll x = y; // 当前批次的起点 ll prod = 1; for (int i = 0; i < step; i++) { y = f(y, c, n); prod = mul_mod(prod, llabs(x - y), n); } d = gcd(prod, n); if (d == n) { // 被掩盖了，回退搜索 y = x; for (int i = 0; i < step; i++) { y = f(y, c, n); d = gcd(llabs(x - y), n); if (d != 1) break; } } if (d != 1) break; } if (d == n) return pollard_rho(n); return d; }几个细节值得说：
- 每批开始时
x = y，把上一批结束的位置作为新起点，这样整个序列始终向前推进，不会重复探测已经扫过的区间。 - 回退循环里，y 从起点 x 重新走 step 步，每一步都求一次 gcd。这里最坏情况会退回到朴素版的每步 gcd，但实际很少触发，因为 d==n 意味着整批乘积恰好是 n 的倍数，概率远低于 d==1。
- 如果回退结束 d 仍然等于 n，说明这批里可能只有模 n 的完全碰撞，没有贡献非平凡因子，重新选 c 再来。
我最初从朴素版换到批量版之后，P4718 的用时肉眼可见降了一个档次，稳定性也更好了。
5. 快速乘与随机数：决定模板稳定的两只黑手
5.1 乘法溢出：不是数学错，是位数错
P4718 的 n 可以接近 4×10^18，这意味着任意两个小于 n 的数相乘，结果可能接近 1.6×10^37，远远超出 64 位整数的表示范围。直接写x*x%mod会在乘法瞬间溢出，结果变成一个随机的大数，后续全错。
所以必须有一个“快速乘”函数，用不溢出的方式计算 a×b mod mod。最直观的方案是__int128——这是 GCC 等编译器提供的 128 位扩展类型，足够装下两个 64 位数的乘积：
ll mul_mod(ll a, ll b, ll mod) { return (ll)((__int128)a * b % mod); }这是我在竞赛环境下的首选，简洁、准确、不容易写错。注意__int128并不是标准 C++ 类型，但主流 OJ 的编译器都支持，放心用。
如果换到不支持__int128的环境，还有一招利用 long double 的近似商法：
ll mul_mod(ll a, ll b, ll mod) { ll q = (ll)((long double)a * b / mod); ll res = a * b - q * mod; if (res < 0) res += mod; if (res >= mod) res -= mod; return res; }原理是先用 long double 算出 a*b/mod 的近似商，再用整数运算修正。long double 在 x86 下有 64 位尾数，对 64 位整数范围内的商精度足够，最后通过两次加减修正误差。这个写法在个别环境下能跑，但可读性和安全性都不如__int128，我一般只在没法用__int128时才考虑。
一个有符号数细节：llabs(x-y)一定要用，因为 x 和 y 的大小关系不确定，差值可能是负数，直接传给 gcd 会出问题。
5.2 随机数：rand()在这个量级早就不够用了
rand()的返回值在大多数实现下只有 2^31 量级，而我们的 n 是 4×10^18 量级，用rand()生成 PR 的初始值，本质上是在一个极小的范围内取样，初始点分布极不均匀，算法很容易陷入长循环。
正确姿势是用mt19937_64，它能生成 64 位随机数：
mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());注意生成 c 和 y 时要避开 0 和 1 这类退化值：
ll c = rng() % (n - 1) + 1; // c ∈ [1, n-1] ll y = rng() % (n - 2) + 2; // y ∈ [2, n-2]如果 n 很小（比如 n=4），这两个区间会非常窄，所以 pollard_rho 开头特判偶数直接返回 2 就显得很重要——它把这一类边界情况直接挡在了随机逻辑之外。
5.3 运算符号与类型边界
在 P4718 的数据范围内，long long可以装下 n，__int128可以装下乘法结果。但如果哪天你把模板扩展用到更大的数，记得检查三个边界：
- 有符号 64 位最大是 2^63-1，如果再往上需要用
unsigned long long； llabs(LLONG_MIN)会溢出，虽然 PR 场景里差值不太可能达到这个量级，但心里要有数；f(x) = (mul_mod(x,x,n)+c)%n中，mul_mod(x,x,n)+c最大约 2n-2，只要 n 不超过 2^63 量级，加法不会溢出。
这些边界不是玄学，是模板题专门用来卡人的地方。
6. 完整模板与递归分解：一份能直接交的代码
6.1 递归分解 factor 的设计
factor 的逻辑前面已经提过，这里展开。判断顺序是：n==1 直接返回，MR 判质数就记录，否则用 PR 拿因子，分别递归。这个递归必须保证最终只剩质因子，所以 PR 返回的 d 就算本身是合数，也会被下一层 factor 继续拆分。
这个过程中最容易忽略的是偶数。虽然 pollard_rho 开头已经特判n % 2 == 0返回 2，但 factor 递归时 n/2 仍可能是偶数，比如 n=4，会连续递归两次 factor(2)。实际上 PR 特判偶数把最麻烦的分支直接短路了，剩下的奇数合数交给 MR 和 PR 处理，逻辑非常干净。
重复质因子也能正确处理：factor(8) 会得到三个 2，最终 max 仍然是 2，不影响“最大质因子”这个答案。
6.2 完整模板代码
以下是我整理后可以直接提交的完整模板，按 P4718 的范围写的：
#include <bits/stdc++.h> using namespace std; using ll = long long; using ull = unsigned long long; mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count()); inline ll mul_mod(ll a, ll b, ll mod) { return (ll)((__int128)a * b % mod); } inline ll f(ll x, ll c, ll mod) { return (mul_mod(x, x, mod) + c) % mod; } ll gcd(ll a, ll b) { while (b) { ll t = a % b; a = b; b = t; } return a; } ll qpow(ll a, ll b, ll mod) { ll res = 1 % mod; while (b) { if (b & 1) res = mul_mod(res, a, mod); a = mul_mod(a, a, mod); b >>= 1; } return res; } bool isPrime(ll n) { if (n < 2) return false; static const ll smallPrimes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}; for (ll p : smallPrimes) { if (n == p) return true; if (n % p == 0) return false; } ll d = n - 1, s = 0; while ((d & 1) == 0) d >>= 1, s++; static const ll bases[] = {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL}; for (ll a : bases) { if (a % n == 0) continue; ll x = qpow(a, d, n); if (x == 1 || x == n - 1) continue; bool ok = false; for (int r = 0; r < s - 1; r++) { x = mul_mod(x, x, n); if (x == n - 1) { ok = true; break; } } if (!ok) return false; } return true; } ll pollard_rho(ll n) { if (n % 2 == 0) return 2; ll c = rng() % (n - 1) + 1; ll y = rng() % (n - 2) + 2; ll d = 1; for (int step = 1; ; step <<= 1) { ll x = y; ll prod = 1; for (int i = 0; i < step; i++) { y = f(y, c, n); prod = mul_mod(prod, llabs(x - y), n); } d = gcd(prod, n); if (d == n) { y = x; for (int i = 0; i < step; i++) { y = f(y, c, n); d = gcd(llabs(x - y), n); if (d != 1) break; } } if (d != 1) break; } if (d == n) return pollard_rho(n); return d; } void factor(ll n, vector<ll>& fac) { if (n == 1) return; if (isPrime(n)) { fac.push_back(n); return; } ll d = pollard_rho(n); factor(d, fac); factor(n / d, fac); } int main() { int T; scanf("%d", &T); while (T--) { ll n; scanf("%lld", &n); vector<ll> fac; factor(n, fac); ll ans = 1; for (ll p : fac) ans = max(ans, p); printf("%lld\n", ans); } return 0; }6.3 我的踩坑记录与自测建议
说几个我实际踩过的坑。
第一，一开始没写mul_mod，直接在快速幂里用a*a%mod，小数据全对，一上大数就 WA。排查了很久才发现是乘法溢出。从那以后我形成习惯：只要模数可能超过 1e9，所有乘法一律走mul_mod，不给自己留侥幸空间。
第二，随机数生成器一开始用了rand()，交上去偶尔会 TLE，像运气不好。后来换成mt19937_64，并排除 0 和 1，算法稳定性明显好了。PR 本质上是个随机算法，种子质量直接决定最坏情况的触发概率，别在这种地方省钱。
第三，忘记特判 n=1。题目数据如果包含 1，fac 会为空，max 函数直接报错。输出前把 ans 初始化为 1 就能兜底。
自测建议：自己写一个暴力分解函数，随机生成 1e5 以内的 n，把模板输出和暴力结果对比，几秒钟就能把绝大多数边界问题暴露出来。再随机生成一批大质数和大合数，重点看会不会死循环或返回错误因子。我每次改完模板都这么跑一遍，比对着样例数据干瞪眼效率高得多。
最后说一个心得：P4718 这类模板题，真正的门槛不在背板子，而在于你有没有把“为什么会这样写”想清楚。不信你可以试试，把第 4 章的批量优化临时改回朴素版，数据不强的时候照样 AC，但一旦碰上精心构造的数据，差距立刻显现。模板可以有默认写法，但每一处优化背后的原因，才是这题真正想让你掌握的东西。
