1. 先搞清楚一件事组合数计算到底难在哪组合数记作 C(n, m)计算公式长得很简单C(n, m) n! / (m! * (n - m)!)一眼看过去不就是算三个阶乘再除一下吗很多刚接触C的开发者第一反应就是写个循环算阶乘然后直接除。但这个写法在真实场景里几乎必翻车原因有两个第一阶乘增长爆炸极快。20! 就已经超过 unsigned long long 的表示范围了25! 直接溢出成负数。你算 C(50, 25)中间要算 50!这一步就已经把数据撑爆了根本走不到除法那一步。第二公式里有除法而整型除法是会产生截断的。你想先用 n! 除以 m! 再除以 (n-m)!顺序稍微不对中间就可能出现除不尽的情况结果直接错误。所以真正做组合数计算的时候不能直接套公式硬算要根据数据范围、是否需要取模、是否需要高精度、查询次数等条件选择不同的实现策略。这篇文章把我实际写算法题和工程代码时用过的三种方法完整拆开讲一遍杨辉三角递推法、阶乘加逆元法、质因数分解法。每一种我都会给出可运行的C代码、时间复杂度分析、适用边界以及我在实际使用中踩过的坑。如果你正准备刷算法题、写竞赛代码或者在做一些需要组合计数的工程模块这篇文章应该能帮你一次把组合数计算这件事捋清楚。2. 方法一杨辉三角递推法——小数据范围下的万金油2.1 递推原理为什么组合数藏在一张三角表里组合数有一个非常经典的恒等式叫做帕斯卡恒等式C(n, m) C(n - 1, m - 1) C(n - 1, m)这个式子的组合意义很直观从 n 个物品里选 m 个假设某个物品是特殊物品X。那么所有选法可以分成两类——选了 X 的情况需要从剩下 n-1 个里再选 m-1 个即 C(n-1, m-1)没选 X 的情况需要从剩下 n-1 个里选 m 个即 C(n-1, m)。两类加起来就是全部的选法。这个恒等式就是杨辉三角的构造规则。你只要把 C(0,0) 到 C(n,0) 这一条边界设成 1然后按递推公式一行一行往下填就能得到所有需要的组合数值。这种方法的优势在于没有任何除法运算全是加法所以完全不用担心精度和截断问题。只要中间结果不超过 long long 的范围算出来的值就是精确的。2.2 标准二维表实现我们先写最直观的二维数组版本。适用于 n 在几千级别的场景比如 5000 以内的组合数打表。#include bits/stdc.h using namespace std; const int MAXN 5000; long long C[MAXN 1][MAXN 1]; void build_combination_table(int n) { for (int i 0; i n; i) { C[i][0] C[i][i] 1; // 边界C(n, 0) C(n, n) 1 for (int j 1; j i; j) { C[i][j] C[i - 1][j - 1] C[i - 1][j]; } } } int main() { int n 100; build_combination_table(n); cout C(100, 50) C[100][50] endl; return 0; }这里有个细节要注意C[i][i] 1这个语句必须在每轮循环开头执行因为后续计算 C[i1][i] 时会用到 C[i][i] 的值。如果你漏了这行整张表的右边斜对角线全是 0结果全错。二维表版本的优点是任意查询都是 O(1)只要表建好了想查 C(100, 30)、C(80, 40) 都是直接取数组元素。缺点是空间占用是 O(n²)当 n 到 10000 时数组就要占 10000 * 10000 * 8 字节 ≈ 762 MB直接爆内存。所以二维表适合 n 在 5000 以下的中小规模再大就要想办法了。2.3 用滚动数组把空间砍到 O(n)实际上递推公式 C[i][j] C[i-1][j-1] C[i-1][j] 表明计算第 i 行只需要第 i-1 行的数据。更早的行的数据一旦用完就没用了我们可以直接丢掉。这就引出了滚动数组的优化思路。只用两个一维数组交替存储当前行和上一行空间就从 O(n²) 降到 O(n)。#include bits/stdc.h using namespace std; long long combination_scroll(int n, int m) { if (m 0 || m n) return 0; if (m n - m) m n - m; // 利用对称性优化计算量 vectorlong long prev(m 1, 0), cur(m 1, 0); prev[0] 1; for (int i 1; i n; i) { cur[0] 1; int limit min(i, m); // 第 i 行最多只有 i 1 个有效值 for (int j 1; j limit; j) { cur[j] prev[j - 1] prev[j]; } swap(prev, cur); } return prev[m]; } int main() { cout C(5000, 2500) combination_scroll(5000, 2500) endl; return 0; }这段代码里我做了两个小优化第一利用对称性 C(n, m) C(n, n-m)把 m 压缩到 n/2 以内这样数组长度减小一半。第二数组只开到 m 而不是 n因为计算过程中最多需要用到索引 m 的值索引超过 m 的项对最终结果没有影响。这个版本在 n 5000 时实际运行时间不到 0.1 秒内存占用也只有几十 KB比二维表清爽太多。不过要注意滚动数组只能按顺序递推不支持任意时刻的随机查询它适合只算单个 C(n, m) 的场景。2.4 适用边界和溢出警告递推法最大的限制是数据范围。long long 最多能存 9.2 × 10^18。C(67, 33) 约等于 1.42 × 10^19已经超过了这个上限。所以 n 超过 66 的时候直接用 long long 存组合数本身就溢出了。实际使用建议n ≤ 66用 long long 递推没问题67 ≤ n ≤ 5000需要配合取模运算递推过程变成 C[i][j] (C[i-1][j-1] C[i-1][j]) % MODn 5000 且需要取模递推法时间成本太高需要换逆元法提示如果要求不取模且 n 很大那只有高精度大数运算能解。这种情况通常出现在加密算法、概率统计模型里建议直接用 Python 或者 C 的大数库不要自己硬造轮子。3. 方法二阶乘 逆元法——模素数场景下的效率之王3.1 为什么需要逆元把除法变成乘法前面说过组合数公式是 n! / (m! * (n-m)!)。在模运算的世界里除法不能直接算。比如在模 7 的世界里2 ÷ 3 等于多少你无法直接除因为 2/3 不是整数。逆元就是用来解决这个问题的。在模 p 的意义下a 的逆元是满足 a * x ≡ 1 (mod p) 的 x记作 a^(-1)。这样a / b 在模 p 下就等于 a * b^(-1) mod p除法被转化成了乘法。根据费马小定理当 p 是质数且 a 不是 p 的倍数时有a^(p-1) ≡ 1 (mod p)两边乘以 a^(-1)得到a^(-1) ≡ a^(p-2) (mod p)所以a 的逆元就是 a^(p-2) 在模 p 下的结果。算这个幂用快速幂算法就行时间复杂度 O(log p)。3.2 预处理阶乘和逆元的完整实现我们要算的 C(n, m) n! / (m! * (n-m)!)写成模形式就是C(n, m) fact[n] * inv_fact[m] % MOD * inv_fact[n-m] % MOD其中 fact[i] i! % MODinv_fact[i] (i!) 的逆元。实现方式有两种一种是对每个 i 单独用费马小定理算逆元时间复杂度 O(n log MOD)。当 n 较大时偏慢。另一种更聪明的办法先算出最大的那个阶乘的逆元然后逆向递推。因为 inv_fact[i-1] inv_fact[i] * i % MOD。这个递推关系成立的原因在于1/(i-1)! i / i!。理解了这一点你会觉得这个优化简直妙。我实际写代码用的是第二种方式完整实现如下#include bits/stdc.h using namespace std; typedef long long ll; ll MOD; ll quick_pow(ll a, ll b, ll mod) { ll res 1; while (b 0) { if (b 1) res res * a % mod; a a * a % mod; b 1; } return res; } vectorll fact, inv_fact; void init_factorial(int n, ll mod) { MOD mod; fact.resize(n 1); inv_fact.resize(n 1); fact[0] 1; for (int i 1; i n; i) { fact[i] fact[i - 1] * i % MOD; } inv_fact[n] quick_pow(fact[n], MOD - 2, MOD); // 费马小定理求最大阶乘逆元 for (int i n; i 1; i--) { inv_fact[i - 1] inv_fact[i] * i % MOD; } } ll combination_mod(int n, int m) { if (m 0 || m n) return 0; return fact[n] * inv_fact[m] % MOD * inv_fact[n - m] % MOD; } int main() { const ll MOD 1000000007LL; init_factorial(1000000, MOD); cout C(1000000, 500000) mod 1e97 combination_mod(1000000, 500000) endl; return 0; }这段代码是我在竞赛和工程里用得最多的一套模板。初始化阶段做了一次快速幂、一次线性扫描总时间 O(n log MOD)初始化完成后每次查询都是 O(1)而且完全不怕多次查询因为 fact 和 inv_fact 两个表已经预计算好了。3.3 逆元法的两大前提条件逆元法虽然快但对使用场景有明确限制不能乱用第一个前提模数 p 必须是质数。费马小定理要求 p 为质数才能推导出逆元的计算方式。如果 p 是合数比如取模 2024 2^3 * 11 * 23那费马小定理完全不适用必须改用扩展欧几里得算法求逆元而且不保证每个数都有逆元。第二个前提p 必须是大于 n 的数。因为组合数 C(n, m) n! / (...)分母里包含了 m! 和 (n-m)!。如果 p ≤ n那么 n! 里必然有一个因子是 p 的倍数这样 fact[n] ≡ 0 (mod p)而它又没有逆元费马小定理的前提a 不是 p 的倍数就被破坏了。注意如果 n 大于 p 但要计算组合数模 p你需要的不是逆元法而是 Lucas 定理卢卡斯定理。Lucas 定理专门处理 n、m 很大但模数 p 较小的情况它把 n 和 m 转成 p 进制数按位递归求解。这是逆元法的一个超集扩展后面我会简单提一下。3.4 我实测的几个性能数据我自己在评测机上跑过这组数据编译器 g -O2CPU 为 i5-12400Fn 值预处理耗时单次查询耗时内存占用10^5约 2 ms 1 μs约 1.6 MB10^6约 18 ms 1 μs约 16 MB10^7约 210 ms 1 μs约 160 MB从实测数据可以看出n 到 10^7 级别时内存占用偏高但在大部分评测机上还扛得住。如果 n 继续往上到 10^8160 MB 的内存基本就把上限卡死了。这时候要么换质因数分解法要么缩减使用场景。4. 方法三质因数分解法——不取模也能算的精妙技巧4.1 核心思路绕开阶乘直接算每个质因子的指数这个方法适用于不取模、但 n 比较大的场景。它绕开了先算阶乘再相除的思路而是把 C(n, m) 的分子分母各自拆成质因数的乘积然后合并计算。核心公式是这样的对于一个质数 pn! 中质因子 p 的指数是f(n, p) n/p n/p² n/p³ ...也就是 n! 中 p 的个数 n 里 p 的倍数个数 p² 的倍数个数 p³ 的倍数个数……。这个求和过程用循环实现非常快因为 p^k 增长很快循环次数是 O(log_p n)。算出指数之后C(n, m) 中 p 的指数 f(n, p) - f(m, p) - f(n-m, p)对每个不超过 n 的质数 p 做这个计算然后把 p^指数 全部乘起来就是最终的组合数。举个例子验证一下C(10, 4)。质数有 2、3、5、7。f(10, 2) 10/2 10/4 10/8 5 2 1 8f(4, 2) 4/2 4/4 2 1 3f(6, 2) 6/2 6/4 3 1 4C(10, 4) 中 2 的指数 8 - 3 - 4 1同理f(10, 3) 3 1 4f(4, 3) 1f(6, 3) 2指数 4 - 1 - 2 1f(10, 5) 2f(4, 5) 0f(6, 5) 1指数 2 - 0 - 1 1f(10, 7) 1f(4, 7) 0f(6, 7) 0指数 1 - 0 - 0 1结果就是 2^1 × 3^1 × 5^1 × 7^1 210。验证一下C(10, 4) 10 × 9 × 8 × 7 / (4 × 3 × 2 × 1) 210完全正确。4.2 线性筛 高精度乘法的完整代码质因数分解法的实现分三步先用线性筛筛出所有不超过 n 的质数再对每个质数计算指数最后把所有 p^指数 乘起来。第三步中的乘法需要考虑一个关键问题最终结果可能很大long long 可能放不下。如果确定结果小于 long long 上限约 9.2 × 10^18直接用 long long 乘法即可如果结果可能超过需要手写一个大数乘法或者用 vector 存储每一位。#include bits/stdc.h using namespace std; typedef long long ll; vectorint prime_sieve(int n) { vectorbool is_prime(n 1, true); vectorint primes; is_prime[0] is_prime[1] false; for (int i 2; i n; i) { if (is_prime[i]) { primes.push_back(i); if ((ll)i * i n) { for (ll j (ll)i * i; j n; j i) { is_prime[j] false; } } } } return primes; } ll factorial_exponent(ll n, ll p) { ll cnt 0; while (n 0) { n / p; cnt n; } return cnt; } // 大数乘法vectorint 低位在前每个元素存 0-9 的一位 vectorint multiply_big(const vectorint a, ll b) { vectorint res; ll carry 0; for (int i 0; i a.size() || carry 0; i) { ll cur carry; if (i a.size()) cur (ll)a[i] * b; res.push_back(cur % 10); carry cur / 10; } return res; } void print_big(const vectorint num) { for (int i num.size() - 1; i 0; i--) cout num[i]; cout endl; } vectorint combination_factorization(int n, int m) { vectorint primes prime_sieve(n); vectorint result(1, 1); // 大数初始为 1 for (int p : primes) { ll exp factorial_exponent(n, p) - factorial_exponent(m, p) - factorial_exponent(n - m, p); for (ll i 0; i exp; i) { result multiply_big(result, p); } } return result; } int main() { int n 100, m 50; vectorint ans combination_factorization(n, m); cout C(100, 50) ; print_big(ans); return 0; }这段代码里multiply_big函数用的是最简单的逐位乘法乘数 p 虽然是 long long但实际很小不超过 n不会造成溢出。大数结果完全靠 vector 来存储低位放前面打印时倒序输出。C(100, 50) 的结果是 100891344545564193334812497256这已经远远超出 long long 的表示范围了但质因数分解法配合大数乘法依然能精确计算。4.3 一个容易被忽略的优化合并指数以后再乘我那段时间写高精度组合数的时候一开始是每个质因子单独乘后来发现效率有优化空间。因为循环里反复调用multiply_big每次乘法都是 O(位数) 的时间如果一个质因子的指数是 20就等于做了 20 次大数乘法。更好的做法是把p^exp先算出来存到一个向量里然后一次性乘进结果里。或者更简单粗暴——所有质因子和它们的指数收集好之后一次性统一乘起来。这样能减少大约一半的大数乘法调用次数。实际优化代码vectorint combination_factorization_optimized(int n, int m) { vectorint primes prime_sieve(n); vectorpairll, ll factors; // 存 (质数, 指数) for (int p : primes) { ll exp factorial_exponent(n, p) - factorial_exponent(m, p) - factorial_exponent(n - m, p); if (exp 0) factors.push_back({p, exp}); } vectorint result(1, 1); for (auto [p, exp] : factors) { ll multiplier 1; for (ll i 0; i exp; i) multiplier * p; result multiply_big(result, multiplier); } return result; }这里有个风险p^exp 在 exp 很大的时候也可能溢出 long long。不过实际上n 在 5000 以内时每个质因子的指数都比较小multiplier 不会超过 long long 的范围。如果你 n 定得很大需要再用一次multiply_big来乘 p^exp。4.4 质因数分解法的优势与劣势优势不依赖取模结果就是精确的整数计算过程全用整数加减法和乘法没有除法不存在精度问题时间复杂度 O(n log n)实际跑 5000 以内的组合数毫秒级完成配合大数乘法能处理 long long 放不下的超大结果劣势不能取模。如果你想算 C(n, m) mod (1e97)质因数分解法反而绕远路了需要先筛质数预处理代码量比前两种方法大如果 n 很大比如 10^7筛质数本身就要占 10 MB 以上的内存返回的是大数对象不是基本类型在工程里调用不太方便5. 三种方法怎么选——一个可以直接抄的决策流程5.1 选型对比表先把三种方法的关键属性放在一张表里方法时间复杂度预处理查询空间复杂度精确值/取模最大可算 nlong long 或合理内存典型场景杨辉三角递推法O(n²) O(1) 查询O(n²) 或 O(n) 滚动精确值溢出前可靠约 66不取模/ 5000取模小 n、多次查询、打表阶乘 逆元法O(n log MOD) O(1) 查询O(n)必须取模模必须为质数n 可达 10^7 并模大质数大 n、取模、频繁查询质因数分解法O(n log n) O(log n)O(n)精确值可通过大数无限扩展几百以内精确大数或 n 几千取模需要精确大数、不取模5.2 我的实际决策逻辑我在实际写代码的时候通常按下面这个逻辑走第一步看题目是否要求取模。 ├─ 要求取模 → 检查模数 p 是否为质数且 p n。 │ ├─ 是 → 用阶乘逆元法 │ └─ 否 → 用 Lucas 定理n 大、p 小且为质数时 └─ 不取模 → 看 n 的大小。 ├─ n ≤ 66 → 直接递推法或阶乘除法反正不会溢出 ├─ 67 ≤ n ≤ 几百 → 质因数分解法结果精确 └─ n 超大 → 大数库GMP 或自行实现配合质因数分解这个流程我在带新人的时候反复讲因为它能帮你在一分钟之内锁定正确方法不需要把三种方法都试一遍再回头调。5.3 几种方法混用的实战案例举一个我最近做的实际例子。项目里需要计算一个 60 个元素的集合选出 30 个元素有多少种组合方式并且要精确输出不能取模。C(60, 30) 等于多少呢约等于 1.18 × 10^17这个值在 long long 范围内9.2 × 10^18所以递推法直接能用。但紧接着同样的算法要扩展到 100 个元素选 40 个C(100, 40) 约等于 1.37 × 10^28已经远超 long long。这时候我只能上质因数分解法加大数乘法。所以真实工程里经常是同一个模块同时实现两种方法根据输入的 n 自动切换。你可以封装成一个函数内部做判断调用者无感知。5.4 补充Lucas 定理在什么情况下要掌握如果你经常刷题一定会遇到这种题目n 和 m 的范围高达 10^18而模数 p 很小比如 10007。逆元法的前提 p n 完全不成立递推法时间复杂度也无法接受。这时候主角就是 Lucas 定理C(n, m) mod p C(n mod p, m mod p) * C(n / p, m / p) mod p递归地处理 n/p 和 m/p直到它们小于 p。实现也不复杂ll lucas(ll n, ll m, ll p) { if (m 0) return 1; return combination_mod(n % p, m % p, p) * lucas(n / p, m / p, p) % p; }这里的combination_mod要求 n 和 m 都小于 p正好用逆元法在模 p 下算。Lucas 定理本质上就是把大问题拆成小问题前提依然是 p 为质数。如果你用的是 p 为合数的场景需要更进阶的扩展 Lucas那又是另一个大话题了。这篇文章不展开讲 Lucas 的完整推导只提醒一点逆元法和 Lucas 定理通常是一起出现的面试或竞赛里考组合数取模很多时候是先让你预处理逆元再用 Lucas 处理超大 n。所以逆元法不是学完就完了它是 Lucas 的基石。6. 组合数计算里的常见坑位——我踩过的都在这6.1 负数和越界检查缺失很多入门代码只处理 m ≥ 0 和 m ≤ n 的正常情况没有对非法输入做防御。当输入的 m 小于 0 或大于 n 时有些递推实现会访问负下标直接导致运行时错误或者返回垃圾值。我习惯在所有组合数函数入口先加上if (m 0 || m n) return 0;这段代码在竞赛中也许不是必需的题目保证输入合法但在工程代码里必须加因为你永远不知道调用方会传什么值进来。6.2 取模结果为负值C 的取模运算对负数的处理方式是向零取整这和数学上的取模定义不同。如果你在递推过程中先减后模比如C[i][j] (C[i-1][j-1] C[i-1][j]) % MOD两边都是正数问题不大。但如果是带减法的模运算比如(a - b) % MODa 小于 b 时结果会变成负数。正确处理方式ll mod_sub(ll a, ll b, ll mod) { return ((a - b) % mod mod) % mod; }先加一个 mod 再取模保证结果始终是非负的。6.3 乘法溢出导致连锁错误逆元法里有一个关键式子fact[n] * inv_fact[m] % MOD * inv_fact[n-m] % MOD。在 C 里long long 乘法可能溢出即使最终结果在 mod 范围内。比如 MOD 1e97fact[n] 和 inv_fact[m] 都在 1e9 左右两者相乘约 1e18刚好在 long long 上限9.2e18附近但已经比较接近了。如果 MOD 再大一点比如 1e13两个 1e13 的数相乘就是 1e26完全溢出。对策有两个一是在乘法过程中每步都取模比如((a * b) % MOD * c) % MOD这个我们已经在写二是引入__int128GCC 支持在乘法可能超出 long long 时使用ll mul_mod(ll a, ll b, ll mod) { return (ll)((__int128)a * b % mod); }6.4 快速幂里 b 的位运算优先级快速幂实现中if (b 1)必须有括号因为的优先级低于如果你写成if (b 1 0)实际上会被解析成if (b (1 0))也就是永远是 false代码直接死循环或者返回错误结果。这个坑我在带实习生时见过不止一次。正确写法if (b 1) { res res * a % mod; }6.5 用递归写 Lucas 时小心栈溢出Lucas 定理的递归深度是 O(log_p n)p 比较小、n 特别大的时候递归深度也可能到几百层。大部分评测机默认栈空间 8 MB几百层普通函数通常没问题但如果你在里面递归调用了逆元初始化函数每层都申请 vector栈就可能爆掉。稳妥的做法是迭代实现ll lucas_iter(ll n, ll m, ll p) { ll res 1; while (n 0 || m 0) { ll ni n % p; ll mi m % p; if (mi ni) return 0; res res * combination_mod(ni, mi, p) % p; n / p; m / p; } return res; }6.6 对拍验证低成本高保障的实测手段不管用哪种方法写完以后我都建议拿 Python 的math.comb小规模对拍一遍。Python 自带高精度整数组合数公式直接算在 n ≤ 1000 时结果一定是精确的。写一个简单的 C 程序输出结果再和 Python 结果对比能快速排除逻辑错误。对拍脚本框架import math n 100 m 50 expected math.comb(n, m) print(fPython 计算 C({n}, {m}) {expected})C 侧输出到文件两个结果比对一致基本可以确认算法正确性。7. 收尾前再分享一点个人经验刷了这么多年算法题写过无数个组合数相关的模块我觉得组合数计算是那种看起来简单、细想全是细节的问题。三种方法不是互相替代的关系而是互补的群体递推法负责小规模频繁查询逆元法负责大模数频繁查询质因数分解法负责精确大数。选错了方法轻则效率低下重则结果根本算不出来。如果你刚开始学我建议从递推法入手因为它最直观、最不容易写错再学逆元法因为它是竞赛里的主力质因数分解法可以作为进阶等真的遇到输出精确整数、范围超 long long的场景再上手也不迟。另外写这类数学函数时随手把输入合法性检查、取模修正、乘法保护性转换写进去。这些代码乍看冗余但在实际工程中能帮你省掉很多调试时间。把核心模板存好需要的时候直接调用比我每次从零推导快得多。