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