Lucas定理与类欧几里得算法:解决组合数前缀和取模问题
1. 从一道“超能”题目说起当组合数学遇上数论最近在整理一些经典的算法竞赛题目时又翻到了这道“P4345-[SHOI2015]超能粒子炮·改”。光看这个标题一股浓浓的科幻感扑面而来但熟悉 OI 题目的朋友都知道这背后往往藏着硬核的数学与算法。这道题的核心正如其关键词所揭示的是Lucas定理与类欧几里得算法的巧妙结合。它不是一个简单的模板套用而是要求解题者深刻理解这两个工具的本质并将它们无缝衔接来解决一个关于组合数前缀和取模的问题。很多人在第一次遇到时会觉得无从下手——明明每个知识点都学过为什么组合起来就这么难今天我们就来彻底拆解这道题不仅告诉你“怎么做”更要讲清楚“为什么这么做”以及我在反复推导和实现中踩过的那些坑。题目通常简化为如下形式给定质数 $p$$p \le 10^6$以及多组询问每组询问给出 $n, k$$n, k \le 10^{18}$要求计算 $$ S(n, k) \sum_{i0}^{k} C_{n}^{i} \mod p $$ 其中 $C_n^i$ 是组合数。$n$ 和 $k$ 的巨大范围直接宣告了任何 $O(n)$ 或 $O(k)$ 的预处理或暴力计算都是不可能的。而质数 $p$ 虽然不大但直接计算组合数再求和即使使用逆元复杂度也无法承受。这正是题目将Lucas定理和类欧推到你面前的原因——前者将大数运算分解到 $p$ 进制下后者则高效处理某种形式的求和。2. 基石深入理解Lucas定理的递归本质要攻克这道题第一步必须牢牢掌握 Lucas 定理。定理的标准形式大家都很熟悉对于质数 $p$有 $C_n^m \equiv C_{n \mod p}^{m \mod p} \cdot C_{\lfloor n/p \rfloor}^{\lfloor m/p \pmod p}$。但很多人只记住了这个公式却忽略了它本质上是一个递归过程。这种递归结构正是我们解决本题求和问题的突破口。我们可以这样形象地理解把 $n$ 和 $m$ 都写成 $p$ 进制数。计算 $C_n^m \mod p$相当于在 $p$ 进制的每一位上独立计算组合数 $C_{n_i}^{m_i}$其中 $n_i, m_i$ 是 $p$ 进制下的对应位然后再把所有位的结果乘起来。这里有一个关键限制如果某一位上 $m_i n_i$那么 $C_{n_i}^{m_i} 0$从而导致整个乘积为 $0$。这个性质在后续分析中至关重要。对于我们的求和式 $S(n, k) \sum_{i0}^{k} C_n^i$直接应用 Lucas 定理似乎不行因为定理处理的是单个组合数。但我们可以利用递归思想。设 $n n_1 \cdot p n_0$$k k_1 \cdot p k_0$其中 $0 \le n_0, k_0 p$。我们把求和下标 $i$ 也按 $p$ 进制分解$i i_1 \cdot p i_0$。那么原求和可以按照 $i_1$即 $i$ 除以 $p$ 的商进行分组$$ S(n, k) \sum_{i0}^{k} C_n^i \equiv \sum_{i_10}^{k_1-1} \left( \sum_{i_00}^{p-1} C_{n_1}^{i_1} \cdot C_{n_0}^{i_0} \right) \sum_{i_00}^{k_0} C_{n_1}^{k_1} \cdot C_{n_0}^{i_0} \pmod{p} $$这个推导需要仔细理解第一项是对所有“完整的” $p$ 块进行求和$i_1$ 从 $0$ 到 $k_1-1$。在每个块内$i_0$ 可以取遍 $[0, p-1]$。第二项是最后一个不完整的块其中 $i_1 k_1$而 $i_0$ 只能取 $[0, k_0]$。根据 Lucas 定理$C_n^i \equiv C_{n_1}^{i_1} \cdot C_{n_0}^{i_0}$。注意这里 $C_{n_1}^{i_1}$ 对于同一个 $i_1$ 是常数可以提到对 $i_0$ 的求和号外面。于是式子可以整理为 $$ S(n, k) \equiv \left( \sum_{i_10}^{k_1-1} C_{n_1}^{i_1} \right) \cdot \left( \sum_{i_00}^{p-1} C_{n_0}^{i_0} \right) C_{n_1}^{k_1} \cdot \left( \sum_{i_00}^{k_0} C_{n_0}^{i_0} \right) \pmod{p} $$观察这个式子我们发现出现了非常熟悉的结构$\sum_{i_10}^{k_1-1} C_{n_1}^{i_1}$ 就是 $S(n_1, k_1-1)$。$\sum_{i_00}^{p-1} C_{n_0}^{i_0}$ 和 $\sum_{i_00}^{k_0} C_{n_0}^{i_0}$ 中$n_0, k_0 p$这些是“小规模”问题。$C_{n_1}^{k_1}$ 也可以用 Lucas 定理递归计算但这里 $n_1, k_1$ 可能仍然很大。这揭示了一个重要的递归关系大规模问题 $S(n, k)$ 可以转化为一个小规模问题 $S(n_1, k_1-1)$、一个小规模子问题与 $n_0, k_0$ 有关、以及一个组合数 $C_{n_1}^{k_1}$ 的计算。这里 $n_1 \lfloor n/p \rfloor$$k_1 \lfloor k/p \rfloor$。注意这里有一个非常关键的边界条件陷阱。当 $k_1 0$ 时第一项求和 $\sum_{i_10}^{-1}$ 是没有意义的。在实际递归中我们需要判断如果 $k_1 0$那么第一项直接为 0整个计算退化到只处理最后一个块即 $S(n, k) \sum_{i_00}^{k_0} C_{n_0}^{i_0}$此时 $n_0, k_0 p$可以直接预处理或暴力计算。3. 递归式的化简与类欧几里得算法的浮现上一节我们得到了一个递归式但它看起来还是有些复杂尤其是其中包含 $\sum_{i_00}^{p-1} C_{n_0}^{i_0}$ 这一项。我们需要进一步化简。这里需要用到组合数的一个经典恒等式 $$ \sum_{i0}^{m} C_n^i 2^n \quad \text{当且仅当} \quad m \ge n $$ 这个等式直观理解是从 $n$ 个元素的所有子集个数就是 $2^n$。当求和上限 $m$ 大于等于 $n$ 时求和就覆盖了所有子集。在我们的式子中$\sum_{i_00}^{p-1} C_{n_0}^{i_0}$由于 $n_0 p$所以 $p-1 \ge n_0$ 恒成立因为 $n_0$ 最大为 $p-1$。因此根据上述恒等式有 $$ \sum_{i_00}^{p-1} C_{n_0}^{i_0} 2^{n_0} $$ 这是一个非常重要的简化它将一个求和式变成了一个简单的幂运算。当然这里的所有运算都是在模 $p$ 意义下的所以 $2^{n_0}$ 也需要模 $p$ 计算。于是我们的递归式可以大大简化 $$ S(n, k) \equiv S(n_1, k_1-1) \cdot 2^{n_0} C_{n_1}^{k_1} \cdot S(n_0, k_0) \pmod{p} $$ 其中$n n_1 \cdot p n_0,\quad 0 \le n_0 p$$k k_1 \cdot p k_0,\quad 0 \le k_0 p$$S(a, b) \sum_{i0}^{b} C_a^i$$n_1 \lfloor n/p \rfloor, \quad k_1 \lfloor k/p \rfloor$现在问题规模从 $(n, k)$ 递归到了 $(n_1, k_1-1)$ 和 $(n_0, k_0)$。$(n_0, k_0)$ 因为都小于 $p$可以直接求解预处理所有 $n, k p$ 的 $S(n, k)$ 或者暴力计算。而 $(n_1, k_1-1)$ 可能仍然很大需要继续递归。这个递归式和欧几里得算法辗转相除的结构非常相似我们不断用 $p$ 去除 $n$ 和 $k$得到商和余数然后对商进行递归。这正是“类欧几里得算法”名称的由来——它处理的对象不是数字的除法而是将 Lucas 定理的递归过程与求和问题相结合形成的一个递归求解过程。4. 算法核心递归实现与细节处理有了清晰的递归式我们就可以着手实现算法了。算法的核心是一个递归函数S(n, k)用于计算 $\sum_{i0}^{k} C_n^i \mod p$。4.1 递归函数的设计与边界首先我们需要明确递归的边界条件。当问题规模小到可以直接计算时递归就应当终止。当 $k 0$ 时求和的下限大于上限结果为 $0$。当 $n p$ 且 $k p$ 时此时 $n$ 和 $k$ 都很小小于质数 $p$我们可以直接计算。通常的做法是预处理出 $0 \le i, j p$ 的所有 $C_i^j \mod p$利用递推公式 $C_i^j C_{i-1}^{j-1} C_{i-1}^{j}$ 和 $C_i^0 1$然后预处理出二维前缀和数组pre[i][j] sum_{t0}^{j} C_i^t mod p。这样当 $n, k p$ 时我们可以 $O(1)$ 地返回pre[n][k]。递归过程如下def S(n, k, p): if k 0: return 0 if n p and k p: return pre[n][k] # 直接返回预处理好的小规模前缀和 # 将 n, k 分解为 p 进制下的高位和低位 n1, n0 divmod(n, p) k1, k0 divmod(k, p) # 递归计算 S(n1, k1-1) part1 S(n1, k1 - 1, p) # 计算 2^{n0} mod p pow2_n0 pow(2, n0, p) # 计算 C(n1, k1) mod p这里需要用到 Lucas 定理或递归 C_n1_k1 Lucas_C(n1, k1, p) # 注意这个函数也是递归的 # 计算 S(n0, k0)此时 n0, k0 p可直接查表 part2 pre[n0][k0] # 合并结果 ans (part1 * pow2_n0) % p ans (ans C_n1_k1 * part2) % p return ans4.2 组合数C(n1, k1)的计算陷阱在上面的递归中我们需要计算 $C_{n_1}^{k_1} \mod p$。这里 $n_1$ 和 $k_1$ 是 $\lfloor n/p \rfloor$ 和 $\lfloor k/p \rfloor$它们可能仍然非常大高达 $10^{18}/p$不能直接计算。我们必须使用Lucas 定理本身来递归计算这个组合数。因此我们需要另一个函数Lucas_C(n, m, p)。它的实现也是标准的 Lucas 定理递归def Lucas_C(n, m, p): if m n: return 0 # 边界条件当 m 0 时C(n, 0) 1 if m 0: return 1 # 当 n, m 都小于 p 时可以直接计算组合数利用预处理的阶乘和逆元 if n p and m p: return C_direct(n, m, p) # 直接计算 C(n, m) mod p # 否则递归计算 n1, n0 divmod(n, p) m1, m0 divmod(m, p) return (Lucas_C(n1, m1, p) * C_direct(n0, m0, p)) % p其中C_direct(n, m, p)用于计算当 $n, m p$ 时的组合数。这需要我们在预处理阶段准备好模 $p$ 下的阶乘数组fac[i]和阶乘逆元数组inv_fac[i]然后通过 $C_n^m \frac{fac[n]}{fac[m] \cdot fac[n-m]}$ 来计算。这里有一个极其关键的细节在S(n, k)的递归中我们调用了Lucas_C(n1, k1, p)。而Lucas_C函数内部当n1, k1 p时它会调用C_direct。这意味着我们预处理的阶乘和逆元必须至少能处理到p-1。同时我们为小规模 $S(n, k)$ 预处理的二维前缀和pre其维度也是 $p \times p$。因此算法的初始化步骤预处理复杂度是 $O(p^2)$这对于 $p \le 10^6$ 来说是可行的但需要留意内存使用pre数组需要 $p^2$ 个整型对于 $p10^6$ 会爆内存。实际上由于 $p$ 是质数且 $p \le 10^6$$O(p^2)$ 的预处理在时间和空间上都是不可接受的。这是我们遇到的第一个大坑。4.3 空间优化避免O(p^2)的预处理直接预处理 $p \times p$ 的二维前缀和pre[i][j]对于 $p10^5$ 级别就已经很吃力了内存约 40GB更别说 $10^6$。我们必须优化。观察递归式 $S(n, k) \equiv S(n_1, k_1-1) \cdot 2^{n_0} C_{n_1}^{k_1} \cdot S(n_0, k_0)$我们发现只有当递归到边界 $n p$ 且 $k p$ 时才需要小规模的 $S(n, k)$ 值。而在边界情况下我们真的需要一个 $p \times p$ 的完整二维数组吗注意到在边界时$n$ 和 $k$ 都小于 $p$但 $n$ 和 $k$ 是随机的。我们能否在需要时快速计算 $S(n, k) \sum_{i0}^{k} C_n^i$其中 $n, k p$由于 $p$ 不大$\le 10^6$而单次计算 $C_n^i$ 是 $O(1)$ 的利用预处理的阶乘和逆元那么暴力计算这个求和的时间复杂度是 $O(k)$最坏是 $O(p)$。在递归树中这样的边界情况会出现很多次如果每次都 $O(p)$ 暴力计算总复杂度可能会退化。一个经典的优化是在边界情况下我们只预处理一维信息。具体来说我们预处理出所有 $i p$ 的 $2^i \mod p$因为 $S(n, p-1) 2^n$以及阶乘和逆元用于计算单个 $C_n^i$。当需要计算 $S(n, k)$$n, k p$时我们无法 $O(1)$ 得到但可以 $O(k)$ 计算。然而在递归中$k$ 是 $k_0$它是 $k \mod p$所以 $k_0 p$。最坏情况下$O(k_0) O(p)$ 的单次计算仍然可能成为瓶颈尤其是在递归较深、调用次数多的时候。实际上通过分析递归的深度和调用次数可以发现总的边界计算次数是 $O(\log_p n)$ 级别的每次计算是 $O(p)$总复杂度约为 $O(p \log_p n)$。对于 $p10^6, n10^{18}$$\log_p n \approx 3$所以 $O(3 \times 10^6)$ 的运算量是可以接受的。因此我们可以放弃预处理二维数组pre改为在边界时实时计算 $S(n, k)$。这样空间复杂度从 $O(p^2)$ 降到了 $O(p)$只需要存储阶乘、逆元和2的幂次这是一个巨大的优化。实操心得在算法竞赛中遇到 $p$ 高达 $10^6$ 时首先要警惕 $O(p^2)$ 的空间。Lucas定理的递归结构通常与 $O(p)$ 的预处理阶乘、逆元相容。当问题要求前缀和时考虑是否能在边界进行 $O(p)$ 的暴力计算并分析递归调用次数来证明复杂度可行。这比盲目开大数组要可靠得多。5. 完整算法流程与复杂度分析让我们整合以上所有步骤梳理出完整的算法流程第一步预处理$O(p)$预处理模 $p$ 下的阶乘数组fac[0..p-1]满足fac[i] i! mod p。预处理阶乘的逆元数组inv_fac[0..p-1]通常使用费马小定理和递推公式inv_fac[i] inv_fac[i1] * (i1) % p。预处理2的幂次数组pow2[0..p-1]满足pow2[i] 2^i mod p。第二步实现辅助函数C_direct(n, m, p): 当 $0 \le m \le n p$ 时直接计算组合数 $C_n^m \mod p$。公式为fac[n] * inv_fac[m] % p * inv_fac[n-m] % p。Lucas_C(n, m, p): 递归计算 $C_n^m \mod p$。边界是m 0返回1或n, m p时调用C_direct。否则分解n1, n0 divmod(n, p); m1, m0 divmod(m, p)递归计算Lucas_C(n1, m1, p) * C_direct(n0, m0, p) % p。S_direct(n, k, p): 当 $n, k p$ 时直接计算 $\sum_{i0}^{k} C_n^i \mod p$。通过循环for i in range(0, k1): ans (ans C_direct(n, i, p)) % p实现。复杂度 $O(k)$。第三步实现主递归函数S(n, k, p)def S(n, k, p): # 边界条件1: 求和上限为负 if k 0: return 0 # 边界条件2: 小规模问题直接计算 if n p and k p: return S_direct(n, k, p) # 分解为 p 进制 n1, n0 divmod(n, p) k1, k0 divmod(k, p) # 递归计算第一部分: S(n1, k1-1) * 2^{n0} part1 S(n1, k1 - 1, p) term1 part1 * pow2[n0] % p # 计算第二部分需要的组合数 C(n1, k1) C_val Lucas_C(n1, k1, p) # 注意这里可能递归调用 Lucas_C # 计算第二部分: C(n1, k1) * S(n0, k0) part2 S_direct(n0, k0, p) # 此时 n0, k0 p term2 C_val * part2 % p # 合并结果 return (term1 term2) % p第四步处理多组询问对于每一组询问的(n, k)调用S(n, k, p)即可得到答案。复杂度分析预处理$O(p)$可接受。单次询问递归函数S和Lucas_C的递归深度均为 $O(\log_p n)$因为每次递归n和k都除以p。在每一层递归中S函数会调用一次Lucas_C而Lucas_C自身也是递归的。但注意S和Lucas_C的递归树是相似的它们处理的参数(n1, k1)是相同的。存在大量的重复子问题因此我们需要使用**记忆化搜索Memoization**来存储已经计算过的S(n, k)和Lucas_C(n, m)的结果否则复杂度会指数爆炸。记忆化状态(n, k)和(n, m)的数量是 $O(\log_p n \times \log_p k)$ 级别的实际上非常少对于 $n,k \le 10^{18}, p \approx 10^6$深度约为3-4层状态数很少。因此加上记忆化后单次询问的复杂度可以认为是 $O(\log_p n)$ 乘以一个较小的常数。边界计算在递归底层S_direct的调用次数等于递归树中叶节点的数量也是 $O(\log_p n)$ 级别每次计算是 $O(p)$。所以总边界计算复杂度为 $O(p \log_p n)$。对于 $p10^6, \log_p n \approx 3$这大约是三百万次运算在合理范围内。避坑指南记忆化搜索是关键中的关键。如果不加记忆化递归树会爆炸因为S(n1, k1-1)和Lucas_C(n1, k1)会反复计算相同的参数。可以使用两个字典或映射分别存储S和Lucas_C的结果。键可以用(n, k)元组。由于参数很大但状态数很少记忆化的效果非常显著。6. 代码实现细节与常见错误理论分析之后我们来看看具体的代码实现中需要注意的细节。我将以 Python 为例便于理解但其中关于取模、递归、记忆化的思想适用于所有语言。6.1 预处理模块的实现def preprocess(p): global fac, inv_fac, pow2 fac [1] * p for i in range(1, p): fac[i] fac[i-1] * i % p inv_fac [1] * p # 费马小定理求 (p-1)! 的逆元 inv_fac[p-1] pow(fac[p-1], p-2, p) for i in range(p-2, -1, -1): inv_fac[i] inv_fac[i1] * (i1) % p pow2 [1] * p for i in range(1, p): pow2[i] pow2[i-1] * 2 % p这里需要注意pow(fac[p-1], p-2, p)利用了费马小定理求逆元。因为p是质数所以a^(p-1) ≡ 1 (mod p)因此a^(p-2)就是a的逆元。6.2 带记忆化的Lucas_C函数from functools import lru_cache lru_cache(maxsizeNone) def Lucas_C(n, m, p): if m n: return 0 if m 0: return 1 if n p and m p: return C_direct(n, m, p) n1, n0 divmod(n, p) m1, m0 divmod(m, p) return (Lucas_C(n1, m1, p) * C_direct(n0, m0, p)) % p def C_direct(n, m, p): if m n: return 0 return fac[n] * inv_fac[m] % p * inv_fac[n-m] % p使用lru_cache装饰器可以自动实现记忆化非常方便。注意边界条件m n时组合数为 0。6.3 带记忆化的主函数Slru_cache(maxsizeNone) def S(n, k, p): if k 0: return 0 if n p and k p: # 边界情况直接计算 res 0 for i in range(k 1): res (res C_direct(n, i, p)) % p return res n1, n0 divmod(n, p) k1, k0 divmod(k, p) part1 S(n1, k1 - 1, p) # 递归计算 S(n1, k1-1) term1 part1 * pow2[n0] % p c_val Lucas_C(n1, k1, p) # 计算 C(n1, k1) # 注意S(n0, k0) 中 n0, k0 p但这里我们直接调用 S_direct 而不是递归 S避免创建不必要的记忆化状态 # 也可以专门写一个 S_direct 函数 part2 0 for i in range(k0 1): part2 (part2 C_direct(n0, i, p)) % p term2 c_val * part2 % p return (term1 term2) % p这里有一个非常重要的优化点在计算S(n0, k0)时因为n0, k0 p我们直接用了暴力循环计算而不是递归调用S(n0, k0, p)。为什么因为如果递归调用虽然n0, k0很小但会被记录到记忆化缓存中。而n0, k0的取值范围是[0, p-1]组合很多这会导致记忆化缓存被大量的小规模状态填满增加空间开销和查询时间。直接暴力计算 $O(k0)$ 的代价在边界是可接受的并且避免了污染记忆化缓存。6.4 多组询问的处理与初始化对于多组询问我们需要在每处理一组新的(p, n, k)时清空记忆化缓存并重新预处理。因为p不同预处理数组和记忆化结果都不同。def solve(): import sys input sys.stdin.read data input().split() T int(data[0]) idx 1 for _ in range(T): n int(data[idx]); idx 1 k int(data[idx]); idx 1 p int(data[idx]); idx 1 # 每次处理新的p需要重新预处理和清空缓存 preprocess(p) Lucas_C.cache_clear() # 清空Lucas_C的记忆化缓存 S.cache_clear() # 清空S的记忆化缓存 ans S(n, k, p) print(ans)常见错误总结忘记清空缓存在不同测试用例之间如果p变化了而缓存没清空会得到错误结果。记忆化键值冲突记忆化函数的参数必须包含模数p或者确保不同p的调用不会共用缓存。上述代码通过在每个用例开始前清空缓存来解决。整数溢出在计算阶乘、乘法时即使每一步都取模中间结果也可能很大在取模前。Python 大整数没问题但在 C 中要小心使用long long并及时取模。递归深度Python 的默认递归深度限制可能不够。对于 $n10^{18}, p10^6$递归深度大约为 $\log_p n \approx 3$远远小于默认的1000所以安全。但如果p很小比如2深度会达到60左右也是安全的。不过为了保险可以设置sys.setrecursionlimit(10000)。时间复杂度估计错误最坏情况下的边界计算S_direct是 $O(p)$如果p很大如 $10^6$且递归树中这样的计算次数很多可能会超时。需要相信复杂度分析$O(p \log_p n)$ 对于 $p10^6$ 和 $n10^{18}$ 是可行的约几百万次运算。7. 举一反三思维扩展与类似问题解决了“超能粒子炮”这道题我们掌握的不只是两个算法的拼凑更是一种递归分解和合并结果的思想。这种思想可以应用到许多其他问题上。扩展1模数非质数的情况如果 $p$ 不是质数Lucas 定理不再适用。此时需要用到扩展Lucas定理 (ExLucas)。ExLucas 的核心思想是将 $p$ 分解为质因数的幂 $p \prod p_i^{e_i}$分别计算答案模每个 $p_i^{e_i}$最后用中国剩余定理合并。在计算模 $p_i^{e_i}$ 时需要处理分母的逆元可能不存在的问题方法是将组合数 $C_n^m$ 中的阶乘所含的 $p_i$ 因子提取出来单独计算。这比本题要复杂得多但递归和分治的思想是相通的。扩展2求和上限k的变化本题求和是 $\sum_{i0}^{k}$。如果求和式变成 $\sum_{il}^{r} C_n^i$我们可以利用前缀和思想转化为 $S(n, r) - S(n, l-1)$。如果求和式带有其他权重比如 $\sum_{i0}^{k} i \cdot C_n^i$能否用类似的方法通常需要更复杂的推导可能会引入新的递归关系。但核心思路依然是利用数位DP按 $p$ 进制分解和组合恒等式进行化简。扩展3多维推广考虑二维求和 $\sum_{i0}^{a} \sum_{j0}^{b} C_{n}^{i} C_{m}^{j}$。这似乎可以分解为两个独立求和的乘积 $(\sum_{i0}^{a} C_n^i) \cdot (\sum_{j0}^{b} C_m^j)$。但如果求和范围有交叉限制问题会变得复杂可能就需要更高级的数学工具。个人体会这道题之所以经典是因为它将 Lucas 定理的递归性质运用到了极致。它告诉我们面对巨大的 $n$ 和 $k$不要害怕利用模数 $p$ 将其分解到 $p$ 进制下大规模问题就会自然崩塌成小规模问题。而“类欧”在这里更像是一种指导思想——模仿欧几里得算法的递归形式去设计和分析算法。在实现时记忆化搜索和精细的边界处理是保证效率的关键。最后永远要对 $O(p^2)$ 的空间复杂度保持警惕在竞赛中这通常是错误方向的信号。