快速傅里叶变换(FFT)

多项式表示

系数表示法：

一个$n$次多项式可以用$n+1$个系数表示出来：$f(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_{n-1}{x}^{n-1}+a_n{x}^n$。

点值表示法：

通过线性代数，高斯消元我们可以知道，一个$n$次多项式可以通过$n+1$个点联立方程组解得：

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​))}。

有离散傅里叶变换$DFT$（把一个多项式从系数表示变成点值表示），$IDFT$（把一个多项式从点值表示变成系数表示），

而$FFT$，就是通过选取某些特殊$x$点，来加速$DFT,IDFT$的一种方法。

点值表示法的多项式相乘：

$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​))}

由此我们想要得到两个多项式相乘的系数，只需要先对两个多项式进行$DFT$，然后对应的点值相乘，再做一次$IDFT$，即可求得系数。

引入复数

两个复数相乘的结果为，模长相乘，辐角相加，证明如下：

有两复数$A(a\cos {\theta }_1,a\sin {\theta }_{1}i),B(b\cos {\theta }_2,b\sin {\theta }_{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)$。

引入$n$次复根，即$x^n=1$，这样的解显然有$n$个，设$w_n^i=e^{\frac{2\pi }{n}i}$，在复平面内即是把一个圆分成了$n$等份。

我们取$n$等分的第一个交所对应的向量$w_n=\cos (\frac{2\pi }{n})+\sin (\frac{2\pi }{n})i$，则其他复根都可用$w_n$的$i$次幂来表示。

快速傅里叶变换

考虑如何分治求解：

$f(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}^7$

按照$x$的次幂，分奇偶，再在右边提出一个$x$，

$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)$

$G(x)=a_0+a_{2}x+a_4{x}^2+a_6{x}^3$

$H(x)=a_1+a_{3}x+a_5{x}^2+a_7{x}^3$

有$f(x)=G(x^2)+xH(x^2)$

由单位复根有

$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(w_n^{k+\frac{n}{2}}))=DFT(G(w_{\frac{n}{2}}^k))-w_n^{k}DFT(H(w_{\frac{n}{2}}^k))$

由此，求出$DFT(G(w_{\frac{n}{2}}^k))$和$DFT(H(w_{\frac{n}{2}}^k))$即可知$DFT(f(w_n^k)),DFT(f(w_n^{k+\frac{n}{2}}))$然后对$G,H$再分别递归求解即可。

快速傅里叶逆变换

把单位复根值代入多项式，得到的是如下结果：

$\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ ⋮\\ y_{n-2}\\ y_{n-1}\end{array}\right]$ = $\left[\begin{array}{cccccc}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)}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 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{array}\right]$ $\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_{n-2}\\ a_{n-1}\end{array}\right]$

经过$DFT$我们已经得到了左边的矩阵，考虑如何变换得到右边的系数矩阵，线性代数我们知道，只要在左边乘上一个中间大矩阵的逆，我们即可得到右边的系数矩阵。

由于这个矩阵的元素非常特殊，他的逆矩阵也有特殊的性质，就是每一项取倒数，再除以$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})$，所以我们只要将这个代入做一次$DFT$，也就是$IDFT$，最后再对整体除以$n$即可得到系数矩阵。

对以上进行证明

$f(x)=\sum _{i=0}^{n-1}{a}_i{x}^i$，$y_i=f(w_n^i)$，构造$A(x)=\sum _{i=0}^{n-1}{y}_i{x}^i$，将$b_i=w_n^{-i}$代入多项式$A(x)$

有$A(b_k)=\sum _{i=0}^{n-1}{y}_i{w}_n^{-ik}=\sum _{i=0}^{n1}{w}_n^{-ik}\sum _{j=0}^{n-1}{a}_j{w}_n^{ij}=\sum _{j=0}^{n-1}{a}_j\sum _{i=0}^{n-1}(w_n^{j-k}{)}^i$

令$S(w_n^a)=\sum _{i=0}^{n-1}(w_n^a{)}^i$

显然有$a=0$，$S(w_n^a)=n$

$a\ne 0$，时我们取$S(w_n^a),w_n^{a}S(w_n^a)$，两者相减，除以一个系数有$S(w_n^a)=\frac{\sum _{i=1}^n(w_n^a{)}^i-\sum _{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}=0$

所以有$S(w_n^a)=[a=1]$

$A(b_k)=a_k\times n$

仔细想想，这个证明，就是把我们$DFT$过程中得到的点值作为系数去做一遍$DFT$，得到的也就是$A(x)$的点值表达式，同时对其除以$n$，也就是$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​}

这个过程中有一个规律，例如$1=001$，倒置后变成了$100$，$4$，也即是最后$a_1$所在的位置。

$r[i]$表示$i$翻转之后的数字，考虑如何从小到大递推得到$r[i]$，有$r[0]=0$，当我们在求$x$时，先考虑除个位数以外的数，就是$r[x>>1]>>1$了，如果个位是$1$则加上$lim>>1$，就有了如下代码

void change(int lim) {for (int i = 0; i < lim; i++) {r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);}
}

真正可用的$FFT$代码

P3803 【模板】多项式乘法（FFT）

#include <bits/stdc++.h>using 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++.h>using 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)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$，称$f(x)$为$g(x)$或者$g(x)$为$f(x)$膜$x^n$意义下的逆元。

下面我们讨论给定$f(x)$，求其逆$f^{-1}(x)$。

倍增求解

假设我们已经求得$f(x)$膜$x^{\lceil \frac{n}{2}}\rceil$下的逆元$f_0^{-1}(x)$，要求$f^{-1}(x)$，即膜$x^n$下的逆元，则

$f(x)f_0^{-1}(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$

显然$f(x)f^{-1}(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$也是成立的

对两边同时乘以$f_0^{-1}(x)$并移项有

$f^{-1}(x)-f_0^{-1}(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$

对两边同时开方得到

$f^{-2}(x)-2f^{-1}{f}_0^{-1}(x)+f_0^{-2}(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$

我们再对两边乘上一个$f(x)$，则有

$f^{-1}(x)-2f_0^{-1}+f(x)f_0^{-2}(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$

再对其进行移项可得

$f^{-1}(x)\equiv f_0^{-1}(x)\left(2-f(x)f_0^{-1}(x)\right)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$

由此我们递归求解即可。

#include <bits/stdc++.h>using 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)$，求$f(x)$，满足$f^2(x)=g(x)$。

假设我们已经得到了$g(x)$，膜$x^{\lceil \frac{n}{2}\rceil }$下的根$f_0(x)$，要求膜$x^n$下的根$f(x)$

有$f_0^2(x)\equiv g(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\lceil \frac{n}{2}\rceil })$

移项再开方有${\left(f_0^2(x)-g(x)\right)}^2\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$

则，${\left(f_0^2(x)+g(x)\right)}^2\equiv 4f_0^2(x)g(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$

$g(x)\equiv {\left(\frac{f_0^2(x)+g(x)}{2f_0(x)}\right)}^2(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$

所以$f(x)\equiv \frac{f_0^2(x)+g(x)}{2f_0(x)}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。

所以有$f(x)\equiv 2^{-1}{f}_0(x)+2^{-1}{f}_0^{-1}(x)g(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$，

对于$g(0)=1$的特殊情况

#include <bits/stdc++.h>using 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++.h>using 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_distribution<int> 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)=\sum _{j=1}^i{f}_{i-j}g(j)$，给定$g(x)$，求$f(x)$，边界条件$f(0)=1$。

假设我们已经算出$[l,mid]$，考虑计算其对$[mid+1,r]$的贡献$w(i)$，

有$w(x)=\sum _{i=l}^{mid}f(i)g(x-i)$，因为$[mid+1,r]$区间还没开始计算，所以$f(i)=0,i\in [mid,x-1]$，则$w(x)=\sum _{i=l}^{x-1}f(i)g(x-i)$

我们设$a(i)=f(i+l),b(i)=g(i+1)$，上式$w(x)=\sum _{i=0}^{x-l-1}a(i)b(x-l-i+1)$

有$w(x)=c(x-l+1)=\sum _{i=0}^{x-l-1}a(i)b(x-l-i+1)$

#include <bits/stdc++.h>using 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)$，任意选取一个点$(x_0,f(x_0))$作为当前我们预估的零点，取他的泰勒展开前两项$g(x)=f(0)+f^{\prime }(x)(x-x_0)$，

解出方程$g(x)=0$，得到$x_1$，然后重复上述操作，最后$x_n$会趋近于我们所要的正解。

考虑如何应用到多项式上

边界条件$n=1$时，$[x^0]g(f(x))=0$，的解单独求出。

假设我们已经求得膜$x^{\lceil \frac{n}{2}\rceil }$下的解$f_0(x)$，要求膜$x^n$下的解$f(x)$，得到该点的泰勒展开：
$$\begin{gathered} \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(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n),i>=2 \\ 有g(f_0(x))+g^{\prime }(f_0(x))(f(x)-f_0(x))\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n) \\ f(x)\equiv f_0(x)-\frac{g(f_0(x))}{g^{\prime }(f_0(x))}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n) \end{gathered}$$

应用

多项式求逆

$$\begin{gathered} 对于给定的h(x)， \\ 有g(f(x))=\frac{1}{f(x)}-h(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n) \\ f(x)\equiv f_0(x)-\frac{\frac{1}{f_0(x)}-h(x)}{-\frac{1}{f_0^2(x)}}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n) \\ f(x)\equiv 2f_0(x)-h(x)f_0^2(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n) \\ f(x)\equiv f_0(x)(2-h(x)f_0(x))(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n) \end{gathered}$$

多项式开根

$$\begin{gathered} 对于给定的h(x) \\ 有g(f(x))=f^2(x)-h(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n) \\ f(x)\equiv f_0(x)-\frac{f_0^2(x)-h(x)}{2f_0(x)}(\mathrm{m}\mathrm{o}\mathrm{d}()x^n) \\ f(x)\equiv \frac{f_0^2(x)+h(x)}{2f_0(x)}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n) \\ f(x)\equiv 2^{-1}{f}_0(x)+2^{-1}{f}_0^{-1}(x)h(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n) \end{gathered}$$

多项式$\exp$

$$\begin{gathered} 对于给定的h(x) \\ 有g(f(x))=\ln f(x)-h(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n) \\ f(x)\equiv f_0(x)-\frac{\ln f_0(x)-h(x)}{\frac{1}{f_0(x)}}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n) \\ f(x)\equiv f_0(x)(1-\ln f_0(x)+h(x))(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n) \end{gathered}$$

多项式对数函数$\ln f(x)$

如果存在解必然有$[x^0]f(x)=1$，

对$\ln f(x)$求导，有$\frac{d\ln f(x)}{dx}\equiv \frac{f^{\prime }(x)}{f(x)}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$，

$dx$乘到右边，再求积分有：
$$\begin{gathered} \int d\ln f(x)\equiv \int \frac{f^{\prime }(x))}{f(x)}dx(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n) \\ \ln f(x)\equiv \int \frac{f^{\prime }(x)}{f(x)}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n) \end{gathered}$$
然后只要对先对$f(x)$求个导，求个逆，最后求一次积分即可，整体复杂度$O(n\log n)$。

#include <bits/stdc++.h>using 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$

牛顿迭代推导式子有$f(x) \equiv f_0(x)(1 - \ln f_0(x) + h(x)) \pmod {x ^ n}\$

#include <bits/stdc++.h>using 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_distribution<int> 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存放多项式对数ln，e存放多项式指数exp，t作为中间转移数组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)$，要求$g(x)\equiv f^k(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$。
$$g(x)={\exp }^{\ln f^k(x)}={\exp }^{k\ln f(x)}$$

#include <bits/stdc++.h>using 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_distribution<int> 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存放多项式对数ln，e存放多项式指数exp，t作为中间转移数组,如果要用到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;
}

多项式除法

给定一个$n$次多项式$F(x)$和$m$次多项式$G(x)$，要求$R(x),Q(x)$，满足$F(x)=R(x)G(x)+Q(x)$。

$R(x)$是一个$n-m$阶多项式，$Q(x)$是一个小于$m$阶的多项式。
$$\begin{gathered} 有F(x)\equiv R(x)G(x)+Q(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n+1}) \\ F(\frac{1}{x})\equiv R(\frac{1}{x})G(\frac{1}{x})+Q(\frac{1}{x})(\mathrm{m}\mathrm{o}\mathrm{d}{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-de{g}_Q}{Q}^{rev}(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n+1}) \\ {F}^{rev}(x)\equiv R^{rev}(x)G^{rev}(x)+Q^{rev}(x)x^{n-de{g}_Q}(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n+1}) \\ 有de{g}_Q<m,n-de{g}_Q>=n-m+1,所以有F^{rev}(x)\equiv R^{rev}(x)G^{rev}(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n-m+1}) \end{gathered}$$
只要多项式求逆，即可得到$R(x)$，然后代入原式求得$Q(x)$。

#include <bits/stdc++.h>using 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_distribution<int> 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存放多项式对数ln，e存放多项式指数exp，t作为中间转移数组,如果要用到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;
}