尧图网站建设 
1. 项目概述从一道题到一个工具箱如果你在刷算法题特别是组合数学相关的题目时遇到了一个模数p不是质数甚至可能是一个合数的情况需要计算组合数C(n, m) mod p那么“扩展卢卡斯定理”就是你绕不开的一道坎。p4720这个题号在很多在线评测系统中指的就是要求实现扩展卢卡斯定理的模板题。它不像普通的卢卡斯定理那样要求模数是质数其核心思想是将非质数模数p分解质因数分别求出组合数对每个质数幂p_i^{k_i}取模的结果最后再用中国剩余定理CRT合并得到最终答案。这个“模板”的价值远不止于解决一道题。它是一套完整的、用于处理“模数为任意正整数时求组合数”的通用解决方案。掌握了它你就拥有了一个强大的数学工具可以应对竞赛、面试乃至一些特定科研场景中更广泛的模运算问题。今天我们就来彻底拆解这个模板不仅告诉你代码怎么写更要讲清楚每一个步骤背后的数学原理和工程化实现的细节让你能真正理解并灵活运用它。2. 核心思路拆解化整为零分而治之面对一个非质数模数p直接计算组合数C(n, m) n! / (m! * (n-m)!)的模运算是行不通的因为分母的逆元可能不存在当分母与模数不互质时。扩展卢卡斯定理的精妙之处在于它采用了一种“化整为零”的策略。2.1 第一步质因数分解模数假设我们需要计算C(n, m) mod p其中p是任意正整数。 首先我们对p进行质因数分解p p1^{k1} * p2^{k2} * ... * pt^{kt}其中pi是质数ki是正整数。为什么这么做因为对于每一个质数幂pi^{ki}我们有可能利用数论知识在一个“局部”的模数下处理阶乘中包含pi因子的特殊情况。而最终的结果可以通过中国剩余定理从这些“局部”解唯一地确定“全局”解mod p。注意这里质因数分解的算法效率很重要。对于p在1e9范围内的题目通常使用试除法即可如果p可能更大则需要更高效的算法如 Pollard-Rho但在竞赛模板题中p的范围通常是可控的。2.2 第二步分别求解模每个质数幂现在问题转化为对于每一个质数幂q pi^{ki}计算C(n, m) mod q。 这是整个算法的核心难点。我们无法直接计算因为n!, m!, (n-m)!可能与q不互质它们包含因子pi。解决思路剔除因子分离计算我们以计算n! mod q为例其中q p^k。分离出所有 p 的因子将n!中所有是p的倍数的数都提取出来。例如n! 1*2*3*...*n。其中p, 2p, 3p, ..., floor(n/p)*p都是p的倍数。把它们全部提取出来一共可以提取出floor(n/p)个p。提取后式子变成了n! p^{floor(n/p)} * (1*2*3*...*floor(n/p)) * [所有不能被 p 整除的数的乘积]。注意(1*2*3*...*floor(n/p))这部分又是floor(n/p)!我们可以递归地处理它处理剩余部分与 p 互质的部分剩下的数是1到n中所有不能被p整除的数。它们的乘积对q取模存在一个重要的周期性规律。考虑数列1, 2, ..., p-1, p1, p2, ..., 2p-1, 2p1, ...。以q 3^2 9为例我们看1到9中与3互质的数1,2,4,5,7,8。它们的乘积1*2*4*5*7*8 mod 9 224 mod 9 8。可以发现1到q中与p互质的数的乘积模q是一个定值记作prod(q)。对于q9prod(9) 8。那么对于1到n中所有与p互质的数我们可以将它们分成若干完整的周期每个周期长度为q和一个不完整的尾巴。完整周期的乘积就是prod(q)^{floor(n/q)} mod q。不完整尾巴即最后不足一个周期的部分需要暴力计算但长度小于q可以接受。将以上两步结合起来我们就得到了计算n! mod q的递归公式通常被称为fact(n, p, q)函数fact(n, p, q) (fact(floor(n/p), p, q) * prod(q)^{floor(n/q)} * tail(n mod q, p, q)) mod q其中tail(r, p, q)计算1到r中所有与p互质的数的乘积模q。2.3 第三步应用公式与合并结果有了计算n! mod q的能力我们就可以计算C(n, m) mod qC(n, m) mod q fact(n, p, q) * inverse(fact(m, p, q), q) * inverse(fact(n-m, p, q), q) * p^{e_n - e_m - e_{n-m}} mod q等等最后一项p^{e_n - e_m - e_{n-m}}是什么别忘了 p 的因子在第二步计算fact函数时我们刻意把所有的p因子都提取出来了。fact(n, p, q)返回的结果是n!剔除所有p因子后剩余部分模q的值。同时我们可以通过递归过程很容易地计算出n!中p因子的总指数e(n) floor(n/p) floor(n/p^2) floor(n/p^3) ...。因此真正的n!对模数q的“贡献”是fact(n, p, q) * p^{e(n)} mod q。 所以组合数C(n, m) n! / (m! * (n-m)!)对模数q的“贡献”是[fact(n, p, q) * p^{e(n)}] / [fact(m, p, q) * p^{e(m)} * fact(n-m, p, q) * p^{e(n-m)}] mod q化简后得到C(n, m) mod q fact(n, p, q) * inverse(fact(m, p, q), q) * inverse(fact(n-m, p, q), q) * p^{e(n) - e(m) - e(n-m)} mod q这里inverse(a, q)表示a在模q意义下的乘法逆元。由于fact函数返回的值是与p互质的我们剔除了所有p因子所以它必然与q互质因为q只包含质因子p逆元一定存在可以用扩展欧几里得算法求解。最后我们得到了t个同余方程C(n, m) ≡ a_i (mod p_i^{k_i}), 对于i 1, 2, ..., t。 其中a_i就是我们上面计算出的C(n, m) mod p_i^{k_i}。2.4 第四步中国剩余定理CRT合并现在我们有了一个同余方程组。中国剩余定理告诉我们如果模数两两互质这里p_i^{k_i}显然互质那么这个方程组在模p下有唯一解。 合并公式为 设M p p1^{k1} * p2^{k2} * ... * pt^{kt}。 设M_i M / p_i^{k_i}。 设t_i是M_i在模p_i^{k_i}意义下的逆元即M_i * t_i ≡ 1 (mod p_i^{k_i})。 则方程组的解为x ≡ Σ (a_i * M_i * t_i) (mod M)。这个x就是我们要的C(n, m) mod p。3. 模板实现与关键代码解析理解了原理我们来看代码实现。一个健壮的扩展卢卡斯模板需要以下几个函数3.1 快速幂与扩展欧几里得这是基础工具。// 快速幂 (a^b % mod) long long pow_mod(long long a, long long b, long long mod) { long long res 1; while (b) { if (b 1) res res * a % mod; a a * a % mod; b 1; } return res; } // 扩展欧几里得求解 ax by gcd(a, b)并返回gcd(a,b) // 同时通过引用返回x, y使得 ax by gcd(a,b) long long exgcd(long long a, long long b, long long x, long long y) { if (b 0) { x 1; y 0; return a; } long long d exgcd(b, a % b, y, x); y - a / b * x; return d; } // 求 a 在模 mod 下的逆元要求 gcd(a, mod) 1 long long inv(long long a, long long mod) { long long x, y; exgcd(a, mod, x, y); return (x % mod mod) % mod; // 调整到 0~mod-1 }3.2 计算 fact(n, p, pk) 和 指数 e(n)这个函数计算n!中剔除所有p因子后剩余部分模pk的值同时通过引用或全局变量记录p因子的总指数。// pk p^k long long fact(long long n, long long p, long long pk) { if (n 0) return 1; long long res 1; // 1. 计算完整周期的乘积 // 周期长度为 pk每个周期中与p互质部分的乘积是固定的我们预处理为 prod_pk // 假设我们已经预处理好了 prod[pk] (1*2*...*pk)中与p互质的数之积 % pk // 在实际代码中prod(pk) 可以递归或循环计算这里为了清晰假设已获得 res res * pow_mod(prod_pk, n / pk, pk) % pk; // 2. 计算最后一个不完整周期 // 暴力计算 1 到 (n % pk) 中与p互质的数的乘积 for (long long i 1; i n % pk; i) { if (i % p ! 0) { // 与p互质 res res * i % pk; } } // 3. 递归处理 floor(n/p)! // 注意递归的部分是 floor(n/p)!它可能还包含p的因子需要继续提取 res res * fact(n / p, p, pk) % pk; return res; } // 计算 n! 中质因子 p 的指数 long long e_func(long long n, long long p) { long long cnt 0; while (n) { cnt n / p; n / p; } return cnt; }在实际模板中prod_pk通常不会全局预处理因为pk随分解变化而是在fact函数内部通过循环计算出来或者用一个辅助函数calc_prod(n, p, pk)来计算1..n中与p互质的数之积模pk。为了效率常见的写法是将fact函数写成迭代形式并同时计算乘积和周期。3.3 计算 C(n, m) mod pk这是对单个质数幂模数的求解。// 计算 C(n, m) % pk其中 pk p^k long long C_mod_pk(long long n, long long m, long long p, long long pk) { if (m n) return 0; // 计算三个阶乘剔除p因子后的值 long long fact_n fact(n, p, pk); long long fact_m fact(m, p, pk); long long fact_nm fact(n - m, p, pk); // 计算p因子的指数差 long long e_n e_func(n, p); long long e_m e_func(m, p); long long e_nm e_func(n - m, p); long long e_diff e_n - e_m - e_nm; // 组合计算结果 long long res fact_n * inv(fact_m, pk) % pk; res res * inv(fact_nm, pk) % pk; res res * pow_mod(p, e_diff, pk) % pk; // 注意如果e_diff为负则结果为0因为分子中p因子不够 return res; }实操心得e_diff可能为负数吗理论上对于组合数C(n, m)其值一定是整数所以p因子的指数差e_n - e_m - e_{n-m}一定0。如果计算出现负数只可能是m n的非法输入我们已经做了判断。但在一些特殊写法或递归中确保指数计算正确很重要。3.4 质因数分解与CRT合并这是主函数exLucas(n, m, p)的逻辑。long long exLucas(long long n, long long m, long long p) { // 1. 质因数分解 p vectorpairlong long, long long factors; // 存储 (质数pi, 幂次ki) long long temp p; for (long long i 2; i * i temp; i) { if (temp % i 0) { long long cnt 0; long long pk 1; while (temp % i 0) { temp / i; cnt; pk * i; } factors.push_back({i, pk}); // 这里存储 (pi, pi^ki) } } if (temp 1) { factors.push_back({temp, temp}); } // 2. 分别求解每个质数幂下的答案 vectorlong long a, mods; for (auto fac : factors) { long long pi fac.first; long long pki fac.second; long long ai C_mod_pk(n, m, pi, pki); a.push_back(ai); mods.push_back(pki); } // 3. CRT 合并 long long res 0; long long M p; // 总模数 for (int i 0; i factors.size(); i) { long long Mi M / mods[i]; long long ti inv(Mi, mods[i]); // Mi 在模 mods[i] 下的逆元 res (res a[i] * Mi % M * ti % M) % M; } return res; }4. 优化技巧与边界处理一个工业级的模板不能只满足功能正确还需要考虑效率和鲁棒性。4.1 优化 fact 函数的实现上面给出的fact函数递归版本清晰但可能有重复计算。更高效的写法是迭代式并预先计算出一个周期内的乘积。long long fact(long long n, long long p, long long pk) { if (n 0) return 1; long long res 1; // 预处理周期乘积 prod_pk: 1..pk 中与p互质的数之积 % pk // 这个计算只需要一次可以放在C_mod_pk函数里传入pk // 假设我们通过一个函数 get_prod(pk, p) 得到了这个值 // 计算完整周期 res pow_mod(prod_pk, n / pk, pk); // 计算不完整尾巴 for (long long i 1; i n % pk; i) { if (i % p ! 0) { res res * i % pk; } } // 递归处理但注意递归的是 n/p而不是 n/pk res res * fact(n / p, p, pk) % pk; return res; }这里的关键是prod_pk的计算。我们可以写一个函数long long calc_prod(long long pk, long long p) { long long res 1; for (long long i 1; i pk; i) { if (i % p ! 0) { res res * i % pk; } } return res; }然后在C_mod_pk中调用long long prod_pk calc_prod(pk, p);并传给fact函数。注意calc_prod的复杂度是O(pk)而pk是p^k。在p较小而k较大时例如p2, k30pk超过10亿这个计算是无法完成的。幸运的是在算法竞赛的数据范围内pk通常不会太大因为p是p的质因子p本身不超过1e9k也不会太大。如果遇到极端情况需要更巧妙的数学方法。4.2 处理大数运算与溢出在整个计算过程中尤其是计算fact和乘法时即使模数pk在long long范围内中间结果a * b % pk也可能溢出因为a和b都是mod pk下的数它们的乘积可能超过2^63-1。有几种解决方案使用__int128如果编译器支持这是最方便的方法。在计算a * b % pk时可以写成(long long)((__int128)a * b % pk)。使用快速乘模仿快速幂实现一个mul_mod函数通过加法来模拟乘法避免溢出。long long mul_mod(long long a, long long b, long long mod) { long long res 0; a % mod; b % mod; while (b) { if (b 1) res (res a) % mod; a (a a) % mod; b 1; } return res; }但需要注意这个O(log b)的复杂度在频繁调用时可能成为瓶颈。使用编译器内置的溢出检查或直接使用Python在竞赛中如果允许使用Python其内置的大整数可以完美规避此问题实现起来更简单。4.3 边界情况与测试m 0或m n组合数C(n, 0) C(n, n) 1。你的模板应该能正确处理。m n组合数为0。需要在入口函数检查。模数p 1根据定义任何数模1都是0。可以特判。p本身就是质数此时扩展卢卡斯定理退化为普通卢卡斯定理如果p较小或直接求逆元计算如果p较大但仍是质数。我们的模板仍然适用但效率可能不如专门的质数模数组合数算法。不过作为一个通用模板正确性优先。测试用例// 测试1: 普通情况 // C(5, 2) mod 6 10 mod 6 4 assert(exLucas(5, 2, 6) 4); // 测试2: 模数为质数 // C(10, 3) mod 7 120 mod 7 1 assert(exLucas(10, 3, 7) 1); // 测试3: 大数情况 (需要确保使用防溢出乘法) // C(100, 50) mod 1013 (质数) 应与直接计算一致 // 可以用小规模验证或对比已知结果 // 测试4: 边界 assert(exLucas(5, 0, 100) 1); assert(exLucas(5, 6, 100) 0);5. 常见问题与调试实录即使理解了原理实现时也难免踩坑。下面是我在实现和教学过程中遇到的一些典型问题。5.1 为什么结果总是0这是最常见的问题。可能的原因有p的指数计算错误在C_mod_pk函数中e_diff e_n - e_m - e_{n-m}。如果这个值大于kpk p^k中的k那么pow_mod(p, e_diff, pk)就会是0导致整个结果为0。但这是正确的因为当组合数本身包含的p因子指数足够高使得整个数能被p^k整除时模p^k的结果就是0。例如C(4, 2) 6模2^24余2模2^38余6但模2^416呢6 mod 16 6不是0。等一下这里需要仔细。 实际上pow_mod(p, e_diff, pk)只有在e_diff k时才会为0。如果组合数C(n, m)是p^k的倍数那么模p^k的结果应该是0。我们的公式p^{e_diff} mod pk确实能反映这一点。所以结果输出0不一定是错误可能是正确答案。你需要用小的样例验证。乘法溢出中间计算a * b % pk时发生了溢出导致结果错误可能偶然为0。务必使用__int128或mul_mod。逆元计算错误inv(fact_m, pk)或inv(fact_nm, pk)计算错误。确保fact_m和fact_nm是与p互质的这正是fact函数的作用并且inv函数能正确返回逆元。可以用小数据测试inv函数。5.2 递归深度过深导致栈溢出或超时fact函数是递归的递归调用是fact(n/p, ...)。当n很大而p很小时例如p2,n1e18递归深度约为log_p(n)对于p2大约是60这在可接受范围内。通常不会栈溢出。但如果你的实现有问题比如递归参数传错可能导致无限递归。5.3 CRT合并结果错误确保CRT的公式正确M是总模数p。M_i M / mods[i]。t_i是M_i在模mods[i]下的逆元不是模M下的逆元。最终累加时每次都要取模Mres (res a[i] * Mi % M * ti % M) % M;可以用简单的例子验证求解x ≡ 2 (mod 3), x ≡ 3 (mod 5)。这里M15,M15,M23。t1 inv(5,3)2,t2inv(3,5)2。解x (2*5*2 3*3*2) mod 15 (2018) mod 15 38 mod 15 8。验证8 mod 3 2,8 mod 5 3正确。5.4 性能瓶颈对于极大的n, m如1e18和p含有小质因子如2、3、5时e_func函数是O(log_p n)很快。fact函数是主要瓶颈。其时间复杂度约为O(pk log_p n)其中O(pk)来自计算prod_pk和尾巴的循环。如果pk很大比如p2, k30, pk~1e9这个循环是无法承受的。优化方向记忆化prod_pk对于相同的(p, pk)prod_pk是固定的。可以全局缓存避免重复计算。更高效的calc_prod存在基于威尔逊定理及其推广的O(pk)计算方法但常数更优。对于pk非常大的情况可能需要利用阶乘的周期性进行更复杂的数学化简但这通常超出了竞赛范围。竞赛题的数据会保证pk在一个可暴力计算的范围内例如pk 1e6左右。6. 模板的最终形态与使用指南结合以上所有讨论这里给出一个相对完整、考虑了溢出处理的扩展卢卡斯模板使用__int128。请注意为了清晰部分细节如prod_pk的缓存未完全展开。#include bits/stdc.h using namespace std; using ll long long; using i128 __int128_t; ll pow_mod(ll a, ll b, ll mod) { ll res 1; a % mod; while (b) { if (b 1) res (ll)((i128)res * a % mod); a (ll)((i128)a * a % mod); b 1; } return res; } ll exgcd(ll a, ll b, ll x, ll y) { if (!b) { x 1; y 0; return a; } ll d exgcd(b, a % b, y, x); y - a / b * x; return d; } ll inv(ll a, ll mod) { ll x, y; exgcd(a, mod, x, y); return (x % mod mod) % mod; } // 计算 n! 中剔除因子 p 后模 pk 的值 ll fact(ll n, ll p, ll pk) { if (n 0) return 1; ll res 1; // 计算周期乘积 prod_pk: 1..pk 中与p互质的数之积 % pk // 这里采用循环计算可优化为记忆化 ll prod_pk 1; for (ll i 1; i pk; i) { if (i % p) prod_pk (ll)((i128)prod_pk * i % pk); } res pow_mod(prod_pk, n / pk, pk); // 计算不完整尾巴 for (ll i 1; i n % pk; i) { if (i % p) res (ll)((i128)res * i % pk); } // 递归处理 res (ll)((i128)res * fact(n / p, p, pk) % pk); return res; } // 计算 n! 中质因子 p 的指数 ll e_func(ll n, ll p) { ll cnt 0; while (n) { cnt n / p; n / p; } return cnt; } // 计算 C(n, m) mod pk ll C_mod_pk(ll n, ll m, ll p, ll pk) { if (m n) return 0; ll f_n fact(n, p, pk); ll f_m fact(m, p, pk); ll f_nm fact(n - m, p, pk); ll e_n e_func(n, p); ll e_m e_func(m, p); ll e_nm e_func(n - m, p); ll e_diff e_n - e_m - e_nm; ll res (ll)((i128)f_n * inv(f_m, pk) % pk); res (ll)((i128)res * inv(f_nm, pk) % pk); res (ll)((i128)res * pow_mod(p, e_diff, pk) % pk); return res; } // 扩展卢卡斯主函数 ll exLucas(ll n, ll m, ll p) { if (m n) return 0; // 质因数分解 p vectorpairll, ll factors; // (质数, 质数幂) ll temp p; for (ll i 2; i * i temp; i) { if (temp % i 0) { ll pk 1; while (temp % i 0) { temp / i; pk * i; } factors.emplace_back(i, pk); } } if (temp 1) { factors.emplace_back(temp, temp); } // 分别求解 vectorll a, mods; for (auto [pi, pki] : factors) { a.push_back(C_mod_pk(n, m, pi, pki)); mods.push_back(pki); } // CRT 合并 ll res 0; ll M p; for (size_t i 0; i factors.size(); i) { ll Mi M / mods[i]; ll ti inv(Mi, mods[i]); res (res (i128)a[i] * Mi % M * ti % M) % M; } return res; } int main() { // 示例 ll n, m, p; cin n m p; cout exLucas(n, m, p) endl; return 0; }使用指南将上述代码保存为模板。注意long long的范围通常n, m可以很大但p一般在1e6以内pk也不会太大否则calc_prod的循环会超时。如果编译器不支持__int128需要实现mul_mod函数替换所有(i128)a * b % mod的运算。对于特别大的pk需要优化prod_pk的计算可能需预处理或使用更快的算法。这个模板解决了模数为合数时求组合数的通用问题是处理数论组合问题的利器。理解其每一步的数学原理比单纯套用代码更重要。当你下次遇到p不是质数的组合数问题时希望这份详细的拆解能让你从容应对。 扩展卢卡斯定理详解：模任意数组合数计算与模板实现
1. 项目概述从一道题到一个工具箱如果你在刷算法题特别是组合数学相关的题目时遇到了一个模数p不是质数甚至可能是一个合数的情况需要计算组合数C(n, m) mod p那么“扩展卢卡斯定理”就是你绕不开的一道坎。p4720这个题号在很多在线评测系统中指的就是要求实现扩展卢卡斯定理的模板题。它不像普通的卢卡斯定理那样要求模数是质数其核心思想是将非质数模数p分解质因数分别求出组合数对每个质数幂p_i^{k_i}取模的结果最后再用中国剩余定理CRT合并得到最终答案。这个“模板”的价值远不止于解决一道题。它是一套完整的、用于处理“模数为任意正整数时求组合数”的通用解决方案。掌握了它你就拥有了一个强大的数学工具可以应对竞赛、面试乃至一些特定科研场景中更广泛的模运算问题。今天我们就来彻底拆解这个模板不仅告诉你代码怎么写更要讲清楚每一个步骤背后的数学原理和工程化实现的细节让你能真正理解并灵活运用它。2. 核心思路拆解化整为零分而治之面对一个非质数模数p直接计算组合数C(n, m) n! / (m! * (n-m)!)的模运算是行不通的因为分母的逆元可能不存在当分母与模数不互质时。扩展卢卡斯定理的精妙之处在于它采用了一种“化整为零”的策略。2.1 第一步质因数分解模数假设我们需要计算C(n, m) mod p其中p是任意正整数。 首先我们对p进行质因数分解p p1^{k1} * p2^{k2} * ... * pt^{kt}其中pi是质数ki是正整数。为什么这么做因为对于每一个质数幂pi^{ki}我们有可能利用数论知识在一个“局部”的模数下处理阶乘中包含pi因子的特殊情况。而最终的结果可以通过中国剩余定理从这些“局部”解唯一地确定“全局”解mod p。注意这里质因数分解的算法效率很重要。对于p在1e9范围内的题目通常使用试除法即可如果p可能更大则需要更高效的算法如 Pollard-Rho但在竞赛模板题中p的范围通常是可控的。2.2 第二步分别求解模每个质数幂现在问题转化为对于每一个质数幂q pi^{ki}计算C(n, m) mod q。 这是整个算法的核心难点。我们无法直接计算因为n!, m!, (n-m)!可能与q不互质它们包含因子pi。解决思路剔除因子分离计算我们以计算n! mod q为例其中q p^k。分离出所有 p 的因子将n!中所有是p的倍数的数都提取出来。例如n! 1*2*3*...*n。其中p, 2p, 3p, ..., floor(n/p)*p都是p的倍数。把它们全部提取出来一共可以提取出floor(n/p)个p。提取后式子变成了n! p^{floor(n/p)} * (1*2*3*...*floor(n/p)) * [所有不能被 p 整除的数的乘积]。注意(1*2*3*...*floor(n/p))这部分又是floor(n/p)!我们可以递归地处理它处理剩余部分与 p 互质的部分剩下的数是1到n中所有不能被p整除的数。它们的乘积对q取模存在一个重要的周期性规律。考虑数列1, 2, ..., p-1, p1, p2, ..., 2p-1, 2p1, ...。以q 3^2 9为例我们看1到9中与3互质的数1,2,4,5,7,8。它们的乘积1*2*4*5*7*8 mod 9 224 mod 9 8。可以发现1到q中与p互质的数的乘积模q是一个定值记作prod(q)。对于q9prod(9) 8。那么对于1到n中所有与p互质的数我们可以将它们分成若干完整的周期每个周期长度为q和一个不完整的尾巴。完整周期的乘积就是prod(q)^{floor(n/q)} mod q。不完整尾巴即最后不足一个周期的部分需要暴力计算但长度小于q可以接受。将以上两步结合起来我们就得到了计算n! mod q的递归公式通常被称为fact(n, p, q)函数fact(n, p, q) (fact(floor(n/p), p, q) * prod(q)^{floor(n/q)} * tail(n mod q, p, q)) mod q其中tail(r, p, q)计算1到r中所有与p互质的数的乘积模q。2.3 第三步应用公式与合并结果有了计算n! mod q的能力我们就可以计算C(n, m) mod qC(n, m) mod q fact(n, p, q) * inverse(fact(m, p, q), q) * inverse(fact(n-m, p, q), q) * p^{e_n - e_m - e_{n-m}} mod q等等最后一项p^{e_n - e_m - e_{n-m}}是什么别忘了 p 的因子在第二步计算fact函数时我们刻意把所有的p因子都提取出来了。fact(n, p, q)返回的结果是n!剔除所有p因子后剩余部分模q的值。同时我们可以通过递归过程很容易地计算出n!中p因子的总指数e(n) floor(n/p) floor(n/p^2) floor(n/p^3) ...。因此真正的n!对模数q的“贡献”是fact(n, p, q) * p^{e(n)} mod q。 所以组合数C(n, m) n! / (m! * (n-m)!)对模数q的“贡献”是[fact(n, p, q) * p^{e(n)}] / [fact(m, p, q) * p^{e(m)} * fact(n-m, p, q) * p^{e(n-m)}] mod q化简后得到C(n, m) mod q fact(n, p, q) * inverse(fact(m, p, q), q) * inverse(fact(n-m, p, q), q) * p^{e(n) - e(m) - e(n-m)} mod q这里inverse(a, q)表示a在模q意义下的乘法逆元。由于fact函数返回的值是与p互质的我们剔除了所有p因子所以它必然与q互质因为q只包含质因子p逆元一定存在可以用扩展欧几里得算法求解。最后我们得到了t个同余方程C(n, m) ≡ a_i (mod p_i^{k_i}), 对于i 1, 2, ..., t。 其中a_i就是我们上面计算出的C(n, m) mod p_i^{k_i}。2.4 第四步中国剩余定理CRT合并现在我们有了一个同余方程组。中国剩余定理告诉我们如果模数两两互质这里p_i^{k_i}显然互质那么这个方程组在模p下有唯一解。 合并公式为 设M p p1^{k1} * p2^{k2} * ... * pt^{kt}。 设M_i M / p_i^{k_i}。 设t_i是M_i在模p_i^{k_i}意义下的逆元即M_i * t_i ≡ 1 (mod p_i^{k_i})。 则方程组的解为x ≡ Σ (a_i * M_i * t_i) (mod M)。这个x就是我们要的C(n, m) mod p。3. 模板实现与关键代码解析理解了原理我们来看代码实现。一个健壮的扩展卢卡斯模板需要以下几个函数3.1 快速幂与扩展欧几里得这是基础工具。// 快速幂 (a^b % mod) long long pow_mod(long long a, long long b, long long mod) { long long res 1; while (b) { if (b 1) res res * a % mod; a a * a % mod; b 1; } return res; } // 扩展欧几里得求解 ax by gcd(a, b)并返回gcd(a,b) // 同时通过引用返回x, y使得 ax by gcd(a,b) long long exgcd(long long a, long long b, long long x, long long y) { if (b 0) { x 1; y 0; return a; } long long d exgcd(b, a % b, y, x); y - a / b * x; return d; } // 求 a 在模 mod 下的逆元要求 gcd(a, mod) 1 long long inv(long long a, long long mod) { long long x, y; exgcd(a, mod, x, y); return (x % mod mod) % mod; // 调整到 0~mod-1 }3.2 计算 fact(n, p, pk) 和 指数 e(n)这个函数计算n!中剔除所有p因子后剩余部分模pk的值同时通过引用或全局变量记录p因子的总指数。// pk p^k long long fact(long long n, long long p, long long pk) { if (n 0) return 1; long long res 1; // 1. 计算完整周期的乘积 // 周期长度为 pk每个周期中与p互质部分的乘积是固定的我们预处理为 prod_pk // 假设我们已经预处理好了 prod[pk] (1*2*...*pk)中与p互质的数之积 % pk // 在实际代码中prod(pk) 可以递归或循环计算这里为了清晰假设已获得 res res * pow_mod(prod_pk, n / pk, pk) % pk; // 2. 计算最后一个不完整周期 // 暴力计算 1 到 (n % pk) 中与p互质的数的乘积 for (long long i 1; i n % pk; i) { if (i % p ! 0) { // 与p互质 res res * i % pk; } } // 3. 递归处理 floor(n/p)! // 注意递归的部分是 floor(n/p)!它可能还包含p的因子需要继续提取 res res * fact(n / p, p, pk) % pk; return res; } // 计算 n! 中质因子 p 的指数 long long e_func(long long n, long long p) { long long cnt 0; while (n) { cnt n / p; n / p; } return cnt; }在实际模板中prod_pk通常不会全局预处理因为pk随分解变化而是在fact函数内部通过循环计算出来或者用一个辅助函数calc_prod(n, p, pk)来计算1..n中与p互质的数之积模pk。为了效率常见的写法是将fact函数写成迭代形式并同时计算乘积和周期。3.3 计算 C(n, m) mod pk这是对单个质数幂模数的求解。// 计算 C(n, m) % pk其中 pk p^k long long C_mod_pk(long long n, long long m, long long p, long long pk) { if (m n) return 0; // 计算三个阶乘剔除p因子后的值 long long fact_n fact(n, p, pk); long long fact_m fact(m, p, pk); long long fact_nm fact(n - m, p, pk); // 计算p因子的指数差 long long e_n e_func(n, p); long long e_m e_func(m, p); long long e_nm e_func(n - m, p); long long e_diff e_n - e_m - e_nm; // 组合计算结果 long long res fact_n * inv(fact_m, pk) % pk; res res * inv(fact_nm, pk) % pk; res res * pow_mod(p, e_diff, pk) % pk; // 注意如果e_diff为负则结果为0因为分子中p因子不够 return res; }实操心得e_diff可能为负数吗理论上对于组合数C(n, m)其值一定是整数所以p因子的指数差e_n - e_m - e_{n-m}一定0。如果计算出现负数只可能是m n的非法输入我们已经做了判断。但在一些特殊写法或递归中确保指数计算正确很重要。3.4 质因数分解与CRT合并这是主函数exLucas(n, m, p)的逻辑。long long exLucas(long long n, long long m, long long p) { // 1. 质因数分解 p vectorpairlong long, long long factors; // 存储 (质数pi, 幂次ki) long long temp p; for (long long i 2; i * i temp; i) { if (temp % i 0) { long long cnt 0; long long pk 1; while (temp % i 0) { temp / i; cnt; pk * i; } factors.push_back({i, pk}); // 这里存储 (pi, pi^ki) } } if (temp 1) { factors.push_back({temp, temp}); } // 2. 分别求解每个质数幂下的答案 vectorlong long a, mods; for (auto fac : factors) { long long pi fac.first; long long pki fac.second; long long ai C_mod_pk(n, m, pi, pki); a.push_back(ai); mods.push_back(pki); } // 3. CRT 合并 long long res 0; long long M p; // 总模数 for (int i 0; i factors.size(); i) { long long Mi M / mods[i]; long long ti inv(Mi, mods[i]); // Mi 在模 mods[i] 下的逆元 res (res a[i] * Mi % M * ti % M) % M; } return res; }4. 优化技巧与边界处理一个工业级的模板不能只满足功能正确还需要考虑效率和鲁棒性。4.1 优化 fact 函数的实现上面给出的fact函数递归版本清晰但可能有重复计算。更高效的写法是迭代式并预先计算出一个周期内的乘积。long long fact(long long n, long long p, long long pk) { if (n 0) return 1; long long res 1; // 预处理周期乘积 prod_pk: 1..pk 中与p互质的数之积 % pk // 这个计算只需要一次可以放在C_mod_pk函数里传入pk // 假设我们通过一个函数 get_prod(pk, p) 得到了这个值 // 计算完整周期 res pow_mod(prod_pk, n / pk, pk); // 计算不完整尾巴 for (long long i 1; i n % pk; i) { if (i % p ! 0) { res res * i % pk; } } // 递归处理但注意递归的是 n/p而不是 n/pk res res * fact(n / p, p, pk) % pk; return res; }这里的关键是prod_pk的计算。我们可以写一个函数long long calc_prod(long long pk, long long p) { long long res 1; for (long long i 1; i pk; i) { if (i % p ! 0) { res res * i % pk; } } return res; }然后在C_mod_pk中调用long long prod_pk calc_prod(pk, p);并传给fact函数。注意calc_prod的复杂度是O(pk)而pk是p^k。在p较小而k较大时例如p2, k30pk超过10亿这个计算是无法完成的。幸运的是在算法竞赛的数据范围内pk通常不会太大因为p是p的质因子p本身不超过1e9k也不会太大。如果遇到极端情况需要更巧妙的数学方法。4.2 处理大数运算与溢出在整个计算过程中尤其是计算fact和乘法时即使模数pk在long long范围内中间结果a * b % pk也可能溢出因为a和b都是mod pk下的数它们的乘积可能超过2^63-1。有几种解决方案使用__int128如果编译器支持这是最方便的方法。在计算a * b % pk时可以写成(long long)((__int128)a * b % pk)。使用快速乘模仿快速幂实现一个mul_mod函数通过加法来模拟乘法避免溢出。long long mul_mod(long long a, long long b, long long mod) { long long res 0; a % mod; b % mod; while (b) { if (b 1) res (res a) % mod; a (a a) % mod; b 1; } return res; }但需要注意这个O(log b)的复杂度在频繁调用时可能成为瓶颈。使用编译器内置的溢出检查或直接使用Python在竞赛中如果允许使用Python其内置的大整数可以完美规避此问题实现起来更简单。4.3 边界情况与测试m 0或m n组合数C(n, 0) C(n, n) 1。你的模板应该能正确处理。m n组合数为0。需要在入口函数检查。模数p 1根据定义任何数模1都是0。可以特判。p本身就是质数此时扩展卢卡斯定理退化为普通卢卡斯定理如果p较小或直接求逆元计算如果p较大但仍是质数。我们的模板仍然适用但效率可能不如专门的质数模数组合数算法。不过作为一个通用模板正确性优先。测试用例// 测试1: 普通情况 // C(5, 2) mod 6 10 mod 6 4 assert(exLucas(5, 2, 6) 4); // 测试2: 模数为质数 // C(10, 3) mod 7 120 mod 7 1 assert(exLucas(10, 3, 7) 1); // 测试3: 大数情况 (需要确保使用防溢出乘法) // C(100, 50) mod 1013 (质数) 应与直接计算一致 // 可以用小规模验证或对比已知结果 // 测试4: 边界 assert(exLucas(5, 0, 100) 1); assert(exLucas(5, 6, 100) 0);5. 常见问题与调试实录即使理解了原理实现时也难免踩坑。下面是我在实现和教学过程中遇到的一些典型问题。5.1 为什么结果总是0这是最常见的问题。可能的原因有p的指数计算错误在C_mod_pk函数中e_diff e_n - e_m - e_{n-m}。如果这个值大于kpk p^k中的k那么pow_mod(p, e_diff, pk)就会是0导致整个结果为0。但这是正确的因为当组合数本身包含的p因子指数足够高使得整个数能被p^k整除时模p^k的结果就是0。例如C(4, 2) 6模2^24余2模2^38余6但模2^416呢6 mod 16 6不是0。等一下这里需要仔细。 实际上pow_mod(p, e_diff, pk)只有在e_diff k时才会为0。如果组合数C(n, m)是p^k的倍数那么模p^k的结果应该是0。我们的公式p^{e_diff} mod pk确实能反映这一点。所以结果输出0不一定是错误可能是正确答案。你需要用小的样例验证。乘法溢出中间计算a * b % pk时发生了溢出导致结果错误可能偶然为0。务必使用__int128或mul_mod。逆元计算错误inv(fact_m, pk)或inv(fact_nm, pk)计算错误。确保fact_m和fact_nm是与p互质的这正是fact函数的作用并且inv函数能正确返回逆元。可以用小数据测试inv函数。5.2 递归深度过深导致栈溢出或超时fact函数是递归的递归调用是fact(n/p, ...)。当n很大而p很小时例如p2,n1e18递归深度约为log_p(n)对于p2大约是60这在可接受范围内。通常不会栈溢出。但如果你的实现有问题比如递归参数传错可能导致无限递归。5.3 CRT合并结果错误确保CRT的公式正确M是总模数p。M_i M / mods[i]。t_i是M_i在模mods[i]下的逆元不是模M下的逆元。最终累加时每次都要取模Mres (res a[i] * Mi % M * ti % M) % M;可以用简单的例子验证求解x ≡ 2 (mod 3), x ≡ 3 (mod 5)。这里M15,M15,M23。t1 inv(5,3)2,t2inv(3,5)2。解x (2*5*2 3*3*2) mod 15 (2018) mod 15 38 mod 15 8。验证8 mod 3 2,8 mod 5 3正确。5.4 性能瓶颈对于极大的n, m如1e18和p含有小质因子如2、3、5时e_func函数是O(log_p n)很快。fact函数是主要瓶颈。其时间复杂度约为O(pk log_p n)其中O(pk)来自计算prod_pk和尾巴的循环。如果pk很大比如p2, k30, pk~1e9这个循环是无法承受的。优化方向记忆化prod_pk对于相同的(p, pk)prod_pk是固定的。可以全局缓存避免重复计算。更高效的calc_prod存在基于威尔逊定理及其推广的O(pk)计算方法但常数更优。对于pk非常大的情况可能需要利用阶乘的周期性进行更复杂的数学化简但这通常超出了竞赛范围。竞赛题的数据会保证pk在一个可暴力计算的范围内例如pk 1e6左右。6. 模板的最终形态与使用指南结合以上所有讨论这里给出一个相对完整、考虑了溢出处理的扩展卢卡斯模板使用__int128。请注意为了清晰部分细节如prod_pk的缓存未完全展开。#include bits/stdc.h using namespace std; using ll long long; using i128 __int128_t; ll pow_mod(ll a, ll b, ll mod) { ll res 1; a % mod; while (b) { if (b 1) res (ll)((i128)res * a % mod); a (ll)((i128)a * a % mod); b 1; } return res; } ll exgcd(ll a, ll b, ll x, ll y) { if (!b) { x 1; y 0; return a; } ll d exgcd(b, a % b, y, x); y - a / b * x; return d; } ll inv(ll a, ll mod) { ll x, y; exgcd(a, mod, x, y); return (x % mod mod) % mod; } // 计算 n! 中剔除因子 p 后模 pk 的值 ll fact(ll n, ll p, ll pk) { if (n 0) return 1; ll res 1; // 计算周期乘积 prod_pk: 1..pk 中与p互质的数之积 % pk // 这里采用循环计算可优化为记忆化 ll prod_pk 1; for (ll i 1; i pk; i) { if (i % p) prod_pk (ll)((i128)prod_pk * i % pk); } res pow_mod(prod_pk, n / pk, pk); // 计算不完整尾巴 for (ll i 1; i n % pk; i) { if (i % p) res (ll)((i128)res * i % pk); } // 递归处理 res (ll)((i128)res * fact(n / p, p, pk) % pk); return res; } // 计算 n! 中质因子 p 的指数 ll e_func(ll n, ll p) { ll cnt 0; while (n) { cnt n / p; n / p; } return cnt; } // 计算 C(n, m) mod pk ll C_mod_pk(ll n, ll m, ll p, ll pk) { if (m n) return 0; ll f_n fact(n, p, pk); ll f_m fact(m, p, pk); ll f_nm fact(n - m, p, pk); ll e_n e_func(n, p); ll e_m e_func(m, p); ll e_nm e_func(n - m, p); ll e_diff e_n - e_m - e_nm; ll res (ll)((i128)f_n * inv(f_m, pk) % pk); res (ll)((i128)res * inv(f_nm, pk) % pk); res (ll)((i128)res * pow_mod(p, e_diff, pk) % pk); return res; } // 扩展卢卡斯主函数 ll exLucas(ll n, ll m, ll p) { if (m n) return 0; // 质因数分解 p vectorpairll, ll factors; // (质数, 质数幂) ll temp p; for (ll i 2; i * i temp; i) { if (temp % i 0) { ll pk 1; while (temp % i 0) { temp / i; pk * i; } factors.emplace_back(i, pk); } } if (temp 1) { factors.emplace_back(temp, temp); } // 分别求解 vectorll a, mods; for (auto [pi, pki] : factors) { a.push_back(C_mod_pk(n, m, pi, pki)); mods.push_back(pki); } // CRT 合并 ll res 0; ll M p; for (size_t i 0; i factors.size(); i) { ll Mi M / mods[i]; ll ti inv(Mi, mods[i]); res (res (i128)a[i] * Mi % M * ti % M) % M; } return res; } int main() { // 示例 ll n, m, p; cin n m p; cout exLucas(n, m, p) endl; return 0; }使用指南将上述代码保存为模板。注意long long的范围通常n, m可以很大但p一般在1e6以内pk也不会太大否则calc_prod的循环会超时。如果编译器不支持__int128需要实现mul_mod函数替换所有(i128)a * b % mod的运算。对于特别大的pk需要优化prod_pk的计算可能需预处理或使用更快的算法。这个模板解决了模数为合数时求组合数的通用问题是处理数论组合问题的利器。理解其每一步的数学原理比单纯套用代码更重要。当你下次遇到p不是质数的组合数问题时希望这份详细的拆解能让你从容应对。