1. 項(xiàng)目緣起從一道“卡脖子”的模數(shù)組合數(shù)說起幾年前我在準(zhǔn)備一場(chǎng)算法競(jìng)賽時(shí)遇到了一道讓我記憶猶新的題目。題目本身描述很簡(jiǎn)單給定一個(gè)巨大的組合數(shù) C(n, m)以及一個(gè)模數(shù) P要求計(jì)算 C(n, m) mod P 的值。這聽起來像是數(shù)論基礎(chǔ)題我信心滿滿地寫下了標(biāo)準(zhǔn)的預(yù)處理階乘和逆元的代碼。然而當(dāng)我提交時(shí)系統(tǒng)返回了一個(gè)大大的“Wrong Answer”。仔細(xì)一看模數(shù) P 的描述心里頓時(shí)涼了半截——這個(gè) P 不是一個(gè)質(zhì)數(shù)甚至不是一個(gè)質(zhì)數(shù)的冪而是一個(gè)任意的正整數(shù)比如 10007 或者 999983 這種質(zhì)數(shù)還好但題目給的可能是 1000、2016 甚至 1000000007 * 1000000009 這種合數(shù)。這就是經(jīng)典的“組合數(shù)取模”問題在模數(shù)非質(zhì)數(shù)時(shí)遇到的困境。我們熟悉的用費(fèi)馬小定理或擴(kuò)展歐幾里得求逆元的方法其前提是模數(shù)為質(zhì)數(shù)這樣才能保證在模意義下每個(gè)非零數(shù)都有乘法逆元。當(dāng)模數(shù)是合數(shù)時(shí)許多數(shù)沒有逆元整個(gè)基于階乘和逆元的遞推公式就失效了。我當(dāng)時(shí)卡在這道題上很久直到后來系統(tǒng)學(xué)習(xí)了“擴(kuò)展盧卡斯定理”Extended Lucas Theorem才豁然開朗。而 P4720正是洛谷上一道專門練習(xí)這個(gè)定理的經(jīng)典模板題。今天我就結(jié)合自己踩坑和實(shí)戰(zhàn)的經(jīng)驗(yàn)把這個(gè)強(qiáng)大工具的原理、實(shí)現(xiàn)細(xì)節(jié)和避坑指南掰開揉碎了講清楚。2. 核心問題拆解為什么普通盧卡斯定理不夠用在深入擴(kuò)展盧卡斯之前我們必須先理解普通盧卡斯定理Lucas Theorem的局限性這樣才能明白我們到底要解決什么問題。2.1 盧卡斯定理的適用場(chǎng)景與限制普通盧卡斯定理表述為對(duì)于質(zhì)數(shù) p有 C(n, m) ≡ C(n mod p, m mod p) * C(n/p, m/p) (mod p)。這是一個(gè)遞歸公式能將大規(guī)模組合數(shù)計(jì)算分解為小規(guī)模組合數(shù)計(jì)算通常配合預(yù)處理小范圍的階乘和逆元來快速求解。它的效率很高代碼也很簡(jiǎn)潔。然而它的核心限制就藏在前提里模數(shù) p 必須是質(zhì)數(shù)。這是因?yàn)樵谶f歸的底層我們需要計(jì)算 C(n‘, m’) mod p這通常通過公式 C n! / (m! * (n-m)!) 來計(jì)算而除法在模運(yùn)算中需要轉(zhuǎn)化為乘以其乘法逆元。逆元存在的充要條件就是該數(shù)與模數(shù)互質(zhì)。當(dāng)模數(shù)是質(zhì)數(shù) p 時(shí)只要分母的階乘不被 p 整除其逆元就一定存在。但一旦模數(shù) p 是合數(shù)分母的階乘很可能與模數(shù)有公因子導(dǎo)致逆元不存在整個(gè)計(jì)算鏈就斷裂了。2.2 合數(shù)模數(shù)帶來的真正挑戰(zhàn)當(dāng)模數(shù) P 是合數(shù)時(shí)直接計(jì)算 C(n, m) mod P 的難點(diǎn)可以歸結(jié)為兩點(diǎn)非互質(zhì)導(dǎo)致的逆元缺失在計(jì)算 n! / (m! * (n-m)!) 時(shí)分母可能與模數(shù) P 不互質(zhì)因此無法直接求逆元進(jìn)行模除。模數(shù)非質(zhì)數(shù)中國剩余定理CRT成為橋梁解決這個(gè)問題的核心思路是將合數(shù)模數(shù) P 質(zhì)因數(shù)分解為 P p1^k1 * p2^k2 * ... * pt^kt。如果我們能分別求出 C(n, m) 對(duì)每個(gè)質(zhì)數(shù)冪模數(shù) pi^ki 的余數(shù) ai即求解一系列同余方程x ≡ a1 (mod p1^k1) x ≡ a2 (mod p2^k2) ... x ≡ at (mod pt^kt) 那么根據(jù)中國剩余定理我們就可以唯一確定出 x 在模 P 意義下的值。所以問題的關(guān)鍵轉(zhuǎn)化為如何計(jì)算 C(n, m) mod p^k其中 p 是質(zhì)數(shù)k是正整數(shù)。這就是擴(kuò)展盧卡斯定理要解決的核心子問題。普通盧卡斯定理處理的是 mod pk1而擴(kuò)展盧卡斯將其推廣到了 mod p^k。3. 擴(kuò)展盧卡斯定理的核心原理剝離p因子與遞歸求解計(jì)算 C(n, m) mod p^k 不能直接用階乘逆元因?yàn)榉帜缚赡馨蜃?p導(dǎo)致與模數(shù) p^k 不互質(zhì)。擴(kuò)展盧卡斯定理的精妙之處在于它通過一種“剝離”技巧將階乘中所有 p 的因子分離出來單獨(dú)處理。3.1 第一步將階乘分解為“與p互質(zhì)部分”和“p的冪次部分”定義函數(shù)F(n, p, pk)用于計(jì)算 n! 中所有與 p 互質(zhì)的因子的乘積再對(duì) pk (即 p^k) 取模。同時(shí)我們記錄下 n! 中 p 這個(gè)質(zhì)因子的總次數(shù)記為G(n, p)。以 n22, p3 為例計(jì)算 22! mod 3^2 22! 1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20 * 21 * 22 我們可以把它重寫為 22! (124578101113141617192022) * (36912151821) 進(jìn)一步把第二組每個(gè)數(shù)中的因子3提出來 22! (124578101113141617192022) * 3^7 * (1234567) 你會(huì)發(fā)現(xiàn)(124578101113141617192022) 這些數(shù)都與3互質(zhì)而 (1234567) 正好是 floor(22/3) 7 的階乘即 7!。于是我們得到一個(gè)遞歸定義 n! ≡ F(n, p, pk) * p^{G(n, p)} * (n/p)! (mod pk) 其中F(n, p, pk)計(jì)算了1到n中所有不被p整除的數(shù)的乘積模 pk。G(n, p) floor(n/p) floor(n/p^2) floor(n/p^3) ...即n!中質(zhì)因子p的個(gè)數(shù)。(n/p)!是遞歸部分。對(duì)于F(n, p, pk)的計(jì)算也有技巧。因?yàn)槟?shù)是 pk而1到pk中與p互質(zhì)的數(shù)會(huì)形成一個(gè)長(zhǎng)度為 φ(pk)pk-p^{k-1} 的循環(huán)節(jié)。我們可以先計(jì)算一個(gè)完整循環(huán)節(jié)內(nèi)所有與p互質(zhì)的數(shù)的乘積模 pk記為prod。那么n 以內(nèi)這樣的完整循環(huán)節(jié)有n / pk個(gè)每個(gè)循環(huán)節(jié)的貢獻(xiàn)是prod^{n/pk} mod pk。最后再乘以剩余的不完整部分即從floor(n/pk)*pk 1到 n 之間且與p互質(zhì)的數(shù)的乘積。3.2 第二步計(jì)算組合數(shù)模 p^k有了上面的分解組合數(shù)可以表示為 C(n, m) n! / (m! * (n-m)!) 將其用 F 和 G 函數(shù)表示 C(n, m) [F(n) * p^{G(n)}] / [F(m) * p^{G(m)} * F(n-m) * p^{G(n-m)}] [F(n) / (F(m) * F(n-m))] * p^{G(n) - G(m) - G(n-m)}我們的目標(biāo)是求 C(n, m) mod p^k。首先計(jì)算指數(shù)部分e G(n) - G(m) - G(n-m)。如果 e k說明組合數(shù)本身包含了至少 p^k 這個(gè)因子那么 C(n, m) mod p^k 0。如果 e k則繼續(xù)。計(jì)算互質(zhì)部分num F(n, p, pk) * inv(F(m, p, pk), pk) * inv(F(n-m, p, pk), pk) mod pk。這里inv(a, pk)表示 a 在模 pk 意義下的逆元。由于 F 函數(shù)計(jì)算的結(jié)果都是與 p 互質(zhì)的所以它們對(duì)模數(shù) pk 的逆元一定存在可以用擴(kuò)展歐幾里得算法求解。最終結(jié)果C(n, m) mod p^k num * p^e mod pk。3.3 第三步中國剩余定理CRT合成最終答案假設(shè)我們對(duì)合數(shù)模數(shù) P 分解后得到了 t 個(gè)方程 x ≡ ans_i (mod pi^ki), i1, 2, ..., t。 其中 ans_i 就是我們用上述方法計(jì)算出來的 C(n, m) mod pi^ki。中國剩余定理的求解過程如下計(jì)算M P。對(duì)于每個(gè) i計(jì)算Mi M / (pi^ki)。計(jì)算Mi在模pi^ki意義下的逆元inv_i因?yàn)?Mi 與 pi^ki 互質(zhì)逆元存在。最終解為x Σ(ans_i * Mi * inv_i) mod M。這一步在算法實(shí)現(xiàn)中通常使用“增量法”合并同余方程每次合并兩個(gè)方程逐步得到最終解比一次性計(jì)算所有逆元更易于編碼。4. 手把手實(shí)現(xiàn)擴(kuò)展盧卡斯模板理解了原理我們來看代碼實(shí)現(xiàn)。我將結(jié)合 P4720 這道模板題的要求給出一個(gè)清晰、健壯且包含詳細(xì)注釋的 C 實(shí)現(xiàn)。代碼會(huì)分為幾個(gè)核心函數(shù)。4.1 基礎(chǔ)工具函數(shù)快速冪與擴(kuò)展歐幾里得這些是數(shù)論算法的基石。// 快速冪計(jì)算 (base^exp) % mod long long qpow(long long base, long long exp, long long mod) { long long res 1 % mod; // 注意 mod1 的情況 base % mod; while (exp) { if (exp 1) res (res * base) % mod; base (base * base) % mod; exp 1; } return res; } // 擴(kuò)展歐幾里得算法求解 ax by gcd(a, b) // 返回 gcd(a, b)并通過引用返回 x, y 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); // 將逆元調(diào)整到 [0, mod) 范圍內(nèi) return (x % mod mod) % mod; }4.2 核心函數(shù) F計(jì)算剔除了p因子的階乘模 p^k這個(gè)函數(shù)對(duì)應(yīng)原理部分的F(n, p, pk)。/** * 計(jì)算 n! 中所有與質(zhì)數(shù) p 互質(zhì)的因子的乘積再對(duì) pk (p^k) 取模。 * param n 階乘的上限 * param p 質(zhì)數(shù) * param pk p^k即當(dāng)前處理的質(zhì)數(shù)冪模數(shù) * return n! 中與 p 互質(zhì)部分的乘積模 pk */ long long factorial_prime(long long n, long long p, long long pk) { if (n 0) return 1; long long res 1; // 1. 處理完整循環(huán)節(jié)周期為 pk每個(gè)周期內(nèi)與p互質(zhì)的數(shù)乘積相同 // 計(jì)算一個(gè)周期內(nèi)的乘積 long long cycle_prod 1; for (long long i 1; i pk; i) { if (i % p ! 0) { // 與p互質(zhì) cycle_prod (cycle_prod * i) % pk; } } // 共有 n/pk 個(gè)完整周期 res qpow(cycle_prod, n / pk, pk); // 2. 處理最后一個(gè)不完整的周期 for (long long i (n / pk) * pk 1; i n; i) { if (i % p ! 0) { res (res * (i % pk)) % pk; // i % pk 防止溢出且結(jié)果等價(jià) } } // 3. 遞歸處理 (n/p)! 中與p互質(zhì)的部分 // 因?yàn)?n! (1*2*...*n) (所有與p互質(zhì)的數(shù)) * p * (所有與p互質(zhì)的數(shù)) * ... // 遞歸部分正是 (n/p)! 中與p互質(zhì)的部分 return res * factorial_prime(n / p, p, pk) % pk; }注意這里有一個(gè)非常關(guān)鍵的優(yōu)化和易錯(cuò)點(diǎn)。在計(jì)算cycle_prod時(shí)我們是在模pk下計(jì)算1到pk之間與p互質(zhì)的數(shù)的乘積。pk可能很大比如 5^10直接循環(huán)pk次在n很大時(shí)會(huì)被多次調(diào)用可能成為性能瓶頸。在實(shí)際的高性能模板中通常會(huì)預(yù)處理這個(gè)值。但為了代碼清晰這里展示了最基本的邏輯。在真正解題時(shí)需要根據(jù)數(shù)據(jù)范圍權(quán)衡是否預(yù)處理。4.3 核心函數(shù) G計(jì)算 n! 中質(zhì)因子 p 的個(gè)數(shù)這個(gè)函數(shù)用于計(jì)算指數(shù)e。/** * 計(jì)算 n! 中質(zhì)因子 p 的個(gè)數(shù)。 * 公式G(n, p) floor(n/p) floor(n/p^2) floor(n/p^3) ... * param n 階乘的上限 * param p 質(zhì)數(shù) * return n! 中質(zhì)因子 p 的個(gè)數(shù) */ long long count_prime_factor(long long n, long long p) { long long cnt 0; while (n) { cnt n / p; n / p; } return cnt; }4.4 核心函數(shù) C_mod_pk計(jì)算組合數(shù)模質(zhì)數(shù)冪這個(gè)函數(shù)整合前兩步計(jì)算 C(n, m) mod p^k。/** * 計(jì)算組合數(shù) C(n, m) 對(duì)質(zhì)數(shù)冪 p^k 取模的結(jié)果。 * param n 組合數(shù)上標(biāo) * param m 組合數(shù)下標(biāo) * param p 質(zhì)數(shù) * param pk p^k * return C(n, m) mod pk */ long long combination_mod_prime_power(long long n, long long m, long long p, long long pk) { if (m n) return 0; if (m 0 || m n) return 1 % pk; // 1. 計(jì)算指數(shù) e G(n) - G(m) - G(n-m) long long e count_prime_factor(n, p) - count_prime_factor(m, p) - count_prime_factor(n - m, p); if (e 0) { // 實(shí)際上e永遠(yuǎn)0這里判斷是為了邏輯清晰也可以判斷 if (e k) return 0; // 如果 e k (即pk中p的冪次)那么結(jié)果模pk為0。這里k可以通過pk和p計(jì)算得到。 // 簡(jiǎn)便寫法如果 e 足夠大使得 p^e 是 pk 的倍數(shù)則直接返回0。 // 更嚴(yán)謹(jǐn)?shù)淖龇ㄊ怯?jì)算 k log(pk) / log(p) 的整數(shù)部分然后比較 e 和 k。 // 下面我們采用另一種方式先計(jì)算互質(zhì)部分如果 e 很大最后乘 p^e 時(shí)再取模。 } else { // 理論上不會(huì)出現(xiàn)負(fù)數(shù)因?yàn)榻M合數(shù)是整數(shù) return 0; } // 2. 計(jì)算互質(zhì)部分F(n) / (F(m) * F(n-m)) long long fn factorial_prime(n, p, pk); long long fm factorial_prime(m, p, pk); long long fnm factorial_prime(n - m, p, pk); long long num fn * inv(fm, pk) % pk * inv(fnm, pk) % pk; // 3. 乘以 p^e long long pe qpow(p, e, pk); // 注意這里是對(duì) pk 取模因?yàn)樽罱K結(jié)果是模 pk // 但是當(dāng) e k 時(shí)p^e mod pk 為 0所以這一步包含了 ek 時(shí)結(jié)果為0的情況。 return num * pe % pk; }踩坑點(diǎn)在計(jì)算num時(shí)一定要先對(duì)fm和fnm分別求逆元然后連乘取模。不能先計(jì)算fm * fnm % pk再求一次逆元因?yàn)槌朔赡芷茐幕ベ|(zhì)性導(dǎo)致逆元不存在。必須保證每個(gè)與pk互質(zhì)的數(shù)單獨(dú)求逆。4.5 核心函數(shù) exLucas主函數(shù)與CRT合并這是對(duì)外的接口處理合數(shù)模數(shù) P。/** * 擴(kuò)展盧卡斯定理主函數(shù)計(jì)算 C(n, m) mod PP 可為任意正整數(shù)。 * param n 組合數(shù)上標(biāo) * param m 組合數(shù)下標(biāo) * param P 模數(shù)任意正整數(shù) * return C(n, m) mod P */ long long exLucas(long long n, long long m, long long P) { if (m n) return 0; if (m 0 || m n) return 1 % P; long long mod P; // 存儲(chǔ) (質(zhì)數(shù), 質(zhì)數(shù)冪, 余數(shù)) 三元組 vectortuplelong long, long long, long long factors; // 1. 對(duì)模數(shù) P 進(jìn)行質(zhì)因數(shù)分解 for (long long i 2; i * i mod; i) { if (mod % i 0) { long long pk 1; while (mod % i 0) { mod / i; pk * i; } // 計(jì)算 C(n, m) mod i^pk long long res combination_mod_prime_power(n, m, i, pk); factors.emplace_back(i, pk, res); } } if (mod 1) { // 處理剩余的大質(zhì)數(shù) factors.emplace_back(mod, mod, combination_mod_prime_power(n, m, mod, mod)); } // 2. 如果只有一個(gè)質(zhì)因數(shù)直接返回結(jié)果 if (factors.size() 1) { return get2(factors[0]); } // 3. 使用中國剩余定理CRT合并所有同余方程 // 增量法合并x ≡ a1 (mod m1), x ≡ a2 (mod m2) // 合并為 x ≡ new_a (mod new_m)其中 new_m m1 * m2 long long a1 get2(factors[0]), m1 get1(factors[0]); for (size_t i 1; i factors.size(); i) { long long a2 get2(factors[i]), m2 get1(factors[i]); // 合并方程x a1 k1*m1 a2 k2*m2 // 即 k1*m1 - k2*m2 a2 - a1 // 令 g gcd(m1, m2)用擴(kuò)展歐幾里得求解 long long k1, k2; long long g exgcd(m1, m2, k1, k2); long long c a2 - a1; if (c % g ! 0) { // 理論上不會(huì)發(fā)生因?yàn)楦髂?shù)兩兩互質(zhì) return -1; // 無解 } long long t m2 / g; // 調(diào)整 k1 為最小非負(fù)特解 k1 (k1 * (c / g) % t t) % t; // 新的余數(shù)和模數(shù) long long new_a (a1 k1 * m1) % (m1 / g * m2); // 注意新模數(shù)是 lcm(m1, m2) m1/g*m2 long long new_m m1 / g * m2; a1 new_a; m1 new_m; } return (a1 % P P) % P; // 確保結(jié)果在 [0, P) 范圍內(nèi) }5. 實(shí)戰(zhàn)測(cè)試與性能優(yōu)化要點(diǎn)將上述代碼整合就可以通過 P4720 這道模板題了。輸入 n, m, P調(diào)用exLucas(n, m, P)即可。但是直接使用上面的代碼可能會(huì)在數(shù)據(jù)較大時(shí)超時(shí)我們需要關(guān)注幾個(gè)性能瓶頸和優(yōu)化點(diǎn)。5.1 性能瓶頸分析factorial_prime函數(shù)中的循環(huán)計(jì)算cycle_prod每次遞歸調(diào)用都會(huì)計(jì)算一次從1到pk的循環(huán)積。如果pk很大比如 10^6 級(jí)別且n也很大導(dǎo)致遞歸深度不淺這個(gè)開銷是巨大的。遞歸調(diào)用factorial_prime遞歸本身有一定開銷但更主要的是重復(fù)計(jì)算。factorial_prime(n/p, p, pk)會(huì)再次計(jì)算cycle_prod。質(zhì)因數(shù)分解對(duì) P 的分解是 O(√P) 的在 P 很大如 10^9時(shí)可以接受但也是常數(shù)開銷。5.2 關(guān)鍵優(yōu)化策略優(yōu)化1預(yù)處理循環(huán)節(jié)乘積這是最重要的優(yōu)化。對(duì)于給定的p和pkcycle_prod是一個(gè)定值。我們可以在計(jì)算combination_mod_prime_power之前先計(jì)算并存儲(chǔ)它避免在遞歸中重復(fù)計(jì)算。// 在 combination_mod_prime_power 函數(shù)內(nèi)部或外部預(yù)處理 long long cycle_prod 1; for (long long i 1; i pk; i) { if (i % p ! 0) { cycle_prod cycle_prod * i % pk; } } // 然后將 cycle_prod 作為參數(shù)傳遞給 factorial_prime或者設(shè)為全局/靜態(tài)變量。 // 修改 factorial_prime 函數(shù)接收這個(gè)預(yù)計(jì)算好的 cycle_prod。 long long factorial_prime(long long n, long long p, long long pk, long long cycle_prod) { if (n 0) return 1; long long res qpow(cycle_prod, n / pk, pk); // ... 剩余部分不變 }優(yōu)化2將遞歸改為迭代factorial_prime的遞歸形式清晰但可以改為迭代形式效率略高且避免了遞歸棧溢出的風(fēng)險(xiǎn)雖然此題一般不會(huì)。long long factorial_prime_iter(long long n, long long p, long long pk, long long cycle_prod) { long long res 1; while (n 0) { res res * qpow(cycle_prod, n / pk, pk) % pk; for (long long i (n / pk) * pk 1; i n; i) { if (i % p ! 0) { res res * (i % pk) % pk; } } n / p; // 關(guān)鍵對(duì)應(yīng)遞歸中的 n/p } return res; }這個(gè)迭代版本模擬了遞歸過程每次循環(huán)處理當(dāng)前n的互質(zhì)部分和剩余部分然后將n更新為n/p直到n為 0。優(yōu)化3使用更快的質(zhì)因數(shù)分解對(duì)于巨大的 P可以使用 Pollard-Rho 算法進(jìn)行質(zhì)因數(shù)分解但這超出了模板題的一般范圍。P4720 的數(shù)據(jù)范圍下試除法足夠。5.3 一個(gè)優(yōu)化后的整合示例核心部分結(jié)合優(yōu)化combination_mod_prime_power函數(shù)可以這樣寫long long combination_mod_prime_power(long long n, long long m, long long p, long long pk) { if (m n) return 0; // 預(yù)處理循環(huán)節(jié)乘積 long long cycle_prod 1; for (long long i 1; i pk; i) { if (i % p ! 0) cycle_prod cycle_prod * i % pk; } auto factorial [](long long x) - long long { long long res 1; long long tx x; while (tx 0) { res res * qpow(cycle_prod, tx / pk, pk) % pk; for (long long i (tx / pk) * pk 1; i tx; i) { if (i % p ! 0) res res * (i % pk) % pk; } tx / p; } return res; }; long long e count_prime_factor(n, p) - count_prime_factor(m, p) - count_prime_factor(n - m, p); // 如果 e 已經(jīng)大于等于 k (即 pk 中 p 的冪次)可以提前返回0。 // 計(jì)算 k: 通過不斷除以 p 得到 long long temp_pk pk, k 0; while (temp_pk % p 0) { temp_pk / p; k; } if (e k) return 0; long long fn factorial(n); long long fm factorial(m); long long fnm factorial(n - m); long long num fn * inv(fm, pk) % pk * inv(fnm, pk) % pk; long long pe qpow(p, e, pk); return num * pe % pk; }6. 邊界條件與常見錯(cuò)誤排查即使理解了原理和代碼在實(shí)際編碼和調(diào)試中依然會(huì)遇到一些隱蔽的坑。6.1 數(shù)據(jù)范圍與溢出處理這是數(shù)論題最經(jīng)典的坑。題目中 n, m 可能高達(dá) 10^18P 在 10^6 以內(nèi)。qpow中的乘法溢出res * base或base * base可能超過long long范圍約 9e18。即使對(duì)mod取模在乘法運(yùn)算前就可能溢出了。必須使用快速乘龜速乘或__int128。// 使用 __int128 的快速冪推薦前提是編譯器支持 long long qpow(long long base, long long exp, long long mod) { __int128 res 1 % mod; __int128 b base % mod; while (exp) { if (exp 1) res (res * b) % mod; b (b * b) % mod; exp 1; } return (long long)res; } // 或者在無法使用 __int128 時(shí)使用快速乘 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; }遞歸/迭代中的變量范圍在factorial_prime的循環(huán)for (long long i (n / pk) * pk 1; i n; i)中(n / pk) * pk的計(jì)算可能溢出。更安全的寫法是long long start (n / pk) * pk; if (start 0) start 0;或者直接利用取模性質(zhì)。6.2 特殊模數(shù)處理模數(shù) P 1根據(jù)定義任何數(shù)模 1 都為 0。在代碼開頭應(yīng)特判。質(zhì)數(shù)冪 pk 可能為 1在質(zhì)因數(shù)分解時(shí)如果 p^k 1這沒有意義。實(shí)際上當(dāng) P 分解時(shí)pk 至少為 p。但若 P1已特判。CRT 合并中的模數(shù)在增量法合并同余方程時(shí)新的模數(shù)是m1 / g * m2即lcm(m1, m2)。務(wù)必注意計(jì)算順序先除后乘避免中間結(jié)果溢出?？梢允褂胈_int128輔助計(jì)算。6.3 調(diào)試技巧當(dāng)結(jié)果錯(cuò)誤時(shí)可以按以下步驟隔離問題測(cè)試小數(shù)據(jù)用小的 n, m 和小的合數(shù) P如 P6, 10進(jìn)行測(cè)試與暴力計(jì)算的結(jié)果對(duì)比。分離測(cè)試子函數(shù)測(cè)試count_prime_factor驗(yàn)證 n! 中因子 p 的個(gè)數(shù)計(jì)算是否正確。測(cè)試factorial_prime選擇小的 n, p, pk手動(dòng)計(jì)算驗(yàn)證。測(cè)試combination_mod_prime_power針對(duì)單個(gè)質(zhì)數(shù)冪模數(shù)測(cè)試。最后測(cè)試exLucas和 CRT 合并。驗(yàn)證質(zhì)因數(shù)分解確保對(duì) P 的分解是正確的。檢查逆元計(jì)算確保在求逆元時(shí)inv函數(shù)傳入的參數(shù)與模數(shù)互質(zhì)。在combination_mod_prime_power中求逆元前可以加斷言assert(gcd(fm, pk)1 gcd(fnm, pk)1)調(diào)試時(shí)。7. 擴(kuò)展盧卡斯的其他應(yīng)用與變體掌握了這個(gè)模板你不僅能解決 P4720還能處理一系列衍生問題。7.1 計(jì)算大組合數(shù)模任意數(shù)這是最直接的應(yīng)用。在一些計(jì)數(shù)問題中模數(shù)可能不是質(zhì)數(shù)比如998244353 * 1000000007這種擴(kuò)展盧卡斯是唯一通用的方法。7.2 處理模數(shù)較小但n,m巨大的情況即使模數(shù) P 是質(zhì)數(shù)如果 n 和 m 巨大遠(yuǎn)超 P盧卡斯定理可以將問題規(guī)?？s小到 P 以內(nèi)。而如果 P 不是質(zhì)數(shù)擴(kuò)展盧卡斯是唯一選擇。它通過遞歸將 n! 分解巧妙地處理了 n 很大的情況。7.3 與多項(xiàng)式、生成函數(shù)結(jié)合在一些更復(fù)雜的組合恒等式證明或求和問題中需要處理模任意數(shù)的二項(xiàng)式系數(shù)。擴(kuò)展盧卡斯提供的C(n, m) mod P的能力可以作為子程序嵌入到更大的算法框架中。7.4 局限性擴(kuò)展盧卡斯定理的時(shí)間復(fù)雜度主要取決于模數(shù) P 的質(zhì)因數(shù)分解和每個(gè)質(zhì)數(shù)冪的大小。設(shè) P 分解為 ∏ pi^ki則時(shí)間復(fù)雜度約為 O(∑ (ki * pi log n))。當(dāng) P 包含大的質(zhì)數(shù)冪如 2^30時(shí)計(jì)算cycle_prod的循環(huán)會(huì)非常慢。在這種情況下算法可能不再適用需要尋找其他數(shù)學(xué)方法或題目給定的特殊約束。在我自己的使用經(jīng)驗(yàn)里擴(kuò)展盧卡斯是一個(gè)“知道原理就能寫但想寫對(duì)、寫快需要很多細(xì)節(jié)打磨”的算法。它不像快速冪或歐拉篩那樣有幾乎固定的短代碼它的實(shí)現(xiàn)長(zhǎng)度和細(xì)節(jié)處理恰恰體現(xiàn)了數(shù)論算法從理論到實(shí)踐的復(fù)雜性。理解F(n, p, pk)那個(gè)遞歸式是第一步而處理好循環(huán)節(jié)、溢出、CRT合并以及各種邊界條件才是能在比賽中穩(wěn)定拿分的關(guān)鍵。建議在理解的基礎(chǔ)上親手實(shí)現(xiàn)并通過 P4720 這道模板題進(jìn)行測(cè)試過程中遇到的每一個(gè)錯(cuò)誤都會(huì)讓你對(duì)模運(yùn)算和遞歸分解有更深的認(rèn)識(shí)。