大阶乘计算进阶：高精度大数与分治乘法实战解析

最近整理一份算法题单时被一道“大阶乘计算-进阶题”卡了一下。题目看着很简单不过就是输出一个正整数 n 的阶乘n 最大能到 100000 甚至更高。当时我的第一反应是“for 循环乘上去不就行了”真正动手才发现这一题考的并不是循环而是三个层面的综合问题大数怎么存、乘法怎么优化、内存和输出怎么规划。对于没接触过这类题的人来说最容易踩的坑就是觉得“n! 就是乘 n 次”却不知道某些语言里连第 21 次乘法都会把 64 位整数直接撑爆。这篇文章把我从踩坑到最终跑通的全过程写下来。内容安排上从基础数据类型极限说起到高精度大数的底层表示再到分治乘法、FFT 加速、素数分解法这几种进阶思路的对比最后针对“精确值、末尾零、位数、取模”四类常见考法给出可以直接用的方案。如果你正在刷题准备笔试面试或者想在工程中处理超大整数运算这篇应该能帮你省不少时间。1. 为什么普通整型存不下阶乘数值增长的恐怖曲线1.1 从 int 到 long long类型上限对比很多新手写阶乘第一版代码通常长这样long long fact(int n) { long long res 1; for (int i 2; i n; i) { res * i; } return res; }这段代码在 n 比较小时没什么问题但只要 n 超过某个阈值结果就开始“放飞自我”。我们先看一组具体数字nn! 的近似值int 能否保存unsigned long long 能否保存double 的表现5120能能120136,227,020,800溢出已经超过 21 亿上限能6227020800202,432,902,008,176,640,000溢出能约 2.43e182.43290200817664e182151,090,942,171,709,440,000溢出溢出已经达到 5.1e195.109094217170944e19170约 7.26e306溢出溢出勉强能表示171约 1.24e309溢出溢出爆掉变成 infint 的上限大概是 21 亿也就是 2^31 - 1。13! 是 62 亿多直接超过这个值。unsigned long long 上限大约是 1.84e1920! 是 2.43e18勉勉强强装得下但到 21! 就变成 5.1e19也超了。这是很多人第一次遇到“大数”概念的瞬间不是代码逻辑写错了而是数据类型本身装不下数。1.2 精度陷阱浮点数能表示却会出错有人可能会说“那用 double 不就行了double 上限很大好像能表示到 10 的 308 次方。”double 确实能表示很大的数但问题是它只能精确表示大约 15 到 17 位有效数字后面的尾巴全是约数。比如 22! 的精确值是 1,124,000,727,777,607,680,000而 double 显示出来是 1.1240007277776077e21从第 17 位往后就跟精确值对不上了。在工程计算里这种误差可能可以接受但“大阶乘计算-进阶题”这类题目通常要求输出精确的每一位数字。一旦要求精确float、double 就全部出局必须引入任意精度整数也就是我们常说的高精度大数。2. 高精度大数表示一切方案的地基2.1 用数组手工实现大整数既然一个整型变量放不下很大的数那就用一组数字来表现它。最直观的想法是用一个数组每个元素存一位十进制数字从低位到高位排列。比如数字 12345 存为[5, 4, 3, 2, 1]第 0 位是个位第 1 位是十位。乘法的时候就模拟我们在纸上列竖式的过程逐位相乘记录进位最后统一处理。一位一位存逻辑最简单但空间和时间效率都很差。更通用的做法是“万进制”也就是每个数组元素存 0 到 9999 之间的一个四位十进制块。这样数组长度直接减少到原来的四分之一而且乘法时在一次循环里能处理四位。为什么是万进制而不是十万进制、百万进制关键在乘法运算的溢出边界。当前结果位乘以当前乘数时单次乘积最大值是 9999 * 9999约等于 1 亿再加进位int 完全接得住如果选更大进制两个大数相乘时中间结果就可能超出 int 安全范围需要被迫改成 long long性能反而下降。下面是一个精简但能直接跑的 C 高精度类#include bits/stdc.h using namespace std; class BigInt { public: vectorint d; // 万进制低位在前 BigInt() { d.clear(); d.push_back(0); } BigInt(long long x) { d.clear(); if (x 0) { d.push_back(0); return; } while (x 0) { d.push_back(x % 10000); x / 10000; } } // 大数乘以一个普通整数 void multiply(int x) { long long carry 0; for (size_t i 0; i d.size(); i) { long long cur (long long)d[i] * x carry; d[i] cur % 10000; carry cur / 10000; } while (carry 0) { d.push_back(carry % 10000); carry / 10000; } } string toString() const { stringstream ss; ss d.back(); for (int i (int)d.size() - 2; i 0; i--) { ss setw(4) setfill(0) d[i]; } return ss.str(); } };用这个类算阶乘就很简单BigInt factorialNaive(int n) { BigInt res(1); for (int i 2; i n; i) { res.multiply(i); } return res; }很多高精度题目的标准解法就是这个版本数据范围小于 10000 时完全够用代码也短。但如果你拿它去算 50000! 或 100000!就会明显感觉程序“卡住”了。2.2 第二版两个大数的乘法如果想进一步优化光会“大数乘小数”还不够还得实现“大数乘大数”。这一步是分治乘法的基础。做法也不复杂把两个大数看成两个整数数组逐位相乘再累加最后统一进位。BigInt multiplyBig(const BigInt a, const BigInt b) { vectorint res(a.d.size() b.d.size(), 0); for (size_t i 0; i a.d.size(); i) { long long carry 0; for (size_t j 0; j b.d.size(); j) { long long cur (long long)res[i j] (long long)a.d[i] * b.d[j] carry; res[i j] cur % 10000; carry cur / 10000; } int k i (int)b.d.size(); while (carry 0) { long long cur (long long)res[k] carry; res[k] cur % 10000; carry cur / 10000; k; } } while (res.size() 1 res.back() 0) { res.pop_back(); } BigInt c; c.d res; return c; }这段代码的复杂度是 O(len(a) * len(b))也就是两个大数位数乘积。后续所有进阶优化的目标就是想尽办法把这个乘法的复杂度压下去。2.3 朴素连乘的复杂度瓶颈到底在哪很多人以为朴素连乘慢是因为乘了 n 次其实真正的问题是“位数也在涨”。我们算 n! 的时候第 i 次乘法的操作对象是一个已经积累了大约 O(i log i) 位的大数。也就是说每乘一次大数的位数都在变长。总操作量大约是n 次乘法 × 平均位数 O(n log n) ≈ O(n² log n)举个直观的例子求 50000! 时50000! 大约有 21 万位。接近尾部的时候每一次乘法都要遍历大概 10 万到 20 万个数组元素而这样的乘法要做 5 万次。总操作量量级在几十亿次在普通 PC 上跑起来要几十秒甚至几分钟完全不能忍。所以“进阶题”的进阶点其实不是要不要用高精度而是怎么在 n 很大的情况下把高精度乘法做快。3. 三种优化思路的横向对比与原理拆解3.1 分治乘法把阶乘分成两半再合并朴素连乘是严格按顺序算1 × 2 × 3 × ... × n。分治乘法的想法很直接不按顺序而是把区间 [1, n] 切成两半先递归算出左半边的乘积和右半边的乘积最后把两个大数乘起来。写成伪代码大概是BigInt solve(int l, int r) { if (l r) return BigInt(l); int mid (l r) / 2; return multiplyBig(solve(l, mid), solve(mid 1, r)); }递归深度只有 O(log n) 层每一层做一次大数乘法。真正的收益在于乘法次数从 n 次降到 log n 次尽管每次乘法是两个大数相乘比“大数乘小数”更重但整体效率还是大幅提升。我实际测试时的体会是分治乘法对 50000! 到 100000! 的优化效果非常明显通常能把运行时间从十几秒压到一两秒量级。这里有一个实现细节不必真的递归到底层否则递归调用开销很大。更稳的做法是先把区间分成若干块每块内用朴素“大数乘小整数”算出来再用大数乘大数合并。常用的分块数是 4 到 16 块或者按线程数分块。3.2 FFT 和 NTT让乘法从 O(n²) 降到 O(n log n)如果 n 继续往上涨比如求 100 万阶乘甚至 1000 万阶乘分治乘法也扛不住。问题出在multiplyBig的 O(len(a) * len(b)) 上。两个 10 万位的数相乘暴力乘法要做 100 亿次基础运算依然很痛。这时就走 FFT 路线两个大整数相乘本质上可以看作两个数字序列的卷积。FFT快速傅里叶变换能在 O(n log n) 时间内算出卷积然后再处理进位。不过 FFT 有一个经典问题它用浮点数运算有精度误差。当大数位数达到几十万位时误差累积可能让最后某一位对不上。所以更严谨的选手会用 NTT快速数论变换也就是在整数模域里做同样的卷积完全避免浮点误差。NTT 写起来比 FFT 复杂不少需要选好模数还得会处理逆变换。竞赛里如果允许引入现成库最省心的方案是用 GMPGNU Multiple Precision Arithmetic Library它内部已经实现了大规模乘法的优化调度。你只需要把算法逻辑写好把乘法交给 GMP 即可。这里有一个比较真实的选型表可以帮你判断在什么场景下用哪种方案场景推荐方案复杂度量级n ≤ 10000 精确值朴素高精度连乘O(n² log n)10000 n ≤ 500000 精确值分治乘法O(n log² n) 量级n 500000 精确值分治 FFT/NTT或直接用 GMPO(n log² n) 但常数更小只求位数斯特林公式O(1)只求末尾零个数统计因子 5 的数量O(log n)只求对 M 取模边乘边模O(n)3.3 素数分解法绕过乘法次数问题的数学路径第三种思路是跳出“连乘”框架改用数论方法。n! 可以分解成质因数的幂次乘积。统计质数 p 在 n! 中出现的次数用的是经典的勒让德公式v_p(n!) floor(n/p) floor(n/p²) floor(n/p³) ...比如v_2(10!) 10/2 10/4 10/8 5 2 1 8所以 10! 中因子 2 的幂次是 8。于是算法变成用筛法得到 n 以内的所有质数对每个质数 p利用上面的公式计算指数 e_p最终结果就是所有 p^e_p 的乘积。这个算法的最大好处是把“连乘 n 次”变成“乘质数次幂的集合”质数的数量级大约是 n / ln n比 n 小得多。遇到需要“乘上自己指数的幂次”时还可以用快速幂优化。不过这个方案在实践里并不总是最快的单个体积小的质数直接乘一次或几次普通整数就行但遇到大质数的高次幂时依然要调用大数乘法。所以素数分解法更多是作为一种漂亮的数学解法存在在需要做取模或者分析结构性规律的题目中非常有用单纯求精确值时竞争力不如分治乘法。3.4 方案选型单纯求精确值怎么选性价比最高我自己总结的实际策略是如果 n ≤ 10000朴素高精度连乘就够代码短不容易出 bug完全没必要上分治。如果 10000 n ≤ 500000优先分治乘法。这个范围内它和 FFT 的性能差距不大但实现难度低很多。如果 n 再往上建议直接考虑系统大数库或者手写 NTT。另外Python 用户有天然优势。Python 内置的 int 本身就是任意精度而且底层实现混合了 Karatsuba、Toom-Cook、FFT 等多种乘法策略。像math.factorial(100000)这种调用Python 内部已经优化得很好有时候比新手手写 C 还快。但这也导致一个现象用 Python 刷简单题非常爽一旦面试要求手写高精度反而容易懵。所以理解底层原理永远不该跳过。4. 从算法到实战大阶乘题目的四类考法4.1 求精确值高精度乘法的主场这是最经典也最硬核的考法直接要求输出 n! 的每一位数字。对于数据范围在 10000 左右的小题直接使用第 2 节里朴素连乘的 BigInt 类就行。真正的坑在输出环节万进制存储时除了最高位需要正常输出其余每一块都必须补零到四位。比如[1, 20]表示的数是 20001也就是 1 在前20 在后输出成1 setw(4) setfill(0) 0020最终结果是10020少了前导零补位就会错。对于 n 达到 100000 以上的大范围分治乘法几乎是必选项。用前面分段分治的思路把区间先切成若干块各块内用朴素 multiply再逐层合并输出结果就能控制在一两秒内。4.2 求末尾零个数数因子就够了这种题不让你输出全值问你“100! 的末尾有多少个 0”。末尾零来自因子 10也就是因子 2 × 因子 5。由于阶乘里因子 2 的数量永远多于因子 5 的数量所以末尾零的个数只取决于因子 5 的个数。做法简单到令人发指int trailingZeros(int n) { int ans 0; while (n 0) { n / 5; ans n; } return ans; }100! 末尾零个数的计算过程是100 / 5 20100 / 25 4100 / 125 0总数是 20 4 24。这种题考察的核心不是高精度而是快速找到数学规律避免真的去算大阶乘。4.3 求位数斯特林公式的经典应用如果题目问“n! 有多少位”也不需要真的把阶乘算出来。对位数有一个现成公式digits(n!) floor(log10(n!)) 1直接算 n! 再取对数是不可能的所以要用斯特林公式近似n! ≈ sqrt(2πn) × (n/e)^n于是log10(n!) ≈ 0.5 × log10(2πn) n × log10(n/e)代码实现int digitsOfFact(int n) { if (n 1) return 1; double PI acos(-1.0); double E exp(1.0); double res 0.5 * log10(2.0 * PI * n) n * log10((double)n / E); return (int)res 1; }例如 n 10这个公式算出来约 6.559取整后加 1结果是 7。10! 是 3628800正好 7 位。n 100000 时算出来是 456574 位这也是我常用来自检验证分治乘法结果是否完整的参考值。4.4 求对 M 取模边乘边模的注意事项有些题不要求全值只要求“n! mod M”。这时完全不需要大数因为模运算满足(a × b) mod M ((a mod M) × (b mod M)) mod M所以从 1 乘到 n每一步都取一次模即可long long modFact(int n, long long mod) { long long res 1; for (int i 2; i n; i) { res (res * i) % mod; } return res; }这里有一个隐藏得很深的坑如果 mod 本身很大比如接近 10^18那么res * i这一步可能直接溢出 long long。解决方式是用__int128做中间乘法long long modFactSafe(int n, long long mod) { long long res 1; for (int i 2; i n; i) { res (long long)((__int128)res * i % mod); } return res; }还有一个竞赛中经常出现的小推论如果 M 是质数且 n ≥ M那么 n! 里包含因子 M所以 n! mod M 0。这个结论能帮你快速排除大多数情况。5. 实测踩坑记录与调试建议5.1 数组越界与进位丢失我排查了一晚上的经典错误第一次写高精度乘法时我最容易犯的错是“只处理了一次进位”。比如万进制下a[i] * b[j]的结果可能高达 9999 × 9999 99980001这中间可能产生 4 位进位。如果代码只在res[ij]上做一次取模和一次进位传递没有循环处理进位最终结果就会缺位。这类 bug 很隐蔽因为算小数字比如 10!、20! 时完全正常但一旦到 1000!结果的某一位就开始错乱。我的排查方法比较“土”但很有效先算 100! 的标准值网上能查到精确值然后把程序输出逐段对照特别是中间靠后的位置。发现错位后在 multiplyBig 内部打印每轮进位的值很快就定位到是进位循环没写完整。提示做高精度题时至少准备一个 100! 和 1000! 的精确值作为回归测试数据每次改完代码先跑这两个用例能避免很多低级错误。5.2 特殊输入问题0!、1!、以及负数阶乘很多实现里阶乘循环是从 2 开始的for (int i 2; i n; i) res.multiply(i);这个写法在 n 0 时也能正确输出 1因为循环一次都不执行res 的初始值正好是 1。但如果你把 res 初始化为 0或者在某处修改过全局状态0! 就很容易变成 0。还有人在读题时忽略了 0! 1 这个约定直接写if (n 0) return 0这种错误一旦出现在提交里就是 WA。负数阶乘在数学上不定义题目一般不会出现但防御性处理还是要做。最好在函数入口处判一下 n 0直接返回错误或者抛出异常避免后续访问数组时出诡异问题。5.3 内存与性能从 1 万到 10 万的实测体感我的测试机是一台普通笔记本数据大概是这样朴素连乘算 10000! 大约一百毫秒级别算 50000! 已经慢到能明显等待数秒算 100000! 基本要十几秒。分治乘法算 100000! 大概压到一秒内算 1000000! 也可以接受但内存占用开始明显上升。这里有个内存问题很容易被忽视一个大数在万进制下100000! 有 456574 位十进制数组长度大约是 11 万多个 int单独看并不大。但分治递归合并时如果每个递归层级都保存临时大数对象总内存会成倍上涨。解决方式是用迭代式分段合并或者及时释放不再使用的临时变量。5.4 调试大阶乘结果的技巧模数哈希校验法大数的精确值动辄几万位甚至几十万位肉眼根本看不出对不对。我推荐的校验方式是“多重模数哈希”用主算法算出 n! 的结果对结果分别取模几个不同的质数比如 1000000007 和 998244353用另一个可信实现比如 Python 的math.factorial算出同样的模值如果所有模值都一致基本可以认为大数乘法逻辑正确。这个技巧在工程里也很有用尤其是你修改了乘法内部实现之后用它能快速确认没有引入隐蔽的进位错误。另外性能调优时可以分段打印时间节点算到 10000、30000、50000、100000 各记录一次。这样能清楚看到瓶颈出现在哪个区间而不用等到整个程序跑完才能判断。最后再分享一个我个人的习惯做这类题前永远先确认数据范围。n 不超过 10000直接朴素高精度省时省力n 到 10 万以上分治乘法基本就是刚需。盲目套最复杂的算法除了增加 bug 风险并不会让你的交付更快。这也算是我刷过几十道高精度题之后最想提醒后来人的一句话。