目錄函數
若是有
\[ A(x) = \sum_{i = 0}^{n} c_i x^i \]優化
那麼因爲求導的加法原則,每一項能夠單獨考慮,則有
\[ A'(x) = \sum_{i = 1}^{n} ic_ix^{i - 1} \]ui
void Der(int *f, int *g, int lf) { For (i, 1, lf) g[i - 1] = 1ll * i * f[i] % Mod; g[lf] = 0; }
若是有
\[ A(x) = \sum_{i = 0}^{n} c_i x^i \]spa
一樣因爲積分的加法原則,每一項能夠單獨考慮,則有(常數項能夠忽略)
\[ \int A(x) \mathrm{d} x = \sum_{i = 1}^{n + 1} \frac{c_{i - 1}}{i} x^{i} \].net
void Int(int *f, int *g, int lf) { g[0] = 0; For (i, 1, lf + 1) g[i] = 1ll * f[i - 1] * inv[i] % Mod; }
定義兩個多項式 \(A(x) = \sum_{i = 0}^{n} a_i x^i, B(x) = \sum_{i = 0}^m b_i x^i\) 。那麼乘積爲 \(C(x) = \sum_{i = 0}^{n} \sum_{j = 0}^m a_ib_jx^{i + j}\) 。code
咱們考慮利用快速傅里葉變換解決,咱們求出 \(\omega_n^k\) 代入的全部的點值,而後用逆變換就能夠求出全部係數了。blog
對於模意義下的咱們利用 \(g^{\frac{p - 1}n}\) 代替 \(\omega\) (其中 \(g\) 爲原根)就能夠了。遞歸
咱們接下來都默認用費馬質數來做爲模數進行 \(\text{NTT}\) 。get
int rev[Maxn], W[Maxn], len; void NTT(int *P, int opt) { Rep (i, len) if (i < rev[i]) swap(P[i], P[rev[i]]); for (int i = 2, p; p = i >> 1, i <= len; i <<= 1) { W[0] = 1; W[1] = fpm(~opt ? g : invg, (Mod - 1) / i); For (k, 2, p - 1) W[k] = 1ll * W[k - 1] * W[1] % Mod; for (int j = 0; j < len; j += i) Rep (k, p) { int u = P[j + k], v = 1ll * P[j + k + p] * W[k] % Mod; P[j + k] = (u + v) % Mod; P[j + k + p] = (u - v + Mod) % Mod; } } if (!~opt) { int invn = fpm(len, Mod - 2); Rep (i, len) P[i] = 1ll * P[i] * invn % Mod; } } void Prepare(int lc) { int cnt = -1; for (len = 1; len <= lc; len <<= 1) ++ cnt; Rep (i, len) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << cnt); } int A[Maxn], B[Maxn], C[Maxn]; void Mult(int *a, int *b, int *c, int la, int lb) { int lc = la + lb; Prepare(lc); Rep (i, len) A[i] = i <= la ? a[i] : 0; NTT(A, 1); Rep (i, len) B[i] = i <= lb ? b[i] : 0; NTT(B, 1); Rep (i, len) C[i] = 1ll * A[i] * B[i] % Mod; NTT(C, -1); For (i, 0, lc) c[i] = C[i]; }
給出函數 \(f\) ,找到一個多項式 \(A\),使得 \(f(A(x)) \equiv 0 \pmod {x^n}\) 。博客
在 \(\pmod x\) 時,能夠求出 \(A(x)\) 的常數項。
已經求出 \(f(A(x)) \equiv 0 \bmod {x^n}\) ,如今求 \(f(B(x)) \equiv 0 \bmod {x^{2n}}\) ,\(B(x)\) 與 \(A(x)\) 前 \(n\) 項係數相同。
在 \(\bmod {x^{2n}}\) 下,在 \(A(x)\) 處展開 \(f(B(x))\) 。
\[ f(B(x)) = \sum_{i \ge 0} \frac{f^{(i)}(A(x))}{i!}(B(x) - A(x))^i \]
不難發如今 \(i \ge 2\) 的時候,最低項已經超過了 \(x^{2n}\) ,因此只有 \(i = 0, 1\) 要考慮,也就是
\[ \begin{aligned} f(B(x)) &\equiv f(A(x)) + f'(A(x))(B(x) - A(x)) &\pmod {x^{2n}}\\ B(x) &\equiv A(x) - \frac{f(A(x))}{f'(A(x))} &\pmod {x^{2n}} \end{aligned} \]
大多數狀況下,多項式複合是一個比較難求解的一個東西,但對於一些特殊的 \(f\) 能夠方便求出。
好像這個很差套牛迭,咱們現推算了。
設 \(F(x)\) 是 \(P(x)\) 在 \(\pmod {x^n}\) 下的逆元,那麼表示出來就是
\[ F(x)P(x) \equiv 1 \pmod {x^n} \]
當 \(n = 1\) 求逆元便可,假設咱們求出在 \(\bmod x^{\frac n 2}\) 意義下的逆元 \(G(x)\) ,那麼有
\[ G(x)P(x) \equiv 1 \pmod {x^\frac n 2}\\ \]
咱們有
\[ \begin{aligned} F(x) - G(x) &\equiv 0 \pmod {x^\frac n 2}\\ F^2(x) - 2F(x)G(x) + G^2(x) &\equiv 0 \pmod {x^n}\\ F(x) &\equiv 2G(x) - P(x)G^2(x) \pmod {x^n} \end{aligned} \]
那麼遞歸求解便可,複雜度 \(\displaystyle \mathcal T(n) = \mathcal T(\frac n2) + \mathcal O(n \log n) = \mathcal O(n \log n)\) 。
因爲這個是大部分多項式操做的一個基礎操做(加減乘除)之一,因此要儘可能加快速度。
能夠先 \(\text{NTT}\) 而後作點值乘法,能少一半常數。(注意一開始須要是 \(2^k \ge len\) 形式)
void Inv(int *f, int *g, int lf) { if (lf == 1) return void(g[0] = fpm(f[0], Mod - 2)); Inv(f, g, lf >> 1); Prepare(lf << 1); Rep (i, len) A[i] = i < lf ? f[i] : 0; NTT(A, 1); Rep (i, len) B[i] = i < lf ? g[i] : 0; NTT(B, 1); Rep (i, len) C[i] = 1ll * A[i] * B[i] % Mod * B[i] % Mod; NTT(C, -1); Rep (i, lf) g[i] = (g[i] * 2ll + Mod - C[i]) % Mod; }
求 \(F^2(x) \equiv P(x) \pmod {x^n}\) ,此時的 \(f(A(x)) = A^2 - P\) 咱們要求它的零點。
那麼咱們套用牛迭的式子就獲得了
\[ \begin{aligned} B(x) &= A(x) - \frac{A^2(x) - P(x)}{2A(x)}\\ &= \frac{A(x)}2 + \frac{P(x)}{2A(x)} \end{aligned} \]
複雜度 \(\displaystyle \mathcal T(n) = \mathcal T(\frac n2) + \mathcal O(n \log n) = \mathcal O(n \log n)\) 。 。。主定理是真的強大
注意常數項非 \(1\) 的時候,一般都會保證存在模意義下的二次剩餘,而後用 \(\text{BSGS}\) 求出它是 \(g^k\) ,而後開根就是 \(g^{k/2}\) 。
void Sqrt(int *f, int *g, int lf) { if (lf == 1) return void(g[0] = getsqrt(f[0])); Sqrt(f, g, lf >> 1); static int tmp[Maxn], tinv[Maxn]; Rep (i, lf) tmp[i] = 2 * g[i] % Mod; Inv(tmp, tinv, lf); Mult(f, tinv, tmp, lf, lf); const int inv2 = fpm(2, Mod - 2); Rep (i, lf) g[i] = (1ll * g[i] * inv2 % Mod + tmp[i]) % Mod; }
求 \(G(x) = \ln F(x)\) ,咱們有
\[ G(x) = \int \frac{F'(x)}{F(x)} \mathrm{d}x \]
那麼咱們求導而後再求逆便可,複雜度 \(\mathcal O(n \log n)\) 。
常數項須要是 \(1\) ,其餘形如 \(\ln x\) 的無理數在模意義下沒法表示。
void Ln(int *f, int *g, int lf) { static int der[Maxn], tmp[Maxn]; Der(f, der, lf); Inv(f, tmp, lf); Mult(der, tmp, tmp, lf, lf); Int(tmp, g, lf); }
能夠套用牛頓迭代和多項式 \(\ln\) 解決。
分治的作法雖然多了個 \(\mathcal O(\log)\) 但好寫許多(並且一般因爲常數優點會比較快,除非你會論文哥那個牛迭優化),原理以下
\[ \begin{aligned} B(x) &= \exp(A(x))\\ B'(x) &= A'(x) \exp(A(x))\\ B'(x) &= A'(x) B(x)\\ B(x) &= \int A'(x)B(x) \mathrm{d}x \end{aligned} \]
咱們考慮把 \(A(x)\) 先求導成 \(A'(x)\) ,令 \(mid = \lfloor \frac{l + r}2 \rfloor\) 而後考慮分治求出 \([l, mid]\) 的 \(B(x)\) ,而後考慮對於 \((mid, r]\) 的貢獻,也就是 \(A_{0 \sim r - l}\) 和 \(B_{l \sim mid}\) 卷而後貢獻上去。
注意 \(A(0) = 0\) 纔在模意義下有含義,此時 \(B(0) = 1\) ,以及積分須要平移一位並乘上 \(\frac{1}{n}\) 。複雜度 \(\mathcal O(n \log^2 n)\) 。
void Exp(int *a, int *b, int l, int r) { if (l == r) return void(b[l] = (l ? 1ll * b[l] * inv[l] % Mod : 1)); int mid = (l + r) >> 1; Exp(a, b, l, mid); static int tmp[Maxn]; Mult(a, b + l, tmp, r - l, mid - l); For (i, mid + 1, r) b[i] = (b[i] + tmp[i - l - 1]) % Mod; Exp(a, b, mid + 1, r); } void Exp(int *f, int *g, int lf) { static int der[Maxn]; Der(f, der, lf); der[lf] = 0; memset(g, 0, sizeof(int) * (lf + 1)); Exp(der, g, 0, lf); }
求 \(G(x) = F^k(x)\) ,咱們取對數有 \(\ln G(x) = k \ln F(x)\) 也就是 \(G(x) = \exp(k \ln F(x))\) 。
可是 \(F(0)\) 可能不爲 \(1\) ,那麼咱們把 \(F(x)\) 除掉 \(F(0)\) 便可。那麼還可能爲 \(0\) 咱們平移便可,也就是除掉 \(x^a\) 。
而後 \(\exp\) 結束記得乘回來它的 \(k\) 次冪便可。
不想判 \(0\) 了。。
void Pow(int *f, int *g, int lf, int power) { static int tf[Maxn], tmp[Maxn], invf; assert(f[0]); invf = fpm(f[0], Mod - 2); Rep (i, lf) tf[i] = 1ll * f[i] * invf % Mod; Ln(tf, tmp, lf); Rep (i, lf) tmp[i] = 1ll * tmp[i] * power % Mod; Exp(tmp, g, lf); invf = fpm(fpm(invf, Mod - 2), power); Rep (i, lf) g[i] = 1ll * g[i] * invf % Mod; }
設 \(G(F(x)) = x\) 求 \(F(x)\) 的一項。此時 \(G(x)\) 多是和 \(F\) 同階的一個多項式,套用牛頓迭代項數會特別多,根本沒法求解。
但要求單項即 \([x^n]F(x)\) 有一個能夠快速計算的式子
\[ [x^n] F(x) = \frac 1n [x^{n - 1}] (\frac{x}{G(x)})^n \]
至於證實能夠參看這篇博客 。
int Lagrange(int *f, int lf, int x) { static int tmp[Maxn], res[Maxn]; Rep (i, lf) res[i] = f[i + 1]; Inv(res, tmp, lf); Power(tmp, res, n); return 1ll * res[x - 1] * fpm(x, Mod - 2) % Mod; }