取模與歐拉降冪的經(jīng)典應(yīng)用)
1. 項目概述一道融合數(shù)論精華的經(jīng)典競賽題看到這個標(biāo)題老競賽黨們DNA估計要動了?!肮糯i文”這個聽起來有點無厘頭的名字背后其實是信息學(xué)奧賽OI和算法競賽圈子里一道相當(dāng)經(jīng)典的題目來自洛谷Luogu題庫的P2480原題出自SDOI2010。這道題之所以讓人印象深刻甚至有點“談虎色變”是因為它完美地融合了數(shù)論中幾個核心且有一定難度的知識點組合數(shù)取模、Lucas定理以及最終的中國剩余定理CRT與歐拉定理的應(yīng)用。它不像單純的動態(tài)規(guī)劃或者圖論題有比較固定的套路可循而是需要你真正理解這些數(shù)論工具的原理并像搭積木一樣將它們靈活地組合起來解決一個復(fù)雜問題。簡單來說題目會給你兩個巨大的整數(shù)N和G要求你計算 G^{sum} mod 999911659 的值。這里的 sum 是什么是所有滿足 d 能整除 N 的 C(N, d) 之和。也就是說你需要計算 N 的所有正因數(shù)的組合數(shù)之和。N 可以非常大最大到 10^9G 也可能很大。直接計算組合數(shù)再求和是絕對不可能的更別提最后還有一個大指數(shù)的模冪運算。這道題的價值在于它提供了一個絕佳的“戰(zhàn)場”讓你系統(tǒng)性地演練如何解決“大組合數(shù)取模”和“大指數(shù)模運算”這兩個在密碼學(xué)、計算數(shù)學(xué)等領(lǐng)域非常實際的問題。如果你能獨立啃下這道題那么你對數(shù)論的理解和應(yīng)用能力會上一個大臺階。接下來我就帶你一步步拆解這道題不僅告訴你“怎么做”更重點講清楚“為什么這么做”以及我在反復(fù)提交和調(diào)試中積累的那些“血淚教訓(xùn)”。2. 核心思路拆解化整為零分而治之面對這種復(fù)合型難題最忌諱的就是一頭扎進去蠻干。我們必須先站在高處看清整個問題的結(jié)構(gòu)然后制定一個清晰的作戰(zhàn)計劃。整個計算過程可以分解為三個核心階段它們環(huán)環(huán)相扣2.1 第一階段目標(biāo)轉(zhuǎn)換與歐拉定理降冪我們的終極目標(biāo)是計算G^sum mod MOD其中MOD 999911659而sum Σ C(N, d) (d|N)。第一步就要注意到一個關(guān)鍵性質(zhì)999911659 是一個質(zhì)數(shù)。這是題目給出的重要條件也是我們所有后續(xù)操作的基石。對于一個質(zhì)數(shù)模數(shù)p我們有費馬小定理歐拉定理在質(zhì)數(shù)下的特例若G不是p的倍數(shù)則G^(p-1) ≡ 1 (mod p)。這里就遇到了第一個陷阱如果G是p的倍數(shù)即G % MOD 0那么結(jié)果直接就是0因為0的任何正數(shù)次冪對MOD取模都是0。這是一個需要特判的邊界情況。如果G不是p的倍數(shù)我們就可以利用歐拉定理進行降冪。因為φ(p) p-1所以對于指數(shù)sum我們可以將其對p-1取模G^sum ≡ G^(sum mod (p-1)) (mod p)這樣我們就把一個可能巨大的指數(shù)sum的問題轉(zhuǎn)化為了計算sum mod (p-1)的問題。而p-1 999911658。注意新的模數(shù)p-1不再是質(zhì)數(shù)了這為下一步埋下了伏筆。關(guān)鍵理解為什么模p-1這是歐拉定理的直接推論a^φ(n) ≡ 1 (mod n)當(dāng)gcd(a, n)1。這里np是質(zhì)數(shù)φ(p)p-1。所以a^k ≡ a^(k mod φ(p)) (mod p)。降冪是處理大指數(shù)模運算的標(biāo)準(zhǔn)操作。2.2 第二階段分解模數(shù)與Lucas定理登場現(xiàn)在問題轉(zhuǎn)化為計算sum_mod Σ C(N, d) mod (p-1)其中p-1 999911658。 直接計算C(N, d)對999911658取模依然極其困難因為模數(shù)不是質(zhì)數(shù)我們無法使用求逆元等簡便方法。這里就是第二個關(guān)鍵技巧分解模數(shù)。我們嘗試分解999911658999911658 2 * 3 * 4679 * 35617我們發(fā)現(xiàn)它恰好是四個互質(zhì)的質(zhì)數(shù)的乘積。這簡直是天賜良機它讓我們可以運用中國剩余定理CRT。中國剩余定理告訴我們要計算X mod MM m1*m2*...*mk且兩兩互質(zhì)等價于分別計算X mod m1,X mod m2, ...,X mod mk得到一組同余方程然后再用CRT合并回X mod M。因此我們只需要分別計算a1 sum mod 2a2 sum mod 3a3 sum mod 4679a4 sum mod 35617然后解同余方程組x ≡ a1 (mod 2) x ≡ a2 (mod 3) x ≡ a3 (mod 4679) x ≡ a4 (mod 35617)解出的x在模M999911658意義下就是我們要的sum_mod。那么如何計算sum mod mi呢這里mi是質(zhì)數(shù)2, 3, 4679, 35617。對于質(zhì)數(shù)模數(shù)下的組合數(shù)計算Lucas定理就派上用場了。Lucas定理對于質(zhì)數(shù)p有C(n, m) mod p C(n mod p, m mod p) * C(n/p, m/p) mod p。 這是一個遞歸過程它將大規(guī)模的組合數(shù)計算分解為對p以下小數(shù)字的組合數(shù)計算極大地簡化了問題。因為n mod p和m mod p都小于p我們可以預(yù)處理出0!到(p-1)!的階乘和階乘逆元然后O(1)地計算小組合數(shù)C(n%p, m%p)。所以對于每一個質(zhì)因子mi我們預(yù)處理該模數(shù)下的階乘數(shù)組fac[0..mi-1]和階乘逆元數(shù)組invfac[0..mi-1]。枚舉N的所有因數(shù)d。對于每個d使用 Lucas 定理計算C(N, d) mod mi。將所有C(N, d) mod mi相加并對mi取模得到ai。2.3 第三階段中國剩余定理合并與最終計算得到四個同余方程的余數(shù)a1, a2, a3, a4后我們使用中國剩余定理求解x。 設(shè)M 2*3*4679*35617 999911658Mi M / mi。 我們需要找到ti使得Mi * ti ≡ 1 (mod mi)即ti是Mi在模mi意義下的逆元。因為mi是質(zhì)數(shù)且Mi與mi互質(zhì)這個逆元一定存在可以用擴展歐幾里得算法或快速冪費馬小定理求出。那么方程組的解為x ≡ Σ(ai * Mi * ti) (mod M)。 計算出的x就是sum_mod sum mod 999911658。最后我們回到第一階段如果G % MOD 0輸出0。否則計算ans fast_pow(G, sum_mod, MOD)即快速冪取模得到最終答案。整個算法的流程圖如下邏輯描述輸入 N, G MOD 999911659 if G % MOD 0: print(0); return M MOD - 1 999911658 分解 M 為質(zhì)因子列表 primes [2, 3, 4679, 35617] 初始化余數(shù)數(shù)組 a[] for each p in primes: sum_p 0 預(yù)處理模 p 下的階乘和逆元階乘 for each divisor d of N: sum_p (sum_p lucas(N, d, p)) % p a[p] sum_p 使用中國剩余定理(CRT)根據(jù) a[] 和 primes[] 求解 sum_mod ans fast_pow(G, sum_mod, MOD) 輸出 ans3. 關(guān)鍵技術(shù)細節(jié)與實現(xiàn)要點思路清晰了但魔鬼藏在細節(jié)里。每個步驟的實現(xiàn)都有需要注意的坑點。3.1 質(zhì)因數(shù)分解與因數(shù)枚舉質(zhì)因數(shù)分解題目給定的MOD-1是固定的所以我們直接硬編碼分解結(jié)果[2, 3, 4679, 35617]即可不需要寫通用的分解函數(shù)。因數(shù)枚舉這是性能的關(guān)鍵點之一。N 最大為 1e9其因數(shù)個數(shù)不會太多1e9以內(nèi)因數(shù)最多的數(shù)其因數(shù)個數(shù)大約在 1300 多個。我們可以用O(sqrt(N))的方法枚舉所有因數(shù)。vectorlong long divisors; for (long long i 1; i * i N; i) { if (N % i 0) { divisors.push_back(i); if (i ! N / i) { // 避免重復(fù)添加平方根 divisors.push_back(N / i); } } }這樣得到的divisors數(shù)組包含了N的所有正因數(shù)。3.2 Lucas定理的高效實現(xiàn)Lucas定理的遞歸實現(xiàn)非常簡潔long long lucas(long long n, long long m, long long p) { if (m 0) return 1; // C(n, m) % p C(n%p, m%p) * lucas(n/p, m/p, p) % p return (C(n % p, m % p, p) * lucas(n / p, m / p, p)) % p; }其中C(n, m, p)函數(shù)用于計算當(dāng)n, m p時的組合數(shù)取模。這要求我們提前預(yù)處理模p下的階乘數(shù)組fac和階乘逆元數(shù)組invfac。預(yù)處理階乘與逆元fac[0] 1; for (int i 1; i p; i) { fac[i] fac[i-1] * i % p; } // 計算階乘逆元利用費馬小定理invfac[i] (i!)^(p-2) mod p invfac[p-1] fast_pow(fac[p-1], p-2, p); for (int i p-2; i 0; --i) { invfac[i] invfac[i1] * (i1) % p; }計算小組合數(shù)long long C(long long n, long long m, long long p) { if (n m) return 0; // 組合數(shù)定義n必須大于等于m // C(n, m) n! / (m! * (n-m)!) return fac[n] * invfac[m] % p * invfac[n - m] % p; }實操心得在實現(xiàn)lucas函數(shù)時一定要先判斷if (m 0)返回 1這是遞歸的基準(zhǔn)情況。另外在C函數(shù)中判斷if (n m)返回 0 非常重要因為在遞歸過程中n%p有可能小于m%p此時的組合數(shù)定義為 0。3.3 中國剩余定理CRT的合并實現(xiàn)我們有方程組x ≡ ai (mod mi)i1 to 4mi兩兩互質(zhì)。CRT的標(biāo)準(zhǔn)求解公式為計算總模數(shù)M m1 * m2 * m3 * m4。對于每個i計算Mi M / mi。計算ti是Mi在模mi意義下的逆元即Mi * ti ≡ 1 (mod mi)。解為x Σ(ai * Mi * ti) mod M。實現(xiàn)代碼long long crt(const vectorlong long a, const vectorlong long m) { long long M 1, x 0; int k a.size(); for (int i 0; i k; i) M * m[i]; for (int i 0; i k; i) { long long Mi M / m[i]; // 求 Mi 模 m[i] 的逆元 ti long long ti inv(Mi, m[i]); // 需要實現(xiàn)擴展歐幾里得求逆元 x (x a[i] * Mi % M * ti % M) % M; } return (x % M M) % M; // 確保返回正數(shù) }其中inv(a, mod)是求a在模mod下的逆元。由于這里的mod即mi都是質(zhì)數(shù)可以用快速冪fast_pow(a, mod-2, mod)來求。但為了通用性通常使用擴展歐幾里得算法。注意事項在累加x的過程中每次乘法后都要對M取模防止中間結(jié)果溢出。即使使用long longai * Mi * ti這三個數(shù)相乘也可能溢出所以需要步步取模。3.4 快速冪算法快速冪是基礎(chǔ)但必須寫對的算法。用于最后的G^sum_mod mod MOD計算以及在求逆元時計算a^(p-2) mod p。long long fast_pow(long long base, long long exp, long long mod) { long long result 1; base % mod; // 先取模防止 base 過大 while (exp 0) { if (exp 1) { result (result * base) % mod; } base (base * base) % mod; exp 1; } return result; }4. 完整代碼實現(xiàn)與逐段解析下面我將結(jié)合代碼詳細講解每個模塊的實現(xiàn)和聯(lián)動。我們假設(shè)輸入為N和G模數(shù)MOD 999911659。#include iostream #include vector #include cmath #include algorithm using namespace std; typedef long long ll; const ll MOD 999911659; // MOD-1 的質(zhì)因子分解 ll primes[4] {2, 3, 4679, 35617}; ll fac[36000]; // 最大質(zhì)因子是35617數(shù)組開大一點 ll invfac[36000]; ll a[4]; // 存儲 sum 對每個質(zhì)因子取模的結(jié)果 vectorll divisors; // 快速冪 ll qpow(ll base, ll exp, ll mod) { ll res 1; base % mod; while (exp) { if (exp 1) res (res * base) % mod; base (base * base) % mod; exp 1; } return res; } // 擴展歐幾里得求逆元 (ax ≡ 1 mod m) ll exgcd(ll a, ll b, ll x, ll y) { if (b 0) { 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; } // 預(yù)處理階乘和階乘逆元 (模 p) void init_fac(ll p) { fac[0] 1; for (ll i 1; i p; i) { fac[i] fac[i-1] * i % p; } // 費馬小定理求階乘逆元 invfac[p-1] qpow(fac[p-1], p-2, p); for (ll i p-2; i 0; --i) { invfac[i] invfac[i1] * (i1) % p; } } // 計算小組合數(shù) C(n, m) % p (要求 n, m p) ll C_small(ll n, ll m, ll p) { if (n m) return 0; return fac[n] * invfac[m] % p * invfac[n - m] % p; } // Lucas定理計算 C(n, m) % p ll lucas(ll n, ll m, ll p) { if (m 0) return 1; return (C_small(n % p, m % p, p) * lucas(n / p, m / p, p)) % p; } // 中國剩余定理合并 ll crt() { ll M MOD - 1; // 即 999911658 ll res 0; for (int i 0; i 4; i) { ll Mi M / primes[i]; ll ti inv(Mi, primes[i]); // 求 Mi 模 primes[i] 的逆元 res (res a[i] * Mi % M * ti % M) % M; } return (res % M M) % M; } int main() { ll N, G; cin N G; // 特判如果 G 是 MOD 的倍數(shù) if (G % MOD 0) { cout 0 endl; return 0; } // 1. 枚舉 N 的所有因數(shù) for (ll i 1; i * i N; i) { if (N % i 0) { divisors.push_back(i); if (i ! N / i) { divisors.push_back(N / i); } } } // 2. 對每個質(zhì)因子 p計算 sum_p Σ lucas(N, d, p) mod p for (int i 0; i 4; i) { ll p primes[i]; init_fac(p); // 預(yù)處理模 p 下的階乘 ll sum_p 0; for (ll d : divisors) { sum_p (sum_p lucas(N, d, p)) % p; } a[i] sum_p; // 記錄余數(shù) } // 3. 用中國剩余定理合并得到 sum_mod sum % (MOD-1) ll sum_mod crt(); // 4. 計算最終答案 G^sum_mod % MOD ll ans qpow(G, sum_mod, MOD); cout ans endl; return 0; }代碼逐段解析全局定義與輸入定義了模數(shù)、質(zhì)因子數(shù)組、階乘數(shù)組、余數(shù)數(shù)組和存儲因數(shù)的向量。讀入N和G。特判首先檢查G % MOD 0這是歐拉定理應(yīng)用的前提。如果為真直接輸出0結(jié)束。枚舉因數(shù)使用O(sqrt(N))的方法找出N的所有正因數(shù)存入divisors向量。核心循環(huán)對每個質(zhì)因子for (int i 0; i 4; i)遍歷四個質(zhì)因子2, 3, 4679, 35617。init_fac(p)針對當(dāng)前質(zhì)數(shù)p預(yù)處理0!到(p-1)!的階乘及其逆元。這是Lucas定理能高效計算的基礎(chǔ)。內(nèi)層循環(huán)for (ll d : divisors)遍歷N的每個因數(shù)d。sum_p (sum_p lucas(N, d, p)) % p使用Lucas定理計算C(N, d) mod p并累加到sum_p中同時保持取模。將最終累加結(jié)果sum_p存入a[i]。這樣我們就得到了四個同余方程的余數(shù)a[0]...a[3]。中國剩余定理合并調(diào)用crt()函數(shù)根據(jù)四個余數(shù)a[i]和模數(shù)primes[i]計算出sum_mod sum % (MOD-1)。最終計算與輸出使用快速冪qpow(G, sum_mod, MOD)計算最終答案并輸出。5. 常見問題與調(diào)試心得這道題在實現(xiàn)和提交時很容易遇到各種問題。下面是我總結(jié)的幾個常見坑點和調(diào)試技巧。5.1 溢出問題這是數(shù)論題最經(jīng)典的錯誤來源。中間結(jié)果溢出即使在long long64位范圍內(nèi)三個long long相乘也可能溢出。例如在crt()函數(shù)中a[i] * Mi * ti就可能超出2^63-1。務(wù)必在每次乘法后立即取模(a[i] * Mi % M * ti) % M??焖賰缰械囊绯鲈趒pow函數(shù)中base * base也可能溢出所以在函數(shù)開頭先base % mod是很好的習(xí)慣。階乘預(yù)處理溢出預(yù)處理fac[i] fac[i-1] * i % p時i最大為p-1而p最大為35617fac[i-1] * i不會超過1e10在long long范圍內(nèi)是安全的。排查技巧當(dāng)你懷疑溢出時可以在關(guān)鍵計算步驟后添加調(diào)試輸出打印中間變量的值看看是否突然變成了負數(shù)這是補碼溢出的典型表現(xiàn)。5.2 邊界條件與特判G % MOD 0這是最重要的特判。如果不加當(dāng)G是MOD倍數(shù)時gcd(G, MOD) ! 1歐拉定理不成立降冪步驟就是錯誤的會導(dǎo)致答案錯誤。N 的因數(shù) 1 和 N枚舉因數(shù)時1和N本身也是因數(shù)必須包含在內(nèi)。組合數(shù)C(N, 1) N,C(N, N) 1它們對sum有貢獻。Lucas 遞歸基準(zhǔn)情況lucas函數(shù)中必須判斷if (m 0) return 1。因為當(dāng)m為0時組合數(shù)C(n, 0) 1。C_small 中的 n m 判斷在遞歸調(diào)用中n % p可能小于m % p此時C(n%p, m%p)應(yīng)為0。這個判斷必須加上。5.3 中國剩余定理CRT的實現(xiàn)細節(jié)逆元的正確性crt()中需要求Mi模primes[i]的逆元ti。確保你的inv函數(shù)能正確工作。由于primes[i]是質(zhì)數(shù)用qpow(Mi, primes[i]-2, primes[i])求逆元更簡潔且不易錯。最終解的模數(shù)CRT求出的解x是模M意義下的。M MOD-1。最后返回(x % M M) % M是為了保證結(jié)果是[0, M-1]范圍內(nèi)的最小非負整數(shù)。模數(shù)一致性在計算Mi M / primes[i]時M必須是所有模數(shù)的乘積即999911658。5.4 性能優(yōu)化點因數(shù)枚舉O(sqrt(N))枚舉已經(jīng)足夠高效。對于N1e9循環(huán)次數(shù)約3e4到4e4次完全可以接受。預(yù)處理復(fù)用對于每個質(zhì)因子p我們都需要預(yù)處理階乘。由于p很小最大35617預(yù)處理是O(p)的總共做4次開銷很小。Lucas遞歸深度由于p最小是2遞歸深度大約是log_p(N)對于N1e9深度最多30層遞歸開銷可以忽略。5.5 調(diào)試與測試建議當(dāng)你寫完后可以用一些小的樣例進行測試。樣例1N2, G3因數(shù)1, 2sum C(2,1)C(2,2)213分別計算 mod 2,3,4679,35617 的余數(shù)都是 3 mod p。CRT合并后 sum_mod 3。答案 3^3 mod 999911659 27。 你可以手動計算驗證。樣例2N4, G2因數(shù)1, 2, 4sum C(4,1)C(4,2)C(4,4)46111計算 mod 2: 11%21; mod 3: 11%32; mod 4679: 11; mod 35617: 11。解同余方程組得到 sum_mod。最終計算 2^sum_mod mod 999911659。如果小樣例對了但提交到OJOnline Judge還是Wrong Answer可以嘗試檢查是否遺漏了G % MOD 0的特判。檢查crt()函數(shù)中乘法取模是否寫全了% M。檢查lucas函數(shù)中的C_small調(diào)用是否包含了n m返回 0 的判斷。用cout打印出中間變量比如四個余數(shù)a[i]、CRT合并后的sum_mod與手算或?qū)ε某绦虻慕Y(jié)果進行對比。這道題綜合性很強幾乎考察了數(shù)論競賽中的所有常用技巧。把它吃透對于理解模運算、組合數(shù)計算、定理的綜合應(yīng)用有極大的好處。在實際編碼中耐心和細心是關(guān)鍵尤其是處理取模和邊界條件時。希望這篇詳細的拆解能幫助你徹底掌握“古代豬文”這道經(jīng)典題目。