网站建设厘金手指下拉15,物流百度推广怎么做网站,南阳网站建设报价,公司要想做个网站这么弄快速傅里叶变换(FFT)
多项式表示
系数表示法#xff1a;
一个nnn次多项式可以用n1n 1n1个系数表示出来#xff1a;f(x)a0a1xa2x2⋯an−1xn−1anxnf(x) a_0 a_1 x a_2 x ^ 2 \dots a_{n - 1} x ^{n- 1} a_n x ^nf(x)a0a1xa2x2⋯an−1xn−1anxn。
点值表示…快速傅里叶变换(FFT)
多项式表示
系数表示法
一个nnn次多项式可以用n1n 1n1个系数表示出来f(x)a0a1xa2x2⋯an−1xn−1anxnf(x) a_0 a_1 x a_2 x ^ 2 \dots a_{n - 1} x ^{n- 1} a_n x ^nf(x)a0a1xa2x2⋯an−1xn−1anxn。
点值表示法
通过线性代数高斯消元我们可以知道一个nnn次多项式可以通过n1n 1n1个点联立方程组解得
f(x){(x0,f(x0),(x1,f(x1)),(x2,f(x2)))…(xn−1,f(xn−1)),(xn,f(xn))}f(x) \{(x_0, f(x_0), (x_1, f(x_1)), (x_2, f(x_2))) \dots (x_{n - 1}, f(x_{n - 1})), (x_n, f(x_n)) \}f(x){(x0,f(x0),(x1,f(x1)),(x2,f(x2)))…(xn−1,f(xn−1)),(xn,f(xn))}。
有离散傅里叶变换DFTDFTDFT把一个多项式从系数表示变成点值表示IDFTIDFTIDFT把一个多项式从点值表示变成系数表示
而FFTFFTFFT就是通过选取某些特殊xxx点来加速DFT,IDFTDFT, IDFTDFT,IDFT的一种方法。
点值表示法的多项式相乘
F(x)f(x)g(x)F(x) f(x) g(x)F(x)f(x)g(x)
F(x){(x0,f(x0)g(x0)),(x1,f(x1)g(x1)),(x2,f(x2)g(x2))…(xn−1,f(xn−1)g(xn−1)),(xn,f(xn)g(xn))}F(x) \{(x_0, f(x_0)g(x_0)), (x_1, f(x_1)g(x_1)), (x_2, f(x_2)g(x_2)) \dots (x_{n - 1}, f(x_{n - 1})g(x_{n - 1})), (x_n, f(x_n)g(x_n)) \}F(x){(x0,f(x0)g(x0)),(x1,f(x1)g(x1)),(x2,f(x2)g(x2))…(xn−1,f(xn−1)g(xn−1)),(xn,f(xn)g(xn))}
由此我们想要得到两个多项式相乘的系数只需要先对两个多项式进行DFTDFTDFT然后对应的点值相乘再做一次IDFTIDFTIDFT即可求得系数。
引入复数
两个复数相乘的结果为模长相乘辐角相加证明如下
有两复数A(acosθ1,asinθ1i),B(bcosθ2,bsinθ2i)A(a \cos \theta_1, a \sin \theta_1 i), B(b \cos \theta_2, b \sin \theta_2i)A(acosθ1,asinθ1i),B(bcosθ2,bsinθ2i)用极角 模长来表示。
两复数相乘有A×B(ab(cosθ1cosθ2−sinθ1sinθ2),ab(sinθ1cosθ2cosθ1sinθ2)i)(abcos(θ1θ2),absin(θ1θ2)i)A \times B (ab(\cos \theta_1 \cos \theta_2 - \sin \theta_1 \sin \theta_2), ab(\sin \theta_1 \cos \theta_2 \cos \theta_1 \sin \theta_2) i) (ab \cos (\theta_1 \theta_2), ab \sin (\theta_1 \theta_2) i)A×B(ab(cosθ1cosθ2−sinθ1sinθ2),ab(sinθ1cosθ2cosθ1sinθ2)i)(abcos(θ1θ2),absin(θ1θ2)i)。
引入nnn次复根即xn1x ^ n 1xn1这样的解显然有nnn个设wnie2πniw_{n} ^{i} e ^{\frac{2 \pi}{n} i}wnien2πi在复平面内即是把一个圆分成了nnn等份。
我们取nnn等分的第一个交所对应的向量wncos(2πn)sin(2πn)iw_n \cos(\frac{2 \pi}{n}) \sin(\frac{2 \pi}{n}) iwncos(n2π)sin(n2π)i则其他复根都可用wnw_nwn的iii次幂来表示。
快速傅里叶变换
考虑如何分治求解
f(x)a0a1xa2x2a3x3a4x4a5x5a6x6a7x7f(x) a_0 a_1 x a_2 x ^ 2 a_3 x ^ 3 a_4 x ^ 4 a_5 x ^ 5 a_6 x ^ 6 a_ 7 x ^ 7f(x)a0a1xa2x2a3x3a4x4a5x5a6x6a7x7
按照xxx的次幂分奇偶再在右边提出一个xxx
f(x)(a0a2x2a4x4a6x6)x(a1a3x2a5x4a7x6)f(x) (a_0 a_2 x ^ 2 a_4 x ^ 4 a_6 x ^ 6) x (a_1 a_3 x ^ 2 a_5 x ^ 4 a_7 x ^ 6)f(x)(a0a2x2a4x4a6x6)x(a1a3x2a5x4a7x6)
G(x)a0a2xa4x2a6x3G(x) a_0 a_2 x a_4 x ^ 2 a_6 x ^ 3G(x)a0a2xa4x2a6x3
H(x)a1a3xa5x2a7x3H(x) a_1 a_3 x a_5 x ^ 2 a_7 x ^ 3H(x)a1a3xa5x2a7x3
有f(x)G(x2)xH(x2)f(x) G(x ^ 2) x H(x ^ 2)f(x)G(x2)xH(x2)
由单位复根有
DFT(f(wnk))DFT(G(wn2k))wnkDFT(H(wn2k))DFT(G(wn2k))wnkDFT(H(wn2k))DFT(f(w _n ^ k)) DFT(G(w _n ^{2k})) w_n ^ k DFT(H(w _n ^{2k})) DFT(G(w _{\frac{n}{2}} ^ k)) w_n ^ k DFT(H(w_{\frac{n}{2}} ^ k))DFT(f(wnk))DFT(G(wn2k))wnkDFT(H(wn2k))DFT(G(w2nk))wnkDFT(H(w2nk))
DFT(f(wnkn2))DFT(G(wn2k))−wnkDFT(H(wn2k))DFT(f(w _n ^{k \frac{n}{2}})) DFT(G(w_{\frac{n}{2}} ^ k)) - w_n ^ k DFT(H(w _{\frac{n}{2}} ^ k))DFT(f(wnk2n))DFT(G(w2nk))−wnkDFT(H(w2nk))
由此求出DFT(G(wn2k))DFT(G(w _{\frac{n}{2}} ^k))DFT(G(w2nk))和DFT(H(wn2k))DFT(H(w_{\frac{n}{2}} ^ k))DFT(H(w2nk))即可知DFT(f(wnk)),DFT(f(wnkn2))DFT(f(w _n ^ k)), DFT(f(w _n ^{k \frac{n}{2}}))DFT(f(wnk)),DFT(f(wnk2n))然后对G,HG, HG,H再分别递归求解即可。
快速傅里叶逆变换
把单位复根值代入多项式得到的是如下结果
[y0y1y2⋮yn−2yn−1]\left[ \begin{matrix} y_0\\ y_1\\ y_2\\ \vdots\\ y_{n - 2}\\ y_{n - 1}\\ \end{matrix} \right]⎣⎢⎢⎢⎢⎢⎢⎢⎡y0y1y2⋮yn−2yn−1⎦⎥⎥⎥⎥⎥⎥⎥⎤ [111⋯111wn1wn2⋯wnn−2wnn−11wn2wn4⋯wn2(n−2)wn2(n−1)⋮⋮⋮⋱⋮⋮1wnn−2wn2(n−2)⋯wn(n−2)(n−2)wn(n−1)(n−2)1wnn−1wn2(n−1)⋯wn(n−2)(n−1)wn(n−1)(n−1)]\left[ \begin{matrix} 1 1 1 \cdots 1 1\\ 1 w_n ^ 1 w_n ^ 2 \cdots w_n ^ {n - 2} w_n ^{n - 1}\\ 1 w_n ^ 2 w_n ^ 4 \cdots w_n ^ {2(n - 2)} w_n ^{2(n - 1)}\\ \vdots \vdots \vdots \ddots \vdots \vdots\\ 1 w_n ^{n - 2} w_n ^ {2(n - 2)} \cdots w_n ^{(n - 2)(n - 2)} w_n ^{(n - 1)(n - 2)}\\ 1 w_n ^{n - 1} w_n ^ {2(n - 1)} \cdots w_n ^{(n - 2)(n - 1)} w_n ^{(n - 1)(n - 1)}\\ \end{matrix} \right]⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡111⋮111wn1wn2⋮wnn−2wnn−11wn2wn4⋮wn2(n−2)wn2(n−1)⋯⋯⋯⋱⋯⋯1wnn−2wn2(n−2)⋮wn(n−2)(n−2)wn(n−2)(n−1)1wnn−1wn2(n−1)⋮wn(n−1)(n−2)wn(n−1)(n−1)⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤ [a0a1a2⋮an−2an−1]\left[ \begin{matrix} a_0\\ a_1\\ a_2\\ \vdots\\ a_{n - 2}\\ a_{n - 1}\\ \end{matrix} \right]⎣⎢⎢⎢⎢⎢⎢⎢⎡a0a1a2⋮an−2an−1⎦⎥⎥⎥⎥⎥⎥⎥⎤
经过DFTDFTDFT我们已经得到了左边的矩阵考虑如何变换得到右边的系数矩阵线性代数我们知道只要在左边乘上一个中间大矩阵的逆我们即可得到右边的系数矩阵。
由于这个矩阵的元素非常特殊他的逆矩阵也有特殊的性质就是每一项取倒数再除以nnn就能得到他的逆矩阵。
每一项取倒数有1wnwn−1e−2πincos(2πn)isin(−2πn)\frac{1}{w_n} w_{n} ^{-1} e ^{-\frac{2 \pi i}{n}} \cos(\frac{2 \pi}{n}) i \sin (- \frac{2 \pi}{n})wn1wn−1e−n2πicos(n2π)isin(−n2π)所以我们只要将这个代入做一次DFTDFTDFT也就是IDFTIDFTIDFT最后再对整体除以nnn即可得到系数矩阵。
对以上进行证明
f(x)∑i0n−1aixif(x) \sum\limits_{i 0} ^{n - 1} a_i x ^ if(x)i0∑n−1aixiyif(wni)y_i f(w_n ^ i)yif(wni)构造A(x)∑i0n−1yixiA(x) \sum\limits_{i 0} ^{n - 1}y_i x ^ iA(x)i0∑n−1yixi将biwn−ib_i w_{n} ^{-i}biwn−i代入多项式A(x)A(x)A(x)
有A(bk)∑i0n−1yiwn−ik∑i0n1wn−ik∑j0n−1ajwnij∑j0n−1aj∑i0n−1(wnj−k)iA(b_k) \sum\limits_{i 0} ^{n - 1} y_i w_n ^{-ik} \sum\limits_{i 0} ^{n 1}w_n ^{-ik} \sum\limits_{j 0} ^{n - 1} a_j w_{n} ^{ij} \sum\limits_{j 0} ^{n - 1} a_j\sum\limits_{i 0} ^{n - 1} (w_{n} ^{j - k}) ^ iA(bk)i0∑n−1yiwn−iki0∑n1wn−ikj0∑n−1ajwnijj0∑n−1aji0∑n−1(wnj−k)i
令S(wna)∑i0n−1(wna)iS(w_{n} ^ a) \sum\limits _{i 0} ^{n - 1} (w_{n} ^{a}) ^ iS(wna)i0∑n−1(wna)i
显然有a0a 0a0S(wna)nS(w_n ^ a) nS(wna)n
a≠0a \neq 0a0时我们取S(wna),wnaS(wna)S(w_n ^ a),w_n ^ a S(w_n ^ a)S(wna),wnaS(wna)两者相减除以一个系数有S(wna)∑i1n(wna)i−∑i0n−1(wna)iwna−1(wna)n−(wna)0wna−10S(w_n ^ a) \frac{\sum\limits_{i 1} ^{n} (w_n ^ a) ^ i - \sum\limits_{i 0} ^{n - 1} (w_n ^ a) ^ i}{w_n ^ a - 1} \frac{(w_n ^ a) ^ n - (w_n ^ a) ^ 0}{w_n ^ a - 1} 0S(wna)wna−1i1∑n(wna)i−i0∑n−1(wna)iwna−1(wna)n−(wna)00
所以有S(wna)[a1]S(w_n ^ a) [a 1]S(wna)[a1]
A(bk)ak×nA(b_k) a_k \times nA(bk)ak×n
仔细想想这个证明就是把我们DFTDFTDFT过程中得到的点值作为系数去做一遍DFTDFTDFT得到的也就是A(x)A(x)A(x)的点值表达式同时对其除以nnn也就是f(x)f(x)f(x)的系数表达式了。
如何优化蝴蝶变换
分治过程中考虑系数如何变换 {a0,a1,a2,a3,a4,a5,a6,a7}{a0,a2,a4,a6}{a1,a3,a5,a7}{a0,a4}{a2,a6}{a1,a5}{a3,a7}{a0}{a4}{a2}{a6}{a1}{a5}{a3}{a7}\{a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7\}\\ \{a_0, a_2, a_4, a_6\}\{a_1, a_3, a_5, a_7\}\\ \{a_0, a_4\}\{a_2, a_6\}\{a_1, a_5\}\{a_3, a_7\}\\ \{a_0\}\{a_4\}\{a_2\}\{a_6\}\{a_1\}\{a_5\}\{a_3\}\{a_7\}\\ {a0,a1,a2,a3,a4,a5,a6,a7}{a0,a2,a4,a6}{a1,a3,a5,a7}{a0,a4}{a2,a6}{a1,a5}{a3,a7}{a0}{a4}{a2}{a6}{a1}{a5}{a3}{a7}
这个过程中有一个规律例如10011 0011001倒置后变成了100100100444也即是最后a1a_1a1所在的位置。
r[i]r[i]r[i]表示iii翻转之后的数字考虑如何从小到大递推得到r[i]r[i]r[i]有r[0]0r[0] 0r[0]0当我们在求xxx时先考虑除个位数以外的数就是r[x1]1r[x 1] 1r[x1]1了如果个位是111则加上lim1lim 1lim1就有了如下代码
void change(int lim) {for (int i 0; i lim; i) {r[i] (i 1) * (lim 1) (r[i 1] 1);}
}真正可用的FFTFFTFFT代码
P3803 【模板】多项式乘法FFT
#include bits/stdc.husing namespace std;struct Complex {double r, i;Complex(double _r 0, double _i 0) : r(_r), i(_i) {}
};Complex operator (const Complex a, const Complex b) {return Complex(a.r b.r, a.i b.i);
}Complex operator - (const Complex a, const Complex b) {return Complex(a.r - b.r, a.i - b.i);
}Complex operator * (const Complex a, const Complex b) {return Complex(a.r * b.r - a.i * b.i, a.r * b.i a.i * b.r);
}Complex operator / (const Complex a, const Complex b) {return Complex((a.r * b.r a.i * b.i) / (b.r * b.r b.i * b.i), (a.i * b.r - a.r * b.i) / (b.r * b.r b.i * b.i));
}typedef long long ll;const int N 5e6 10;int r[N], n, m;Complex a[N], b[N];void get_r(int lim) {for (int i 0; i lim; i) {r[i] (i 1) * (lim 1) (r[i 1] 1);}
}void FFT(Complex *f, int lim, int rev) {for (int i 0; i lim; i) {if (i r[i]) {swap(f[i], f[r[i]]);}}const double pi acos(-1.0);for (int mid 1; mid lim; mid 1) {Complex wn Complex(cos(pi / mid), rev * sin(pi / mid));for (int len mid 1, cur 0; cur lim; cur len) {Complex w Complex(1, 0);for (int k 0; k mid; k, w w * wn) {Complex x f[cur k], y w * f[cur mid k];f[cur k] x y, f[cur mid k] x - y;}}}if (rev -1) {for (int i 0; i lim; i) {a[i].r / lim;}}
}int main() {// freopen(in.txt, r, stdin);// freopen(out.txt, w, stdout);scanf(%d %d, n, m);n 1, m 1;for (int i 0; i n; i) {scanf(%lf, a[i].r);}for (int i 0; i m; i) {scanf(%lf, b[i].r);}int lim 1;while (lim n m) {lim 1;}get_r(lim);FFT(a, lim, 1);FFT(b, lim, 1);for (int i 0; i lim; i) {a[i] a[i] * b[i];}FFT(a, lim, -1);for (int i 0; i n m - 1; i) {printf(%lld , ll(a[i].r 0.5));}puts();return 0;
}快速数论变换(NTT)
#include bits/stdc.husing namespace std;const int N 5e6 10, mod 998244353;int a[N], b[N], r[N], n, m;int quick_pow(int a, int n) {int ans 1;while (n) {if (n 1) {ans 1ll * ans * n % mod;}a 1ll * a * a % mod;n 1;}return ans;
}void get_r(int lim) {for (int i 0; i lim; i) {r[i] (i 1) * (lim 1) (r[i 1] 1);}
}void NTT(int *f, int lim, int rev) {for (int i 0; i lim; i) {if (i r[i]) {swap(f[i], f[r[i]]);}}for (int mid 1; mid lim; mid 1) {int wn quick_pow(3, (mod - 1) / (mid 1));for (int len mid 1, cur 0; cur lim; cur len) {int w 1;for (int k 0; k mid; k, w 1ll * w * wn % mod) {int x f[cur k], y 1ll * w * f[cur mid k] % mod;f[cur k] (x y) % mod, f[cur mid k] (x - y mod) % mod;}}}if (rev -1) {int inv quick_pow(lim, mod - 2);reverse(f 1, f lim);for (int i 0; i lim; i) {f[i] 1ll * f[i] * inv % mod;}}
}int main() {// freopen(in.txt, r, stdin);// freopen(out.txt, w, stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);return 0;
}多项式求逆
f(x)g(x)≡1(modxn)f(x) g(x) \equiv 1 \pmod {x ^ n}f(x)g(x)≡1(modxn)称f(x)f(x)f(x)为g(x)g(x)g(x)或者g(x)g(x)g(x)为f(x)f(x)f(x)膜xnx ^ nxn意义下的逆元。
下面我们讨论给定f(x)f(x)f(x)求其逆f−1(x)f ^{-1}(x)f−1(x)。
倍增求解
假设我们已经求得f(x)f(x)f(x)膜x⌈n2⌉x ^{\lceil \frac{n}{2}} \rceilx⌈2n⌉下的逆元f0−1(x)f_0 ^{-1} (x)f0−1(x)要求f−1(x)f ^{-1}(x)f−1(x)即膜xnx ^{n}xn下的逆元则
f(x)f0−1(x)≡1(modx⌈n2⌉)f(x) f_0 ^{-1}(x) \equiv 1 \pmod{x ^{\lceil\frac{n}{2}\rceil} }f(x)f0−1(x)≡1(modx⌈2n⌉)
显然f(x)f−1(x)≡1(modx⌈n2⌉)f(x) f^{-1}(x) \equiv 1 \pmod {x ^{ \lceil\frac{n}{2}\rceil}}f(x)f−1(x)≡1(modx⌈2n⌉)也是成立的
对两边同时乘以f0−1(x)f_0 ^{-1}(x)f0−1(x)并移项有
f−1(x)−f0−1(x)≡0(modx⌈n2⌉)f ^{-1}(x) - f_0 ^{-1}(x) \equiv 0 \pmod{x ^{\lceil\frac{n}{2}\rceil}}f−1(x)−f0−1(x)≡0(modx⌈2n⌉)
对两边同时开方得到
f−2(x)−2f−1f0−1(x)f0−2(x)≡0(modxn)f ^{-2}(x) - 2 f^{-1} f_0 ^{-1}(x) f_0 ^{-2}(x) \equiv 0 \pmod {x ^n}f−2(x)−2f−1f0−1(x)f0−2(x)≡0(modxn)
我们再对两边乘上一个f(x)f(x)f(x)则有
f−1(x)−2f0−1f(x)f0−2(x)≡0(modxn)f ^{-1}(x) - 2 f_0 ^{-1} f(x) f_0 ^{-2}(x) \equiv 0 \pmod{x ^n}f−1(x)−2f0−1f(x)f0−2(x)≡0(modxn)
再对其进行移项可得
f−1(x)≡f0−1(x)(2−f(x)f0−1(x))(modxn)f ^{-1}(x) \equiv f_0 ^{-1}(x)\left( 2 - f(x) f_0 ^{-1}(x) \right) \pmod {x ^n}f−1(x)≡f0−1(x)(2−f(x)f0−1(x))(modxn)
由此我们递归求解即可。
#include bits/stdc.husing namespace std;const int N 1e6 10, mod 998244353;int a[N], b[N], c[N], r[N], n;int quick_pow(int a, int n) {int ans 1;while (n) {if (n 1) {ans 1ll * ans * a % mod;}a 1ll * a * a % mod;n 1;}return ans;
}void get_r(int lim) {for (int i 0; i lim; i) {r[i] ((i 1) * (lim 1)) (r[i 1] 1);}
}void NTT(int *f, int lim, int rev) {for (int i 0; i lim; i) {if (i r[i]) {swap(f[i], f[r[i]]);}}for (int mid 1; mid lim; mid 1) {int wn quick_pow(3, (mod - 1) / (mid 1));for (int len mid 1, cur 0; cur lim; cur len) {int w 1;for (int k 0; k mid; k, w 1ll * w * wn % mod) {int x f[cur k], y 1ll * w * f[cur mid k] % mod;f[cur k] (x y) % mod, f[cur mid k] (x - y mod) % mod;}}}if (rev -1) {int inv quick_pow(lim, mod - 2);reverse(f 1, f lim);for (int i 0; i lim; i) {f[i] 1ll * inv * f[i] % mod;}}
}void polyinv(int *a, int *b, int n) {if (n 1) {b[0] quick_pow(a[0], mod - 2);return ;}polyinv(a, b, n 1 1);int lim 1;while (lim 2 * n) {lim 1;}get_r(lim);for (int i 0; i n; i) {c[i] a[i];}for (int i n; i lim; i) {c[i] 0;}NTT(b, lim, 1);NTT(c, lim, 1);for (int i 0; i lim; i) {int cur (2 - 1ll * c[i] * b[i] % mod mod) % mod;b[i] 1ll * b[i] * cur % mod;}NTT(b, lim, -1);for (int i n; i lim; i) {b[i] 0;}
}int main() {// freopen(in.txt, r, stdin);// freopen(out.txt, w, stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);scanf(%d, n);for (int i 0; i n; i) {scanf(%d, a[i]);}polyinv(a, b, n);for (int i 0; i n; i) {printf(%d%c, b[i], i 1 n ? \n : );}return 0;
}多项式开根
给定多项式g(x)g(x)g(x)求f(x)f(x)f(x)满足f2(x)g(x)f ^ 2(x) g(x)f2(x)g(x)。
假设我们已经得到了g(x)g(x)g(x)膜x⌈n2⌉x ^{\lceil \frac{n}{2} \rceil}x⌈2n⌉下的根f0(x)f_0 (x)f0(x)要求膜xnx ^ nxn下的根f(x)f(x)f(x)
有f02(x)≡g(x)(modx⌈n2⌉)f_0 ^2(x) \equiv g(x) \pmod {x ^{\lceil \frac{n}{2} \rceil}}f02(x)≡g(x)(modx⌈2n⌉)
移项再开方有(f02(x)−g(x))2≡0(modxn)\left(f_0 ^2(x) - g(x) \right) ^ 2 \equiv 0 \pmod {x ^ n}(f02(x)−g(x))2≡0(modxn)
则(f02(x)g(x))2≡4f02(x)g(x)(modxn)\left( f_0 ^ 2(x) g(x) \right) ^ 2 \equiv 4 f_0 ^ 2(x) g(x) \pmod {x ^ n}(f02(x)g(x))2≡4f02(x)g(x)(modxn)
g(x)≡(f02(x)g(x)2f0(x))2(modxn)g(x) \equiv \left(\frac{f_0 ^ 2(x) g(x)}{2f_0 (x)} \right) ^ 2 \pmod {x ^ n}g(x)≡(2f0(x)f02(x)g(x))2(modxn)
所以f(x)≡f02(x)g(x)2f0(x)(modxn)f(x) \equiv \frac{f_0 ^ 2(x) g(x)}{2f_0(x)} \pmod {x ^ n}f(x)≡2f0(x)f02(x)g(x)(modxn)。
所以有f(x)≡2−1f0(x)2−1f0−1(x)g(x)(modxn)f(x) \equiv 2 ^{-1} f_0 (x) 2 ^{-1} f_0 ^{-1}(x) g(x) \pmod {x ^ n}f(x)≡2−1f0(x)2−1f0−1(x)g(x)(modxn)
对于g(0)1g(0) 1g(0)1的特殊情况
#include bits/stdc.husing namespace std;typedef long long ll;const int N 5e6 10, mod 998244353, inv2 mod 1 1;int a[N], b[N], c[N], d[N], r[N];int quick_pow(int a, int n) {int ans 1;while (n) {if (n 1) {ans 1ll * ans * a % mod;}a 1ll * a * a % mod;n 1;}return ans;
}void get_r(int lim) {for (int i 0; i lim; i) {r[i] (i 1) * (lim 1) (r[i 1] 1);}
}void NTT(int *f, int lim, int rev) {for (int i 0; i lim; i) {if (i r[i]) {swap(f[i], f[r[i]]);}}for (int mid 1; mid lim; mid 1) {int wn quick_pow(3, (mod - 1) / (mid 1));for (int len mid 1, cur 0; cur lim; cur len) {int w 1;for (int k 0; k mid; k, w 1ll * w * wn % mod) {int x f[cur k], y 1ll * w * f[cur mid k] % mod;f[cur k] (x y) % mod, f[cur mid k] (x - y mod) % mod;}}}if (rev -1) {int inv quick_pow(lim, mod - 2);reverse(f 1, f lim);for (int i 0; i lim; i) {f[i] 1ll * f[i] * inv % mod;}}
}void polyinv1(int *a, int *b, int n) {if (n 1) {b[0] quick_pow(a[0], mod - 2);return ;}polyinv1(a, b, n 1 1);int lim 1;while (lim 2 * n) {lim 1;}get_r(lim);for (int i 0; i n; i) {c[i] a[i];}for (int i n; i lim; i) {c[i] 0;}NTT(b, lim, 1);NTT(c, lim, 1);for (int i 0; i lim; i) {int cur (2 - 1ll * c[i] * b[i] % mod mod) % mod;b[i] 1ll * b[i] * cur % mod;}NTT(b, lim, -1);for (int i n; i lim; i) {b[i] 0;}
}void polysqrt(int *a, int *b, int n) {if (n 1) {b[0] 1;return ;}polysqrt(a, b, n 1 1);int lim 1;while (lim 2 * n) {lim 1;}for (int i 0; i lim; i) {d[i] 0;}polyinv1(b, d, n);for (int i 0; i n; i) {c[i] a[i];}for (int i n; i lim; i) {c[i] 0;}get_r(lim);NTT(b, lim, 1);NTT(c, lim, 1);NTT(d, lim, 1);for (int i 0; i lim; i) {b[i] (1ll * inv2 * b[i] % mod 1ll * inv2 * d[i] % mod * c[i] % mod) % mod; }NTT(b, lim, -1);for (int i n; i lim; i) {b[i] 0;}
}int main() {// freopen(in.txt, r, stdin);// freopen(out.txt, w, stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);int n;scanf(%d, n);for (int i 0; i n; i) {scanf(%d, a[i]);}polysqrt(a, b, n);for (int i 0; i n; i) {printf(%d%c, b[i], i 1 n ? \n : );}return 0;
}二次剩余解一般情况
#include bits/stdc.husing namespace std;typedef long long ll;const int N 5e6 10, mod 998244353, inv2 mod 1 1;int a[N], b[N], c[N], d[N], r[N];namespace Quadratic_residue {struct Complex {int r, i;Complex(int _r 0, int _i 0) : r(_r), i(_i) {}};int I2;Complex operator * (const Complex a, Complex b) {return Complex((1ll * a.r * b.r % mod 1ll * a.i * b.i % mod * I2 % mod) % mod, (1ll * a.r * b.i % mod 1ll * a.i * b.r % mod) % mod);}Complex quick_pow(Complex a, int n) {Complex ans Complex(1, 0);while (n) {if (n 1) {ans ans * a;}a a * a;n 1;}return ans;}int get_residue(int n) {mt19937 e(233);if (n 0) {return 0;}if(quick_pow(n, (mod - 1) 1).r mod - 1) {return -1;}uniform_int_distributionint r(0, mod - 1);int a r(e);while(quick_pow((1ll * a * a % mod - n mod) % mod, (mod - 1) 1).r 1) {a r(e);}I2 (1ll * a * a % mod - n mod) % mod;int x quick_pow(Complex(a, 1), (mod 1) 1).r, y mod - x;if(x y) swap(x, y);return x;}
}int quick_pow(int a, int n) {int ans 1;while (n) {if (n 1) {ans 1ll * ans * a % mod;}a 1ll * a * a % mod;n 1;}return ans;
}void get_r(int lim) {for (int i 0; i lim; i) {r[i] (i 1) * (lim 1) (r[i 1] 1);}
}void NTT(int *f, int lim, int rev) {for (int i 0; i lim; i) {if (i r[i]) {swap(f[i], f[r[i]]);}}for (int mid 1; mid lim; mid 1) {int wn quick_pow(3, (mod - 1) / (mid 1));for (int len mid 1, cur 0; cur lim; cur len) {int w 1;for (int k 0; k mid; k, w 1ll * w * wn % mod) {int x f[cur k], y 1ll * w * f[cur mid k] % mod;f[cur k] (x y) % mod, f[cur mid k] (x - y mod) % mod;}}}if (rev -1) {int inv quick_pow(lim, mod - 2);reverse(f 1, f lim);for (int i 0; i lim; i) {f[i] 1ll * f[i] * inv % mod;}}
}void polyinv1(int *a, int *b, int n) {if (n 1) {b[0] quick_pow(a[0], mod - 2);return ;}polyinv1(a, b, n 1 1);int lim 1;while (lim 2 * n) {lim 1;}get_r(lim);for (int i 0; i n; i) {c[i] a[i];}for (int i n; i lim; i) {c[i] 0;}NTT(b, lim, 1);NTT(c, lim, 1);for (int i 0; i lim; i) {int cur (2 - 1ll * c[i] * b[i] % mod mod) % mod;b[i] 1ll * b[i] * cur % mod;}NTT(b, lim, -1);for (int i n; i lim; i) {b[i] 0;}
}void polysqrt(int *a, int *b, int n) {if (n 1) {b[0] Quadratic_residue::get_residue(a[0]);return ;}polysqrt(a, b, n 1 1);int lim 1;while (lim 2 * n) {lim 1;}for (int i 0; i lim; i) {d[i] 0;}polyinv1(b, d, n);for (int i 0; i n; i) {c[i] a[i];}for (int i n; i lim; i) {c[i] 0;}get_r(lim);NTT(b, lim, 1);NTT(c, lim, 1);NTT(d, lim, 1);for (int i 0; i lim; i) {b[i] (1ll * inv2 * b[i] % mod 1ll * inv2 * d[i] % mod * c[i] % mod) % mod; }NTT(b, lim, -1);for (int i n; i lim; i) {b[i] 0;}
}int main() {// freopen(in.txt, r, stdin);// freopen(out.txt, w, stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);int n;scanf(%d, n);for (int i 0; i n; i) {scanf(%d, a[i]);}polysqrt(a, b, n);for (int i 0; i n; i) {printf(%d%c, b[i], i 1 n ? \n : );}return 0;
}分治FFT
考虑计算这么一个式子f(i)∑j1ifi−jg(j)f(i) \sum\limits_{j 1} ^{i} f_{i - j}g(j)f(i)j1∑ifi−jg(j)给定g(x)g(x)g(x)求f(x)f(x)f(x)边界条件f(0)1f(0) 1f(0)1。
假设我们已经算出[l,mid][l, mid][l,mid]考虑计算其对[mid1,r][mid 1, r][mid1,r]的贡献w(i)w(i)w(i)
有w(x)∑ilmidf(i)g(x−i)w(x) \sum\limits_{i l} ^{mid} f(i) g(x - i)w(x)il∑midf(i)g(x−i)因为[mid1,r][mid 1, r][mid1,r]区间还没开始计算所以f(i)0,i∈[mid,x−1]f(i) 0, i \in [mid, x - 1]f(i)0,i∈[mid,x−1]则w(x)∑ilx−1f(i)g(x−i)w(x) \sum\limits_{i l} ^{x - 1} f(i) g(x - i)w(x)il∑x−1f(i)g(x−i)
我们设a(i)f(il),b(i)g(i1)a(i) f(i l), b(i) g(i 1)a(i)f(il),b(i)g(i1)上式w(x)∑i0x−l−1a(i)b(x−l−i1)w(x) \sum\limits_{i 0} ^{x - l - 1} a(i)b(x - l - i 1)w(x)i0∑x−l−1a(i)b(x−l−i1)
有w(x)c(x−l1)∑i0x−l−1a(i)b(x−l−i1)w(x) c(x - l 1) \sum\limits_{i 0} ^{x - l - 1}a(i) b(x - l - i 1)w(x)c(x−l1)i0∑x−l−1a(i)b(x−l−i1)
#include bits/stdc.husing namespace std;typedef long long ll;const int N 5e6 10, mod 998244353, inv2 mod 1 1;int a[N], b[N], c[N], d[N], r[N];int quick_pow(int a, int n) {int ans 1;while (n) {if (n 1) {ans 1ll * ans * a % mod;}a 1ll * a * a % mod;n 1;}return ans;
}void get_r(int lim) {for (int i 0; i lim; i) {r[i] (i 1) * (lim 1) (r[i 1] 1);}
}void NTT(int *f, int lim, int rev) {for (int i 0; i lim; i) {if (i r[i]) {swap(f[i], f[r[i]]);}}for (int mid 1; mid lim; mid 1) {int wn quick_pow(3, (mod - 1) / (mid 1));for (int len mid 1, cur 0; cur lim; cur len) {int w 1;for (int k 0; k mid; k, w 1ll * w * wn % mod) {int x f[cur k], y 1ll * w * f[cur mid k] % mod;f[cur k] (x y) % mod, f[cur mid k] (x - y mod) % mod;}}}if (rev -1) {int inv quick_pow(lim, mod - 2);reverse(f 1, f lim);for (int i 0; i lim; i) {f[i] 1ll * f[i] * inv % mod;}}
}void DACFFT(int l, int r) {if (l r) {return ;}int mid l r 1;DACFFT(l, mid);for (int i 0; i mid - l; i) {c[i] b[i l];}for (int i 0; i r - l; i) {d[i] a[i 1];}int lim 1;while (lim r - l mid - l 1) {lim 1;}get_r(lim);NTT(c, lim, 1);NTT(d, lim, 1);for (int i 0; i lim; i) {c[i] 1ll * c[i] * d[i] % mod;}NTT(c, lim, -1);for (int i mid - l; i r - l - 1; i) {b[i l 1] (b[i l 1] c[i]) % mod;}for (int i 0; i lim; i) {c[i] d[i] 0;}DACFFT(mid 1, r);
}int main() {// freopen(in.txt, r, stdin);// freopen(out.txt, w, stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);int n;scanf(%d, n);for (int i 1; i n; i) {scanf(%d, a[i]);}b[0] 1;DACFFT(0, n - 1);for (int i 0; i n; i) {printf(%d%c, b[i], i n ? \n : );}return 0;
}牛顿迭代
迭代求函数零点
对于f(x)f(x)f(x)任意选取一个点(x0,f(x0))(x_0, f(x_0))(x0,f(x0))作为当前我们预估的零点取他的泰勒展开前两项g(x)f(0)f′(x)(x−x0)g(x) f(0) f(x)(x - x_0)g(x)f(0)f′(x)(x−x0)
解出方程g(x)0g(x) 0g(x)0得到x1x_1x1然后重复上述操作最后xnx_nxn会趋近于我们所要的正解。
考虑如何应用到多项式上
边界条件n1n 1n1时[x0]g(f(x))0[x ^0]g(f(x)) 0[x0]g(f(x))0的解单独求出。
假设我们已经求得膜x⌈n2⌉x ^{\lceil\frac{n}{2} \rceil}x⌈2n⌉下的解f0(x)f_0(x)f0(x)要求膜xnx ^nxn下的解f(x)f(x)f(x)得到该点的泰勒展开 ∑i0∞g(i)(f0(x))i!(f(x)−f0(x))nf(x)−f0(x)前面的项已经被截了,所以最低次幂是大于⌈n2⌉的有(f(x)−f0(x))i≡0(modxn),i2有g(f0(x))g′(f0(x))(f(x)−f0(x))≡0(modxn)f(x)≡f0(x)−g(f0(x))g′(f0(x))(modxn)\sum_{i 0} ^{\infty} \frac{g ^{(i)}(f_0(x))}{i !}(f(x) - f_0(x)) ^n\\ f(x) - f_0(x)前面的项已经被截了,所以最低次幂是大于\lceil \frac{n}{2} \rceil的\\ 有(f(x) - f_0 (x)) ^ i \equiv 0 \pmod {x ^ n}, i 2\\ 有g(f_0(x)) g(f_0(x))(f(x) - f_0(x)) \equiv 0 \pmod {x ^ n}\\ f(x) \equiv f_0(x) - \frac{g(f_0(x))}{g(f_0(x))} \pmod {x ^ n}\\ i0∑∞i!g(i)(f0(x))(f(x)−f0(x))nf(x)−f0(x)前面的项已经被截了,所以最低次幂是大于⌈2n⌉的有(f(x)−f0(x))i≡0(modxn),i2有g(f0(x))g′(f0(x))(f(x)−f0(x))≡0(modxn)f(x)≡f0(x)−g′(f0(x))g(f0(x))(modxn)
应用
多项式求逆
对于给定的h(x)有g(f(x))1f(x)−h(x)≡0(modxn)f(x)≡f0(x)−1f0(x)−h(x)−1f02(x)(modxn)f(x)≡2f0(x)−h(x)f02(x)(modxn)f(x)≡f0(x)(2−h(x)f0(x))(modxn)对于给定的h(x)\\ 有g(f(x)) \frac{1}{f(x)} - h(x) \equiv 0 \pmod {x ^ n}\\ f(x) \equiv f_0(x) - \frac{\frac{1}{f_0(x)} - h(x)}{-\frac{1}{f_0 ^ 2(x)}} \pmod{x ^n}\\ f(x) \equiv 2f_0(x) - h(x) f_0 ^ 2(x) \pmod {x ^ n}\\ f(x) \equiv f_0(x)(2 - h(x)f_0(x)) \pmod {x ^ n}\\ 对于给定的h(x)有g(f(x))f(x)1−h(x)≡0(modxn)f(x)≡f0(x)−−f02(x)1f0(x)1−h(x)(modxn)f(x)≡2f0(x)−h(x)f02(x)(modxn)f(x)≡f0(x)(2−h(x)f0(x))(modxn)
多项式开根
对于给定的h(x)有g(f(x))f2(x)−h(x)≡0(modxn)f(x)≡f0(x)−f02(x)−h(x)2f0(x)(mod()xn)f(x)≡f02(x)h(x)2f0(x)(modxn)f(x)≡2−1f0(x)2−1f0−1(x)h(x)(modxn)对于给定的h(x)\\ 有g(f(x)) f ^ 2(x) - h(x) \equiv 0 \pmod {x ^ n}\\ f(x) \equiv f_0(x) - \frac{f_0 ^2(x) - h(x)}{2f_0(x)} \pmod (x ^n)\\ f(x) \equiv \frac{f_0 ^ 2(x) h(x)}{2f_0(x)} \pmod {x ^ n}\\ f(x) \equiv 2 ^{-1} f_0 (x) 2 ^{-1} f_0 ^{-1}(x) h(x) \pmod {x ^n}\\ 对于给定的h(x)有g(f(x))f2(x)−h(x)≡0(modxn)f(x)≡f0(x)−2f0(x)f02(x)−h(x)(mod()xn)f(x)≡2f0(x)f02(x)h(x)(modxn)f(x)≡2−1f0(x)2−1f0−1(x)h(x)(modxn)
多项式exp\expexp
对于给定的h(x)有g(f(x))lnf(x)−h(x)≡0(modxn)f(x)≡f0(x)−lnf0(x)−h(x)1f0(x)(modxn)f(x)≡f0(x)(1−lnf0(x)h(x))(modxn)对于给定的h(x)\\ 有g(f(x)) \ln f(x) - h(x) \equiv 0 \pmod {x ^ n}\\ f(x) \equiv f_0(x) - \frac{\ln f_0(x) - h(x)}{\frac{1}{f_0(x)}} \pmod {x ^ n}\\ f(x) \equiv f_0(x)(1 - \ln f_0(x) h(x)) \pmod {x ^ n}\\ 对于给定的h(x)有g(f(x))lnf(x)−h(x)≡0(modxn)f(x)≡f0(x)−f0(x)1lnf0(x)−h(x)(modxn)f(x)≡f0(x)(1−lnf0(x)h(x))(modxn)
多项式对数函数lnf(x)\ln f(x)lnf(x)
如果存在解必然有[x0]f(x)1[ x ^ 0]f(x) 1[x0]f(x)1
对lnf(x)\ln f(x)lnf(x)求导有dlnf(x)dx≡f′(x)f(x)(modxn)\frac{d \ln f(x)}{dx} \equiv \frac{f(x)}{f(x)} \pmod {x ^ n}dxdlnf(x)≡f(x)f′(x)(modxn)
dxdxdx乘到右边再求积分有 ∫dlnf(x)≡∫f′(x))f(x)dx(modxn)lnf(x)≡∫f′(x)f(x)(modxn)\int d \ln f(x) \equiv \int \frac{f(x))}{f(x)} dx \pmod {x ^ n}\\ \ln f(x) \equiv \int \frac{f(x)}{f(x)} \pmod {x ^ n}\\ ∫dlnf(x)≡∫f(x)f′(x))dx(modxn)lnf(x)≡∫f(x)f′(x)(modxn) 然后只要对先对f(x)f(x)f(x)求个导求个逆最后求一次积分即可整体复杂度O(nlogn)O(n \log n)O(nlogn)。
#include bits/stdc.husing namespace std;typedef long long ll;const int N 5e6 10, mod 998244353, inv2 mod 1 1;int a[N], b[N], c[N], d[N], r[N], inv[N];int quick_pow(int a, int n) {int ans 1;while (n) {if (n 1) {ans 1ll * ans * a % mod;}a 1ll * a * a % mod;n 1;}return ans;
}void get_r(int lim) {for (int i 0; i lim; i) {r[i] (i 1) * (lim 1) (r[i 1] 1);}
}void get_inv(int n) {inv[1] 1;for (int i 2; i n; i) {inv[i] 1ll * (mod - mod / i) * inv[mod % i] % mod;}
}void NTT(int *f, int lim, int rev) {for (int i 0; i lim; i) {if (i r[i]) {swap(f[i], f[r[i]]);}}for (int mid 1; mid lim; mid 1) {int wn quick_pow(3, (mod - 1) / (mid 1));for (int len mid 1, cur 0; cur lim; cur len) {int w 1;for (int k 0; k mid; k, w 1ll * w * wn % mod) {int x f[cur k], y 1ll * w * f[cur mid k] % mod;f[cur k] (x y) % mod, f[cur mid k] (x - y mod) % mod;}}}if (rev -1) {int inv quick_pow(lim, mod - 2);reverse(f 1, f lim);for (int i 0; i lim; i) {f[i] 1ll * f[i] * inv % mod;}}
}void polyinv(int *a, int *b, int n) {if (n 1) {b[0] quick_pow(a[0], mod - 2);return ;}polyinv(a, b, n 1 1);int lim 1;while (lim 2 * n) {lim 1;}get_r(lim);for (int i 0; i n; i) {c[i] a[i];}for (int i n; i lim; i) {c[i] 0;}NTT(b, lim, 1);NTT(c, lim, 1);for (int i 0; i lim; i) {int cur (2 - 1ll * c[i] * b[i] % mod mod) % mod;b[i] 1ll * b[i] * cur % mod;}NTT(b, lim, -1);for (int i n; i lim; i) {b[i] 0;}
}void derivative(int *a, int *b, int n) {for (int i 0; i n; i) {b[i] 1ll * a[i 1] * (i 1) % mod;}
}void integrate(int *a, int n) {get_inv(n);for (int i n - 1; i 1; i--) {a[i] 1ll * a[i - 1] * inv[i] % mod;}a[0] 0;
}void polyln(int *a, int *b, int n) {derivative(a, b, n);polyinv(a, d, n);int lim 1;while (lim 2 * n) {lim 1;}get_r(lim);NTT(b, lim, 1);NTT(d, lim, 1);for (int i 0; i lim; i) {b[i] 1ll * b[i] * d[i] % mod;}NTT(b, lim, -1);integrate(b, n);
}int main() {// freopen(in.txt, r, stdin);// freopen(out.txt, w, stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);int n;scanf(%d, n);for (int i 0; i n; i) {scanf(%d, a[i]);}polyln(a, b, n);for (int i 0; i n; i) {printf(%d%c, b[i], i 1 n ? \n : );}return 0;
}多项式exp\expexp
牛顿迭代推导式子有$f(x) \equiv f_0(x)(1 - \ln f_0(x) h(x)) \pmod {x ^ n}\$
#include bits/stdc.husing namespace std;const int mod 998244353, inv2 mod 1 1;namespace Quadratic_residue {struct Complex {int r, i;Complex(int _r 0, int _i 0) : r(_r), i(_i) {}};int I2;Complex operator * (const Complex a, Complex b) {return Complex((1ll * a.r * b.r % mod 1ll * a.i * b.i % mod * I2 % mod) % mod, (1ll * a.r * b.i % mod 1ll * a.i * b.r % mod) % mod);}Complex quick_pow(Complex a, int n) {Complex ans Complex(1, 0);while (n) {if (n 1) {ans ans * a;}a a * a;n 1;}return ans;}int get_residue(int n) {mt19937 e(233);if (n 0) {return 0;}if(quick_pow(n, (mod - 1) 1).r mod - 1) {return -1;}uniform_int_distributionint r(0, mod - 1);int a r(e);while(quick_pow((1ll * a * a % mod - n mod) % mod, (mod - 1) 1).r 1) {a r(e);}I2 (1ll * a * a % mod - n mod) % mod;int x quick_pow(Complex(a, 1), (mod 1) 1).r, y mod - x;if(x y) swap(x, y);return x;}
}const int N 1e6 10;int r[N], inv[N], a[N], b[N], c[N], d[N], e[N], t[N], n;//a是输入数组b存放多项式逆c存放多项式开根d存放多项式对数lne存放多项式指数expt作为中间转移数组int quick_pow(int a, int n) {int ans 1;while (n) {if (n 1) {ans 1ll * a * ans % mod;}a 1ll * a * a % mod;n 1;}return ans;
}void get_r(int lim) {for (int i 0; i lim; i) {r[i] (i 1) * (lim 1) (r[i 1] 1);}
}void get_inv(int n) {inv[1] 1;for (int i 2; i n; i) {inv[i] 1ll * (mod - mod / i) * inv[mod % i] % mod;}
}void NTT(int *f, int lim, int rev) {for (int i 0; i lim; i) {if (i r[i]) {swap(f[i], f[r[i]]);}}for (int mid 1; mid lim; mid 1) {int wn quick_pow(3, (mod - 1) / (mid 1));for (int len mid 1, cur 0; cur lim; cur len) {int w 1;for (int k 0; k mid; k, w 1ll * w * wn % mod) {int x f[cur k], y 1ll * w * f[cur mid k] % mod;f[cur k] (x y) % mod, f[cur mid k] (x - y mod) % mod;}}}if (rev -1) {int inv quick_pow(lim, mod - 2);reverse(f 1, f lim);for (int i 0; i lim; i) {f[i] 1ll * f[i] * inv % mod;}}
}void polyinv(int *f, int *g, int n) {if (n 1) {g[0] quick_pow(f[0], mod - 2);return ;}polyinv(f, g, n 1 1);for (int i 0; i n; i) {t[i] f[i];}int lim 1;while (lim 2 * n) {lim 1;}get_r(lim);NTT(t, lim, 1);NTT(g, lim, 1);for (int i 0; i lim; i) {int cur (2 - 1ll * g[i] * t[i] % mod mod) % mod;g[i] 1ll * g[i] * cur % mod;t[i] 0;}NTT(g, lim, -1);for (int i n; i lim; i) {g[i] 0;}
}void polysqrt(int *f, int *g, int n) {if (n 1) {g[0] Quadratic_residue::get_residue(f[0]);return ;}polysqrt(f, g, n 1 1);polyinv(g, b, n);int lim 1;while (lim 2 * n) {lim 1;}get_r(lim);for (int i 0; i n; i) {t[i] f[i];}NTT(g, lim, 1);NTT(b, lim, 1);NTT(t, lim, 1);for (int i 0; i lim; i) {g[i] (1ll * inv2 * g[i] % mod 1ll * inv2 * b[i] % mod * t[i] % mod) % mod;b[i] t[i] 0;}NTT(g, lim, -1);for (int i n; i lim; i) {g[i] 0;}
}void derivative(int *a, int *b, int n) {for (int i 0; i n; i) {b[i] 1ll * a[i 1] * (i 1) % mod;}
}void integrate(int *a, int n) {for (int i n - 1; i 1; i--) {a[i] 1ll * a[i - 1] * inv[i] % mod;}a[0] 0;
}void polyln(int *f, int *g, int n) {polyinv(f, b, n);derivative(f, g, n);int lim 1;while (lim 2 * n) {lim 1;}get_r(lim);NTT(g, lim, 1);NTT(b, lim, 1);for (int i 0; i lim; i) {g[i] 1ll * g[i] * b[i] % mod;b[i] 0;}NTT(g, lim, -1);for (int i n; i lim; i) {g[i] 0;}integrate(g, n);
}void polyexp(int *f, int *g, int n) {if (n 1) {g[0] 1;return ;}polyexp(f, g, n 1 1);int lim 1;while (lim 2 * n) {lim 1;}polyln(g, d, n);for (int i 0; i n; i) {t[i] (f[i] - d[i] mod) % mod;}t[0] (t[0] 1) % mod;get_r(lim);NTT(g, lim, 1);NTT(t, lim, 1);for (int i 0; i lim; i) {g[i] 1ll * g[i] * t[i] % mod;t[i] d[i] 0;}NTT(g, lim, -1);for (int i n; i lim; i) {g[i] 0;}
}int main() {// freopen(in.txt, r, stdin);// freopen(out.txt, w, stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);scanf(%d, n);for (int i 0; i n; i) {scanf(%d, a[i]);}get_inv(4 * n);polyexp(a, e, n);for (int i 0; i n; i) {printf(%d%c, e[i], i 1 n ? \n : );}return 0;
}多项式快速幂
给定f(x)f(x)f(x)要求g(x)≡fk(x)(modxn)g(x) \equiv f ^{k}(x) \pmod{x ^ n}g(x)≡fk(x)(modxn)。 g(x)explnfk(x)expklnf(x)g(x) \exp ^{\ln f^{k}(x)} \exp ^{k \ln f(x)}\\ g(x)explnfk(x)expklnf(x)
#include bits/stdc.husing namespace std;const int mod 998244353, inv2 mod 1 1;namespace Quadratic_residue {struct Complex {int r, i;Complex(int _r 0, int _i 0) : r(_r), i(_i) {}};int I2;Complex operator * (const Complex a, Complex b) {return Complex((1ll * a.r * b.r % mod 1ll * a.i * b.i % mod * I2 % mod) % mod, (1ll * a.r * b.i % mod 1ll * a.i * b.r % mod) % mod);}Complex quick_pow(Complex a, int n) {Complex ans Complex(1, 0);while (n) {if (n 1) {ans ans * a;}a a * a;n 1;}return ans;}int get_residue(int n) {mt19937 e(233);if (n 0) {return 0;}if(quick_pow(n, (mod - 1) 1).r mod - 1) {return -1;}uniform_int_distributionint r(0, mod - 1);int a r(e);while(quick_pow((1ll * a * a % mod - n mod) % mod, (mod - 1) 1).r 1) {a r(e);}I2 (1ll * a * a % mod - n mod) % mod;int x quick_pow(Complex(a, 1), (mod 1) 1).r, y mod - x;if(x y) swap(x, y);return x;}
}const int N 1e6 10;int r[N], inv[N], a[N], b[N], c[N], d[N], e[N], t[N], n;int quick_pow(int a, int n) {int ans 1;while (n) {if (n 1) {ans 1ll * a * ans % mod;}a 1ll * a * a % mod;n 1;}return ans;
}void get_r(int lim) {for (int i 0; i lim; i) {r[i] (i 1) * (lim 1) (r[i 1] 1);}
}void get_inv(int n) {inv[1] 1;for (int i 2; i n; i) {inv[i] 1ll * (mod - mod / i) * inv[mod % i] % mod;}
}void NTT(int *f, int lim, int rev) {for (int i 0; i lim; i) {if (i r[i]) {swap(f[i], f[r[i]]);}}for (int mid 1; mid lim; mid 1) {int wn quick_pow(3, (mod - 1) / (mid 1));for (int len mid 1, cur 0; cur lim; cur len) {int w 1;for (int k 0; k mid; k, w 1ll * w * wn % mod) {int x f[cur k], y 1ll * w * f[cur mid k] % mod;f[cur k] (x y) % mod, f[cur mid k] (x - y mod) % mod;}}}if (rev -1) {int inv quick_pow(lim, mod - 2);reverse(f 1, f lim);for (int i 0; i lim; i) {f[i] 1ll * f[i] * inv % mod;}}
}void polyinv(int *f, int *g, int n) {if (n 1) {g[0] quick_pow(f[0], mod - 2);return ;}polyinv(f, g, n 1 1);for (int i 0; i n; i) {t[i] f[i];}int lim 1;while (lim 2 * n) {lim 1;}get_r(lim);NTT(t, lim, 1);NTT(g, lim, 1);for (int i 0; i lim; i) {int cur (2 - 1ll * g[i] * t[i] % mod mod) % mod;g[i] 1ll * g[i] * cur % mod;t[i] 0;}NTT(g, lim, -1);for (int i n; i lim; i) {g[i] 0;}
}void polysqrt(int *f, int *g, int n) {if (n 1) {g[0] Quadratic_residue::get_residue(f[0]);return ;}polysqrt(f, g, n 1 1);polyinv(g, b, n);int lim 1;while (lim 2 * n) {lim 1;}get_r(lim);for (int i 0; i n; i) {t[i] f[i];}NTT(g, lim, 1);NTT(b, lim, 1);NTT(t, lim, 1);for (int i 0; i lim; i) {g[i] (1ll * inv2 * g[i] % mod 1ll * inv2 * b[i] % mod * t[i] % mod) % mod;b[i] t[i] 0;}NTT(g, lim, -1);for (int i n; i lim; i) {g[i] 0;}
}void derivative(int *a, int *b, int n) {for (int i 0; i n; i) {b[i] 1ll * a[i 1] * (i 1) % mod;}
}void integrate(int *a, int n) {for (int i n - 1; i 1; i--) {a[i] 1ll * a[i - 1] * inv[i] % mod;}a[0] 0;
}void polyln(int *f, int *g, int n) {polyinv(f, b, n);derivative(f, g, n);int lim 1;while (lim 2 * n) {lim 1;}get_r(lim);NTT(g, lim, 1);NTT(b, lim, 1);for (int i 0; i lim; i) {g[i] 1ll * g[i] * b[i] % mod;b[i] 0;}NTT(g, lim, -1);for (int i n; i lim; i) {g[i] 0;}integrate(g, n);
}void polyexp(int *f, int *g, int n) {if (n 1) {g[0] 1;return ;}polyexp(f, g, n 1 1);int lim 1;while (lim 2 * n) {lim 1;}polyln(g, d, n);for (int i 0; i n; i) {t[i] (f[i] - d[i] mod) % mod;}t[0] (t[0] 1) % mod;get_r(lim);NTT(g, lim, 1);NTT(t, lim, 1);for (int i 0; i lim; i) {g[i] 1ll * g[i] * t[i] % mod;t[i] d[i] 0;}NTT(g, lim, -1);for (int i n; i lim; i) {g[i] 0;}
}/*a是输入数组b存放多项式逆c存放多项式开根d存放多项式对数lne存放多项式指数expt作为中间转移数组,如果要用到polyinv得提前调用get_inv(n)先预先得到我们想要得到的逆元范围。
*/char str[N];int main() {// freopen(in.txt, r, stdin);// freopen(out.txt, w, stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);scanf(%d %s, n, str 1);for (int i 0; i n; i) {scanf(%d, a[i]);}int k 0;for (int i 1; str[i]; i) {k (1ll * k * 10 (str[i] - 0)) % mod;}get_inv(4 * n);polyln(a, d, n);for (int i 0; i n; i) {a[i] 1ll * d[i] * k % mod;d[i] 0;}polyexp(a, e, n);for (int i 0; i n; i) {printf(%d%c, e[i], i 1 n ? \n : );}return 0;
}多项式除法
给定一个nnn次多项式F(x)F(x)F(x)和mmm次多项式G(x)G(x)G(x)要求R(x),Q(x)R(x), Q(x)R(x),Q(x)满足F(x)R(x)G(x)Q(x)F(x) R(x)G(x) Q(x)F(x)R(x)G(x)Q(x)。
R(x)R(x)R(x)是一个n−mn - mn−m阶多项式Q(x)Q(x)Q(x)是一个小于mmm阶的多项式。 有F(x)≡R(x)G(x)Q(x)(modxn1)F(1x)≡R(1x)G(1x)Q(1x)(modxn1)同时乘上一个xn,Frev(x)≡(xmRrev(x))(xn−mGrev(x))xn−degQQrev(x)(modxn1)Frev(x)≡Rrev(x)Grev(x)Qrev(x)xn−degQ(modxn1)有degQm,n−degQn−m1,所以有Frev(x)≡Rrev(x)Grev(x)(modxn−m1)有F(x) \equiv R(x) G(x) Q(x) \pmod{x ^ {n 1}}\\ F(\frac{1}{x}) \equiv R(\frac{1}{x})G(\frac{1}{x}) Q(\frac{1}{x}) \pmod {x ^{n 1}}\\ 同时乘上一个x ^ n, F^{rev}(x) \equiv \left(x ^m R ^{rev}(x)\right) \left(x ^{n - m}G ^{rev}(x)\right) x ^{n - deg_Q} Q ^{rev} (x) \pmod{x ^{n 1}} \\ F^{rev}(x) \equiv R^{rev}(x) G^{rev}(x) Q^{rev}(x) x ^{n - deg_Q}\pmod{x ^{n 1}}\\ 有deg_Q m, n - deg_Q n - m 1,所以有F^{rev}(x) \equiv R^{rev}(x) G^{rev}(x) \pmod{x ^{n - m 1}}\\ 有F(x)≡R(x)G(x)Q(x)(modxn1)F(x1)≡R(x1)G(x1)Q(x1)(modxn1)同时乘上一个xn,Frev(x)≡(xmRrev(x))(xn−mGrev(x))xn−degQQrev(x)(modxn1)Frev(x)≡Rrev(x)Grev(x)Qrev(x)xn−degQ(modxn1)有degQm,n−degQn−m1,所以有Frev(x)≡Rrev(x)Grev(x)(modxn−m1) 只要多项式求逆即可得到R(x)R(x)R(x)然后代入原式求得Q(x)Q(x)Q(x)。
#include bits/stdc.husing namespace std;const int mod 998244353, inv2 mod 1 1;namespace Quadratic_residue {struct Complex {int r, i;Complex(int _r 0, int _i 0) : r(_r), i(_i) {}};int I2;Complex operator * (const Complex a, Complex b) {return Complex((1ll * a.r * b.r % mod 1ll * a.i * b.i % mod * I2 % mod) % mod, (1ll * a.r * b.i % mod 1ll * a.i * b.r % mod) % mod);}Complex quick_pow(Complex a, int n) {Complex ans Complex(1, 0);while (n) {if (n 1) {ans ans * a;}a a * a;n 1;}return ans;}int get_residue(int n) {mt19937 e(233);if (n 0) {return 0;}if(quick_pow(n, (mod - 1) 1).r mod - 1) {return -1;}uniform_int_distributionint r(0, mod - 1);int a r(e);while(quick_pow((1ll * a * a % mod - n mod) % mod, (mod - 1) 1).r 1) {a r(e);}I2 (1ll * a * a % mod - n mod) % mod;int x quick_pow(Complex(a, 1), (mod 1) 1).r, y mod - x;if(x y) swap(x, y);return x;}
}const int N 6e5 10;int r[N], inv[N], b[N], c[N], d[N], e[N], t[N];int quick_pow(int a, int n) {int ans 1;while (n) {if (n 1) {ans 1ll * a * ans % mod;}a 1ll * a * a % mod;n 1;}return ans;
}void get_r(int lim) {for (int i 0; i lim; i) {r[i] (i 1) * (lim 1) (r[i 1] 1);}
}void get_inv(int n) {inv[1] 1;for (int i 2; i n; i) {inv[i] 1ll * (mod - mod / i) * inv[mod % i] % mod;}
}void NTT(int *f, int lim, int rev) {for (int i 0; i lim; i) {if (i r[i]) {swap(f[i], f[r[i]]);}}for (int mid 1; mid lim; mid 1) {int wn quick_pow(3, (mod - 1) / (mid 1));for (int len mid 1, cur 0; cur lim; cur len) {int w 1;for (int k 0; k mid; k, w 1ll * w * wn % mod) {int x f[cur k], y 1ll * w * f[cur mid k] % mod;f[cur k] (x y) % mod, f[cur mid k] (x - y mod) % mod;}}}if (rev -1) {int inv quick_pow(lim, mod - 2);reverse(f 1, f lim);for (int i 0; i lim; i) {f[i] 1ll * f[i] * inv % mod;}}
}void polyinv(int *f, int *g, int n) {if (n 1) {g[0] quick_pow(f[0], mod - 2);return ;}polyinv(f, g, n 1 1);for (int i 0; i n; i) {t[i] f[i];}int lim 1;while (lim 2 * n) {lim 1;}get_r(lim);NTT(t, lim, 1);NTT(g, lim, 1);for (int i 0; i lim; i) {int cur (2 - 1ll * g[i] * t[i] % mod mod) % mod;g[i] 1ll * g[i] * cur % mod;t[i] 0;}NTT(g, lim, -1);for (int i n; i lim; i) {g[i] 0;}
}void polysqrt(int *f, int *g, int n) {if (n 1) {g[0] Quadratic_residue::get_residue(f[0]);return ;}polysqrt(f, g, n 1 1);polyinv(g, b, n);int lim 1;while (lim 2 * n) {lim 1;}get_r(lim);for (int i 0; i n; i) {t[i] f[i];}NTT(g, lim, 1);NTT(b, lim, 1);NTT(t, lim, 1);for (int i 0; i lim; i) {g[i] (1ll * inv2 * g[i] % mod 1ll * inv2 * b[i] % mod * t[i] % mod) % mod;b[i] t[i] 0;}NTT(g, lim, -1);for (int i n; i lim; i) {g[i] 0;}
}void derivative(int *a, int *b, int n) {for (int i 0; i n; i) {b[i] 1ll * a[i 1] * (i 1) % mod;}
}void integrate(int *a, int n) {for (int i n - 1; i 1; i--) {a[i] 1ll * a[i - 1] * inv[i] % mod;}a[0] 0;
}void polyln(int *f, int *g, int n) {polyinv(f, b, n);derivative(f, g, n);int lim 1;while (lim 2 * n) {lim 1;}get_r(lim);NTT(g, lim, 1);NTT(b, lim, 1);for (int i 0; i lim; i) {g[i] 1ll * g[i] * b[i] % mod;b[i] 0;}NTT(g, lim, -1);for (int i n; i lim; i) {g[i] 0;}integrate(g, n);
}void polyexp(int *f, int *g, int n) {if (n 1) {g[0] 1;return ;}polyexp(f, g, n 1 1);int lim 1;while (lim 2 * n) {lim 1;}polyln(g, d, n);for (int i 0; i n; i) {t[i] (f[i] - d[i] mod) % mod;}t[0] (t[0] 1) % mod;get_r(lim);NTT(g, lim, 1);NTT(t, lim, 1);for (int i 0; i lim; i) {g[i] 1ll * g[i] * t[i] % mod;t[i] d[i] 0;}NTT(g, lim, -1);for (int i n; i lim; i) {g[i] 0;}
}/*b存放多项式逆c存放多项式开根d存放多项式对数lne存放多项式指数expt作为中间转移数组,如果要用到polyinv得提前调用get_inv(n)先预先得到我们想要得到的逆元范围。
*/int f[N], fr[N], g[N], gr[N], rr[N], n, m;int main() {// freopen(in.txt, r, stdin);// freopen(out.txt, w, stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);scanf(%d %d, n, m);for (int i 0; i n; i) {scanf(%d, f[i]);fr[n - i] f[i];}for (int i 0; i m; i) {scanf(%d, g[i]);gr[m - i] g[i];}for (int i n - m 1; i n; i) {fr[i] gr[i] 0;}polyinv(gr, b, n - m 1);for (int i 0; i n - m 1; i) {gr[i] b[i];b[i] 0;}int lim 1;while (lim 2 * (n - m 1)) {lim 1;}get_r(lim);NTT(fr, lim, 1);NTT(gr, lim, 1);for (int i 0; i lim; i) {fr[i] 1ll * fr[i] * gr[i] % mod;}NTT(fr, lim, -1);for (int i 0; i n - m; i) {rr[i] fr[n - m - i];}for (int i 0; rr[i]; i) {printf(%d , rr[i]);}puts();lim 1;while (lim 2 * n) {lim 1;}get_r(lim);NTT(rr, lim, 1);NTT(g, lim, 1);for (int i 0; i lim; i) {g[i] 1ll * g[i] * rr[i] % mod;}NTT(g, lim, -1);for (int i 0; i m; i) {f[i] (f[i] - g[i] mod) % mod;}for (int i 0; i m; i) {printf(%d , f[i]);}puts();return 0;
}
