ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

C++组合数计算全攻略:公式法、递推法与质因数分解实战

C++组合数计算全攻略:公式法、递推法与质因数分解实战 直接聊组合数吧。用C算组合数这件事看起来就是个公式套用C(n, m) n! / (m! * (n-m)!)一行代码的事但真等你把n调到20以上用long long一跑就懵了要么溢出要么慢得像蜗牛要么答案直接错得离谱。我见过太多人在这个基础问题上栽跟头所以这篇就把组合数的3种主流实现方式掰开揉碎讲清楚——什么时候用哪种、每种背后是什么原理、边界在哪、踩过哪些坑一次性给你说明白。这篇文章适合正在学C的初学者、准备算法竞赛的选手以及工作中临时要处理组合数计算的开发者。不需要你有太高深的基础只要能看懂数组和循环就能跟着代码走一遍顺便把背后的数学原理和工程取舍也弄懂。1. 为什么组合数看起来简单真写起来却容易翻车1.1 先弄清要解决什么问题组合数在数学上定义为从n个不同元素中选取m个元素的方案数记作C(n, m)或nCm计算公式是C(n, m) n! / (m! * (n-m)!)这个公式谁都会背但放到C里就有三个麻烦第一阶乘膨胀速度极快。n20时n!大约是2.43e18已经逼近long long的上限9.22e18n21时n!就直接爆了。也就是说你还没算到最终的组合数中间过程就先溢出了。第二直接除不尽。如果你按公式从左到右老老实实算“n!除以m!再除以(n-m)!”那中间结果天生就是极大的整数很可能早在除法发生之前就已经溢出。第三组合数的取值本身也可能超出基础类型范围。比如C(100, 50)大约等于1.0089e29这已经不是long long能装得下的了得用高精度。所以组合数计算根本不是“套公式”的活而是“如何在避免中间爆炸的前提下正确得到最终结果”的工程问题。方法不同适用边界完全不同这也是这篇要讲三种方法的核心原因。1.2 三种实现思路的定位与选型结论我把这三种方法提前摆出来方便你脑子里先有个地图方法一乘除化简法公式法优化版。通过边乘边约分、及时约分来避免中间结果过大优点是代码直观、容易理解适合n不超过60、结果在64位整数范围内的场合。方法二递推法帕斯卡公式/杨辉三角。利用C(n, m) C(n-1, m-1) C(n-1, m)做动态规划优点是可以在取模条件下使用数值稳定适合需要批量计算组合数、且n在几千甚至几万级别的场景。方法三质因数分解法高精度实现。把组合数分解成质因数的乘积再用高精度乘法把结果算出来优点是彻底解决溢出问题n能到几千甚至更大适合需要精确结果的场景。怎么选我的经验是先看结果范围再看是否有取模需求最后看是否只算一次。如果只算一次C(50, 25)这种方法一最省事如果要算一堆组合数做成表方法二最好使如果n到了100以上还要求精确值老老实实上方法三别硬撑。2. 方法一乘除化简公式最直接也最限制明确2.1 从数学公式到第一版代码先写最直观的版本注意看它哪里会出问题#include iostream long long factorial(int n) { long long res 1; for (int i 2; i n; i) res * i; return res; } long long combinationBad(int n, int m) { if (m 0 || m n) return 0; return factorial(n) / (factorial(m) * factorial(n - m)); }这段代码在n20以内还能跑n21直接溢出。而且哪怕n和m本身没让最终结果出界中间的factorial(n)也已经爆了。阶乘是超级膨胀的中间形态能避免就应该尽量避免。优化思路是把公式改写一下让中间值尽量小。利用组合数的对称性先把m改成n-m中较小的那个再通过公式C(n, m) n * (n-1) * ... * (n-m1) / m!这种写法把计算量从三次阶乘简化成一次连乘加一次阶乘。但问题依然存在分子的连乘依然可能很大而且分子除以分母时不一定整除到哪里才能除。最简单粗暴的做法是算完再除但分子在n较大时依然会爆。2.2 乘除交叉进行让溢出概率直线下降第二个优化版本就是边乘边除每乘一个分子因子就尝试除以分母因子尽量控制中间值大小long long combinationGood(int n, int m) { if (m 0 || m n) return 0; if (m n - m) m n - m; long long result 1; for (int i 1; i m; i) { result result * (n - m i) / i; } return result; }这段代码写得非常紧凑逻辑也很经典循环i从1到m分子从n-m1依次乘到n分母从1依次乘到m乘一个除一个。为什么可以保证整除因为连续m个整数相乘后一定能被m!整除这是数学上早就证明的性质。实际执行时每步除法都能整除不会出现小数。但别高兴太早这种交叉乘除只是让中间结果的规模接近最终结果而不是无限制地小。比如算C(100, 50)最终结果约1e29中间结果在最后几步会涨到trillion级别最终还是超过long long上限。所以这个方法的使用边界大致是最终结果不超过long long最大值也就是C(60, 30)约为1.18e17可以安全通过C(67, 33)约为1.42e19就已经超了。再补充一个更稳的变体每步先做gcd约分再相乘能把溢出风险再压一档long long combinationStable(int n, int m) { if (m 0 || m n) return 0; if (m n - m) m n - m; long long result 1; for (int i 1; i m; i) { long long numerator n - m i; long long denominator i; long long g std::gcd(result, denominator); result / g; denominator / g; g std::gcd(numerator, denominator); numerator / g; denominator / g; result result * numerator / denominator; } return result; }每一步先把当前结果与要除的数约分再把要乘的数和剩余分母约分这能让中间值始终贴近最终结果。代价是多调用几次gcd但换来的是安全范围更大实际跑下来非常稳。这个变体我强烈推荐作为“公式法”的默认实现。3. 方法二递推法算法竞赛里的常青树3.1 帕斯卡公式与二维表构建递推法的核心是帕斯卡恒等式C(n, m) C(n-1, m-1) C(n-1, m)边界条件是C(n, 0) C(n, n) 1。这个式子用动态规划填一张二维表就能高效算出一堆组合数时间复杂度O(nm)空间复杂度O(nm)。#include vector long long combinationDP(int n, int m) { if (m 0 || m n) return 0; std::vectorstd::vectorlong long C(n 1, std::vectorlong long(m 1, 0)); for (int i 0; i n; i) { C[i][0] 1; for (int j 1; j std::min(i, m); j) { C[i][j] C[i-1][j-1] C[i-1][j]; } } return C[n][m]; }填表过程其实就是把杨辉三角按行算一遍每个格子只依赖上一行的两个格子思路很直白。这个方法的优势在于只要每一格的值不溢出整张表就算得出来。它不经过阶乘中间值比方法一更友好在n65左右、结果不超过long long范围时非常可靠。3.2 滚动数组压缩空间配合取模发挥最大值二维表在n变大的时候空间有点浪费因为算第i行只需要第i-1行的数据。改成滚动数组后空间降到O(m)long long combinationRolling(int n, int m) { if (m 0 || m n) return 0; if (m n - m) m n - m; std::vectorlong long dp(m 1, 0); dp[0] 1; for (int i 1; i n; i) { // 关键点必须倒着更新确保用的是上一行的旧值 for (int j std::min(i, m); j 1; --j) { dp[j] dp[j] dp[j-1]; } } return dp[m]; }这里有一个特别容易踩的坑内层循环必须倒着来。如果正着更新dp[j-1]已经被这一轮循环改过了用的就不是上一行的值算出来的就是错的结果。我第一次写的时候就是正着循环结果C(5, 2)算出来10看着像对其实只是数字碰巧对换几个参数就露馅。递推法最强大的场景是配合取模运算。普通整数加法容易溢出但取模后就可以放心大胆地用。经典的应用是在模素数p下计算组合数const long long MOD 1000000007LL; long long combinationMod(int n, int m) { if (m 0 || m n) return 0; std::vectorlong long dp(m 1, 0); dp[0] 1; for (int i 1; i n; i) { for (int j std::min(i, m); j 1; --j) { dp[j] (dp[j] dp[j-1]) % MOD; } } return dp[m]; }这个模板在n达到几千、几万时依然秒出结果因为总计算量是n*m取模运算稍微慢一点但完全可以接受。比赛里常见的组合数取模、多项式系数、概率DP等题目基本都是这个套路。4. 方法三质因数分解高精度乘法治大数顽疾4.1 核心原理用质数乘积表示组合数当n超过67精确值超过long long范围时前面两种方法都无能为力但问题还得解。这时候就用到了算术基本定理任何一个正整数都能唯一分解成质因数幂的乘积。组合数C(n, m)当然也是整数所以也能写成C(n, m) p1^e1 * p2^e2 * ... * pk^ek问题就变成怎么求每个质数p在C(n, m)里的指数e这里用到一个经典结论——质数p在n!中的指数是e(n!) floor(n/p) floor(n/p^2) floor(n/p^3) ...这个公式也叫勒让德定理。直观理解就是从1乘到n先数一遍有多少个数是p的倍数每个贡献1个p再看有多少个数是p^2的倍数再贡献1个p以此类推。于是e(C(n, m)) e(n!) - e(m!) - e((n-m)!)把每个质数的指数算出来最后把所有质数的指数幂乘起来就得到了精确结果。这里的“乘”是大数乘法因为结果本身已经超出基础数据类型必须用高精度手段。4.2 完整代码实现与细节我先给出一份完整实现然后逐段解释#include iostream #include vector #include string #include algorithm // 欧拉筛返回[1, n]范围内的所有质数 std::vectorint sievePrimes(int n) { std::vectorbool isPrime(n 1, true); std::vectorint primes; for (int i 2; i n; i) { if (isPrime[i]) { primes.push_back(i); if ((long long)i * i n) { for (long long j (long long)i * i; j n; j i) { isPrime[j] false; } } } } return primes; } // 计算质数p在n!中的指数 long long exponentInFactorial(int n, int p) { long long exp 0; while (n) { n / p; exp n; } return exp; } // 高精度乘法把num这个大数用vectorint逆序存储乘以x void multiplyByInt(std::vectorint num, int x) { int carry 0; for (size_t i 0; i num.size(); i) { int cur num[i] * x carry; num[i] cur % 10; carry cur / 10; } while (carry) { num.push_back(carry % 10); carry / 10; } } // 组合数精确值返回字符串 std::string combinationExact(int n, int m) { if (m 0 || m n) return 0; if (m n - m) m n - m; std::vectorint primes sievePrimes(n); std::vectorlong long exponents; for (int p : primes) { long long e exponentInFactorial(n, p) - exponentInFactorial(m, p) - exponentInFactorial(n - m, p); if (e 0) exponents.push_back(e); else exponents.push_back(0); } std::vectorint result(1, 1); for (size_t i 0; i primes.size(); i) { if (exponents[i] 0) continue; for (long long j 0; j exponents[i]; j) { multiplyByInt(result, primes[i]); } } std::string s; for (auto it result.rbegin(); it ! result.rend(); it) { s.push_back(char(0 *it)); } return s; }逐个说几个容易错的地方筛质数时注意标记循环的边界。for (long long j (long long)i * i; j n; j i)里的(long long)i * i先转长整型防止i*i在int范围内溢出。n到几千几万时没感觉但万一n到了一百万以上这个细节能保命。指数计算里我在函数开头对m做了对称化处理把m取成较小值。这本身不影响最终结果但能减少后面高精度乘法的次数算是一种小优化。高精度乘法这里我用的是一位一位存int的vector从低位开始存这样处理进位最方便。每个质因子的指数可能很大比如C(1000, 500)里质数2的指数有几百甚至上千所以循环里反复乘同一个质数没问题但整体乘法次数会比较多n越大越耗时。4.3 一个更高效的高精度加速思路上面这个版本在n1000以内跑得飞快但n到几千时重复乘以同一个质数几百次的代价就开始显现。想再快的话可以对每个质数做“快速幂”然后高精度乘以一个“很大的数”。但C标准库没有现成的高精度大数乘法自己写FFT快速傅里叶变换又有点过度设计。如果你只是平时做项目或刷题这个版本已经完全够用。我在实际测试里n2000、m1000的组合数大概几百位这个程序能在几十毫秒内算完完全能接受。真要算C(100000, 50000)这种几千位的大数那得专门优化高精度乘法建议直接用Python或GMP库别跟自己过不去。5. 三种方法横向对比怎么选才不后悔5.1 复杂度与溢出风险对照表我把三种方法的关键参数整理成一张表方便你直接对照方法时间复杂度额外空间结果范围主要优点主要限制乘除化简法O(m)O(1)≤ C(67,33)左右约1e19以内代码短、好理解、适合单次计算中间值接近结果仍可能溢出gcd稳定法O(m log n)O(1)同上更稳中间值更小多几次gcd调用范围没本质提升递推法O(n*m)O(m)与滚动数组存储类型一致可批量计算、可配合取模单次查询不如公式法快质因数分解法O(n log log n 结果位数)O(n)无上限精确大数彻底解决溢出实现复杂度高需要筛质数和手写高精度需要多说一句的是递推法的时间复杂度O(n*m)看起来挺高但如果n是几千、m也是几千乘积就是百万量级现代CPU毫秒级搞定根本不叫事。反而是方法一的O(m)只在单次计算时有优势真要算C(1,0)到C(n,m)一整张表递推法一次填表全部解决。5.2 我实际测试过的数据和选型经验我在本地用release模式跑过几组典型数据贴出来给你参考耗时数据因机器而异但比例是可靠的计算目标方法一乘除化简方法二递推方法三质因数分解C(30, 15)微秒级微秒级微秒级C(60, 30)微秒级但需要gcd版本微秒级微秒级C(100, 50)直接溢出溢出毫秒级结果30位C(1000, 500)失败失败约10毫秒结果约300位C(5000, 2500)失败失败约几百毫秒结果约1500位从表里能直观看到方法一和方法二本质是同一条船上的都受限于long long的存储范围。如果结果本身能装进long long这两个方法随便挑如果装不下那就只有方法三能打。选型上我再给几个实战建议只算一次且n60直接上gcd稳定版乘除化简法代码短不易错。需要在模p下算一堆组合数递推法是王者滚动数组加取模内存占用小速度也够快。p是素数还能配Lucas定理走得更远这里先不展开。要求精确大数结果质因数分解法。不建议再尝试用double算大组合数然后四舍五入——浮点误差在30位数字面前根本无法接受我见过有人用double算C(100, 50)得到1.008913445455642e29看着像模像样但最后几位全是错的。6. 常见问题与排查技巧实录6.1 我踩过的溢出、越界和边界坑问题1m n 没处理这个最基础但真有人会漏。公式里C(n,m)当mn时数学上是0但代码里如果不先做判断factorial(n-m)会去算负数的阶乘直接死循环或返回垃圾值。所有实现都要在函数开头加一句if (m 0 || m n) return 0;问题2忘了用对称性优化C(n, m) C(n, n-m)。这个性质不只是为了省几轮循环更重要的是能显著缩小中间值。比如C(100, 98)如果不做对称化公式法要连乘98项中间结果接近最终结果C(100, 2)4950的阶乘级中间体做了对称化只算C(100, 2)分子只有两项中间体最多100*999900完全不是一个量级。问题3递推法内层循环正着写前面强调过滚动数组必须倒着更新。正着写的后果是数据被污染而且某些参数下结果恰好是错的。排查时可以在算完后用对称性C(n,m)C(n,n-m)交叉验证一遍或者拿小n手算核对。问题4把long long当万能药到n21就爆很多人最初以为“用long long就安全了”。实际上long long上限约9.22e18而C(100, 50)约1e29差了10个数量级。我之前带新人的时候有人用long long去算C(100, 50)得到负数也不意外这是有符号整型的溢出特征。遇到这种情况要么检查是不是该换方法三要么确认自己能否接受取模。问题5质因数分解法中指数为负当你的指数计算代码写成exponentInFactorial(n,p) - exponentInFactorial(m,p) - exponentInFactorial(n-m,p)时理论上不可能为负因为C(n,m)一定是整数。但如果你在调用前没做mn的检查或者m在对称化前用了原始值导致n-m为负指数就可能算成负数。所以指数算出来小于等于0时直接跳过比较稳妥。6.2 组合数取模场景的补充模板最后再补充一个竞赛里高频出现的变体模p为素数、n特别大比如n达到1e18这时不能用递推法需要用Lucas定理。Lucas定理的内容是C(n, m) % p C(n/p, m/p) * C(n%p, m%p) % p递归处理即可long long modPow(long long a, long long b, long long p) { long long res 1; while (b) { if (b 1) res res * a % p; a a * a % p; b 1; } return res; } long long lucas(long long n, long long m, long long p) { if (m 0) return 1; long long ni n % p, mi m % p; if (mi ni) return 0; // 这里可以用预处理的阶乘和逆元快速计算C(ni, mi) % p // 注意需要保证 p 是素数 return lucas(n / p, m / p, p) * combSmall(ni, mi, p) % p; }这个扩展不是必须掌握的但当你在编程题里遇到超大n的取模组合数时它就是正解。它的原理是把大n分解成p进制下的若干位逐位套用小范围内的组合数取模用小环境替代大环境既绕开了溢出又保证了效率。回到最初那句话组合数计算真正考验的不是公式记忆而是对“中间结果膨胀”的警惕和对场景的判断。我的习惯是任何一段计算代码先问自己三个问题——结果会不会溢出中间过程会不会溢出需不需要取模把这三个问题想清楚方法自然而然就选对了。希望这篇实战笔记能帮你少踩几个坑也欢迎你拿自己的测试数据来验证这几种方法的边界——理论说得再多都不如亲手跑一遍来得直观。
返回列表