说实话第一眼看到 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 只保证给你一个“非平凡因子”不保证它是质数。比如 n1009×1013PR 第一次很可能直接命中 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如果 c0 且 x0序列恒为 0整个算法直接卡死即使只取 c0当 x 比较特殊时序列也可能迅速退化。随机数生成对 PR 影响很大后面专门讲。dn 意味着什么它说明 |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那么把这批差值全部乘起来得到 prodprod 同样含有 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 保证了单批批量不会太大同时整体采样范围迅速扩大。相比朴素版每轮 gcdbatch 方法把 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但实际很少触发因为 dn 意味着整批乘积恰好是 n 的倍数概率远低于 d1。如果回退结束 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 很小比如 n4这两个区间会非常窄所以 pollard_rho 开头特判偶数直接返回 2 就显得很重要——它把这一类边界情况直接挡在了随机逻辑之外。5.3 运算符号与类型边界在 P4718 的数据范围内long long可以装下 n__int128可以装下乘法结果。但如果哪天你把模板扩展用到更大的数记得检查三个边界有符号 64 位最大是 2^63-1如果再往上需要用unsigned long longllabs(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 的逻辑前面已经提过这里展开。判断顺序是n1 直接返回MR 判质数就记录否则用 PR 拿因子分别递归。这个递归必须保证最终只剩质因子所以 PR 返回的 d 就算本身是合数也会被下一层 factor 继续拆分。这个过程中最容易忽略的是偶数。虽然 pollard_rho 开头已经特判n % 2 0返回 2但 factor 递归时 n/2 仍可能是偶数比如 n4会连续递归两次 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, vectorll 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); vectorll 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 本质上是个随机算法种子质量直接决定最坏情况的触发概率别在这种地方省钱。第三忘记特判 n1。题目数据如果包含 1fac 会为空max 函数直接报错。输出前把 ans 初始化为 1 就能兜底。自测建议自己写一个暴力分解函数随机生成 1e5 以内的 n把模板输出和暴力结果对比几秒钟就能把绝大多数边界问题暴露出来。再随机生成一批大质数和大合数重点看会不会死循环或返回错误因子。我每次改完模板都这么跑一遍比对着样例数据干瞪眼效率高得多。最后说一个心得P4718 这类模板题真正的门槛不在背板子而在于你有没有把“为什么会这样写”想清楚。不信你可以试试把第 4 章的批量优化临时改回朴素版数据不强的时候照样 AC但一旦碰上精心构造的数据差距立刻显现。模板可以有默认写法但每一处优化背后的原因才是这题真正想让你掌握的东西。