多項式求逆元
多項式求逆元，即已知多項式$A(x)$，我們需要找到一個多項式$A^{-1}(x)$
使得
$$A(x)A^{-1}(x)\equiv 1\pmod {x^n}$$
我們稱多項式$A^{-1}(x)$為多項式$A(x)$的逆元
在這里${x^n}$是一個數，模${x^n}$的意義就是將$\ge n$的項都忽略掉
至於我們為什么要模$x^n$，因為通過計算不難發現：除了僅有常數項的多項式的逆元為一個常數之外，其余多項式的逆元均有無窮多項
算法
這里介紹一種比較常用的$O(nlogn)$倍增算法，實際上許多與多項式有關的操作都需要用的倍增算法
假設我們已經求出了多項式$A(x)$在模$x^{\frac{n}{2}}$意義下(就是一半)的逆元，設其為$B(x)$
即
$$A(x)B(x)\equiv 1\pmod {x^{\lceil\frac n2\rceil}}$$
根據取模運算的性質，$A^{-1}(x)$前$\frac{n}{2}$項與$B(x)$是相同的
$$A(x)A^{-1}(x)\equiv 1\pmod {x^{\lceil\frac n2\rceil}}$$
合並兩個式子可以得到
$$A(x)(B(x)-A^{-1}(x))\equiv 0\pmod{x^{\lceil\frac n2\rceil}}$$
這里的$A(x)$不為零
$$B(x)-A^{-1}(x)\equiv 0\pmod{x^{\lceil\frac n2\rceil}}$$
然后兩邊平方后拆開
$$B^2(x)+A^{-2}(x)-2B(x)A^{-1}(x)\equiv 0\pmod {x^n}$$
兩邊同乘$A(x)$
$$B^2(x)A(x)+A^{-1}(x)-2B(x)\equiv 0\pmod {x^n}$$
移項
$$A^{-1}(x)\equiv B(x)(2-B(x)*A(x))\pmod {x^n}$$
這樣我們就得到了$A(x)$和$B(x)$的關系
利用NTT計算復雜度為$O(nlogn)$
代碼實現
一份常數還不錯的代碼
// luogu-judger-enable-o2 #include<cstdio> #include<cstring> #include<algorithm> #define getchar() (p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 1<<21, stdin), p1 == p2) ? EOF : *p1++) #define swap(x,y) (x ^= y, y ^= x, x ^= y) #define mul(a, b) 1ll * a * b % P #define add(a, b) (a + b >= P ? a + b - P : a + b) #define dec(a, b) (a - b < 0 ? a - b + P : a - b) #define rg register const int MAXN = (1 << 21) + 10, P = 998244353, G = 3, Gi = 332748118; char buf[1<<21], *p1 = buf, *p2 = buf; inline int read() { char c = getchar(); int x = 0, f = 1; while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();} while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar(); return x * f; } char obuf[1<<24], *O=obuf; void print(int x) { if(x > 9) print(x / 10); *O++= x % 10 + '0'; } inline int fastpow(int a, int k) { int base = 1; while(k) { if(k & 1) base = mul(a, base); a = mul(a, a); k >>= 1; } return base % P; } int N, r[MAXN], X[MAXN], Y[MAXN], A[MAXN], B[MAXN], Og[MAXN]; inline void NTT(int *A, int type, int len) { int limit = 1, L = 0; while(limit < len) limit <<= 1, L++; for(rg int i = 0; i < limit; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1)); for(rg int i = 0; i < limit; i++) if(i < r[i]) swap(A[i], A[r[i]]); for(rg int mid = 1; mid < limit; mid <<= 1) { int R = mid << 1; int W = fastpow(G, (P - 1) / R); Og[0] = 1; for(rg int j = 1; j < mid; j++) Og[j] = mul(Og[j - 1], W); for(rg int j = 0; j < limit; j += R) { for(rg int k = 0; k < mid; k++) { const int x = A[j + k], y = mul(Og[k], A[j + k + mid]); A[j + k] = add(x, y), A[j + k + mid] = dec(x, y); } } } if(type == -1) { std::reverse(&A[1], &A[limit]); for(int i = 0, inv = fastpow(len , P - 2); i < limit; i++) A[i] = 1ll * A[i] * inv % P; } } void Inv(int *a, int *b, int len) {// a要求的多項式 b逆元 len：要求的逆元的長度 if(len == 1) { b[0] = fastpow(a[0], P - 2); return ; } Inv(a, b, len >> 1); for(rg int i = 0; i < len; i++) A[i] = a[i], B[i] = b[i]; NTT(A, 1, len << 1); NTT(B, 1, len << 1); for(rg int i = 0; i < (len << 1); i++) A[i] = mul(mul(A[i], B[i]), B[i]) ; NTT(A, -1, len << 1); for(rg int i = 0; i < len; i++) b[i] = (1ll * (b[i] << 1) % P + P - A[i] ) % P; } int main() { #ifdef WIN32 freopen("a.in","r",stdin); #endif N = read(); for(int i = 0; i < N; i++) X[i] = (read() + P) % P; int Len; for(Len = 1; Len < N; Len <<= 1); Inv(X, Y, Len); for(int i = 0; i < N; i++) print(Y[i]), *O++ = ' '; fwrite(obuf, O-obuf, 1 , stdout); return 0; }
多項式除法
給定多項式$A(x)$,$B(x)$
我們需要找到多項式$D(x)$，$R(x)$，使得
$$A(x) = D(x)B(x) + R(x)$$
在這里$A(x)$為$N$次多項式，$B(x)$為$M$次多項式
那么$D(x)$為$N-M+1$次多項式，$R(x)$的次數$<M$，這樣我們可以保證解的唯一性
算法
我們考慮如何解決上面的問題
首先$R(x)$具體的值是不用考慮的，因為我們求出$D(x)$后可以把$D(x)$帶入從而求得$R(x)$
另外，根據我們求逆元的經驗，沒有模的多項式除法是有無窮多項的
這個其實也好解決，我們強制讓所有多項式$\pmod {x^{n - m +1}}$
那么接下來我們只需要解決余數，也就是$R(x)$的問題了
考慮到我們的模數為$x^{n - m +1}$，也就是說我們會丟棄所有多項式的前$n-m+1$項
但是$R(x)$是一個$M-1$次多項式，直接模肯定是消不掉的，我們考慮能不能讓它的系數乘上$x^{n-m+1}$還能保證要求的多項式跟原來多項式意義相同
這里，我們定義翻轉操作
$$A^R(x) = x^n A(\frac{1}{x}) $$
也就是將多項式的系數進行翻轉
下面是神仙推導
$$A(x)=B(x)D(x)+R(x)$$
將$x$用$\frac{1}{x}$替代，兩邊同乘$x^N$
$$\begin{eqnarray*} x^n A(\frac{1}{x}) &=& x^{n - m}D(\frac{1}{x}) x^mB(\frac{1}{x}) + x^{n - m + 1}x^{m - 1}R(\frac{1}{x}) \\ A^R(x) &=& D^R(x)B^R(x) + x^{n - m + 1}R^R(x) \end{eqnarray*}$$
然后通過$\pmod {x ^{n - m +1}}$就可以把$R(x)$消掉啦
得到
$$A^R(x)\equiv B^R(x)D^R(x)\pmod{x^{n-m+1}}$$
$$A^R(x) * B^{-R}(x)\equiv D^R(x)\pmod{x^{n-m+1}}$$
這樣的話我們求出$B(x)$的逆元，然后用NTT乘起來就好了
多項式取模
實際就是上面的$R(x)$
用多項式除法做就可以
代碼實現
// luogu-judger-enable-o2 // luogu-judger-enable-o2 #include<cstdio> #include<cstring> #include<algorithm> #define getchar() (p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 1<<21, stdin), p1 == p2) ? EOF : *p1++) #define swap(x,y) (x ^= y, y ^= x, x ^= y) #define mul(a, b) 1ll * a * b % P #define add(a, b) (a + b >= P ? a + b - P : a + b) #define dec(a, b) (a - b < 0 ? a - b + P : a - b) #define rg register using namespace std; const int MAXN = 1e6, P = 998244353, Gi = 3; char buf[1<<21], *p1 = buf, *p2 = buf; inline int read() { char c = getchar(); int x = 0, f = 1; while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();} while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar(); return x * f; } char obuf[1<<24], *O=obuf; void print(int x) { if(x > 9) print(x / 10); *O++= x % 10 + '0'; } inline int fastpow(int a, int k) { int base = 1; while(k) { if(k & 1) base = mul(a, base); a = mul(a, a); k >>= 1; } return base % P; } int N, r[MAXN], A[MAXN], B[MAXN], Og[MAXN], F[MAXN], Q[MAXN], G[MAXN], R[MAXN], Ginv[MAXN], Atmp[MAXN], Btmp[MAXN]; int GetLen(int x) { int len; for(len = 1; len <= x; len <<= 1); return len; } inline void NTT(int *A, int type, int len) { int limit = 1, L = 0; while(limit < len) limit <<= 1, L++; for(rg int i = 0; i < limit; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1)); for(rg int i = 0; i < limit; i++) if(i < r[i]) swap(A[i], A[r[i]]); for(rg int mid = 1; mid < limit; mid <<= 1) { int R = mid << 1; int W = fastpow(Gi, (P - 1) / R); Og[0] = 1; for(rg int j = 1; j < mid; j++) Og[j] = mul(Og[j - 1], W); for(rg int j = 0; j < limit; j += R) { for(rg int k = 0; k < mid; k++) { const int x = A[j + k], y = mul(Og[k], A[j + k + mid]); A[j + k] = add(x, y), A[j + k + mid] = dec(x, y); } } } if(type == -1) { std::reverse(&A[1], &A[limit]); for(int i = 0, inv = fastpow(len , P - 2); i < limit; i++) A[i] = 1ll * A[i] * inv % P; } } void Inv(int *a, int *b, int len) {// a要求的多項式 b逆元 len：要求的逆元的長度 if(len == 1) { b[0] = fastpow(a[0], P - 2); return ; } Inv(a, b, len >> 1); for(rg int i = 0; i < len; i++) A[i] = a[i], B[i] = b[i]; NTT(A, 1, len << 1); NTT(B, 1, len << 1); for(rg int i = 0; i < (len << 1); i++) A[i] = mul(mul(A[i], B[i]), B[i]) ; NTT(A, -1, len << 1); for(rg int i = 0; i < len; i++) b[i] = (1ll * (b[i] << 1) % P + P - A[i] ) % P; for(rg int i = 0; i < (len << 1); i++) A[i] = B[i] = 0; } void Mul(int *a, int *b, int *c, int N, int M) { int len = GetLen(max(N, M)) << 1; for(int i = 0; i <= N; i++) Atmp[i] = a[i]; for(int i = 0; i <= M; i++) Btmp[i] = b[i]; NTT(Atmp, 1, len); NTT(Btmp, 1, len); for(int i = 0; i <= len; i++) c[i] = 1ll * Atmp[i] * Btmp[i] % P, Atmp[i] = Btmp[i] = 0; NTT(c, -1, len); } int main() { #ifdef WIN32 freopen("a.in","r",stdin); #endif int N = read(), M = read(); for(int i = 0; i <= N; i++) F[i] = read(); for(int i = 0; i <= M; i++) G[i] = read(); reverse(F, F + N + 1); reverse(G, G + M + 1); //for(int i = 0; i <= N - M; i++) Ginv[i] = G[i];//tag Inv(G, Ginv, GetLen(N - M)); Mul(F, Ginv, Q, N - M, N - M); reverse(Q, Q + N - M + 1); for(int i = 0; i <= N - M; i++) print(Q[i]), *O++ = ' '; *O++ = '\n'; reverse(F, F + N + 1); reverse(G, G + M + 1); Mul(Q, G, R, N - M, M); for(int i = 0; i < M; i++) print(dec(F[i], R[i])), *O++ = ' '; fwrite(obuf, O-obuf, 1 , stdout); return 0; }
泰勒展開
泰勒公式是一個用函數在某點的信息描述其附近取值的公式。如果函數足夠平滑的話，在已知函數在某一點的各階導數值的情況之下，泰勒公式可以用這些導數值做系數構建一個多項式來近似函數在這一點的鄰域中的值。泰勒公式還給出了這個多項式和實際的函數值之間的偏差。
一句話應用：構造多項式來逼近函數
$$g(x)=g(0)+\frac{f^{1}(0)}{1!}x+\frac{f^{2}(0)}{2!}x^{2}+……+\frac{f^{n}(0)}{n!}x^{n}$$
$f^n$表示對$f$進行$n$次求導
這里的“多項式”我們可以直觀的理解為一種特殊的“函數”
普通牛頓迭代法
用途：求函數$f(x)$的零點
首先任取一個點$x_0$
然后對$f(x)$泰勒展開
$$f(x)=f(x_0)+f'(x_0)(x-x_0)+\frac {f''(x_0)(x-x_0)^2}2+\cdots$$
我們只保留線性部分
$$f(x)=f(x_0)+f'(x_0)(x-x0)=0$$
這樣的話
$$x=x0-\frac{f(x_0)}{f'(x_0)}$$
然后用$x$替換$x_0$，索然我們在保留線性部分時丟掉了很多信息，但是我們進行多次迭代仍然可以得到最終解
多項式牛頓迭代法
用途：
已知多項式$F(x)$和函數$G(x)$，求
$$G(F(x)) \equiv 0 \pmod {x^n}$$
算法
仍然考慮倍增
當$n=1$時，$F(x)$僅有一個常數項，
上面的式子可以化為
$$G(x) \equiv 0$$
我們可以直接利用牛頓迭代求解
按照多項式求逆的思路，假設我們已經求得
$$G(F_0(x)) \equiv 0 \pmod {x^{\lceil \frac{n}{2} \rceil}}$$
然后在$F_0(x)$處進行泰勒展開
$$\begin{eqnarray*} G(F(x)) & = & G(F_0(x)) \\ & + & \frac{G'(F_0(x))}{1!}\left(F(x) - F_0(x)\right) \\ & + & \frac{G''(F_0(x))}{2!}\left(F(x) - F_0(x)\right)^2 \\ & + & \cdots \end{eqnarray*}$$
在這里我們不難發現一個問題$F_0(x)$的前$\frac{n}{2}$項與$F(x)$是相同的
那么$F(x)-F_0(x)$的前$\frac{n}{2}$項必然是$0$
那么對於$(F(x)-F_0)^2$，前$n$項必然是$0$
同理當指數大於$2$時，前$n$項必然也是$0$
那么我們如果只保留整數部分，得到的應該是准確解！
$$G(F(x)) \equiv G(F_0(x)) + G'(F_0(x))\left(F(x) - F_0(x)\right) \pmod {x^n}$$
由於$$G(F(x)) \equiv 0 \pmod {x^n}$$
然后就做完了
$$F(x) \equiv F_0(x) - \frac{G(F_0(x))}{G'(F_0(x))} \pmod {x^n}$$
多項式開根
利用牛頓迭代法可以快速的推出多項式開根的做法
多項式開根即已知多項式$A(x)$，求多項式$B(x)$，滿足
$B^2(x) \equiv A(x) \pmod{x^n}$
設$F(x)$滿足
$F^2(x) - A(x) \equiv 0 \pmod{x^n}$
設$G(F(x)) = F(x)^2-A(x)$
我們要求的也就是$G(F(x))$的零點
根據求導公式，得到$G'(F(x)) = 2F(x)$(在這里我們認為$A(x)$為常數)
$$\begin{eqnarray*} F(x) & \equiv & F_0(x) - \frac{F_0^2(x) - A(x)}{2F_0(x)} \\ & \equiv & \frac{F_0^2(x) + A(x)}{2F_0(x)} \pmod {x^n} \end{eqnarray*}$$
$F_0$是在$\pmod {x^{\frac{n}{2}}}$意義下的結果
這里還有一個問題，就是當多項式僅有一個常數項的時候怎么處理
我們實際是要計算
$$\begin{equation}\label{quadratic}x^2 \equiv a \pmod p~~(p \nmid a)\end{equation}$$
這玩意兒的學名叫做二次剩余
可以看這里https://blog.csdn.net/a_crazy_czy/article/details/51959546
多項式對數與多項式$exp$
麥克勞林級數
在泰勒公式中，取$x=0$，得到的級數
$$\sum ^{\infty }_{n=0}\dfrac {f^{n}\left( 0\right) }{n!}x^{n}$$
稱為麥克勞林級數。
對數的計算
其實我並不知道這玩意兒的真正意義是什么，不過有大佬給出了它的定義，那就默認按定義來吧
多項式的對數可以認為是一個多項式和麥克勞林級數的復合
即對於多項式$A(x)$，有
$$\ln (1 - A(x)) = -\sum_{i \geq 1} \frac{A^i(x)}{i}$$
在計算時直接對$lnA(x)$微分后再積分
$$lnA(x)=\int (lnA(x))'=\int\frac{A'(x)}{A(x)}$$
這里用到了復合函數求導的鏈式法則
$f(g(x))=f’(g(x))g’(x)$
這樣在計算的時候可以先求導再積分，這兩者的復雜度都是$O(n)$的
加上多項式求逆的復雜度，總復雜度為$O(nlogn)$
多項式exp
多項式的定義為，給定多項式$A(x)$
$$\exp(A(x)) = e^{A(x)} = \sum_{i \geq 0} \frac{A^i(x)}{i!}$$
計算方法：
設$F(x) = e^{A(x)}$
$$lnF(x) = A(x)$$
設$G(F(x)) = lnF(x) - A(x) = 0$
這樣就轉換成了求多項式零點的問題，直接上牛頓迭代
$$\begin{eqnarray*} F(x) & \equiv & F_0(x) - \frac{G(F_0(x))}{G'(F_0(x))} \\ & \equiv & F_0(x) \left (1 - \ln F_0(x) + A(x) \right ) \pmod {x^n}\end{eqnarray*}$$
復雜度$O(nlogn)$
多項式冪函數
給定多項式$A(x)$和$k$，求
$$A^k(x)\pmod{x^n}$$
用上面的兩種方法轉換
$$A^k(x)\equiv e^{klnA(x)}\pmod{x^n}$$
時間復雜度$O(nlogn)$，常數巨大！
多項式三角函數
不會
$$e^{iA(x)}=cos(A(x))+isin(A(x))$$
多項式的快速差值與多點求值
此部分參(抄)考(襲)自http://blog.miskcoo.com/2015/05/polynomial-multipoint-eval-and-interpolation
多點求值
給出一個多項式$A(x)$以及$n$個點$x_0, x_1, x_2, \dots, x_{n-1}$
求出$A(x_0), A(x_1), A(x_2),\dots,x_{n-1}$，即每個點在多項式中的取值
快速差值
給出一個集合$$X = \{(x_i, y_i) : 0 \leq i \leq n\}$$
求一個多項式$A(x)$，滿足
$$\forall (x, y) \in X, A(x) = y$$
算法
因為這兩個問題的特殊性，因此在計算過程中可能會用到彼此，大家直接略過就好
多點求值
首先把要求的值分成兩半
$$\begin{eqnarray*} X^{[0]} &=& \{x_0, x_1, \cdots, x_{\lfloor \frac{n}{2} \rfloor}\} \\ X^{[1]} &=& \{x_{\lfloor \frac{n}{2} \rfloor+1},x_{\lfloor \frac{n}{2} \rfloor+2}, \cdots, x_{n-1}\} \end{eqnarray*}$$
我們拿前一半來講， 后一半同理
記用$x^[0]$插值得到的$\frac{n}{2}$次多項式為$A^[0]$
構造多項式
$$\begin{eqnarray*} P^{[0]}(x) &=& \prod_{i=0}^{\lfloor \frac{n}{2} \rfloor} (x-x_i) \end{eqnarray*}$$
對於該多項式來說，當$x \in x^[0]$時，$P^{[0]}(x) = 0$，那么$A(x)$可以表示為
$$A(x) = D(x)P^{[0]}(x) + A^{[0]}(x)$$
這個式子等價於
$$A^{[0]}(x) \equiv A(x) \pmod {P^{[0]}(x)}$$
用多項式取模算就可以
時間復雜度：$O(nlog^n)$
快速插值
$O(N^2)$：
$${A(x)=\sum_{i=1}^n{\frac{\prod_{{j}\neq{i}}{({x}-{x_j})}}{\prod_{{j}\neq{i}}{({x_i}-{x_j})}}{y_i}}}$$
首先按下標分類
$$\begin{eqnarray*} X^{[0]} &=& \{(x_i, y_i) : 0 \leq i \leq \lfloor \frac{n}{2} \rfloor \} \\ X^{[1]} &=& \{(x_i, y_i) : \lfloor \frac{n}{2} \rfloor < i \leq n\} \end{eqnarray*}$$
現在假設已經求出了$X^{[0]}$的插值多項式$A^{[0]}(x)$，
設
$$P(x) = \prod_{i=0}^{\lfloor \frac{n}{2} \rfloor} (x-x_i)$$
$$A(x) = A^{[1]}(x)P(x) + A^{[0]}(x)$$
其中 $A^{[1]}(x)$ 是一個未知的多項式，由於 $\forall (x, y) \in X^{[0]}, P(x) = 0, A^{[0]}(x) = y$
這樣的話就滿足 $X^{[0]}$的點都在$A(x)$上，問題就變成要將$X^{[1]}$ 內的點插值，使得
$\forall (x_i, y_i) \in X^{[1]}, y_i = A^{[1]}(x_i)P(x_i) + A^{[0]}(x_i)$化簡之后得到
$$A^{[1]}(x_i) = \frac{A^{[0]}(x_i)-y_i}{P(x_i)}$$
這樣就得到了新的待插值點，利用同樣的方法求出插值出$A^{[1]}$然后合並就可以了
由於每一次都要多點求值求出$X^{[1]}$內$P(x)$和$A^{[0]}(x)$的值，最終復雜度是
$T(n) = 2T(\frac{n}{2}) + \mathcal O(n \log^2 n) = \mathcal O(n \log^3 n)$