算法詳解:從公式推導到NTT實現與調試)
1. 從一道模板題說起多項式對數函數ln到底是什么如果你在洛谷、Codeforces或者任何一個算法競賽社區混跡過一段時間大概率會刷到過“P4725 【模板】多項式對數函數多項式 ln”這道題。它就像算法競賽選手在多項式領域的一個“成人禮”標志著從只會做加減乘除的“小學生”進階到開始觸及多項式更深刻運算的“中學生”。但很多人在第一次接觸時都會懵多項式還能取對數這玩意兒有什么用難道是把1 2x 3x^2丟進計算器按ln鍵嗎顯然不是。這里的“多項式對數函數”是一個形式冪級數上的形式運算。它解決的核心問題是給定一個常數項為1的多項式或形式冪級數A(x)求另一個多項式B(x)使得在形式冪級數的意義下exp(B(x)) A(x)。這里 exp 是指數函數。換句話說B(x) 就是 A(x) 的“形式對數”。這個運算在組合數學、生成函數、多項式算法中有著極其重要的地位。比如當你用生成函數刻畫一個組合結構時對其取 ln 往往對應著將連通分量拆解出來在多項式牛頓迭代求解中ln 和 exp 是一對關鍵的基礎算子。這道模板題之所以經典是因為它完美地將多項式求導、積分、求逆、乘法這幾個基礎操作串聯了起來形成了一個完整的算法鏈條。網上能找到的題解和代碼很多但大多只給出了“怎么做”的步驟和代碼對于“為什么這么做”、“每一步背后的數學原理是什么”、“實現時有哪些一踩就炸的坑”卻語焉不詳。我這篇文章就想結合我多次實現和調試的經驗把這些隱藏在水面下的東西徹底講透。我們不止要會套模板更要理解這個模板的每一顆螺絲釘是怎么擰上去的。2. 核心公式推導為什么求ln變成了求導、求逆和積分幾乎所有教程都會直接甩給你這個公式 若 A(x) 1 a_1 x a_2 x^2 ...且 A(0)1則ln(A(x)) ∫ [A(x) / A(x)] dx這個公式是整套算法的基石。我們來一步步拆解它看它到底是怎么來的。2.1 從形式微分的定義出發首先我們得認同對形式冪級數也可以定義“導數”。這很直觀對于多項式 A(x) ∑_{i0}^{n} a_i x^i其形式導數 A(x) 就是 ∑_{i1}^{n} i * a_i x^{i-1}。就是把每一項的指數拿下來當系數然后指數減一。現在考慮我們想求的 B(x) ln(A(x))。這里 ln 是一個形式運算。我們對這個等式兩邊同時關于 x 求形式導數利用鏈式法則左邊B(x) 右邊d/dx [ln(A(x))] A(x) / A(x) 這里直接類比了實數域上 ln(f(x)) 的導數為 f(x)/f(x)于是我們得到了一個關鍵等式B(x) A(x) / A(x)。2.2 從微分到積分得到了 B(x) 的導數那么 B(x) 本身自然就是對其積分B(x) ∫ B(x) dx ∫ [A(x) / A(x)] dx注意這里積分會有一個積分常數 C。因為是不定積分。那么 C 是多少我們需要利用初始條件A(0) 1。我們希望 B(x) 也是一個形式冪級數并且通常定義 ln(1) 0。所以 B(0) ln(A(0)) ln(1) 0。另一方面我們對 ∫ [A(x) / A(x)] dx 求出的結果其常數項就是積分產生的常數 C。為了讓 B(0) 0我們必須令 C 0。所以最終公式里我們直接寫為定積分形式從 0 積到 x或者理解為取不定積分后忽略常數項因為常數項為0。在實現時我們做不定積分然后手動將結果的常數項設為0即可。2.3 公式的可行性分析這個公式將 ln 運算轉化為了三個我們已知能做的操作求導 (A(x))O(n) 復雜度極其簡單。求逆 (1 / A(x))這里需要計算 A(x) 的乘法逆元。這需要用到多項式求逆算法通常使用牛頓迭代法復雜度 O(n log n)。積分 (∫ ... dx)O(n) 復雜度是求導的逆過程同樣簡單。所以整個多項式 ln 的算法復雜度就卡在了多項式求逆這一步為 O(n log n)。這也就是為什么多項式求逆是多項式全家桶里更基礎的一個模板。注意這個公式成立有一個絕對的前提A(x) 的常數項必須為 1。為什么 從數學上看ln(A(x)) 要想展開成形式冪級數必須在 x0 處有定義且 ln(A(0)) 需要是一個有限值我們取0。如果 A(0)0ln(0) 無定義如果 A(0) 是其他非1常數 c那么 ln(A(x)) ln(c) ln(1 (A(x)-c)/c)。這里 ln(c) 是一個實數常數但我們的多項式是在某個模數如998244353的有限域上運算的ln(c) 在這個域里可能沒有定義除非 c 是模數的原根相關。為了保證運算純粹在模意義下進行且結果是一個多項式常數項為0最方便且通用的約定就是要求 A(0)1。這樣 ln(1)0一切都很干凈。3. 算法步驟拆解與零基礎實現指南理解了公式我們來把算法步驟徹底細化。假設我們有多項式 A(x)其次數為 n-1通常我們處理長度為 n 的數組下標 0 到 n-1 對應次數 0 到 n-1 的系數且滿足 A[0] 1。3.1 第一步計算 A(x) 的導數 A(x)這一步是熱身。設 A(x) a0 a1x a2x^2 ... a_{n-1}x^{n-1}。 那么 A(x) a1 2a2x 3a3*x^2 ... (n-1)*a_{n-1}*x^{n-2}。在代碼中這就是一個簡單的循環// 假設系數存儲在數組 a 中長度為 n for (int i 1; i n; i) { da[i-1] 1LL * a[i] * i % mod; // da 存儲導數系數 } // da 的有效長度變為 n-1注意邊界導數的次數比原多項式低一次。3.2 第二步計算 A(x) 的乘法逆元 B(x) 1 / A(x)這是整個算法的核心和性能瓶頸。我們需要求一個多項式 B(x)使得 A(x) * B(x) ≡ 1 (mod x^n)。這里mod x^n的意思是我們只關心乘積的前 n 項0 到 n-1 次更高次的項可以忽略。求逆通常使用牛頓迭代法。其思想是假設我們已經求出了在模 x^{ceil(m/2)} 意義下的逆元 B_0(x)如何快速得到在模 x^m 意義下的逆元 B(x)推導過程略涉及泰勒展開結論是迭代公式為B(x) ≡ B_0(x) * (2 - A(x) * B_0(x)) (mod x^m)實際操作時我們采用遞歸或迭代倍增的方式初始條件當 n1 時A(x) 只有一個常數項 a0。由前提 a0 1所以在模 x^1 意義下其逆元就是 1。假設我們已經求出在模 x^{ceil(n/2)} 意義下的逆元 B_0(x)。目標計算模 x^n 意義下的逆元 B(x)。根據公式我們需要計算T(x) A(x) * B_0(x) (mod x^n)// 注意這里模數要提升到 x^nT(x) (2 - T(x)) (mod x^n)// 對 T(x) 的每一項做 2 - t_i 運算B(x) T(x) * B_0(x) (mod x^n)這個過程需要多項式乘法。利用 NTT快速數論變換可以將乘法優化到 O(n log n)。由于這部分是獨立模板代碼較長。其關鍵點在于每次迭代時對于 A(x) 我們只需要前 m 項當前目標長度對于 B_0(x) 我們知道它在模 x^{m/2} 下是精確的。計算A(x) * B_0(x)時結果長度會增長但我們只取前 m 項。然后進行2 - T(x)的系數運算最后再乘一次 B_0(x) 并取前 m 項。3.3 第三步計算 C(x) A(x) * B(x)現在我們有了導數da長度為 n-1和逆元b長度為 n。將它們相乘。注意da的次數是 n-2b的次數是 n-1它們的乘積次數最高為 (n-2)(n-1)2n-3。但我們最終只需要前 n-1 項因為下一步積分后我們要得到 n 項結果。所以我們可以只計算到長度至少為 n-1 的卷積。設dc da * b。我們取dc的前 n-1 項。注意dc[0]對應的是A(x)*B(x)的常數項。3.4 第四步對 C(x) 積分得到最終結果積分是導數的逆運算。如果C(x) c0 c1*x c2*x^2 ...那么它的積分∫ C(x) dx C c0*x (c1/2)*x^2 (c2/3)*x^3 ...其中 C 是積分常數。我們已經知道結果的常數項必須為 0。所以我們計算 對于 i 從 0 到 n-2res[i1] dc[i] * inv(i1) % mod其中inv(i1)是 i1 在模 mod 下的乘法逆元需要預處理。 而res[0] 0。這樣得到的res就是ln(A(x))的前 n 項系數。3.5 完整流程圖示與復雜度分析輸入: A(x), 滿足 A[0]1, 次數界 n 輸出: B(x) ln(A(x)) mod x^n 1. 求導: DA(x) derivative(A(x)) // O(n) 2. 求逆: IA(x) inverse(A(x), n) // O(n log n) 使用牛頓迭代NTT 3. 乘法: C(x) DA(x) * IA(x) mod x^{n-1} // O(n log n) NTT乘法結果取前n-1項 4. 積分: B(x) integral(C(x)) // O(n) 常數項設為0 5. 返回 B(x)總時間復雜度由兩次 O(n log n) 的操作主導即求逆和乘法。空間上需要一些臨時數組進行變換和計算。4. 實戰代碼剖析從模塊構建到邊界處理光說不練假把式。下面我結合一個典型的基于 NTT模數 998244353原根為 3的實現來逐塊解析代碼并指出那些容易寫錯、調試到崩潰的細節。4.1 基礎工具函數快速冪與逆元const int mod 998244353, g 3; // 原根 int qpow(int a, int b) { int res 1; while (b) { if (b 1) res 1LL * res * a % mod; a 1LL * a * a % mod; b 1; } return res; }qpow是標準快速冪。inv函數通常直接調用qpow(a, mod-2)但頻繁調用時建議預處理 1~n 的逆元。4.2 核心NTT 與多項式乘法這是所有多項式操作的基礎。代碼較長但結構固定。關鍵點在于rev數組的蝴蝶變換以及三層循環的迭代實現。這里我強調幾個易錯點長度與界限NTT 要求長度是 2 的冪。每次進行多項式乘法前必須計算lim 1, bit 0; while (lim n m) lim 1, bit;。然后初始化 rev 數組。結果清零對于長度為lim的數組一定要確保lim范圍內的數據是有效的或者在使用前清空。特別是多次調用時舊數據可能殘留。逆變換后的縮放NTT 逆變換后每一項需要乘以lim的逆元。int invlim qpow(lim, mod-2); for (int i0; ilim; i) a[i]1LL*a[i]*invlim%mod;4.3 多項式求逆的實現細節這是最難寫對的部分。我給出一個相對清晰的迭代版本框架void poly_inv(int *a, int *b, int n) { // 計算 b(x)使得 a(x)*b(x) ≡ 1 (mod x^n) static int tmp[N]; // 臨時數組需要足夠大如4倍n b[0] qpow(a[0], mod-2); // 初始條件常數項逆元 for (int len 2; (len 1) n; len 1) { // 當前目標是求出模 x^len 下的逆元 int lim len 1; // 乘法需要長度 // 將 a 的前 len 項拷貝到 tmp并做 NTT for (int i 0; i len; i) tmp[i] i n ? a[i] : 0; for (int i len; i lim; i) tmp[i] b[i] 0; ntt(tmp, lim, 1); ntt(b, lim, 1); // 根據公式 B B0 * (2 - A * B0) 計算 for (int i 0; i lim; i) { b[i] 1LL * b[i] * (2 - 1LL * tmp[i] * b[i] % mod mod) % mod; } ntt(b, lim, -1); // 重要b 中 len 之后的項是無意義的必須清零防止影響下一輪 for (int i len; i lim; i) b[i] 0; } // 最后確保 b 只有前 n 項有效后面清零如果傳入的n不是2的冪 for (int i n; i lim; i) b[i] 0; }踩坑實錄1清零清零清零這是多項式題最經典的錯誤。在牛頓迭代的每一輪結束后b數組在len之后的系數必須手動設為0。因為 NTT 逆變換后這些位置可能留有上一輪或計算過程中的垃圾值。下一輪循環時我們會把整個b數組包括后面的垃圾值做 NTT這些垃圾值會污染整個頻域導致結果完全錯誤。這個 bug 非常隱蔽因為小數據時可能因為長度不夠碰不到垃圾內存而僥幸正確大數據一定掛。4.4 多項式 ln 的完整實現集成了求導、求逆、乘法和積分。void poly_derivative(int *a, int *da, int n) { for (int i 1; i n; i) da[i-1] 1LL * a[i] * i % mod; da[n-1] 0; // 導數長度減一最后一位可置0 } void poly_integral(int *a, int *ia, int n) { ia[0] 0; // 常數項為0 // 預處理1~n的逆元 inv[i] for (int i 1; i n; i) ia[i] 1LL * a[i-1] * inv[i] % mod; } void poly_ln(int *a, int *res, int n) { // 前提檢查a[0] 必須為 1 assert(a[0] 1); static int da[N], ia[N], tmp[N]; // 1. 求導 poly_derivative(a, da, n); // da 長度為 n-1 // 2. 求逆 poly_inv(a, ia, n); // ia 是 A(x) 的逆長度為 n // 3. 乘法da * ia int lim 1, bit 0; while (lim (n-1) n) lim 1, bit; // 結果需要前 n-1 項 for (int i 0; i lim; i) tmp[i] (i n-1) ? da[i] : 0; for (int i 0; i lim; i) { // 注意ia 只有前 n 項有效后面在 poly_inv 中已清零 // 但這里為了安全可以只拷貝前 n 項后面置0 if (i n) ia[i] ia[i]; else if (i lim) ia[i] 0; // 確保 ia 在 lim 長度內有效 } // 這里需要調用一個標準的 NTT 乘法函數輸入 tmp 和 ia結果存回 tmp ntt_mul(tmp, ia, lim); // 假設這個函數處理了 NTT 變換、點乘、逆變換和縮放 // 現在 tmp 的前 n-1 項是 da * ia 的結果 // 4. 積分 poly_integral(tmp, res, n); // 積分后長度變為 n }踩坑實錄2長度對齊與數組越界在調用poly_inv時我們傳入的長度是n它會計算模x^n的逆。但注意poly_inv內部可能按2的冪分配內存。我們傳給poly_ln的數組a其有效長度就是n。但在求導后da有效長度是n-1。進行乘法da * ia時da長度n-1ia長度n卷積結果長度至少需要(n-1)(n-1)2n-2才能保證前n-1項精確。我們設置的lim必須大于等于這個值。同時要確保傳入 NTT 乘法的數組在lim長度內都有定義要么是有效系數要么是0。任何未初始化的值都會導致錯誤。5. 調試技巧與常見問題排查就算你完全理解了算法第一遍代碼也幾乎不可能一次 AC。以下是我總結的排查清單5.1 結果完全不對/隨機數檢查 NTT 的正確性這是根源。寫一個簡單的測試比如計算(12x) * (13x)看結果是不是1 5x 6x^2。確保正變換、點乘、逆變換、縮放每一步都正確。檢查數組清零如 4.3 節所述在poly_inv的每一輪迭代后必須清零b數組len之后的部分。在poly_ln中調用 NTT 前確保tmp和ia在lim范圍內的數據是干凈的。檢查長度計算lim是否足夠大while (lim n m) lim 1中的n和m是否正確對于da * ian n-1da的長度m nia的長度所以lim需要至少2n-1的下一個2的冪。5.2 結果前幾項對后面錯檢查求逆的邊界poly_inv函數是否保證了結果嚴格只有前n項有效在倍增過程中我們計算的是模x^len的逆但len可能大于n。函數最后需要把n之后的系數清零。檢查積分用的逆元表inv[i]是否預處理正確inv[i]是i在模mod下的逆元通常用線性遞推inv[i] mod - 1LL * (mod/i) * inv[mod%i] % mod來求。確保inv[1] 1。5.3 常數項不為0檢查輸入確認輸入多項式a[0]是否真的為 1。題目可能不保證需要自己先判斷或處理。檢查積分函數poly_integral是否將res[0]設為了 05.4 性能問題TLENTT 的蝴蝶變換rev數組是否預處理每次乘法都重新計算會超時。不必要的拷貝在poly_inv和poly_ln中盡量減少大數組的memcpy操作。使用指針和就地計算。乘法優化對于da * ia我們只需要前n-1項。可以使用“半在線卷積”的思路進行優化但模板題通常不需要標準的 NTT 乘法即可通過。5.5 一個實用的調試方法對拍寫一個暴力版本的poly_ln用于小數據范圍比如 n 10的驗證。 暴力版本可以模擬形式冪級數的運算先預處理逆元然后根據定義ln(A(x)) ∑_{k1} (-1)^{k-1} * (A(x)-1)^k / k。因為 A(x)-1 的常數項為0所以這個級數在模 x^n 意義下是有限的只需要算到 kn-1 即可。 用這個暴力程序去驗證你的 NTT 優化版本在小數據n5,6,7...下的結果是否一致。這是定位問題最有效的方式。6. 從模板到應用ln 在生成函數中的意義搞懂了實現我們再來聊聊它到底有什么用這樣下次遇到問題你才能想到用它。6.1 組合意義的連接集合與連通分量這是 ln 最經典的應用。假設我們有一個組合類 A比如所有的圖其指數生成函數EGF為 A(x)。那么A(x)的 expB(x) exp(A(x))通常代表了由 A 中的“連通”對象任意組合而成的“所有”對象。例如A(x) 是連通圖的 EGF那么 B(x) 就是所有圖的 EGF。反過來A(x)的 lnC(x) ln(B(x))就代表了從“所有”對象中提取出“連通”分量。例如已知所有圖的 EGF B(x)那么 ln(B(x)) 就是連通圖的 EGF。在很多計數問題中我們更容易求出所有方案的生成函數而想要得到連通方案的生成函數就需要對其取 ln。6.2 多項式牛頓迭代中的角色牛頓迭代是求解多項式方程F(G(x)) 0的強大工具。例如求exp、求sqrt開根、求復合逆函數等。 在推導這些迭代式時ln和exp經常作為一對互逆的運算出現。例如求G(x) exp(F(x))可以轉化為方程ln(G(x)) - F(x) 0然后應用牛頓迭代。此時poly_ln就成了迭代過程中必須調用的子程序。6.3 形式微分與形式積分的工具ln的公式本身完美結合了求導、求逆和積分。這使得它成為學習多項式形式運算的一個優秀案例。掌握了它你就掌握了處理形式冪級數的一整套基本工具鏈。7. 總結與擴展思考實現一個poly_ln就像搭樂高。你需要準備好“求導”、“求逆”其內部又需要“NTT乘法”、“積分”這幾個基礎模塊然后按照公式∫ (A / A) dx把它們正確地拼接起來。其中多項式求逆是最復雜、最容易出錯的一環務必理解其牛頓迭代的倍增思想并牢記迭代后清零的紀律。在競賽中poly_ln很少單獨出題它往往是更大問題的一塊拼圖。比如你需要先對某個生成函數取 ln進行一些操作再取 exp。因此將它寫對、寫熟封裝成一個可靠的函數是進軍更高級多項式算法如指數函數、三角函數、快速冪、復合逆的必經之路。最后關于常數項不為1的情況理論上可以通過提取公因式解決若 A(0) c ≠ 0則 ln(A(x)) ln(c) ln(A(x)/c)。但 ln(c) 在模意義下需要離散對數來求解這超出了普通多項式模板的范圍。所以模板題和常見應用都默認常數項為1。寫多項式代碼是對耐心和細心的雙重考驗。一個符號的錯誤、一次忘記的清零都可能導致調試數小時。但一旦你徹底征服了它那種對復雜算法了如指掌、對每一行代碼都充滿自信的感覺是無與倫比的。希望這篇超詳細的拆解能幫你少走些彎路真正把這塊硬骨頭啃下來。