- bitworld 的博客
【数论】卢卡斯定理
- 2025-4-26 10:09:29 @
卢卡斯定理
想象一下,当组合数的 和 变得非常巨大,而你只需要它对一个(通常不大的)素数 取模的结果时,直接计算阶乘和它们的逆元可能会因为 或 大于等于 而导致分母出现 的倍数,使得模 意义下的逆元不存在,或者计算量过大。这时,卢卡斯定理就像一位先知,为我们指明了一条优雅而高效的道路。
第一章:组合数
1.1 组合
组合数学是数学的一个迷人分支,它研究有限集合的排列、组合以及计数问题。"组合"这个词本身就充满了魅力,它关乎“选择”。从 个不同的元素中,不考虑顺序地选出 个元素,有多少种不同的选法?这就是我们常说的组合数,记作 或 。
它的计算公式大家应该很熟悉了:
$$\binom{n}{k} = \frac{n!}{k!(n-k)!} \quad (0 \le k \le n) $$其中 是阶乘, 特殊定义为 。如果 或 ,则 。
组合数无处不在。概率论中的抽样问题、二项式展开 的系数、图论中的路径计数、动态规划的状态转移…… 它们是构建许多数学模型和算法的基础砖块。理解和计算组合数,是每一位算法学习者的基本功。
1.2 模运算
在计算机的世界里,我们经常和“溢出”作斗争。当计算结果大到无法用标准数据类型(如 int
或 long long
)存储时,问题就来了。组合数,特别是涉及阶乘的计算,数值增长极快。例如, 就已经远远超出了 64 位整数的表示范围。
为了应对这个问题,算法竞赛题目常常要求我们将最终答案对一个特定的数 取模(Modulo)。模运算,你可以理解为一种“循环”的算术。想象一个只有 个刻度(0 到 )的钟表,任何整数最终都会落在这个钟表的某个刻度上。我们说 和 模 同余,记作 ,意思是 和 除以 的余数相同。
模运算有很好的性质,加、减、乘法都可以直接在模意义下进行:
- $(a + b) \pmod{M} = ((a \pmod{M}) + (b \pmod{M})) \pmod{M}$
- $(a - b) \pmod{M} = ((a \pmod{M}) - (b \pmod{M}) + M) \pmod{M}$ (加 是为了保证结果非负)
- $(a \times b) \pmod{M} = ((a \pmod{M}) \times (b \pmod{M})) \pmod{M}$
这些性质使得我们可以在计算过程中不断取模,将中间结果和最终结果都控制在 范围内,有效避免溢出。
1.3 挑战:巨大的组合数与有限的模
现在,我们面临的问题是计算 。看起来似乎很简单,结合组合数公式和模运算性质,我们是不是可以直接计算 , , ,然后进行“除法”呢?
这里潜藏着一个关键问题:模运算下的除法并不像我们熟悉的实数除法那样直接。计算 不能简单地变成 。我们需要引入“模逆元”的概念。
1.4 关键工具:模逆元与费马小定理
在模 的世界里,我们寻找 的“除数”的替代品,即一个数 ,使得 。这个 就被称为 关于模 的模逆元。如果找到了 ,那么模意义下的除法就可以转化为乘法:
$$(A / B) \pmod{M} \equiv (A \times B^{-1}) \pmod{M} $$模逆元并非总是存在。它存在的充要条件是 与 互质,即它们的最大公约数 。
幸运的是,当模数 是一个素数(Prime Number),我们设为 时,情况变得简单。对于任何不能被 整除的整数 (即 ), 和 必然互质,因此模 的逆元 一定存在。
如何找到这个逆元呢?费马小定理 (Fermat's Little Theorem) 闪亮登场!它告诉我们:如果 是一个素数,且整数 不能被 整除(即 ),那么:
将这个式子两边同乘以 (我们知道它存在),得到:
太棒了!这意味着,要计算 模素数 的逆元,我们只需要计算 即可。而计算一个数的幂模一个数,可以使用快速幂(也叫二进制幂)算法在 的时间复杂度内高效完成。
至此,计算 (其中 是素数)的“朴素”方法似乎已经成型:
- 计算 。
- 计算 。
- 计算 。
- 如果 且 ,则使用快速幂计算 和 。
- 最终结果是 。
但这种方法真的万无一失吗?下一章,我们将深入探讨它的局限性,从而引出卢卡斯定理的必要性。
第二章:为何朴素方法失效?
我们在第一章结尾提出了一个计算 ( 为素数)的看似可行的方法:计算分子 ,计算分母 ,然后利用费马小定理求分母的模逆元,最后将分子与分母的逆元相乘。这个流程在很多情况下是有效的,尤其是当 和 相对于 不太大时。
然而,当 或 变得非常大时,这个“朴素”方法会遇到一个致命的障碍。
2.1 朴素计算:阶乘与逆元
让我们回顾一下这个方法。计算 。 我们需要计算三个阶乘模 :, , 。 然后需要计算 和 的模 逆元。 最后是乘法 $(n! \pmod p) \times (k!)^{-1} \pmod p \times ((n-k)! )^{-1} \pmod p$。
计算阶乘模 看起来没问题,我们可以迭代计算:fact[i] = (fact[i-1] * i) % p
。即使 很大,只要 本身还在可接受范围内(比如 ),计算 似乎不是主要瓶颈,因为我们只需要关心 模 的行为。等等,这里有个误区!计算 并不是只关心 ,而是需要 。如果 很大,直接循环计算 次是不可行的。
但更大的问题出在求逆元上。
2.2 陷阱:当模数 潜入分母
费马小定理告诉我们,计算 的前提是 。 在我们的组合数公式 中,分母是 。 我们需要计算 和 。
那么,什么时候 或者 呢? 答案很简单:当阶乘的计算包含了因子 时。 也就是说,如果 ,那么 $k! = 1 \times 2 \times \dots \times p \times \dots \times k$ 中必然包含因子 ,所以 。 同理,如果 ,那么 。
一旦分母中的某一项模 等于 0,它的模 逆元就不存在了!我们无法使用费马小定理,整个基于“乘以逆元”的计算流程就崩溃了。
核心障碍:当 时,朴素方法计算 可能会因为分母 或 包含 因子而失效。
例如,计算 。 。这里 且 且 。 $7! \pmod 5 = (1 \cdot 2 \cdot 3 \cdot 4 \cdot \underline{0} \cdot 1 \cdot 2) \pmod 5 = 0$。 。 。 。 $4!^{-1} \equiv 4^{-1} \equiv 4^{5-2} \equiv 4^3 = 64 \equiv 4 \pmod 5$。 结果是 。 实际 $C(7, 3) = \frac{7 \times 6 \times 5}{3 \times 2 \times 1} = 35 \equiv 0 \pmod 5$。这次碰巧可以。
再看 。 。, , 。 。 。 。 分母 模 5 为 0,逆元不存在!朴素方法失败。
2.3 瓶颈:巨大的 与 的比较
在算法竞赛中, 的范围可以达到 甚至 ,而模数 可能相对较小,比如 或 这种常见的素数。只要 ,朴素方法就面临理论上的困难(逆元可能不存在)。
即使逆元存在(比如 且 ,但这要求 ,限制性太强),计算 本身也可能是个问题。虽然我们可以发现 $n! \pmod p = (1 \times \dots \times (p-1)) \times (p+1 \times \dots \times (2p-1)) \times \dots \pmod p$,这其中涉及到 Wilson 定理 以及一些周期性,但处理起来依然复杂,特别是当 巨大时。
我们需要一种不依赖于直接计算巨大阶乘,并且能巧妙绕开“分母为零”陷阱的方法。
2.4 我们需要新方法
总结一下,当 相对于素数模 很大时(特别是 ),计算 的朴素方法存在两个主要问题:
- 可行性问题: 当 或 时,分母 或 模 为 0,导致无法求解模逆元。
- 效率问题: 即使分母模 不为 0,当 非常大时,直接计算 也是低效或不可行的。
这强烈的需求推动我们去寻找更高级的工具。幸运的是,法国数学家爱德华·卢卡斯 (Édouard Lucas) 在 19 世纪就为我们提供了这样一个强大的定理。它通过一种精妙的方式,将大组合数模素数的问题分解为一系列小组合数模素数的问题,完美地解决了上述困难。
下一章,我们将正式揭开卢卡斯定理的神秘面纱。
第三章:卢卡斯定理
面对计算大组合数模素数 的困境,卢卡斯定理应运而生。它提供了一个极其优雅且强大的工具,让我们能够将看似棘手的问题分解,从而轻松求解。
3.1 定理叙述:来自卢卡斯的智慧
卢卡斯定理 (Lucas's Theorem): 设 是一个素数,对于非负整数 和 ,我们将 和 分别表示为 进制数:
$$n = n_m p^m + n_{m-1} p^{m-1} + \dots + n_1 p + n_0 $$$$k = k_m p^m + k_{m-1} p^{m-1} + \dots + k_1 p + k_0 $$其中 是 和 在 进制下的第 位数字(不足的位数用 0 补齐,使得两者具有相同的最高次幂 )。 则有如下同余关系:
$$\binom{n}{k} \equiv \prod_{i=0}^{m} \binom{n_i}{k_i} \pmod{p} $$换句话说,计算 的结果,等于将 和 写成 进制后,对每一位 分别计算其对应的组合数 ,然后将这些结果全部乘起来,最后再对 取模。
3.2 进制视角:拆解与对应
这个定理最核心的洞见在于它建立了整体组合数 与其 进制表示下各位数字对应的“局部”组合数 之间的联系。
如何得到 和 的 进制表示呢?这和我们熟悉的十进制转二进制类似,就是不断地用 除以 取余数。
- ... 以此类推。
对 也进行同样的操作得到 。
例如,计算 。 。 先将 和 转换为 13 进制: ,所以 。记作 。 ,所以 。记作 。 (这里最高次幂 )
根据卢卡斯定理:
$$\binom{100}{30} \equiv \binom{n_1}{k_1} \times \binom{n_0}{k_0} \pmod{13} $$$$\binom{100}{30} \equiv \binom{7}{2} \times \binom{9}{4} \pmod{13} $$现在,问题被转化为了计算两个规模小得多的组合数模 13: 和 。
3.3 核心思想:大问题化小,各个击破
卢卡斯定理的威力在于“分治”思想。它把一个可能涉及非常大的 和 的组合数计算 ,分解成了计算一系列 的乘积。
注意到 和 都是 和 的 进制位,因此它们满足 。这意味着:
- 计算规模变小: 我们只需要计算形如 的组合数,其中 且 。
- 朴素方法可行: 对于 ,由于 ,分母 和 都不可能包含因子 (除非它们是 ),因此它们模 都不为 0。这意味着我们可以安全地使用第一章提到的朴素方法(计算阶乘模 ,然后用费马小定理求逆元)来计算每一个 。
卢卡斯定理巧妙地绕过了第二章提到的所有障碍!它不需要直接处理巨大的 ,也不怕分母出现 的倍数。它将问题分解到可以在模 下安全计算的小单元上。
3.4 特殊情况: 或
我们知道组合数的定义要求 。如果给定的 ,那么 。卢卡斯定理是否也适用于这种情况呢?
答案是肯定的。如果 ,那么在它们的 进制表示中,必然存在某个最高位 使得 (否则 不会大于 )。根据组合数的定义,如果 ,则 。由于卢卡斯定理的结果是所有 的乘积,只要其中有一项是 0,整个乘积就等于 0。 所以,如果 ,卢卡斯定理 会自然地得到 0,与 一致。
更一般地,即使 ,如果在 进制表示的某一位 上出现了 ,那么 ,导致整个乘积 也等于 0。
推论: 的一个必要条件是,在 进制下, 的每一位数字都不超过 对应的位数字(即对所有 都有 )。 这有时也被描述为:当且仅当在 进制表示下, 的每一位都“小于等于” 的对应位时, 才可能不被 整除。这在二进制 () 时尤其直观: 当且仅当 的二进制表示是 的二进制表示的子集(即 )。
现在我们理解了卢卡斯定理的内容和核心思想。下一章,我们将通过具体的例子来演练如何使用它。
第四章:实例剖析
4.1 小试牛刀:计算
- 问题: 计算 。
- 参数: , , 。
- 步骤 1: 进制转换
- 转为 5 进制:。所以 。即 。
- 转为 5 进制:。所以 。即 。
- 步骤 2:应用卢卡斯定理$$\binom{10}{3} \equiv \binom{n_1}{k_1} \times \binom{n_0}{k_0} \pmod 5 $$$$\binom{10}{3} \equiv \binom{2}{0} \times \binom{0}{3} \pmod 5 $$
- 步骤 3:计算小组合数模
- 。
- 。因为 ,根据组合数定义,。
- 步骤 4:合并结果
- 验证: $\binom{10}{3} = \frac{10 \times 9 \times 8}{3 \times 2 \times 1} = 120$。 余 0。所以 。结果正确!注意这里我们绕开了第二章中 导致 的问题。
4.2 步步为营:计算
- 问题: 计算 。
- 参数: , , 。
- 步骤 1: 进制转换
- 转为 7 进制:。所以 。即 。
- 转为 7 进制:。所以 。即 。
- 步骤 2:应用卢卡斯定理$$\binom{12}{5} \equiv \binom{n_1}{k_1} \times \binom{n_0}{k_0} \pmod 7 $$$$\binom{12}{5} \equiv \binom{1}{0} \times \binom{5}{5} \pmod 7 $$
- 步骤 3:计算小组合数模
- 。
- 。
- 步骤 4:合并结果
- 验证: $\binom{12}{5} = \frac{12 \times 11 \times 10 \times 9 \times 8}{5 \times 4 \times 3 \times 2 \times 1} = 11 \times 9 \times 8 = 792$。 余 1。所以 。结果正确!
4.3 稍显复杂:计算
- 问题: 计算 。
- 参数: , , 。
- 步骤 1: 进制转换 (在 3.2 节已完成)
- ()
- ()
- 步骤 2:应用卢卡斯定理$$\binom{100}{30} \equiv \binom{7}{2} \times \binom{9}{4} \pmod{13} $$
- 步骤 3:计算小组合数模
- 计算 : 。 。
- 计算 : $\binom{9}{4} = \frac{9 \times 8 \times 7 \times 6}{4 \times 3 \times 2 \times 1} = 9 \times 2 \times 7 = 126$。 现在需要计算 。。所以 。 (或者,我们可以在计算过程中取模: $\binom{9}{4} = \frac{9 \times 8 \times 7 \times 6}{4 \times 3 \times 2 \times 1} \pmod{13}$ 分子: $(9 \times 8 \times 7 \times 6) \pmod{13} = (72 \times 42) \pmod{13}$ 分子 $\equiv (7 \times 3) \pmod{13} = 21 \equiv 8 \pmod{13}$。 分母: 。 现在需要计算 。 求 使用费马小定理: 。 。 。 : 。 所以 。 因此 。 (验证:$11 \times 6 = 66 = 5 \times 13 + 1 \equiv 1 \pmod{13}$) $\binom{9}{4} \equiv 8 \times 6 \pmod{13} = 48 = 3 \times 13 + 9 \equiv 9 \pmod{13}$。结果一致。)
- 步骤 4:合并结果$$\binom{100}{30} \equiv \binom{7}{2} \times \binom{9}{4} \pmod{13} \equiv 8 \times 9 \pmod{13} $$
- 结论: 除以 13 的余数是 7。
4.4 从实践中领悟定理的运作
通过这几个例子,我们可以清晰地看到卢卡斯定理是如何运作的:
- 进制转换是基础: 正确地将 和 转换为模数 的 进制表示是第一步。
- 分解是核心: 定理将原问题分解为一系列对应位上的小组合数计算。
- 小组合数是计算单元: 每个 的计算是独立的,并且因为 ,可以用基础方法(如阶乘+模逆元)安全计算。
- 乘积是最终答案: 将所有小组合数的模 结果乘起来,最后再模 ,得到最终答案。
- 内置检查: 如果在任何一位 上出现 ,对应 ,直接导致最终结果为 0,这与组合数的性质一致。
现在我们不仅知道了卢卡斯定理是什么,也知道了如何使用它。但一个真正好奇的头脑会问:为什么这个定理是正确的?下一章,我们将深入探索其背后的数学证明。
第五章:卢卡斯定理的证明
理解了卢卡斯定理是什么以及如何使用它之后,探索它为何成立是提升我们数学直觉和深刻理解的关键一步。证明卢卡斯定理有多种方法,其中一种非常优美且常用的是利用生成函数(特别是二项式定理)和模运算的性质。
5.1 思路引导:多项式系数的联系
我们知道二项式定理:
$$(1+x)^n = \binom{n}{0} + \binom{n}{1}x + \binom{n}{2}x^2 + \dots + \binom{n}{n}x^n = \sum_{k=0}^{n} \binom{n}{k} x^k $$这表明组合数 正是多项式 展开后 项的系数。如果我们能在模 的意义下找到 的另一种表达形式,并且这种形式能够很自然地将其 的系数与 联系起来,那么定理就得证了。
我们的目标就是要在模 的世界里研究 。
5.2 关键引理:二项式定理与
让我们先看一个特殊情况:。根据二项式定理:
$$(1+x)^p = \binom{p}{0} + \binom{p}{1}x + \binom{p}{2}x^2 + \dots + \binom{p}{p-1}x^{p-1} + \binom{p}{p}x^p $$现在我们来考察这些系数 在模 下的行为:
- 对于 ,我们有 $\binom{p}{k} = \frac{p!}{k!(p-k)!} = \frac{p \times (p-1)!}{k!(p-k)!}$。 因为 ,所以 和 中的所有因子都小于 。由于 是素数,这意味着 和 都不能被 整除。 因此,在 中,分子有一个因子 ,而分母没有因子 。这意味着 一定是 的倍数。 所以,对于 ,我们有 。
将这些系数代回到 的展开式中,并在模 下观察:
$$(1+x)^p \equiv \binom{p}{0} + \sum_{k=1}^{p-1} \binom{p}{k}x^k + \binom{p}{p}x^p \pmod p $$$$(1+x)^p \equiv 1 + \sum_{k=1}^{p-1} (0 \times x^k) + 1 \times x^p \pmod p $$这是一个非常关键且优美的结果!
5.3 “冰封王座”?不,是“新生之梦”:
上述结果 其实是一个更普遍性质的特例,这个性质有时被(非正式地)称为“新生之梦”(Freshman's Dream),因为它看起来像是初学者可能会犯的错误 ,但在模素数 且指数为 时,它竟然是真的!
定理(同态性质/新生之梦): 若 是素数,则对于任意整数 和 (或者更一般地,对于任何交换环中元素),有:
证明: 使用二项式定理展开 。正如我们在 5.2 节分析 时看到的,只有当 和 时,。 所以,$(a+b)^p \equiv \binom{p}{0} a^p b^0 + \binom{p}{p} a^0 b^p \pmod p = 1 \cdot a^p + 1 \cdot b^p = a^p + b^p \pmod p$。 证毕。
利用这个性质,我们可以通过归纳法得到: $(1+x)^{p^j} \equiv ( (1+x)^{p^{j-1}} )^p \equiv (1+x^{p^{j-1}})^p \equiv 1^p + (x^{p^{j-1}})^p \equiv 1 + x^{p^j} \pmod p$。 即对于任意非负整数 ,都有 。
5.4 核心证明:利用 的 进制展开
现在我们准备好了。将 写成 进制:。 考虑 :
$$(1+x)^n = (1+x)^{n_0} \times (1+x)^{n_1 p} \times \dots \times (1+x)^{n_m p^m} $$利用 $(a \times b) \pmod p = (a \pmod p \times b \pmod p) \pmod p$,以及 ,我们可以在模 下处理每一项:
$$(1+x)^{n_j p^j} = ((1+x)^{p^j})^{n_j} \equiv (1 + x^{p^j})^{n_j} \pmod p \quad (\text{使用 } (1+x)^{p^j} \equiv 1 + x^{p^j} \pmod p) $$所以,
$$(1+x)^n \equiv (1+x)^{n_0} \times (1 + x^p)^{n_1} \times (1 + x^{p^2})^{n_2} \times \dots \times (1 + x^{p^m})^{n_m} \pmod p $$让我们称这个乘积为 。
现在,我们想知道 中 项的系数模 是多少。这等价于求 中 项的系数模 。 将 也写成 进制:。 我们如何在 中得到一个 项呢? 是 个因式的乘积,第 个因式是 。 根据二项式定理,$(1+x^{p^j})^{n_j} = \sum_{i_j=0}^{n_j} \binom{n_j}{i_j} (x^{p^j})^{i_j} = \sum_{i_j=0}^{n_j} \binom{n_j}{i_j} x^{i_j p^j}$。
要从 $P(x) = \prod_{j=0}^{m} \left( \sum_{i_j=0}^{n_j} \binom{n_j}{i_j} x^{i_j p^j} \right)$ 中得到 项,我们需要从第 个因式的展开式中选择一项 (其系数为 ),然后将它们相乘。 也就是说,我们需要选择 。 这样选出的项的乘积是:
$$\left( \binom{n_0}{k_0} x^{k_0 p^0} \right) \times \left( \binom{n_1}{k_1} x^{k_1 p^1} \right) \times \dots \times \left( \binom{n_m}{k_m} x^{k_m p^m} \right) $$$$= \left( \prod_{j=0}^{m} \binom{n_j}{k_j} \right) x^{k_0 p^0 + k_1 p^1 + \dots + k_m p^m} $$$$= \left( \prod_{j=0}^{m} \binom{n_j}{k_j} \right) x^k $$这个过程是否唯一?是的。因为 的 进制表示是唯一的,要凑出 的指数 ,只能是从第 个因式 的展开中恰好选择 这一项(如果 ,则 ,选不出这一项,系数为0)。 因此, 中 项的系数恰好是 。
5.5 严谨性:系数比较的逻辑
我们已经建立了两个等式:
- (原始定义)
- $(1+x)^n \equiv P(x) = \sum_{i=?} \left( \prod_{j=0}^{m} \binom{n_j}{i_j \text{ s.t. } \sum i_j p^j = i} \right) x^i \pmod p$ (推导结果)
比较这两个多项式在模 下 项的系数,我们得到:
$$\binom{n}{k} \equiv \prod_{j=0}^{m} \binom{n_j}{k_j} \pmod p $$这正是卢卡斯定理!证明完毕。
这个证明非常巧妙地运用了二项式定理和模 下 这一关键性质,将问题转化到了 进制表示上,并自然地导出了结果。
现在我们不仅会用卢卡斯定理,还理解了它为什么是正确的。下一步,就是将它打造成我们代码库中的利器。
第六章:C++ 实现与注意事项
掌握了理论和证明,现在是时候将卢卡斯定理转化为实实在在的代码了。我们将使用 C++ 来实现,并讨论一些实现中的细节和注意事项。
6.1 基础组件:快速幂与模逆元
卢卡斯定理的核心计算依赖于求解一系列 。而计算这个小组合数通常需要用到模逆元,模逆元又需要快速幂。所以,我们需要先准备好这两个基础函数。
#include <iostream> // C++ 标准输入输出库,但按要求只用这个可能不全
#include <vector> // 如果需要存储阶乘等
// 定义长整型,方便处理可能较大的 n, k
typedef long long ll;
// 函数:快速幂 power(b, e, m)
// 功能:计算 (b^e) % m
// 时间复杂度: O(log e)
ll power(ll b, ll e, ll m) {
ll res = 1;
b %= m; // 底数先取模
while (e > 0) {
if (e % 2 == 1) res = (res * b) % m; // 如果指数是奇数
b = (b * b) % m; // 底数平方
e /= 2; // 指数除以2
}
return res;
}
// 函数:模逆元 inv(a, m)
// 功能:计算 a 在模 m 下的乘法逆元 (要求 m 是素数)
// 使用费马小定理: a^(m-2) % m
// 时间复杂度: O(log m)
ll inv(ll a, ll m) {
// a %= m; // power 函数内部会处理 a % m
// if (a == 0) return 0; // 0 没有逆元,或者按需处理错误
return power(a, m - 2, m);
}
6.2 基础组件:计算小范围组合数
现在我们需要一个函数来计算 ,其中 且 。由于 较小(小于 ),我们可以直接使用组合数公式 。
有几种实现方式:
方法一:直接计算,每次求逆元 这种方法最简单,不需要预处理。
// 函数:组合数 C(n_prime, k_prime, p)
// 功能:计算 C(n', k') % p (要求 p 是素数, n' < p)
// 直接使用公式 C(n, k) = n! / (k! * (n-k)!)
// 时间复杂度: O(k' + log p) 或 O(n'-k' + log p),即 O(min(k', n'-k') + log p)
ll C_small(ll n_prime, ll k_prime, ll p) {
if (k_prime < 0 || k_prime > n_prime) {
return 0; // 无效输入或 k > n
}
if (k_prime == 0 || k_prime == n_prime) {
return 1; // C(n, 0) = C(n, n) = 1
}
// 优化: C(n, k) = C(n, n-k),选择较小的计算
if (k_prime > n_prime / 2) {
k_prime = n_prime - k_prime;
}
// 计算分子: n' * (n'-1) * ... * (n'-k'+1) mod p
ll numerator = 1;
for (ll i = 0; i < k_prime; ++i) {
numerator = (numerator * (n_prime - i)) % p;
}
// 计算分母: k'! mod p
ll denominator = 1;
for (ll i = 1; i <= k_prime; ++i) {
denominator = (denominator * i) % p;
}
// 计算分母的模逆元
ll inv_denominator = inv(denominator, p);
// 结果 = 分子 * 分母逆元 mod p
return (numerator * inv_denominator) % p;
}
方法二:预处理阶乘和阶乘逆元
如果我们需要在同一个模数 下多次计算小组合数(比如卢卡斯定理中会调用多次 C_small
),预处理阶乘和它们的逆元会更高效。我们只需要预处理到 即可。
const int MAX_P = 100005; // 假设 p 不会超过这个值,根据题目调整
ll fact[MAX_P];
ll invFact[MAX_P];
// 函数:预处理阶乘及其逆元
// 时间复杂度: O(p + log p) 或 O(p) 如果线性求逆元
void precompute(ll p) {
fact[0] = 1;
invFact[0] = 1; // 0! 的逆元通常定义为 1
for (int i = 1; i < p; ++i) {
fact[i] = (fact[i - 1] * i) % p;
}
// 计算 (p-1)! 的逆元,然后反向推导
invFact[p - 1] = inv(fact[p - 1], p);
for (int i = p - 2; i >= 1; --i) {
invFact[i] = (invFact[i + 1] * (i + 1)) % p;
}
// 或者,对每个 fact[i] 单独求逆元,复杂度 O(p log p)
// for (int i = 1; i < p; ++i) {
// invFact[i] = inv(fact[i], p);
// }
}
// 函数:使用预处理数据计算组合数 C(n', k', p)
// 时间复杂度: O(1) 在预处理之后
ll C_small_precomputed(ll n_prime, ll k_prime, ll p) {
if (k_prime < 0 || k_prime > n_prime) {
return 0;
}
if (n_prime >= p || k_prime >= p) {
// 这里理论上不应该发生,因为调用 C_small 的前提是 n', k' < p
// 可以加个错误处理或断言
cerr << "Error: C_small called with arguments >= p" << endl;
return -1; // 表示错误
}
// C(n', k') = n'! / (k'! * (n'-k')!)
// = fact[n'] * invFact[k'] * invFact[n'-k'] % p
ll res = fact[n_prime];
res = (res * invFact[k_prime]) % p;
res = (res * invFact[n_prime - k_prime]) % p;
return res;
}
选择哪种方法取决于具体场景。如果 很小且固定,或者需要大量计算 ,预处理是明智的。如果 较大或者只计算少数几次,直接计算可能更方便。
6.3 核心实现:Lucas(n, k, p)
函数
现在我们可以组装卢卡斯定理的核心函数了。它会迭代地处理 和 的 进制位。
// 函数:卢卡斯定理 Lucas(n, k, p)
// 功能:计算 C(n, k) % p (要求 p 是素数)
// 依赖 C_small 或 C_small_precomputed 函数
// 时间复杂度: O(log_p n * Time(C_small))
// 如果使用预处理的 C_small_precomputed (O(1)),则总时间 O(p + log_p n) 或 O(p log p + log_p n)
// 如果使用未预处理的 C_small (O(p)), 则总时间 O(log_p n * p)
ll Lucas(ll n, ll k, ll p) {
if (k < 0 || k > n) { // 基础情况:k 无效或大于 n
return 0;
}
if (k == 0 || k == n) { // 基础情况:C(n, 0) = C(n, n) = 1
return 1;
}
// 优化:如果 k > n/2, 用 C(n, n-k) 计算可能更快,但这会改变 p 进制位,所以不能直接用
// Lucas 定理本身已经处理了最优分解
ll result = 1; // 初始化结果
// 迭代处理 n 和 k 的 p 进制位
while (n > 0 || k > 0) {
ll ni = n % p; // 当前 n 的最低位 (p 进制)
ll ki = k % p; // 当前 k 的最低位 (p 进制)
// 核心检查:如果 ki > ni,则 C(n, k) = 0 mod p
if (ki > ni) {
return 0;
}
// 计算当前位的组合数 C(ni, ki) % p
// 这里选择调用哪个 C_small 版本
// ll term = C_small(ni, ki, p);
ll term = C_small_precomputed(ni, ki, p); // 假设已预处理
// 累乘结果
result = (result * term) % p;
// 移除最低位,继续处理更高位
n /= p;
k /= p;
}
return result; // 返回最终结果
}
注意: 在使用 Lucas
函数前,如果选择了 C_small_precomputed
,需要先调用 precompute(p)
。
6.4 效率考量:预处理阶乘与逆元
- 预处理的优势: 当你需要计算多个 (其中 可能不同,但 相同),或者像卢卡斯定理这样内部需要多次计算小组合数时,预处理阶乘和逆元可以将每次小组合数计算的时间复杂度从 或 降低到 。总的预处理时间是 (如果用线性逆元)或 (如果用快速幂求逆元)。
- 不预处理的情况: 如果只计算一次 ,且 相对较小,那么不预处理,每次用 的时间计算 也是可以接受的。总复杂度大约是 或更粗略地记为 。
选择哪种策略取决于 的大小、查询次数以及具体题目的时间限制。
6.5 边界处理与常见错误
- 负数 或 : 确保在函数开始就处理这些情况,直接返回 0。
- : 在卢卡斯迭代过程中,一旦发现某个 进制位 ,必须立刻返回 0,因为 会导致最终乘积为 0。
- 模运算: 确保所有中间乘法结果都及时取模 ,防止溢出。
- long long: 可能很大,需要使用
long long
类型。中间计算(如快速幂、乘法)也应使用long long
。 - 模数 必须是素数: 卢卡斯定理及其基于费马小定理的实现都依赖于 是素数。如果模数不是素数,需要使用其他方法(如扩展卢卡斯定理处理素数幂,或中国剩余定理处理合数)。
第七章:应用与问题求解
卢卡斯定理是组合计数模素数问题的有力武器。它不仅能直接解决形如“求 ”的问题,更常常作为解决复杂组合问题的关键子程序。
7.1 直接应用:求
最直接的应用就是题目明确要求计算 ,其中 可能很大(例如 ),而 是一个素数(通常不大,比如 或 )。这种情况下,直接套用我们实现的 Lucas(n, k, p)
函数即可。
例题场景: 假设题目要求计算 。 这里 , (是素数)。 巨大,但 相对较小。直接用卢卡斯定理:
- 预处理 到 的阶乘和逆元模 。时间复杂度 或 。
- 调用
Lucas(10^18, 10^9, 10007)
。 这个函数会进行 次迭代。 $\log_{10007} 10^{18} = \frac{\log_{10} 10^{18}}{\log_{10} 10007} \approx \frac{18}{4} = 4.5$。大约需要 5 次迭代。 每次迭代内部调用C_small_precomputed
,时间复杂度 。 所以卢卡斯计算本身非常快,主要时间在预处理。
总时间复杂度约为 ,对于 来说非常高效。
7.2 间接应用:组合计数问题中的模运算
很多组合计数问题,其最终的答案可能是一个复杂的表达式,其中包含了组合数。如果题目要求结果模一个素数 ,并且表达式中的组合数 满足 很大而 相对小,那么卢卡斯定理就派上用场了。
常见模式:
- 计数 DP + 组合数: 某个 DP 状态转移方程中包含 ,而 可能很大。
- 容斥原理 + 组合数: 容斥计算的项涉及 。
- 生成函数 + 组合数: 从生成函数推导出的系数表达式包含组合数。
- 直接组合分析: 通过组合意义直接推导出涉及 的公式。
在这些场景下,当需要计算 时,如果 ,就应该考虑使用卢卡斯定理作为计算子模块。
7.3 典型例题:SDOI2010 古代猪文 (P2480)
这道题是卢卡斯定理与数论知识(尤其是中国剩余定理 CRT 和费马小定理/欧拉定理)结合的经典范例。
题目要求: 计算 。
分析:
- 目标形式: 我们需要计算 ,其中 是一个大素数,指数 。
- 指数处理(费马小定理): 根据费马小定理, (前提是 ;如果 ,结果直接是 0)。 所以,核心任务是计算指数 ,即 $E' = \left( \sum_{d|n} \binom{n}{d} \right) \pmod{M-1}$。
- 模数 的挑战: 。这个数不是素数!我们不能直接在模 下计算组合数。 。其中 都是素数。(注: 需要验证是否为素数,它确实是。题目给的模数通常设计成这样。)
- 中国剩余定理 (CRT) 登场: 我们无法直接计算 。但是,如果我们能分别计算出 模 的结果,即:
- $E_1 = E \pmod 2 = \left( \sum_{d|n} \binom{n}{d} \right) \pmod 2$
- $E_2 = E \pmod 3 = \left( \sum_{d|n} \binom{n}{d} \right) \pmod 3$
- $E_3 = E \pmod {p_3} = \left( \sum_{d|n} \binom{n}{d} \right) \pmod {p_3}$ 那么,我们就可以利用中国剩余定理,根据 , , 这三个同余方程组,唯一确定 在模 下的值 。
- 卢卡斯定理的应用: 现在的问题转化为计算 分别模 , , 的结果。 对于每一个素数模 (),我们需要: a. 找到 的所有约数 。(可以通过质因数分解 然后 DFS 生成所有约数。) b. 对于每个约数 ,计算 。这里,由于 可能远大于 ,就需要使用卢卡斯定理! c. 将所有 的结果加起来,再模 ,得到 。
- 最终计算: a. 计算 。 b. 使用 CRT 解出 。 c. 使用快速幂计算 。
这个例子完美展示了卢卡斯定理如何作为一个关键模块,嵌入在更复杂的数论问题框架中。 它解决了在大组合数计算中遇到素数模的问题。
7.4 思维拓展:与其他算法(如CRT)的结合
正如“古代猪文”所示,卢卡斯定理经常需要和其他数论工具结合使用。
- 与中国剩余定理 (CRT): 当最终模数是合数 时,可以将其分解质因数 。如果能分别计算出答案模每个 的结果,就可以用 CRT 合并。如果 (即 无平方因子),则模 的计算可能用到卢卡斯定理。如果 ,则需要扩展卢卡斯定理(见第八章)。
- 与容斥原理: 某些计数问题用容斥原理展开后,每一项可能是带符号的组合数 。如果需要模素数 且 ,计算每一项时使用卢卡斯定理。
- 与动态规划: 如果 DP 转移涉及组合数 ,其中 来自子问题且可能很大,转移时调用卢卡斯定理。
理解卢卡斯定理不仅是学会一个孤立的算法,更是掌握了一个处理特定类型模运算的强大组件,能够灵活地嵌入到各种算法框架中,解决更广泛的问题。
下一章,我们将总结卢卡斯定理的价值,并展望其局限与扩展。
第八章:超越边界 - 扩展、总结与展望
我们已经一起走过了卢卡斯定理的发现、理解、证明、实现和应用的全过程。在结束这段旅程之前,让我们再次回顾它的核心价值,认识它的局限,并稍微眺望一下更远方的风景。
8.1 卢卡斯定理的核心价值回顾
卢卡斯定理之所以在算法竞赛和组合数论中占有重要地位,主要在于它提供了计算大组合数模素数 的高效方法:
- 解决了可行性问题: 它巧妙地避免了当 时,直接计算 因分母包含 因子而无法求逆元的问题。
- 实现了分治降维: 通过将 转换为 进制,它把一个可能涉及非常大数字的计算,分解为一系列 的小组合数计算。这些小组合数可以用基本方法(阶乘+模逆元)安全高效地求解。
- 效率高: 结合预处理阶乘和逆元,计算 的时间复杂度主要由预处理的 或 以及对数次迭代 构成,对于 很大但 相对较小的情况非常高效。
- 应用广泛: 不仅能直接求解 ,更能作为子程序嵌入到更复杂的计数问题中,与其他数论工具(如 CRT)或算法(如 DP、容斥)结合使用。
掌握卢卡斯定理,意味着你拥有了一把解锁特定类型组合计数难题的钥匙。
8.2 局限性:素数模数的限制
卢卡斯定理最主要的局限性在于它只适用于模数 是素数的情况。其证明的核心依赖于 以及 (for ) 这些只对素数 成立的性质。
如果题目要求的模数 不是素数,怎么办?
- 情况一: 可以分解为互不相同的素数 的乘积。 即 。这时可以使用中国剩余定理 (CRT)。分别计算 (这里可以用卢卡斯定理),得到 个同余方程。然后用 CRT 合并这些结果,得到最终的 。
- 情况二: 包含素数的幂次因子,即 ,其中至少有一个 。 这时,我们仍然可以用 CRT。问题转化为计算 ,其中 是素数,。当 时,用卢卡斯定理。当 时,标准的卢卡斯定理不再适用。我们需要更强大的工具——扩展卢卡斯定理 (Extended Lucas's Theorem)。
8.3 扩展卢卡斯定理 (模素数幂)
计算 的问题比模素数 更复杂。其核心思想是:
- 将 写出来。
- 对于分子 和分母 ,分别计算它们包含素数因子 的次数(这可以通过 Legendre 定理计算: 中 的次数是 )。设 中 的次数分别是 。那么 中 的次数就是 。
- 将 中所有的因子 都“提取”出来。例如,,其中 是 中所有不含 因子的部分相乘。同样得到 和 。
- 那么 $\binom{n}{k} = \frac{p^{v_n} N'}{p^{v_k} K' \times p^{v_{n-k}} R'} = p^{v_n - v_k - v_{n-k}} \times \frac{N'}{K' R'} = p^v \times \frac{N'}{K' R'}$。
- 关键在于计算 。这需要计算 分别模 的值,以及 和 模 的逆元。计算 (即 去掉所有 因子后模 )通常涉及对 到 中不被 整除的数按模 的周期性进行分组计算。求逆元则需要用到扩展欧几里得算法,因为 与 互质。
- 最后将 乘回去,得到 的结果。
扩展卢卡斯定理的过程比标准卢卡斯定理复杂得多,涉及更多的数论知识(Legendre 定理、模 下的阶乘计算、扩展欧几里得求逆元)。但它将我们处理组合数模运算的能力从素数模扩展到了素数幂模,再结合 CRT,就能处理任意整数模。
经典例题
P3807 【模板】卢卡斯定理/Lucas 定理
题目背景
这是一道模板题。
题目描述
给定整数 的值,求出 的值。
输入数据保证 为质数。
注: 表示组合数。
输入格式
本题有多组数据。
第一行一个整数 ,表示数据组数。
对于每组数据:
一行,三个整数 。
输出格式
对于每组数据,输出一行,一个整数,表示所求的值。
输入输出样例 #1
输入 #1
2
1 2 5
2 1 5
输出 #1
3
3
说明/提示
对于 的数据,,。
参考代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll; // 使用 typedef 定义 ll 为 long long
const ll maxn = 1e5 + 15;
ll n, m, p, T, fac[maxn];
// 快速幂求逆元,使用费马小定理
ll pow(ll b, ll q) {
ll ans = 1;
for(; q; q >>= 1, b = b * b % p) {
if(q & 1) ans = ans * b % p;
}
return ans % p;
}
// 计算组合数 C(n, m) % p
ll C(ll n, ll m) {
if (m > n) return 0;
return fac[n] * pow(fac[m] * fac[n - m] % p, p - 2) % p;
}
// 卢卡斯定理
ll Lucas(ll n, ll m) {
if (m == 0) return 1;
return (Lucas(n / p, m / p) * C(n % p, m % p)) % p;
}
int main() {
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0); // 关闭同步流,提高效率
cin >> T;
while (T--) {
cin >> n >> m >> p;
fac[0] = 1;
for (int i = 1; i <= p; i++) {
fac[i] = fac[i - 1] * i % p; // 预处理阶乘
}
cout << Lucas(n + m, m) % p << '\n'; // 输出结果
}
return 0;
}
P4720 【模板】扩展卢卡斯定理/exLucas
题目背景
这是一道模板题。
题目描述
求
其中 为组合数。
输入格式
一行三个整数 ,含义由题所述。
输出格式
一行一个整数,表示答案。
输入输出样例 #1
输入 #1
5 3 3
输出 #1
1
输入输出样例 #2
输入 #2
666 233 123456
输出 #2
61728
说明/提示
对于 的数据,,,不保证 是质数。
参考代码
#include <bits/stdc++.h>
using namespace std;
// 使用 long long 处理大整数
using ll = long long;
// 快速幂 (a^b % md)
// 参数: 底数a, 指数b, 模数md
// 返回: (a^b) % md
ll qp(ll a, ll b, ll md) {
ll res = 1; // 结果初始化为1
a %= md; // 底数先取模
while (b > 0) {
if (b & 1) res = (__int128)res * a % md; // 如果b是奇数,累乘当前a
// 使用 __int128 防止 res * a 溢出 long long
a = (__int128)a * a % md; // a自乘并取模
b >>= 1; // b右移一位 (相当于b /= 2)
}
return res;
}
// 扩展欧几里得算法
// 求解 ax + by = gcd(a, b),并返回 gcd(a, b)
// x 和 y 是通过引用传递的参数,用于接收结果
ll exgcd(ll a, ll b, ll &x, ll &y) {
if (b == 0) { // 递归基:当 b=0 时,gcd(a, 0) = a
x = 1; // 此时 ax + 0y = a,一组解是 x=1, y=0
y = 0;
return a; // 返回最大公约数 a
}
ll x1, y1; // 存储下一层递归的解
ll d = exgcd(b, a % b, x1, y1); // 递归调用
// 从下一层的解 (x1, y1) 推导当前层的解 (x, y)
// ax + by = b*x1 + (a % b)*y1
// = b*x1 + (a - floor(a/b)*b)*y1
// = a*y1 + b*(x1 - floor(a/b)*y1)
x = y1;
y = x1 - (a / b) * y1;
return d; // 返回最大公约数
}
// 模逆元
// 计算 a 在模 md 下的乘法逆元 a^-1
// 要求 a 和 md 互质
ll inv(ll a, ll md) {
ll x, y;
ll g = exgcd(a, md, x, y); // 调用 exgcd 求解 ax + md*y = gcd(a, md)
// assert(g == 1); // 对于求逆元,必须保证 gcd(a, md) == 1
// x 是 ax = 1 (mod md) 的一个解,可能为负
return (x % md + md) % md; // 将 x 调整到 [0, md-1] 范围内
}
// 勒让德定理 (Legendre's Formula)
// 计算 n! 中质因子 pi 的幂次 (v_pi(n!))
ll leg(ll n, ll pi) {
ll cnt = 0; // 计数器
while (n) {
cnt += n / pi; // 加上 floor(n/pi)
n /= pi; // n 更新为 n/pi
}
return cnt; // 返回总幂次
}
// 计算 n! / pi^leg(n, pi) mod pk
// 即计算 n! 除去所有 pi 因子后,对 pk 取模的结果
ll fac(ll n, ll pi, ll pk) {
if (n == 0) return 1; // 0! = 1
ll res = 1;
// 计算 [1, pk] 范围内所有与 pi 互质的数的乘积模 pk
for (ll i = 1; i <= pk; ++i) {
if (i % pi != 0) { // 只乘与 pi 互质的数
res = (__int128)res * i % pk;
}
}
// 这个乘积块会重复 n / pk 次
res = qp(res, n / pk, pk);
// 计算剩余部分 [1, n % pk] 中与 pi 互质的数的乘积
for (ll i = 1; i <= n % pk; ++i) {
if (i % pi != 0) {
res = (__int128)res * i % pk;
}
}
// 递归计算 fac(n / pi, pi, pk)
// 因为 n! = (pi * 2pi * ... * floor(n/pi)pi) * (其他数)
// = pi^floor(n/pi) * floor(n/pi)! * (其他数)
// 我们已经处理了 (其他数) 部分,现在需要乘上 floor(n/pi)! 除去 pi 因子的部分
return (__int128)res * fac(n / pi, pi, pk) % pk;
}
// 计算 C(n, m) mod pk, 其中 pk = pi^k 是一个质数幂
ll calC(ll n, ll m, ll pi, ll pk) {
if (m < 0 || m > n) return 0; // 组合数定义,m 不能小于 0 或大于 n
// 计算 n!, m!, (n-m)! 除去 pi 因子后的部分模 pk
ll fn = fac(n, pi, pk);
ll fm = fac(m, pi, pk);
ll fnm = fac(n - m, pi, pk);
// 计算 C(n, m) 中 pi 因子的最终幂次
// v_pi(C(n, m)) = v_pi(n!) - v_pi(m!) - v_pi((n-m)!)
ll vp = leg(n, pi) - leg(m, pi) - leg(n - m, pi);
// 计算结果
// C(n,m) = [n!/pi^leg(n)] / [m!/pi^leg(m)] / [(n-m)!/pi^leg(n-m)] * pi^vp (mod pk)
ll res = (__int128)fn * inv(fm, pk) % pk; // 乘以 m! 部分的逆元
res = (__int128)res * inv(fnm, pk) % pk; // 乘以 (n-m)! 部分的逆元
res = (__int128)res * qp(pi, vp, pk) % pk; // 乘以 pi^vp 部分
return res;
}
// 中国剩余定理 (CRT)
// 求解线性同余方程组 x = a[i] (mod Pk[i])
// 要求所有 Pk[i] 两两互质 (题目中质数幂次天然满足)
// a: 存储余数的向量 (C(n, m) mod pk_i)
// Pk: 存储模数的向量 (pk_i = pi_i ^ ki_i)
ll crt(const vector<ll>& a, const vector<ll>& Pk) {
ll M = 1; // 所有模数的乘积 M = p
for (ll mod_i : Pk) {
M *= mod_i; // 计算 M,注意 p <= 10^6,所以 M <= 10^6,不会溢出 ll
}
__int128 ans = 0; // 最终结果,使用 __int128 防止中间累加溢出
for (size_t i = 0; i < a.size(); ++i) {
ll Mi = M / Pk[i]; // Mi = M / pk_i
ll invMi = inv(Mi, Pk[i]); // 计算 Mi 在模 Pk[i] 下的逆元
// CRT 的合并公式: ans = sum(a_i * Mi * inv(Mi, pk_i)) mod M
// 使用 __int128 防止 a[i] * Mi * invMi 溢出 ll
ans = (ans + (__int128)a[i] * Mi % M * invMi % M) % M;
}
return (ll)((ans + M) % M); // 保证结果非负并返回 long long
}
int main() {
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
ll n, m, p;
cin >> n >> m >> p;
vector<ll> a; // 存储每个 C(n, m) mod pk_i 的结果
vector<ll> Pk; // 存储 p 的每个质因子幂次 pk_i
ll tmp = p; // 临时变量用于分解质因数 p
for (ll i = 2; i * i <= tmp; ++i) { // 试除法分解质因数
if (tmp % i == 0) { // 如果 i 是 tmp 的因子
ll pi = i; // 当前质因子
ll pk = 1; // 当前质因子的幂次 (pi^k)
while (tmp % i == 0) { // 计算 pi 的最高次幂 pk
pk *= i;
tmp /= i;
}
Pk.push_back(pk); // 存储模数 pk_i
// 计算 C(n, m) mod pk,并存入结果向量 a
a.push_back(calC(n, m, pi, pk));
}
}
// 如果分解结束后 tmp 仍然大于 1,说明它本身是一个大于 sqrt(p) 的质因子
if (tmp > 1) {
ll pi = tmp; // 这个质因子就是 tmp 本身
ll pk = tmp; // 幂次就是 tmp (指数为 1)
Pk.push_back(pk);
a.push_back(calC(n, m, pi, pk));
}
// 所有模数下的结果都已计算完毕,使用 CRT 合并
cout << crt(a, Pk) << endl; /
return 0;
}
P2480 [SDOI2010] 古代猪文
题目背景
“在那山的那边海的那边有一群小肥猪。他们活泼又聪明,他们调皮又灵敏。他们自由自在生活在那绿色的大草坪,他们善良勇敢相互都关心……”
——选自猪王国民歌
很久很久以前,在山的那边海的那边的某片风水宝地曾经存在过一个猪王国。猪王国地理位置偏僻,实施的是适应当时社会的自给自足的庄园经济,很少与外界联系,商贸活动就更少了。因此也很少有其他动物知道这样一个王国。
猪王国虽然不大,但是土地肥沃,屋舍俨然。如果一定要拿什么与之相比的话,那就只能是东晋陶渊明笔下的大家想象中的桃花源了。猪王勤政爱民,猪民安居乐业,邻里和睦相处,国家秩序井然,经济欣欣向荣,社会和谐稳定。和谐的社会带给猪民们对工作火红的热情和对未来的粉色的憧憬。
小猪iPig是猪王国的一个很普通的公民。小猪今年10岁了,在大肥猪学校上小学三年级。和大多数猪一样,他不是很聪明,因此经常遇到很多或者稀奇古怪或者旁人看来轻而易举的事情令他大伤脑筋。小猪后来参加了全猪信息学奥林匹克竞赛(Pig Olympiad in Informatics, POI),取得了不错的名次,最终保送进入了猪王国大学(Pig Kingdom University, PKU)深造。
现在的小猪已经能用计算机解决简单的问题了,比如能用P++语言编写程序计算出A + B的值。这个“成就”已经成为了他津津乐道的话题。当然,不明真相的同学们也开始对他刮目相看啦~
小猪的故事就将从此展开,伴随大家两天时间,希望大家能够喜欢小猪。
题目描述
猪王国的文明源远流长,博大精深。
iPig 在大肥猪学校图书馆中查阅资料,得知远古时期猪文文字总个数为 。当然,一种语言如果字数很多,字典也相应会很大。当时的猪王国国王考虑到如果修一本字典,规模有可能远远超过康熙字典,花费的猪力、物力将难以估量。故考虑再三没有进行这一项劳猪伤财之举。当然,猪王国的文字后来随着历史变迁逐渐进行了简化,去掉了一些不常用的字。
iPig 打算研究古时某个朝代的猪文文字。根据相关文献记载,那个朝代流传的猪文文字恰好为远古时期的 ,其中 是 的一个正约数(可以是 或 )。不过具体是哪 ,以及 是多少,由于历史过于久远,已经无从考证了。
iPig 觉得只要符合文献,每一种 都是有可能的。他打算考虑到所有可能的 。显然当 等于某个定值时,该朝的猪文文字个数为 。然而从 个文字中保留下 个的情况也是相当多的。iPig 预计,如果所有可能的 的所有情况数加起来为 的话,那么他研究古代文字的代价将会是 。
现在他想知道猪王国研究古代文字的代价是多少。由于 iPig 觉得这个数字可能是天文数字,所以你只需要告诉他答案除以 的余数就可以了。
输入格式
一行两个正整数 。
输出格式
输出一行一个整数表示答案。
输入输出样例 #1
输入 #1
4 2
输出 #1
2048
说明/提示
- 对于 的数据,;
- 对于 的数据,;
- 对于 的数据,;
- 对于 的数据,。
参考代码
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int P = 999911659; // 最终模数 (题目给定,是质数)
const int Pm1 = P - 1; // 指数需要对其取模 (合数)
// Pm1 = 999911658 = 2 * 3 * 4679 * 35617
const int MAXF = 35617 + 5; // 预处理阶乘最大范围 (大于Pm1的最大质因子)
ll fac[MAXF]; // 阶乘数组 fac[i] = i! mod p
ll ifac[MAXF]; // 阶乘逆元数组 ifac[i] = (i!)^-1 mod p
// 函数名: qp (Quick Power)
// 功能: 计算 (a^b) % md
ll qp(ll a, ll b, ll md) {
ll res = 1;
a %= md;
while (b > 0) {
if (b & 1) res = (__int128)res * a % md; // 防溢出
a = (__int128)a * a % md; // 防溢出
b >>= 1;
}
return res;
}
// 函数名: exgcd (Extended Euclidean Algorithm)
// 功能: 求解 ax + by = gcd(a, b)
ll exgcd(ll a, ll b, ll &x, ll &y) {
if (b == 0) {
x = 1; y = 0;
return a;
}
ll x1, y1;
ll d = exgcd(b, a % b, x1, y1);
x = y1;
y = x1 - (a / b) * y1;
return d;
}
// 函数名: inv (Modular Inverse)
// 功能: 计算 a 在模 md 下的逆元
ll inv(ll a, ll md) {
ll x, y;
exgcd(a, md, x, y);
return (x % md + md) % md; // 保证结果为正
}
// 函数名: pre (Precomputation)
// 功能: 预处理阶乘和阶乘逆元 mod md (md <= MAXF)
void pre(ll md) {
fac[0] = 1;
ifac[0] = 1;
for (int i = 1; i < md; ++i) {
fac[i] = (__int128)fac[i - 1] * i % md;
}
// 倒推计算逆元: (i!) = ((i+1)! / (i+1)) => (i!)^-1 = ((i+1)!)^-1 * (i+1)
ifac[md - 1] = inv(fac[md - 1], md); // 先求最后一个的逆元
for (int i = md - 2; i >= 1; --i) {
ifac[i] = (__int128)ifac[i + 1] * (i + 1) % md;
}
}
// 函数名: C (Combinations)
// 功能: 计算小组合数 C(n, m) mod md (要求 n, m < md,且 fac/ifac 已预处理)
ll C(ll n, ll m, ll md) {
if (m < 0 || m > n) return 0; // 处理边界和非法情况
if (m == 0 || m == n) return 1; // C(n, 0) = C(n, n) = 1
if (m > n / 2) m = n - m; // 优化: C(n, m) = C(n, n-m)
// 使用预处理的阶乘和逆元计算
// C(n, m) = n! / (m! * (n-m)!) mod md
ll res = fac[n];
res = (__int128)res * ifac[m] % md;
res = (__int128)res * ifac[n - m] % md;
return res;
}
// 函数名: lucas (Lucas' Theorem)
// 功能: 计算 C(n, k) mod md (md 必须是质数)
ll lucas(ll n, ll k, ll md) {
if (k < 0 || k > n) return 0; // 处理边界
if (k == 0) return 1; // C(n, 0) = 1
ll res = 1;
// 循环处理 n 和 k 的 md 进制表示的每一位
while (k > 0) {
ll ni = n % md; // n 在 md 进制下的当前位
ll ki = k % md; // k 在 md 进制下的当前位
// 根据卢卡斯定理,如果某一位 ni < ki, 则 C(ni, ki) = 0,最终结果为 0
if (ni < ki) return 0;
// 累乘每一位的组合数结果 C(ni, ki) mod md
res = (__int128)res * C(ni, ki, md) % md;
// 进入下一位 (更高位)
n /= md;
k /= md;
}
return res;
}
// 函数名: crt (Chinese Remainder Theorem)
// 功能: 求解线性同余方程组 x = rem[i] (mod mods[i])
// 要求: mods[i] 两两互质
ll crt(const vector<ll>& rem, const vector<ll>& mods) {
ll M = 1; // 所有模数之积 (即 Pm1)
for (ll mi : mods) M *= mi;
__int128 ans = 0; // 使用 __int128 防止中间结果溢出
for (size_t i = 0; i < rem.size(); ++i) {
ll Mi = M / mods[i]; // M_i = M / m_i
ll invMi = inv(Mi, mods[i]); // t_i = Mi^-1 (mod m_i)
// CRT 合并公式: ans = sum(rem_i * M_i * t_i) mod M
ans = (ans + (__int128)rem[i] * Mi % M * invMi % M) % M;
}
// 保证结果非负
return (ll)((ans % M + M) % M);
}
int main() {
ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
ll n, g;
cin >> n >> g;
// 特判: 如果 g 是 P 的倍数, g^p mod P = 0 (p > 0)
if (g % P == 0) {
cout << 0 << endl;
return 0;
}
vector<ll> divs; // 存储 n 的所有约数 k
// 找 n 的所有约数,复杂度 O(sqrt(n))
for (ll i = 1; i * i <= n; ++i) {
if (n % i == 0) {
divs.push_back(i);
if (i * i != n) { // 防止对完全平方数重复添加其平方根
divs.push_back(n / i);
}
}
}
// P-1 = 999911658 的质因子 (都是一次幂)
vector<ll> mods = {2, 3, 4679, 35617};
vector<ll> rem(mods.size()); // 存储 p mod mods[i] 的结果
// 对 P-1 的每个质因子 pi,计算 p mod pi
for (size_t i = 0; i < mods.size(); ++i) {
ll md = mods[i]; // 当前计算用的模数 (质数)
pre(md); // 预处理阶乘和逆元 mod md
ll pi = 0; // 用于累加 sum_{k|n} C(n, k) mod md
for (ll d : divs) { // 遍历 n 的所有约数 k (这里用变量 d)
ll cnk = lucas(n, d, md); // 使用卢卡斯定理计算 C(n, d) mod md
pi = (pi + cnk) % md; // 累加结果
}
rem[i] = pi; // 存储 sum(C(n, k)) mod md 的结果
}
// 使用 CRT 合并 rem[i] (p mod qi) 和 mods[i] (qi)
// 得到指数 p' = p mod (P-1)
ll exp = crt(rem, mods);
// 计算最终答案 g^p' mod P
ll ans = qp(g, exp, P);
cout << ans << endl; // 输出答案
return 0; // 程序正常结束
}