1. 项目背景与整体设计思路1.1 P4718到底在考什么P4718这道模板题本质上只要求做两件事第一判断给定的整数n是不是质数第二如果n是合数找出它的最大质因子。听起来很简单对吧但真正做起来你会发现数据规模直接拉到10^18级别甚至更大。这个量级下最直觉的试除法根本走不通要试除到sqrt(n)也就是大约10^9次循环光是单组数据就已经不可接受更何况题面往往是多组询问。所以这道题考查的不是你会不会暴力试除而是你手头有没有一套真正能处理超大整数的因数分解工具。在我个人做数论题的实战里最稳定的组合就是Miller-Rabin素性测试配合Pollard-Rho随机分解。前者负责快速判定质数后者负责把合数拆开。两者配合起来10^18这个量级的数基本可以在几十毫秒内完成分解这也是P4718这道模板题最核心的考点。1.2 为什么是Miller-Rabin Pollard-Rho先看其他方案。试除法前面已经排除了。Fermat分解虽然对某些特殊形式的合数比如两个相近的因子很有效但对一般情况完全看运气。Pollard p-1算法只在p-1由小质因子构成时表现好适用范围太窄。椭圆曲线分解法ECM理论很强但实现复杂度高模板代码动辄几百行而且常数不小。Pollard-Rho的优势在于实现简单核心循环不到30行没有苛刻的适用前提期望复杂度是O(n^{1/4})量级。Miller-Rabin的好处更直接在64位整数范围内只需要固定几个基底就能做到完全确定性判定速度快到可以忽略不计。两个算法结合正好覆盖了先判素、再分解的全部需求。算法适用场景实现难度缺点试除法小整数极低碰上10^18直接报废Fermat分解因子接近sqrt(n)低普通合数效率极差Pollard p-1p-1光滑的数中对一般大质因子束手无策ECM通用大数分解高实现复杂常数大Miller-Rabin Pollard-Rho通用大整数中随机化理论上是概率算法1.3 整体流程怎么串起来整个方案的执行链条是这样的输入n先用Miller-Rabin判断它是否质数。如果是直接输出Prime如果不是就调用Pollard-Rho找出一个非平凡因子d然后递归地对d和n/d继续分解直到两边都变成质数。过程中用一个哈希表或map统计每个质因子的出现次数最后扫一遍取最大的质因子输出即可。这个流程有一个好处Pollard-Rho只负责拆出任意一个因子而不直接追求完整分解。拆出来之后递归层数通常很浅因为每次都能把规模缩小整体跑起来非常稳。我把这套逻辑直接封装成模板存了下来之后在好几道数论题里都复用到了。2. 前置算法Miller-Rabin素性测试2.1 从费马小定理到强伪素数检测Miller-Rabin的理论根基是费马小定理如果n是质数那么对于任意与n互质的a一定有a^(n-1) ≡ 1 (mod n)但反过来不成立。有些合数对某些a也满足这个式子比如所有的Carmichael数。Miller-Rabin在费马检测的基础上加了一道二次探测思路是这样的假设n是奇质数那么n-1一定是偶数。把n-1写成d * 2^s的形式其中d是奇数。根据费马小定理a^(n-1) (a^d)^(2^s) ≡ 1 (mod n)如果某个中间结果x满足x^2 ≡ 1 (mod n)那么x必须是1或n-1因为质数模意义下1的平方根只有这两个。所以算法流程变成先算x a^d然后连续平方s次。如果在平方过程中x从某个不等于±1的数突然变成了1说明n一定是合数这就是强伪素数检测的原理。我把这个过程在下面用小规模的例子验证过。比如n 561这是最小的Carmichael数费马检测会傻傻地认为它是质数对a2成立但Miller-Rabin拿a2走一遍二次探测就能把它揪出来。这正是Miller-Rabin比朴素费马检测强的地方。2.2 基底怎么选才稳妥Miller-Rabin不是绝对判定而是概率判定。但对64位以内的整数有一组经过验证的确定性基底只要这些都测一遍结果就是100%准确的{2, 325, 9375, 28178, 450775, 9780504, 1795265022}这组基底覆盖了2^64以内的所有整数。也就是说在P4718的数据范围内用它就是确定性的不存在误判概率。实际代码里我会先把n对2、3、5、7这些小质数取模排除掉能被小质数整除的情况再进入Miller-Rabin主流程。一方面提前过滤掉大量平凡情况另一方面可以顺手处理n为这些质数本身的case。注意一点基底对n取模后如果等于0要跳过这个基底因为此时a是n的倍数直接用反而会干扰后续判定。这个细节很容易被忽略但实际测试数据里真有可能出现。2.3 乘法溢出最阴间的坑Miller-Rabin的核心操作是快速幂取模而快速幂里必然要算a * b % mod。在64位整数范围内两个10^18量级的数相乘结果可能到10^36远超64位范围。很多人栽在这里死活分不出来到底是算法错了还是溢出错了。我的建议是直接用__int128。在C里把两个long long转成__int128做乘法再取模编译器和运行时都支持效率也够。但这里有个非常容易踩的坑__int128虽然有128个比特但有符号表示时正上限是2^127 - 1。两个有符号long long相乘最大值约是(2^63 - 1)^2 ≈ 2^126恰好没超过上限所以只要坚持用long long配合__int128就不会溢出。如果图省事用了unsigned long long两个最大值相乘接近2^128会把有符号__int128直接撑爆结果错误而且极难排查。我当年在这个问题上浪费过一整个下午。下面是Miller-Rabin的标准实现我加了详细注释using ll long long; ll gcd(ll a, ll b) { while (b) { ll t a % b; a b; b t; } return a; } ll mul_mod(ll a, ll b, ll mod) { return (ll)((__int128)a * b % mod); } ll pow_mod(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 is_prime(ll n) { if (n 2) return false; static ll small[] {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}; for (ll p : small) { if (n % p 0) return n p; } // 把 n-1 拆成 d * 2^s ll d n - 1; int s 0; while ((d 1) 0) { d 1; s; } // 64位范围内确定性基底 static ll bases[] {2, 325, 9375, 28178, 450775, 9780504, 1795265022}; for (ll a : bases) { if (a % n 0) continue; ll x pow_mod(a, d, n); if (x 1 || x n - 1) continue; bool ok false; for (int r 1; r s; r) { x mul_mod(x, x, n); if (x n - 1) { ok true; break; } } if (!ok) return false; } return true; }中间那个二次探测循环的理解方式s代表还能平方的次数。初始x a^d如果已经是±1直接认为这个基底通过否则就平方s-1次如果在某次平方后变成n-1说明它可以继续安全地平方到最终结果1这个基底也通过。如果平方到最后一刻都没有出现n-1说明这一路平方下来必然经过了非±1的平方等于1的违禁情况判定为合数。3. Pollard-Rho算法核心3.1 生日悖论怎么帮我们找因子Pollard-Rho的思路和常规分解完全不同。它不试图系统性地枚举因子而是利用了一个概率现象生日悖论。假设n有一个较小的质因子p我们不知道p是多少但知道它存在。现在在[0, n-1]范围内随机生成一组数由于这些数对p取模后等价于在[0, p-1]范围内随机取值根据生日悖论大约只需要O(sqrt(p))个随机数就会出现两个数对p取模相等。一旦出现这种情况它们的差就是p的倍数那么gcd(|x - y|, n)大概率能把这个p提取出来。这个过程和D罩杯悖论是同一个数学内核23个人里有两人生日相同的概率超过50%而不是很多人直觉想的23/365。放在因子分解里p只有10^6量级时sqrt(p)大约1000很快就能碰撞。即使p是10^9量级sqrt(p)不过是31623次随机尝试依然可行。3.2 伪随机函数与Floyd判圈Pollard-Rho的巧妙之处在于它不用真的存下一大堆随机数而是用一个伪随机函数迭代生成序列f(x) (x^2 c) mod n从某个x0出发不断套f得到一条序列。这条序列在模p意义下必然进入循环因为状态只有p种。进入循环后就会出现两个位置x和y使得x ≡ y (mod p)但x ≠ y (mod n)。此时|x - y|包含因子p。找这个碰撞点不需要存储历史用Floyd判圈算法即可一个指针每次走一步另一个指针每次走两步两者最终会在圈内相遇。相遇时计算gcd(|x - y|, n)如果结果是1就继续迭代如果结果是n说明运气不好x和y在模n意义下也相等了换一个c重来如果结果在1和n之间就找到了一个非平凡因子。这里说句实话单看描述会觉得这个算法有点玄学我第一次看的时候也半信半疑。但实际跑起来它对绝大多数数据都在极短的时间内出结果只有极小概率会反复撞上d n导致重试。这也是随机化算法的典型特点理论上限不保证但实践中非常可靠。3.3 批量gcd性能的关键优化初版Pollard-Rho很容易写出每迭代一步就调一次gcd的版本。但gcd本身有对数级复杂度而且迭代次数是sqrt(p)量级频繁调用会让常数爆炸。优化思路非常直白把连续多次迭代的|x - y|累乘起来对n取模每隔一定步数统一做一次gcd。因为如果其中某个差值已经包含了因子p这些差值的乘积也必然包含因子pgcd一下同样能提取出来。这个技巧把gcd的调用次数压缩了上百倍实测提速非常明显。具体到参数我用的是每128步做一次gcd。这个数不是拍脑袋定的太小了gcd次数多太大了乘积的模运算压力上升而且一旦乘积变成0要回退一整段。128是很多实现里常用的折中值。还有一种倍增写法每次让快指针冲出去2^k步然后统一做一次gcd效果类似。两种我都试过128步批量gcd的版本更容易调试所以最后模板里保留了这个方案。Pollard-Rho核心代码如下ll pollard_rho(ll n) { if (n % 2 0) return 2; if (n % 3 0) return 3; if (n % 5 0) return 5; // 实际运行时用随机种子初始化 while (true) { ll c (ll)mt() % (n - 1) 1; ll x (ll)mt() % (n - 2) 2; // 在[2, n-1]选起点 ll y x; ll d 1; ll prod 1; int step 0; auto f [](ll v) { return (mul_mod(v, v, n) c) % n; }; while (d 1) { x f(x); y f(f(y)); ll delta llabs(x - y); if (delta 0) break; // 模n下也撞了换参重来 prod mul_mod(prod, delta, n); step; if (step 128) { d gcd(prod, n); prod 1; step 0; } } if (step 0 d 1) d gcd(prod, n); if (d 1 d n) return d; } }随机数部分我用了mt19937而不是cstdlib的rand()。原因有两个rand()在多数环境下只有2^31左右的周期对于高维随机需求有点勉强更重要的是rand()的分布质量在Pollard-Rho这种场景下确实出现过连续多次选到不好参数的情况。mt19937配合种子整体稳定性明显好一截。4. 完整模板与P4718应用4.1 可以直接抄的完整代码把Miller-Rabin和Pollard-Rho拼在一起再加上递归分解逻辑就是一套完整的模板。我整理的版本贴在这里已经压掉了所有跟核心无关的修饰可以直接往OJ上交当然平时还是要自己敲一两遍避免比赛时手生#include bits/stdc.h using namespace std; using ll long long; mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count()); ll gcd(ll a, ll b) { while (b) { ll t a % b; a b; b t; } return a; } ll mul_mod(ll a, ll b, ll mod) { return (ll)((__int128)a * b % mod); } ll pow_mod(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 is_prime(ll n) { if (n 2) return false; static ll small[] {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}; for (ll p : small) { if (n % p 0) return n p; } ll d n - 1; int s 0; while ((d 1) 0) { d 1; s; } static ll bases[] {2, 325, 9375, 28178, 450775, 9780504, 1795265022}; for (ll a : bases) { if (a % n 0) continue; ll x pow_mod(a, d, n); if (x 1 || x n - 1) continue; bool ok false; for (int r 1; r s; 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; if (n % 3 0) return 3; if (n % 5 0) return 5; while (true) { ll c (ll)(rng() % (n - 1)) 1; ll x (ll)(rng() % (n - 2)) 2; ll y x; ll d 1; ll prod 1; int step 0; auto f [](ll v) { return (mul_mod(v, v, n) c) % n; }; while (d 1) { x f(x); y f(f(y)); ll delta llabs(x - y); if (delta 0) break; prod mul_mod(prod, delta, n); step; if (step 128) { d gcd(prod, n); prod 1; step 0; } } if (step 0 d 1) d gcd(prod, n); if (d 1 d n) return d; } } void factorize(ll n, mapll, int cnt) { if (n 1) return; if (is_prime(n)) { cnt[n]; return; } ll d pollard_rho(n); factorize(d, cnt); factorize(n / d, cnt); } int main() { int T; scanf(%d, T); while (T--) { ll n; scanf(%lld, n); if (is_prime(n)) { printf(Prime\n); continue; } mapll, int cnt; factorize(n, cnt); ll ans 0; for (auto kv : cnt) ans max(ans, kv.first); printf(%lld\n, ans); } return 0; }4.2 主流程的执行过程详解以n 1234567890123456789为例这个数是合数主函数先调用is_primeMiller-Rabin用一组基底走一遍很快返回false。然后进入factorize第一次调用pollard_rho可能返回一个因子比如3或者更大的因子。假设返回了3那么递归拆3质数直接记录和n/3下一步再拆n/3如果它还是合数继续调用pollard_rho拆出一个因子如此往复。这里有个非常关键的性质递归过程中每一次分解出的两个数p和n/p规模至少有一个缩小到原来的sqrt以下。所以虽然理论最坏情况可能退化但实际递归深度就是O(log n)量级不可能出现无穷递归。map统计质因子时key从小到大排序最后遍历取最大值作为最大质因子这就是P4718要求的结果。4.3 期望复杂度到底是多少Pollard-Rho期望找最小因子p需要O(sqrt(p))次迭代。每次迭代里有两次f调用、一次乘法和每128步一次的gcd。如果n的最小质因子是n^{1/2}量级那么总迭代大约O(n^{1/4})再加上gcd的log n整体期望就是O(n^{1/4} log n)。对10^18量级的数来说n^{1/4}大约31623配合上128步一次gcd实际毫秒级就能分解完。这也是为什么这套组合在竞赛圈成为大数分解事实标准的原因。需要强调一下这个复杂度是期望意义下的。极端数据有可能让某一次pollard_rho反复返回n但概率极低而且通过换c重试基本都能绕过去。我做过的实测中连续分解几千个10^18量级随机合数没有一次卡死超过一秒。5. 常见问题与踩坑记录5.1 死循环和长时间跑不出来这是Pollard-Rho新手最容易遇到的问题程序在某组数据上卡住不动或者明显超过时限。我排查下来最常见的两个原因一是随机参数选得不好。c如果选得让伪随机序列很快就进入一个很短的圈且x和y在这个圈上始终不产生有效的gcd结果就会反复空转。解决方法是给pollard_rho外层套一个while(true)一旦检测到d n或者迭代了足够多次数没有进展就重新生成c和x0。我的模板里已经这么做了。二是delta等于0的情况。当x和y在模n意义下完全相同时llabs(x - y) 0后续prod乘上0直接变0gcd(0, n) n然后程序会判断d 1 d n不成立进入下一轮重试。这个逻辑是自洽的但如果你写的代码在delta 0时没有跳出而是继续乘0就会陷入无限循环。所以delta的那个if判断不能省。5.2 乘法溢出与类型混用前面提过__int128配合long long是安全的但实际代码里很容易手滑引入unsigned类型。比如rng()本身返回unsigned long long如果直接拿它和n运算可能把类型带偏。我的做法是先强转成long long再计算保证后续运算都走mul_mod里的__int128逻辑。另外模运算里不要轻易用减法。比如计算(x - y) % n如果直接写x - y得到负数后再取模行为可能不符合直觉。最好用llabs取绝对值。因为x和y都在[0, n-1]里差值绝对值不超过n-1直接作为乘积因子没问题。5.3 边界数据的处理n 2、3、5这些小数必须特判is_prime里的小质数表已经覆盖了。n 1的情况也要提前处理is_prime返回false但如果直接走pollard_rho会出问题因为pollard_rho里假设n至少是个大于5的合数。我的factorize开头有if (n 1) return这就是防这一手。还有一个隐蔽问题当n是一个大质数的平方比如p^2pollard_rho有时会慢一些。因为伪随机序列在模p意义下碰撞后gcd结果很大概率是p但如果恰好得到n也就是p^2就会重来。这种情况我在极限数据里遇到过不算频繁但如果你要处理p^2这种特殊构造的hack数据建议把pollard_rho里的如果d n就换c重试逻辑调严一点或者增加一个重试次数超过阈值就强制换参数的兜底。5.4 常见问题速查表问题现象可能原因解决方案程序卡住不结束随机参数c不好序列长期空转外层while循环换c重试返回值一直是nx和y在模n下也碰撞重新生成c和起点x0结果错误但小数据正确乘法溢出全部走__int128乘法Miller-Rabin误判小质数基底越界或漏判先用小质数表过滤分解大平方数极慢因子p太小但碰撞后gcd得n增加重试上限或强制换参6. 实测表现与扩展思考6.1 和朴素分解的差距有多大我用随机生成的10^18量级合数做了一组对比。本地环境普通桌面CPU试除法在10^9级别循环已经卡到无法接受而Miller-Rabin Pollard-Rho的完整分解对大多数数耗时在10毫秒以内。十几组数据全部跑完总耗时不到0.2秒。差距是数量级的。6.2 算法还能怎么扩展模板里用的是long long理论范围是到9.22 * 10^18。如果你想处理更大范围的数比如到10^30可以把数据类型全部换成__int128Miller-Rabin的基底需要扩展或改用更多随机基底。Pollard-Rho的核心循环不受影响但mul_mod要改成更通用的模乘实现比如用long double快速乘因为两个__int128相乘在乘法过程中会溢出到256位无法直接取模。另外这套分解框架不只是用在做题上。大整数分解在信息安全领域是很多协议安全性的基石比如RSA私钥的安全性就建立在大整数分解困难这个假设上。虽然实际破解远不能靠Pollard-Rho做到但理解这个算法能让你对为什么RSA密钥要选大质数且两个质因子尽量接近有更直观的认识。密码学里那些关于密钥长度不断演进的分析本质上都是围绕大数分解的算法复杂度展开的。最后分享一个我个人的使用习惯在需要多次调用pollard_rho的程序里我不会在每次调用时重新初始化随机种子而是全局只用一次seed让后续调用自然取得不同的随机序列。这样做的好处是既保证了随机性又避免了多次seed导致前后参数重复的问题。另外如果你打算把这套模板迁移到多线程环境跑批量分解每个线程需要独立的随机数生成器实例千万不要共享同一个rng对象否则并发访问会出现随机数序列错乱直接表现为分解结果异常。