1. 从一道“超纲”的算法题说起最近在刷AcWing的算法提高课做到第887题“求组合数 III”时我愣了一下。题目要求计算组合数 C(a, b) mod p但这里的 p 不是一个普通的数而是一个质数并且 a 和 b 的范围可以大到 10^18p 的范围在 1 到 10^5 之间。如果你尝试用之前学过的递推法O(n^2)或者预处理阶乘逆元法O(n log n)会立刻发现此路不通——因为 a 和 b 的规模远超你能预处理的范围内存和时间都不允许。这就像给你一把螺丝刀却让你去拧一颗需要液压扳手才能动的大螺栓工具完全不对路。这道题的核心是Lucas定理。我第一次接触这个定理是在数论专题里它完美地解决了“大组合数取模小质数”的经典难题。但如果你只停留在背诵定理公式C(a, b) ≡ C(a mod p, b mod p) * C(a/p, b/p) (mod p)的层面那么在实际编码和应对竞赛或面试时你可能会吃大亏。因为定理只是给出了一个递归分解的思路真正高效的实现尤其是处理边界条件和进行必要的优化才是区分“会”与“精通”的关键。网上很多教程只给个模板却很少讲清楚为什么递归终点那样写、为什么需要那些特判、以及如何优化到极致。今天我就结合这道AcWing 887把Lucas定理从原理到实现再到一个经过实战检验的“优化版”实现彻底拆解清楚。无论你是正在备赛的选手还是对算法感兴趣的程序员这篇内容都能让你不仅“抄”到一个好模板更能理解其每一行代码背后的逻辑。2. Lucas定理化“大”为“小”的魔法要理解优化必须先吃透原理。Lucas定理之所以能处理 a, b 很大的情况其精髓在于“降维打击”。2.1 定理的直观理解我们先把组合数 C(a, b) 写成更标准的数学形式C(n, m) n! / (m! * (n-m)!)。当 n 和 m 很大但模数 p 是一个相对较小的质数时直接计算阶乘再取模是不可行的因为 n! 本身可能都算不出来。Lucas定理提供了一种递归计算方法C(n, m) ≡ C(n mod p, m mod p) * C(n/p, m/p) (mod p)这个公式在做什么它实际上是在利用 p 进制表示。我们把 n 和 m 写成 p 进制数 n n_k * p^k n_{k-1} * p^{k-1} ... n_1 * p n_0 m m_k * p^k m_{k-1} * p^{k-1} ... m_1 * p m_0那么Lucas定理告诉我们模 p 意义下的组合数 C(n, m)等于它们在 p 进制下每一位对应的组合数 C(n_i, m_i) 的乘积再对 p 取模。其中C(n/p, m/p) 就是处理更高位的递归过程。为什么是质数 p这是定理成立的核心前提。因为组合数的计算涉及除法阶乘的除法而在模运算中除法需要转化为乘以逆元。逆元存在的充要条件就是除数和模数互质。当模数 p 是质数且我们处理的数阶乘中的每一项都小于 p 时这些数都与 p 互质从而它们的逆元一定存在整个计算过程在模 p 域内才是封闭且良定义的。2.2 一个手工演算的例子假设我们要计算 C(13, 4) mod 3。p 3 是质数。应用 Lucas 定理n 13, m 4。n mod p 13 % 3 1m mod p 4 % 3 1n/p 13 / 3 4 (整数除法)m/p 4 / 3 1所以 C(13, 4) ≡ C(1, 1) * C(4, 1) (mod 3)继续递归计算 C(4, 1) mod 3n 4, m 1。n mod p 4 % 3 1m mod p 1 % 3 1n/p 4 / 3 1m/p 1 / 3 0所以 C(4, 1) ≡ C(1, 1) * C(1, 0) (mod 3)我们知道基础组合数C(1,1)1, C(1,0)1。因此C(4,1) ≡ 1 * 1 ≡ 1 (mod 3)。回溯C(13, 4) ≡ C(1,1) * 1 ≡ 1 * 1 ≡ 1 (mod 3)。你可以验证一下C(13,4)715715 % 3 1结果正确。这个例子展示了如何将计算 C(13,4) 转化为计算更小的 C(1,1), C(4,1), C(1,0)。当问题规模n, m远大于 p 时这种转化的效率优势是指数级的。3. 基础实现与隐藏的“坑点”理解了定理我们可以写出第一版代码。这里以C为例但逻辑是通用的。// 快速幂计算 a^k % p long long qmi(long long a, long long k, long long p) { long long res 1; while (k) { if (k 1) res res * a % p; a a * a % p; k 1; } return res; } // 通过预处理阶乘和阶乘逆元计算 C(a, b) % p其中 a, b p long long C(long long a, long long b, long long p) { if (b a) return 0; // 组合数定义b不能大于a long long res 1; // 根据公式 C(a, b) a! / (b! * (a-b)!) // 在模p下即 a! * inv(b!) * inv((a-b)!) for (long long i 1, j a; i b; i, j--) { res res * j % p; res res * qmi(i, p-2, p) % p; // 费马小定理求逆元 } return res; } // Lucas定理递归实现 long long lucas(long long a, long long b, long long p) { if (a p b p) return C(a, b, p); // 递归终点 return C(a % p, b % p, p) * lucas(a / p, b / p, p) % p; }这个版本很直观但它在效率和健壮性上都有问题是典型的“教科书式”实现离竞赛或工程要求有差距。坑点1C函数中的低效计算C函数用循环连乘再逐个求逆元时间复杂度是 O(b)。当 b 接近 p最大可达1e5时单次计算就是1e5量级。而lucas递归深度是 O(log_p a)最坏情况下可能调用多次C函数导致总时间开销变大。在AcWing 887这种有多组测试数据的题目中很容易超时。坑点2未考虑b a的递归传递注意基础C函数中我们判断了if (b a) return 0。但是在lucas函数中递归调用是lucas(a/p, b/p, p)。如果a p但b p那么在递归终点a p b p这个条件可能永远不成立因为b可能一直大于等于p不对仔细想b在每次递归时也会除以p所以最终b也会变小。但是存在一种情况在某一层递归a % p b % p。根据组合数定义C(a%p, b%p)应该为0。我们的C函数能正确处理。然而更关键的是递归终点条件。如果b a根据定义C(a, b)0。这个情况可能在递归的某一层出现。我们的终点条件if (a p b p)在a p但b可能仍然很大因为b来自上层可能bp时就会直接调用C(a, b, p)而此时的b可能大于a因为a已经小于p而b可能大于等于p。所以递归终点条件应该更严格。坑点3没有利用 p 是小质数的特点进行预处理既然 p 最大只有1e5我们完全可以预处理出所有阶乘fact[i]和阶乘的逆元infact[i]模 p 的结果。这样计算C(a, b, p)时其中 a, b p就可以用fact[a] * infact[b] % p * infact[a-b] % p在 O(1) 时间内得到结果。这比循环乘法和快速幂求逆元快得多。4. 优化版实现每一步的抉择与理由针对上述坑点我们重构代码。目标是单次查询C(a, b) % p的时间复杂度为 O(log_p a)并且健壮。4.1 优化1预处理阶乘与阶乘逆元这是性能提升的关键。我们开两个数组fact和infact长度至少为p因为 a, b p。fact[i] i! % pinfact[i] (i!)^{-1} % p即 i! 的模 p 逆元。如何求infact[i]有两种主流方法对每个fact[i]用快速幂求逆元复杂度 O(p log p)。更高效的方法先求出infact[p-1] fact[p-1]^{p-2}费马小定理然后利用关系infact[i-1] infact[i] * i % p倒着递推求出所有infact[i]。因为(i!)^{-1} ≡ ((i1)!)^{-1} * (i1) (mod p)。这样总复杂度是 O(p)。我们采用第二种方法因为它更快。// 预处理阶乘和阶乘逆元p是模数 vectorlong long fact, infact; void init(int p) { fact.resize(p 1); infact.resize(p 1); fact[0] infact[0] 1; for (int i 1; i p; i) { fact[i] fact[i - 1] * i % p; } // 求 fact[p-1] 的逆元 infact[p - 1] qmi(fact[p - 1], p - 2, p); // 倒推阶乘逆元 for (int i p - 2; i 0; i--) { infact[i] infact[i 1] * (i 1) % p; } // infact[0] 已经是1 }4.2 优化2改写 C 函数与 Lucas 函数现在C(a, b, p)函数中的 a, b 已经满足 a, b p这是由 Lucas 递归保证的我们稍后修改递归终点。我们可以安全地使用预处理数组。// 计算 C(a, b) % p, 要求 a, b p long long C(long long a, long long b, long long p) { if (b a) return 0; // 组合数定义 // 公式C(a, b) a! / (b! * (a-b)!) // 模 p 下 fact[a] * infact[b] * infact[a-b] % p return fact[a] * infact[b] % p * infact[a - b] % p; }接下来是关键的lucas函数。我们需要修正递归终点条件。思考一下递归应该什么时候结束当b 0时根据定义C(a, 0) 1这是一个自然的终点。另外当a p且b p时我们可以直接查表计算。但是更严谨且包含所有情况的条件是当a p且b p时因为此时 a, b 已经落在我们预处理数组的范围内可以直接用 O(1) 的C函数计算。如果b 0但a还很大a p我们还需要继续递归吗不需要因为C(a, 0) 1是恒成立的与 a 的大小无关。所以我们可以把b 0作为一个优先判断的终点。long long lucas(long long a, long long b, long long p) { if (b 0) return 1; // C(a, 0) 1 // 当 a, b 都小于 p 时直接使用预处理结果计算 if (a p b p) return C(a, b, p); // 否则应用 Lucas 定理递归 return C(a % p, b % p, p) * lucas(a / p, b / p, p) % p; }这里有一个非常重要的细节递归调用时我们使用a / p和b / p。在C中对于long long类型的/操作是向零取整的整数除法这正是我们需要的它对应了 p 进制表示中移除最低位的过程。4.3 整合与边界处理将以上整合并考虑多组测试数据。注意预处理init(p)只需要在每组测试数据内部且 p 不变的情况下执行一次。如果题目有多组数据且 p 不同那么每组数据都需要重新初始化。#include iostream #include vector using namespace std; long long qmi(long long a, long long k, long long p) { long long res 1; while (k) { if (k 1) res res * a % p; a a * a % p; k 1; } return res; } vectorlong long fact, infact; void init(int p) { fact.resize(p 1); infact.resize(p 1); fact[0] infact[0] 1; for (int i 1; i p; i) { fact[i] fact[i - 1] * i % p; } infact[p - 1] qmi(fact[p - 1], p - 2, p); for (int i p - 2; i 0; i--) { infact[i] infact[i 1] * (i 1) % p; } } long long C(long long a, long long b, long long p) { if (b a) return 0; // 注意这里 a 和 b 是 long long 类型但实际值一定小于 p所以可以安全转型为 int 用作数组下标 return fact[a] * infact[b] % p * infact[a - b] % p; } long long lucas(long long a, long long b, long long p) { if (b 0) return 1; if (a p b p) return C(a, b, p); return C(a % p, b % p, p) * lucas(a / p, b / p, p) % p; } int main() { int n; cin n; while (n--) { long long a, b, p; cin a b p; init(p); // 对每个不同的 p 进行预处理 cout lucas(a, b, p) endl; } return 0; }5. 针对AcWing 887的最终适配与细节打磨上面的代码已经是一个健壮的优化版 Lucas 实现了。但对于 AcWing 887 这道具体的题目我们还可以根据它的输入输出格式和潜在的数据特点做最后的调整和细节强化。5.1 输入输出的效率题目没有明确说数据量有多大但算法竞赛中多组测试数据是常态。使用cin/cout可能会比scanf/printf慢一些尤其是在关闭同步流的情况下。为了保险起见我们可以使用scanf和printf或者给cin/cout加上加速语句ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);5.2 预处理数组的存储与复用在我们的主循环中对于每一组(a, b, p)我们都调用了init(p)。如果多组数据的模数p相同我们就会重复初始化造成浪费。但题目没有保证p相同所以我们无法复用。这里就保持每轮重新初始化。数组fact和infact我们定义在全局但在init函数中会被resize。注意vector的resize操作本身有一定开销但鉴于 p 1e5这个开销可以接受。5.3 关于“优化版”的再思考真的最优吗我们称其为“优化版”主要是针对那个 O(b) 计算小组合数的朴素方法。我们的优化在于用 O(p) 的预处理换来 O(1) 的查询。这对于需要多次计算不同 (a, b) 组合数但模数 p 固定的场景是巨大的提升。然而在 AcWing 887 的语境下每组数据独立且我们每组数据只查询一次lucas(a, b, p)。那么预处理的 O(p) 成本是否划算我们来粗略估算一下朴素版循环求 C单次lucas递归深度约log_p a假设为 L。每次递归中的C(a%p, b%p, p)计算成本是 O(p)最坏。总成本 O(L * p)。由于 a 可达 1e18p 最小为 2L 最大约 60。p 最大 1e5所以最坏复杂度约为 60 * 1e5 6e6 次运算。对于单组数据尚可但多组数据比如 100 组就可能到 6e8有风险。预处理版每组数据预处理成本 O(p)查询成本 O(L)。总成本 O(p L)。最坏情况 p1e5L60成本约 1e5。对于多组数据是 O(n * p)。如果 n100就是 1e7 次运算比朴素版最坏情况好。实际上在 p 较大接近1e5时预处理版优势明显。在 p 很小时两者差异不大但预处理版代码更清晰且C函数是 O(1) 的常数更小。所以预处理版是更稳定、更通用的选择。5.4 一个易错点long long 的乘法溢出注意lucas函数的返回语句return C(a % p, b % p, p) * lucas(a / p, b / p, p) % p;这里两个部分相乘即使各自都取模过 p它们的乘积仍然可能超过long long的范围大约 9e18导致乘法溢出。例如p1e5两个模 p 后的数最大都是 1e5-1乘积约 1e10在 long long 范围内。但为了绝对安全可以使用慢速乘龟速乘或者在乘法时转为__int128如果编译器支持。在竞赛环境中通常p 1e5乘积小于 1e10用long long是安全的。但这是一个好的编程习惯点。// 使用 __int128 防止中间结果溢出如果环境支持 return (long long)((__int128)C(a % p, b % p, p) * lucas(a / p, b / p, p) % p);5.5 最终代码呈现结合所有考虑AcWing 887 的最终优化版代码如下#include iostream #include vector using namespace std; typedef long long LL; LL qmi(LL a, LL k, LL p) { LL res 1; while (k) { if (k 1) res res * a % p; a a * a % p; k 1; } return res; } vectorLL fact, infact; void init(LL p) { fact.resize(p 1); infact.resize(p 1); fact[0] infact[0] 1; for (int i 1; i p; i) { fact[i] fact[i - 1] * i % p; } infact[p - 1] qmi(fact[p - 1], p - 2, p); for (int i p - 2; i 0; i--) { infact[i] infact[i 1] * (i 1) % p; } } LL C(LL a, LL b, LL p) { if (b a) return 0; // a, b 此时一定小于 p直接作为下标访问数组 return fact[a] * infact[b] % p * infact[a - b] % p; } LL lucas(LL a, LL b, LL p) { if (b 0) return 1; if (a p b p) return C(a, b, p); // 使用显式类型转换确保中间乘法在LL范围内或使用__int128 return C(a % p, b % p, p) * lucas(a / p, b / p, p) % p; } int main() { ios::sync_with_stdio(false); cin.tie(0); int n; cin n; while (n--) { LL a, b, p; cin a b p; init(p); cout lucas(a, b, p) \n; } return 0; }6. 测试与验证如何确保代码正确写完代码尤其是涉及数论和递归的代码必须进行测试。除了题目给的样例我们应该自己构造一些边界和特殊情况的测试用例。测试用例设计基础验证C(5, 2) mod 7。手工计算 10 mod 7 3。C(13, 4) mod 3。我们前面算过结果是1。边界情况b 0:C(100, 0) mod 11应该为 1。b a:C(5, 10) mod 13应该为 0。a 或 b 等于 0:C(0, 0) mod 5定义为 1。p 很小C(10, 4) mod 2。10 mod 20, 4 mod 20, C(0,0)1; 10/25, 4/22, 递归计算 C(5,2) mod 2。5 mod 21, 2 mod 20, C(1,0)1; 5/22, 2/21, 递归 C(2,1) mod 2。2 mod 20, 1 mod 21, C(0,1)0。所以最终结果 1100。验证C(10,4)210, 210 mod 20。正确。p 较大a,b 也很大C(1e18, 1e17) mod 99991。需要程序跑可以和其他可靠方法如Python的大数计算后取模对比。递归深度测试a非常大p很小如2测试递归是否正常终止不会栈溢出。对于a1e18, p2递归深度约 log2(1e18) ≈ 60完全安全。调试技巧可以在lucas和C函数中加入一些打印语句提交时记得删除观察递归过程和中间结果确保与手工推导一致。LL lucas(LL a, LL b, LL p) { // cout lucas: a a , b b , p p endl; if (b 0) return 1; if (a p b p) { LL res C(a, b, p); // cout Direct C( a , b ) res endl; return res; } LL part1 C(a % p, b % p, p); LL part2 lucas(a / p, b / p, p); LL res part1 * part2 % p; // cout Recursive: C( a%p , b%p ) part1 , lucas( a/p , b/p ) part2 , res res endl; return res; }7. 举一反三Lucas定理的变体与相关题目掌握了标准的 Lucas 定理你可以解决一大类“大组合数取模小质数”的问题。但还有一些相关的变体和技巧值得了解。7.1 当 p 不是质数时怎么办Lucas 定理要求 p 是质数。如果 p 是合数我们需要用扩展 Lucas 定理。其核心思想是将合数 p 分解质因数p p1^k1 * p2^k2 * ... * pm^km。然后分别计算C(n, m) mod pi^ki最后用中国剩余定理合并结果。计算C(n, m) mod pi^ki时需要处理阶乘中 pi 的因子因为分母可能有 pi 的因子导致无法求逆元这通常通过剔除因子再计算的方法解决。这比标准 Lucas 复杂得多属于数论进阶内容。7.2 如何计算组合数模任意数不一定是质数除了扩展 Lucas对于模数不是质数但不大比如 1e6的情况我们可以尝试对组合数公式n! / (m! * (n-m)!)中的阶乘进行质因数分解将除法转化为质因子指数的加减最后用快速幂乘起来。这种方法适用于模数可以不是质数但 n, m 不能太大因为要分解 1~n 的阶乘。7.3 相关题目推荐AcWing 886. 求组合数 II模数1e97是质数但 n, m 范围较大1e5需要用预处理阶乘逆元的方法可以看作是 Lucas 定理中p很大大于 n, m时的特例此时直接套用预处理公式即可不需要递归。洛谷 P3807 【模板】卢卡斯定理标准 Lucas 定理模板题。洛谷 P4720 【模板】扩展卢卡斯定理学习扩展 Lucas 的好题目。一些涉及组合数取模的计数问题在许多动态规划或排列组合计数问题中如果结果需要对一个质数取模并且需要计算的组合数下标很大那么 Lucas 定理可能就是破题的关键。回过头看 AcWing 887它确实是一道很好的模板题迫使你理解并实现 Lucas 定理。而“优化版”的实现不仅仅是把 O(b) 的 C 函数优化成 O(1)更重要的是对递归终点、边界条件、溢出处理的周全考虑。把这些细节都琢磨透下次再遇到“大数取模”的问题你心里就有底了。编程实现数论算法往往公式推导只占一半功夫把公式严谨、高效、无漏洞地转化成代码才是真正的挑战。
![[光学原理与应用-521]：对光的错误理解与纠偏](http://pic.xiahunao.cn/yaotu/[光学原理与应用-521]：对光的错误理解与纠偏)
