第一次看到 P2398 GCD SUM 這個題我的第一反應(yīng)是這不就是兩個 for 套一個 gcd 嗎題目要求的東西非常直白給定一個正整數(shù) n計算所有 1≤i≤n、1≤j≤n 的數(shù)對 (i,j) 的 gcd 之和。但刷題這么多年我早就養(yǎng)成了一個習(xí)慣先別急著寫循環(huán)先看一眼數(shù)據(jù)范圍再決定策略。這道題的 n 上限是 100000如果直接枚舉 i 和 j數(shù)對總量是 10^10就算 gcd 寫得再快也扛不住。所以這道題真正的考點不是“你會不會求 gcd”而是“你能不能把求和式子化簡成一個可以預(yù)處理、可以分塊計算的形式”。這篇文章我就把完整的推導(dǎo)過程、代碼實現(xiàn)和幾個容易踩的坑一起寫出來適合剛接觸數(shù)論題、想搞懂歐拉函數(shù)和數(shù)論分塊怎么配合使用的朋友。1. 先看清題目到底在求哪個東西1.1 暴力三重循環(huán)為什么一上來就該被否掉題面濃縮起來就一行ans ∑_{i1}^{n} ∑_{j1}^{n} gcd(i, j)很多新手看到這個式子的直覺就是枚舉 i、j然后調(diào) __gcd 累加。代碼確實短短到十行以內(nèi)就能寫完long long ans 0; for (int i 1; i n; i) for (int j 1; j n; j) ans gcd(i, j);但這種寫法在 n100000 時完全跑不動。n2 已經(jīng)是 10^10gcd 內(nèi)部還要做若干次取模運算總計算量輕松突破 10^11。評測環(huán)境里這屬于必超時的方案代碼再短也沒有意義。競賽里遇到“求和式子里帶著 gcd、lcm 這類數(shù)論函數(shù)”的題目我的經(jīng)驗是立刻換思路不去枚舉數(shù)對而是枚舉函數(shù)值。gcd(i,j) 的取值只可能是 1 到 n 之間的整數(shù)那我們就按 gcd 的取值 d 來分桶數(shù)清楚每個桶里有多少個數(shù)對最后統(tǒng)一加權(quán)求和。1.2 枚舉 d把求和對象從“數(shù)對”換成“gcd 值”設(shè) cnt[d] 表示滿足 gcd(i,j)d 的有序數(shù)對 (i,j) 數(shù)量。因為每一對數(shù)對都會恰好落入某一個 d 對應(yīng)的桶里所以答案可以寫成ans ∑_{d1}^{n} d * cnt[d]這個變形把問題拆成了兩塊外面這層 d 的求和很好處理后面甚至能用等差數(shù)列合并真正麻煩的是 cnt[d] 怎么求。如果只是求“gcd(i,j) 是 d 的倍數(shù)”的數(shù)對個數(shù)那太簡單了i 和 j 都得是 d 的倍數(shù)數(shù)量就是 ?n/d?2。但“恰好等于 d”要嚴(yán)格得多它要求 i/d 和 j/d 不能再有公共因子。這一步處理好了整道題就通了。2. 核心推導(dǎo)用 φ 函數(shù)統(tǒng)計互質(zhì)數(shù)對2.1 從“恰好等于 d”到“a、b 互質(zhì)”我們先把“gcd(i,j)d”翻譯成更可操作的條件。如果最大公約數(shù)正好是 d那么 i 和 j 都是 d 的倍數(shù)可以寫成i d * a j d * b并且 gcd(a,b) 必須等于 1。反過來也成立只要 a、b 互質(zhì)那么 da 和 db 的最大公約數(shù)就恰好是 d不會多也不會少。因為 i、j 都不超過 n所以 a、b 都在 1 到 ?n/d? 之間。令 m?n/d?那么 cnt[d] 就變成了一個純粹的計數(shù)問題在 a,b ∈ [1,m] 中滿足 gcd(a,b)1 的有序數(shù)對有多少個“有序”這個詞要劃重點因為題目里的 (i,j) 和 (j,i) 是不同的數(shù)對后面計數(shù)公式里那個“減 1”就是因為這個。2.2 互質(zhì)有序?qū)τ嫈?shù)公式2×φ前綴和?1先回憶歐拉函數(shù) φ(x) 的定義1 到 x 中與 x 互質(zhì)的正整數(shù)個數(shù)。它天然適合回答“如果固定 b x那么有多少個 a≤b 滿足 gcd(a,b)1”答案就是 φ(x)。如果我們把 b 從 1 枚舉到 m把 φ(b) 累加起來得到的是所有滿足 a≤b 且互質(zhì)的有序?qū)?shù)量。這個數(shù)量我記為S(m) ∑_{b1}^{m} φ(b)但這只統(tǒng)計了“a 在左邊、b 在右邊”的一半情況。對于任意 ab 且互質(zhì)的數(shù)對題目里還應(yīng)該有 (b,a) 這一組它們互為對稱。除了 (1,1) 這種對角線上的特例其他每個互質(zhì)對都能對稱出一組新的有序?qū)?。所以最終的有序互質(zhì)對總數(shù)是cnt[d] 2 * S(m) - 1減掉的 1 就是 (1,1) 被重復(fù)計算的那一次。把 φ(1)1 也納入 S(m)公式就變得非常干凈2 倍前綴和再整體減 1。2.3 手算驗證公式不會漏也不會重這個公式第一次見的人容易心里打鼓我建議自己手算一遍。拿 m3 舉例a、b 都在 {1,2,3} 里互質(zhì)的有序?qū)σ还彩窍旅?7 組(1,1)、(1,2)、(2,1)、(1,3)、(3,1)、(2,3)、(3,2)而 φ(1)1φ(2)1φ(3)2所以 2S(3)-1 2(112)-1 7完全一致。再看 m4φ(4)2S(4)62*6-111。把所有 16 組有序?qū)α谐鰜砣サ?(2,2)、(3,3)、(4,4) 以及 (2,4)、(4,2)剩下的確實是 11 組。算到這里答案表達(dá)式已經(jīng)浮出水面了ans ∑_{d1}^{n} d * ( 2*S(?n/d?) - 1 )現(xiàn)在的問題就變成了怎么快速得到 S(m)也就是 φ 的前綴和。2.4 用樣例數(shù)據(jù)把公式跑一遍拿 n4 手算驗證。φ(1)1φ(2)1φ(3)2φ(4)2所以 S(1)1S(2)2S(3)4S(4)6。d?n/d?2*S(?n/d?)-1d * 該項141111223631134114累加得到 1163424和題目樣例輸出一致。公式從第一道坎就站住了后面可以放心往下寫代碼。3. 預(yù)處理歐拉函數(shù)線性篩而不是逐個分解質(zhì)因數(shù)3.1 積性函數(shù)與三個關(guān)鍵遞推公式里需要 S(m)也就是 φ 的前綴和。n 最大才 100000第一反應(yīng)可能是對每個數(shù)暴力分解質(zhì)因數(shù)再套 φ(n)n*∏(1-1/p) 的公式。那樣做總復(fù)雜度大概是 O(n√n)在這個數(shù)據(jù)范圍也能跑但不夠優(yōu)雅而且一旦題目把 n 加到 10^7這個方案立刻崩盤。正確做法是線性篩。歐拉函數(shù)是積性函數(shù)只要 gcd(a,b)1就有 φ(a*b)φ(a)*φ(b)。線性篩求 φ 依賴下面三條遞推若 p 是素數(shù)那么 φ(p)p-1若當(dāng)前數(shù) i 與素數(shù) p 互質(zhì)那么 φ(i*p)φ(i)φ(p)φ(i)(p-1)若 p 整除 i說明 p 已經(jīng)出現(xiàn)在 i 的質(zhì)因數(shù)分解中此時 φ(i*p)φ(i)*p。第 3 條初看有點反直覺p 已經(jīng)存在了為什么只是乘 p 而不是乘 p-1因為 p 作為一個質(zhì)因子早就在 i 里了i*p 和 i 的質(zhì)因數(shù)集合完全相同互質(zhì)數(shù)的“占比”沒有變只是數(shù)本身變大了 p 倍所以 φ 直接乘 p。3.2 線性篩求 φ 的完整寫法一個標(biāo)準(zhǔn)的篩法實現(xiàn)如下const int MAXN 100000 5; int phi[MAXN]; long long prefPhi[MAXN]; bool notPrime[MAXN]; vectorint primes; void init(int n) { phi[1] 1; for (int i 2; i n; i) { if (!notPrime[i]) { primes.push_back(i); phi[i] i - 1; } for (int p : primes) { long long v 1LL * i * p; if (v n) break; notPrime[v] true; if (i % p 0) { phi[v] phi[i] * p; break; } else { phi[v] phi[i] * (p - 1); } } } for (int i 1; i n; i) { prefPhi[i] prefPhi[i - 1] phi[i]; } }這里有幾個細(xì)節(jié)值得專門說一下。phi[1] 必須手動設(shè)置為 1。這不是可選項而是公式成立的前提。前文所有涉及 S(m) 的地方如果 φ(1) 是 0那么 2*S(m)-1 從一開始就全錯了。當(dāng) i%p0 時一定要 break。歐拉篩之所以是線性的靠的就是“每個合數(shù)只被它的最小質(zhì)因子篩掉一次”。如果不 break同一個合數(shù)會被多個質(zhì)因子重復(fù)標(biāo)記性能會退化到接近埃氏篩甚至更差。v 用 long long 和 n 比較這個習(xí)慣能防止 i*p 在 int 里溢出。本題 n 只有 1e5 還看不出問題但如果哪天把這份代碼搬到 n1e9 的題目里這個細(xì)節(jié)能直接救命。3.3 為什么不直接對每個數(shù)求 φ有人會覺得這道題 n 才 1e5我暴力分解每個數(shù)求 φ 也很快啊何必費勁寫線性篩這話有一定道理在只求 AC 的前提下暴力確實能過。但刷題不能只盯著單題歐拉函數(shù)、莫比烏斯函數(shù)這類積性函數(shù)的線性篩預(yù)處理是數(shù)論題里復(fù)用率極高的工具。我見過太多人平時不練篩法一到需要預(yù)處理 φ 的題就臨時寫暴力最后在小數(shù)據(jù)上 AC換個大數(shù)據(jù)就被卡死。把線性篩寫成肌肉記憶遇到 n 很大的同類題目能直接無縫遷移這個投資很值。另外線性篩求 φ 的代碼量真的很小和暴力分解相比并沒有復(fù)雜多少。與其留一個“數(shù)據(jù)一變大就要重構(gòu)”的隱患不如一次寫到位。4. 數(shù)論分塊把最后的求和壓到 O(√n)4.1 觀察 ?n/d? 的分段特性公式已經(jīng)到手直接 for d1 到 n 累加復(fù)雜度 O(n)對于 1e5 的數(shù)據(jù)完全能過。但我想多說一步因為這道題的答案公式天生適合數(shù)論分塊學(xué)會之后能應(yīng)付更大數(shù)據(jù)范圍。觀察 ?n/d? 這個值的變化。它隨 d 增大而遞減但不是每一步都變而是一段一段地跳。比如 n10 時d1234、56~10?n/d?105321值相同的 d 會聚成連續(xù)的區(qū)間整個 [1,n] 只被分成大約 2√n 段。理解這個結(jié)論可以換個角度當(dāng) d≤√n 時d 本身最多只有 √n 種取法當(dāng) d√n 時?n/d? √n商最多也只有 √n 種取法。兩邊一夾段數(shù)自然就是 O(√n)。既然每一段里的 ?n/d? 都相同那么 (2*S(?n/d?)-1) 在這一整段里都是常數(shù)完全可以提到求和號外面把 d 的貢獻(xiàn)一次性算掉。4.2 分塊的寫法與代碼細(xì)節(jié)假設(shè)當(dāng)前處理到的左端點是 l令 m n/l那么滿足 ?n/d? m 的最遠(yuǎn)右端點就是 r n/m。這一段的貢獻(xiàn)為∑_{dl}^{r} d * (2S(m)-1) (2S(m)-1) * ∑_{dl}^{r} d后面那個 d 的求和是等差數(shù)列直接用 (lr)*(r-l1)/2 計算。完整 AC 代碼如下#include bits/stdc.h using namespace std; int main() { int n; cin n; vectorint phi(n 1); vectorlong long pref(n 1); vectorbool notPrime(n 1, false); vectorint primes; primes.reserve(n / 10); phi[1] 1; for (int i 2; i n; i) { if (!notPrime[i]) { primes.push_back(i); phi[i] i - 1; } for (int p : primes) { long long v 1LL * i * p; if (v n) break; notPrime[v] true; if (i % p 0) { phi[v] phi[i] * p; break; } phi[v] phi[i] * (p - 1); } } for (int i 1; i n; i) pref[i] pref[i - 1] phi[i]; long long ans 0; for (long long l 1, r; l n; l r 1) { long long m n / l; r n / m; long long sumD (l r) * (r - l 1) / 2; ans sumD * (2 * pref[(int)m] - 1); } cout ans \n; return 0; }分塊最容易寫錯的就是右端點 rn/m。我第一次學(xué)這個技巧的時候總想著 rn/(m1) 之類更“精確”的式子結(jié)果把自己繞暈了。記住一個固定套路m 由左端點算出來右端點一定等于 n/m別猶豫。至于為什么 l、r、m、sumD 都用 long long是因為雖然 n 只有 1e5但 (lr)*(r-l1) 這個中間量在更大范圍內(nèi)會超過 int。寫題時寧可全開 long long 也不用去賭“這次數(shù)據(jù)沒那么大”。4.3 三種算法的復(fù)雜度對比方案復(fù)雜度n1e5 時的情況枚舉 i、j 再求 gcdO(n2 log n)完全不可行線性篩 φ d 從 1 到 n 直接循環(huán)O(n)可以通過線性篩 φ 數(shù)論分塊O(n √n)可以通過且常數(shù)極小如果只求 ACO(n) 的樸素求和已經(jīng)足夠。但我還是推薦把分塊寫上去原因很簡單答案表達(dá)式天然就是“常數(shù)段求和”的形狀分塊寫起來不比直接循環(huán)復(fù)雜遇到 n 變成大數(shù)的同類型題還能繼續(xù)用。5. 驗證、易錯點和另一種推導(dǎo)路線5.1 邊界數(shù)據(jù)與一致性檢查代碼寫完別急著交先手測幾個小數(shù)據(jù)n1只有 (1,1) 一組答案就是 1。分塊時 [1,1] 這一段的 m1pref[1]12*1-11sumD1答案 1正確。n2四對 gcd 分別是 1、1、1、2總和 5。公式d1 時 1*(2S(2)-1)3d2 時 2(2*S(1)-1)2總和 5正確。n4前面手算過24與樣例一致。更穩(wěn)妥的做法是寫一個暴力枚舉函數(shù)和正解對拍幾百組隨機(jī)小數(shù)據(jù)。這種方法能一次性覆蓋大量邊界情況比肉眼檢查放心得多。5.2 三個容易踩的坑第一個坑答案和中間量必須用 long long。n100000 時最終答案在 10^10 級別int 直接溢出。這沒什么技巧凡是可能變大的量一律 long long。第二個坑φ(1) 忘記初始化為 1。篩法循環(huán)通常從 2 開始phi 數(shù)組默認(rèn)全是 0如果少了 phi[1]1 這一步pref[1] 就會變成 0后面每一段的 (2*S-1) 都會系統(tǒng)性偏小。這已經(jīng)不是“差一點點”的問題而是公式整體失效。第三個坑分塊右端點寫錯。典型癥狀是小數(shù)據(jù)答案沒問題n 變大后答案偏偏少一段。右端點如果算小了會漏掉區(qū)間尾部的數(shù)對算大了又會把下一段的值混進(jìn)來導(dǎo)致后續(xù) l、r 全部錯位。所以務(wù)必記住 rn/m 這個固定寫法。5.3 莫比烏斯反演視角殊途同歸如果熟悉莫比烏斯反演這道題還有另一條推導(dǎo)路線。設(shè) F(k) 表示 gcd(i,j) 是 k 的倍數(shù)的數(shù)對個數(shù)顯然 F(k)?n/k?2。再設(shè) f(k) 表示 gcd(i,j) 恰好等于 k 的數(shù)對個數(shù)由莫比烏斯反演f(k) ∑_{t1}^{?n/k?} μ(t) * F(k*t)代入答案式并換元 Tk*t可以得到ans ∑_{T1}^{n} ?n/T?2 * ∑_{k|T} k*μ(T/k)內(nèi)層 ∑_{k|T} k*μ(T/k) 正是 φ(T) 的經(jīng)典表達(dá)。于是最終得到另一個同樣簡潔的公式ans ∑_{T1}^{n} φ(T) * ?n/T?2這個形式和前面推出來的 2*S(?n/d?)-1 看起來完全不同但數(shù)值上完全一致。兩條路線最后都落在歐拉函數(shù)上這件事本身就說明這道題的數(shù)學(xué)本質(zhì)就是對互質(zhì)數(shù)對計數(shù)φ 是繞不開的主角。我個人寫這道題更偏愛第一種推導(dǎo)因為每一步都有明確的計數(shù)含義代碼里只需要 φ 前綴和不用額外處理 μ。但如果你想檢驗自己的反演功底第二種推導(dǎo)是很好的練習(xí)題也能作為題解里的“彩蛋”寫給讀者看。這道題隔一段時間重新做一次我都會提醒自己同一個原則看到雙重求和先想換序看到 gcd 先想按值分桶。把“枚舉數(shù)對”轉(zhuǎn)化成“枚舉函數(shù)值”這一步想通了后面線性篩、數(shù)論分塊都只是常規(guī)動作。