组合数取模算法详解：从费马小定理到线性逆元预处理

1. 项目概述为什么我们需要一个“组合数取模的板子”在算法竞赛、密码学或者一些需要处理大规模数据的工程场景里我们经常会遇到一个看似简单却暗藏玄机的问题计算组合数 C(n, m)。比如从10个人里选3个组成小组有多少种选法答案是 C(10, 3) 120。当 n 和 m 比较小时我们直接用阶乘公式 C(n, m) n! / (m! * (n-m)!) 手算或者编程计算都没问题。但现实往往是骨感的。一旦 n 和 m 的值变大比如 n1e5, m5e4直接计算阶乘会得到一个天文数字远远超出任何编程语言中基本数据类型的表示范围。更棘手的是很多问题要求我们对结果取一个模通常是质数如 1e97以得到一个确定范围内的整数。这时直接计算再取模的路就行不通了因为中间过程的阶乘值就已经溢出了。这就是“组合数取模的板子”存在的意义。它不是一个有形的电路板而是算法领域的一个行话指的是一段封装好的、高效的、可靠的代码模块Code Template/Snippet专门用来解决“在模运算下快速计算组合数”这一高频需求。有了这个“板子”就像木匠有了趁手的工具遇到相关问题可以直接调用省去了重复推导和调试的时间把精力集中在问题本身的逻辑上。对于参赛选手这能节省宝贵的比赛时间对于开发者这能提升代码的健壮性和可维护性。2. 核心思路与方案选型从朴素到高效的演进计算模意义下的组合数核心挑战在于处理除法。因为在模运算中a / b 不能直接计算我们需要找到 b 的“模逆元”使得乘法代替除法。整个方案的演进围绕着如何高效地预处理和利用这些逆元展开。2.1 方案一费马小定理求逆元适合 n 较小查询极少这是最直观的方法。根据费马小定理若模数 MOD 是质数且 a 不是 MOD 的倍数则 a^(MOD-2) ≡ a^(-1) (mod MOD)。因此我们可以先预处理出所有阶乘 fact[i] i! % MOD然后对于每次查询 C(n, m)计算公式为C(n, m) fact[n] * inv(fact[m]) % MOD * inv(fact[n-m]) % MOD其中inv(x)表示 x 的模逆元通过快速幂计算pow(x, MOD-2, MOD)得到。为什么这么选思路直接代码简单易于理解。它清晰地展示了“阶乘”和“逆元”这两个核心概念。它的劣势是什么效率。每次查询都需要进行两次快速幂求逆元单次查询时间复杂度是 O(log MOD)。当查询次数 Q 很大时例如 Q1e5总复杂度 O(Q log MOD) 可能成为瓶颈。此外如果 n 很大比如 1e6预处理阶乘数组 O(n) 是必要的但这个方法在查询上的开销较大。2.2 方案二线性预处理阶乘逆元标准板子最常用这是竞赛和工程中最主流、最推荐的方案。核心优化点在于我们不止预处理阶乘数组fact[i]还同步预处理出阶乘的逆元数组inv_fact[i]。如何线性地预处理阶乘逆元这里有一个巧妙的递推方法首先用费马小定理求出最大阶乘fact[n_max]的逆元inv_fact[n_max]。因为inv_fact[n] (n!)^{-1} mod MOD。然后利用关系式inv_fact[i] inv_fact[i1] * (i1) % MOD反向递推。这是因为(i!)^{-1} ≡ ((i1)!)^{-1} * (i1) (mod MOD)。这样一来初始化复杂度为 O(n_max)之后每次查询 C(n, m) 只需要三次模乘C(n, m) fact[n] * inv_fact[m] % MOD * inv_fact[n-m] % MOD查询复杂度是 O(1)。为什么这是最优选它在预处理和查询之间取得了完美的平衡。对于 n 在 1e6 量级、查询次数 Q 在 1e6 量级的情景这个方案的总时间复杂度是 O(n_max Q)完全可接受。代码实现固定逻辑清晰是当之无愧的“板子”首选。2.3 方案三Lucas 定理用于模数较小或非质数模数当模数 MOD 不是很大比如小于 1e5的时候无论 n 和 m 多大都可以用 Lucas 定理将其分解为多个小规模组合数问题的乘积。 Lucas 定理C(n, m) % p C(n%p, m%p) * C(n/p, m/p) % p其中 p 是模数质数。 这是一个递归过程递归终点是 m0 或 nm。为什么需要它方案二要求预处理的阶乘数组大小至少为 n_max。如果题目中 n 可以非常大比如 1e18但模数 p 很小比如 131我们根本无法预处理 1e18 的阶乘。Lucas 定理将问题规模缩小到模数 p 以内此时我们可以用方案一或二先预处理出 0 到 p-1 的阶乘及其逆元然后递归求解。它的局限性是什么模数 p 必须是质数。对于非质数模数需要更复杂的扩展 Lucas 定理。此外由于是递归常数比 O(1) 查询的方案二要大。实操心得95% 以上的组合数取模问题模数都是质数且足够大如 1e97, 998244353n 在 1e6 以内。因此方案二线性预处理阶乘逆元是必须熟练掌握并作为默认板子的。方案一用于理解原理方案三用于应对特定“模数小、n 极大”的题目。3. 核心细节解析与实现要点我们将深入拆解最核心的方案二并给出一个工业强度的板子实现。这里假设模数 MOD 为质数且n_max为我们需要处理的最大 n 值。3.1 数据定义与初始化#include bits/stdc.h using namespace std; typedef long long ll; const int MAXN 1e6 5; // 根据题目要求调整最大范围 const ll MOD 1e9 7; // 常用质数模数 ll fact[MAXN]; // 阶乘数组 fact[i] i! % MOD ll inv_fact[MAXN]; // 阶乘的逆元数组 inv_fact[i] (i!)^{-1} % MOD // 快速幂函数用于求逆元 ll qpow(ll a, ll b) { ll res 1; while (b) { if (b 1) res res * a % MOD; a a * a % MOD; b 1; } return res; } // 初始化函数必须在程序开始调用一次 void init_comb() { fact[0] 1; for (int i 1; i MAXN; i) { fact[i] fact[i - 1] * i % MOD; } // 预处理最大阶乘的逆元 inv_fact[MAXN - 1] qpow(fact[MAXN - 1], MOD - 2); // 线性递推阶乘逆元 for (int i MAXN - 2; i 0; --i) { inv_fact[i] inv_fact[i 1] * (i 1) % MOD; } }关键细节与原理fact[0] 1这是数学定义0的阶乘为1。也是递推的起点。递推求阶乘fact[i] fact[i-1] * i % MOD。这是最基础的线性预处理。求最大阶乘逆元inv_fact[MAXN-1] qpow(fact[MAXN-1], MOD-2)。这里利用了费马小定理fact[MAXN-1]^(MOD-2)就是它的模逆元。注意MAXN-1必须是小于 MOD 的否则fact[MAXN-1]可能为0如果它是MOD的倍数导致求逆元失败。因为 MOD 通常是1e97量级的大质数而MAXN是1e6所以这个条件是满足的。线性递推逆元inv_fact[i] inv_fact[i1] * (i1) % MOD。这是整个算法的精华。推导如下我们知道fact[i1] fact[i] * (i1)。两边同时取逆元inv_fact[i1] inv_fact[i] * inv(i1)。所以inv_fact[i] inv_fact[i1] * (i1)。注意这里的(i1)是在模 MOD 意义下的而inv_fact[i1]我们已经有了从后往前算。这个递推避免了为每个 i 都调用一次快速幂将 O(n log MOD) 优化为 O(n)。3.2 组合数计算函数// 计算组合数 C(n, m) 前提是 0 m n MAXN-1 且 init_comb() 已被调用 ll comb(int n, int m) { if (m 0 || m n) return 0; // 非法输入返回0 return fact[n] * inv_fact[m] % MOD * inv_fact[n - m] % MOD; }为什么这么写边界检查if (m 0 || m n) return 0;这是必要的健壮性处理。根据组合数定义m 不在 [0, n] 范围内时结果为0。计算公式直接套用公式C(n, m) n! / (m! * (n-m)!)用乘法逆元代替除法即n! * (m!)^{-1} * ((n-m)!)^{-1} mod MOD。取模顺序连续乘法取模每次乘法后都取模防止中间结果溢出long long范围在64位系统上两个大约1e9的数相乘可能接近1e18仍在long long范围内但三个数连续乘就可能溢出。fact[n] * inv_fact[m] % MOD先算一次取模再乘下一个逆元并取模是安全的写法。3.3 板子的变体动态调整大小有时题目中 n 的最大值并不固定或者我们想节省内存。我们可以将数组改为vector并根据第一次查询的 n 动态初始化到所需大小。vectorll fact, inv_fact; ll MOD 1e97; void init_comb(int n) { // n 为需要预处理的最大阶乘下标 fact.resize(n 1, 1); inv_fact.resize(n 1, 1); for (int i 1; i n; i) { fact[i] fact[i - 1] * i % MOD; } inv_fact[n] qpow(fact[n], MOD - 2); for (int i n - 1; i 0; --i) { inv_fact[i] inv_fact[i 1] * (i 1) % MOD; } } ll comb(int n, int m) { if (m 0 || m n) return 0; // 如果查询的n超过了当前预处理范围则扩展预处理 if (n (int)fact.size()) { init_comb(n); // 重新初始化到n注意这可能不是线程安全的 } return fact[n] * inv_fact[m] % MOD * inv_fact[n - m] % MOD; }注意事项动态版本在多次查询且 n 递增的场景下可能因多次init_comb导致效率降低。通常在竞赛中我们会根据数据范围一次性分配足够大的静态数组如MAXN1e65这样更简单高效。动态版本更适合于内存敏感或 n 上限不确定的工程场景。4. 完整板子代码与使用示例下面给出一个整合了边界处理、常用功能且注释清晰的终极板子。/** * 组合数取模板子 (Modular Combination Template) * 特性 * 1. 模数 MOD 必须为质数。 * 2. 预处理 O(N)查询 O(1)。 * 3. 包含组合数 C(n, m) 和排列数 A(n, m) 计算。 */ #include bits/stdc.h using namespace std; using ll long long; const int MAXN 1e6 10; // 预处理的最高阶乘 n const ll MOD 1e9 7; // 常用质数模数 ll fact[MAXN], inv_fact[MAXN]; // 快速幂 ll qpow(ll a, ll b) { ll res 1; a % MOD; while (b) { if (b 1) res res * a % MOD; a a * a % MOD; b 1; } return res; } // 初始化阶乘和阶乘逆元 void init_comb() { fact[0] 1; for (int i 1; i MAXN; i) { fact[i] fact[i - 1] * i % MOD; } inv_fact[MAXN - 1] qpow(fact[MAXN - 1], MOD - 2); for (int i MAXN - 2; i 0; --i) { inv_fact[i] inv_fact[i 1] * (i 1) % MOD; } } // 组合数 C(n, m) 即 n 选 m ll C(int n, int m) { if (m 0 || m n) return 0; return fact[n] * inv_fact[m] % MOD * inv_fact[n - m] % MOD; } // 排列数 A(n, m) 即 n 选 m 并排列 ll A(int n, int m) { if (m 0 || m n) return 0; return fact[n] * inv_fact[n - m] % MOD; } int main() { // 初始化整个程序运行期间一次即可 init_comb(); // 示例用法 int n 10, m 3; cout C( n , m ) C(n, m) endl; // 输出 120 cout A( n , m ) A(n, m) endl; // 输出 720 // 处理多个查询 int q; cin q; while (q--) { int a, b; cin a b; cout C(a, b) \n; } return 0; }5. 常见问题、边界情况与排查技巧即使有了板子在实际使用中还是会遇到各种坑。这里记录一些典型问题和解决方法。5.1 模数 MOD 不是质数怎么办这是最常遇到的问题之一。我们的板子基于费马小定理前提是 MOD 为质数。如果 MOD 不是质数比如 MOD1000000000上述方法失效。解决方案分解模数如果 MOD 可以分解为若干个互质的质数幂的乘积如 MOD p1^k1 * p2^k2 * ...可以使用中国剩余定理CRT分别计算 C(n, m) mod p_i^k_i最后合并答案。使用扩展 Lucas 定理这是一个更通用的算法可以直接处理 MOD 为任意正整数的情况。但其原理和实现比普通 Lucas 定理复杂得多通常只在必要时才使用。审视问题在竞赛中非质数模数非常罕见。首先 double-check 题目确认模数是否真的是非质数。有时可能是自己看错了。踩坑记录我曾在一个题目中误将模数998244353质数打成了998244354偶数非质数导致调试了很久。务必仔细核对模数5.2 n 或 m 可能为 0 或负数我们的板子中C(n, m)函数已经做了边界检查if (m 0 || m n) return 0;。这符合组合数的数学定义。C(n, 0) 1C(n, n) 1。我们的计算fact[n] * inv_fact[0] * inv_fact[n]因为inv_fact[0] 1从递推式中可得所以结果正确为1。C(n, m)当 mn 时数学上定义为0我们返回0也是合理的。5.3 预处理数组大小 MAXN 设置不当这是运行时错误的主要来源。如果程序查询的n超过了MAXN-1就会访问fact数组越界。排查与解决仔细阅读题目数据范围题目通常会说明n的最大值。将MAXN设置为这个最大值加一个裕量如 5 或 10。使用动态版本如果不确定范围或者想节省内存可以采用第3.3节的动态vector版本。运行时检查可以在C(n, m)函数开头加入断言assert(n MAXN);来在调试阶段发现问题。5.4 多组测试数据下的初始化错误做法在每组测试数据中都调用init_comb()。正确做法init_comb()只需要在程序开始时调用一次。因为阶乘和逆元只与模数 MOD 和预处理上限 MAXN 有关与输入数据无关。重复初始化会造成大量重复计算可能导致超时。5.5 模数 MOD 为 1 的特殊情况虽然极少见但如果 MOD1那么所有数模1都是0。我们的快速幂求逆元qpow(a, MOD-2)会变成qpow(a, -1)导致错误。当然MOD1 时答案恒为0可以特判。通常的质数模数1e97, 998244353不会出现此问题。5.6 组合数结果的实际意义与溢出即使取了模也要时刻记得C(n, m)本身可能非常大。我们的板子计算的是C(n, m) % MOD。如果有题目需要比较组合数的大小不取模这个板子是无能为力的需要高精度计算或其他数学方法。6. 性能优化与扩展功能一个成熟的板子还可以包含一些优化和便利函数。6.1 使用 constexpr 或编译器优化对于固定的 MOD 和 MAXNC11 以后可以使用constexpr在编译期计算阶乘表但这通常只适用于较小的 MAXN几百以内对于 1e6 的量级编译期计算不现实。更实用的优化是开启编译器优化-O2。6.2 添加记忆化或查询缓存在极端情况下如果查询的 (n, m) 对大量重复可以添加一个mappairint, int, ll来缓存结果。但考虑到 O(1) 的查询已经很快而 map 操作有 log 开销通常不需要。除非查询模式极其特殊且重复率极高。6.3 扩展功能二项式系数相关计算杨辉三角预处理当 n 非常小如 2000时可以直接用动态规划计算杨辉三角C[n][m] C[n-1][m-1] C[n-1][m]并对 MOD 取模。这种方法不需要逆元甚至 MOD 不是质数也可以。const int N 2000; ll C[N5][N5]; void init_pascal() { for (int i0; iN; i) { C[i][0] C[i][i] 1; for (int j1; ji; j) { C[i][j] (C[i-1][j-1] C[i-1][j]) % MOD; } } }组合恒等式有些问题需要用到C(n, m) C(n, n-m)来减少计算或者用到C(n, k) C(n-1, k-1) * n / k的递推需要用到逆元。我们的板子已经足够高效这些恒等式更多用于数学推导。6.4 多模数支持在密码学或一些特殊数学问题中可能需要用多个不同的质数模数分别计算最后用中国剩余定理合并。可以为每个模数维护一套独立的fact和inv_fact数组。结构上可以封装一个Combination类将 MOD 和 MAXN 作为模板参数或构造函数参数。templatell MOD, int MAXN struct Combination { ll fact[MAXN], inv_fact[MAXN]; void init() { /* 初始化代码同上 */ } ll C(int n, int m) { /* 组合数代码同上 */ } }; // 使用 Combination1000000007, 1000005 comb1; Combination998244353, 1000005 comb2;这个“组合数取模的板子”从原理到实现从基础用法到边界处理几乎涵盖了你会遇到的所有情况。核心就是“线性预处理阶乘和阶乘逆元”这个经典算法。把它理解透彻、代码写熟遇到相关问题就能像套公式一样快速解决。在实际比赛或项目中我通常会把这个板子代码单独存成一个头文件如comb.hpp需要时直接#include并调用init_comb()和C(n, m)这能极大提升编码效率和准确性。最后再强调一次安全使用的前提确保模数 MOD 是一个质数并且预处理范围 MAXN 大于你将用到的最大 n。