多項式求逆元,即已知多項式$A(x)$,咱們須要找到一個多項式$A^{-1}(x)$算法
使得函數
$$A(x)A^{-1}(x)\equiv 1\pmod {x^n}$$ui
咱們稱多項式$A^{-1}(x)$爲多項式$A(x)$的逆元spa
在這裏${x^n}$是一個數,模${x^n}$的意義就是將$\ge n$的項都忽略掉.net
至於咱們爲何要模$x^n$,由於經過計算不難發現:除了僅有常數項的多項式的逆元爲一個常數以外,其他多項式的逆元均有無窮多項code
這裏介紹一種比較經常使用的$O(nlogn)$倍增算法,實際上許多與多項式有關的操做都須要用的倍增算法blog
假設咱們已經求出了多項式$A(x)$在模$x^{\frac{n}{2}}$意義下(就是一半)的逆元,設其爲$B(x)$ip
即get
$$A(x)B(x)\equiv 1\pmod {x^{\lceil\frac n2\rceil}}$$string
根據取模運算的性質,$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
在泰勒公式中,取$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)$
多項式的定義爲,給定多項式$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)$