多项式 & 形式幂级数
域的定义：
如果一个包含至少两个元素的集合 \(F\) 及两个运算 \(+\) 和 \(\cdot\) 满足：
- \((F,+)\) 是交换群，如果这个群的单位元为 \(0\)，那我们记 \(F^{\times}=F∖\{0\}\)；
- \((F^{\times},\cdot)\) 也是交换群，一般这个群的单位元被称为 \(1\)。
那我们就称 \((F,+,\cdot)\) 构成一个域。在不引起歧义的情况下，域的运算符默认是 \(+\) 和 \(\cdot\)。
比如我们常见的 \(\mathbb{Q},\mathbb{R},\mathbb{C}\) 都在常规的加法和乘法下构成一个域。
如果 \(p\) 是一个素数，那么模 \(p\) 的完全剩余系 \(\{0,1,\cdots,p-1\}\) 就在模 \(p\) 意义下的加法和乘法构成一个域，这个域常被记作 \(\mathbb{F}_p\) 或 \(\mathbb{Z}_p\)。
多项式和形式幂级数的定义：
如果 \((F,+,\cdot)\) 是一个域，那么 \(F\) 上的多项式指的是形如 \(a_0x^0+a_1x^1+a_2x^2+\cdots+a_nx^n\) 的式子，其中 \(a_0,a_1,\cdots,a_n\in F\)。这里的 \(x\) 只是一个标识符，没有实际含义。
为方便，我们记 \(a_0x^0+a_1x^1+a_2x^2+\cdots+a_nx^n\) 为 \(\displaystyle\sum_{k=0}^n a_kx^k\)，记 \(F\) 上的所有多项式为 \(F[x]\)。
如果 \(f(x)=a_0x^0+a_1x^1+a_2x^2+\cdots+a_nx^n\)，定义最大的使 \(a_k\) 不为 \(0\) 的 \(k\) 为 \(f(x)\) 的次数，记做 \(\operatorname{deg} f\)。如果所有的 \(a_k\) 都是 \(0\)，则称 \(f(x)\) 为 \(0\) 多项式，记作 \(f(x)=0\)，此时补充定义 \(\operatorname{deg} f=-\infty\)。
域上多项式也可以带入具体的 \(x\)，也可以让多项式相加、相减、相乘，只不过运算是在 \(F\) 上进行的。
推广后的多项式仍然保持原来的一些良好性质，比如加法具有结合和交换律，乘法对加减法有分配律。
当多项式的次数为无穷时就是形式幂级数。
部分形式幂级数有其封闭形式，即使有无穷多项，也可以写成一个较为简单的形式。
下面是一些常用的形式幂级数：
接下来，我们分两条路走，一条是形式幂级数的数学运用，一条是多项式算法。
形式幂级数的数学运用：
普通生成函数 / OGF：
对于一个序列 \(a_0,a_1,\cdots\)，它所对应的形式幂级数被称为 \(\{a_k\}_{k=0}^{\infty}\) 的普通生成函数，生成函数一般指普通生成函数。
普通生成函数求递推数列的通项公式：
最典型的，我们取斐波那契数列，\(f_0=0,f_1=1,f_i=f_{i-1}+f_{i-2}(i\ge 2)\)。
我们写出斐波那契数列的普通生成函数。
现在我们要求出 \(F(x)\) 的封闭形式。
先转换通项公式，把边界条件弄掉，我们可以将通项公式改写为 \(f_i=f_{i-1}+f_{i-2}+[i=1]\)，再令 \(f_i=0(i<0)\)。
然后推式子：
解得 \(F(x)=\dfrac{x}{1-x-x^2}\)。
这是一个有理函数，我们要求解它的系数的通项公式。
注意到有理函数 \(\dfrac{a}{(1-\rho x)^m}\) 的系数很简单，\(x^k\) 项系数为 \(a{\rho}^k{k+m-1\choose k}\)。
如果我们可以将 \(F(x)\) 表示为这样的有理函数之和，就可以快速计算通项公式。
所以我们对 \(F(x)\) 做部分分式分解，由于分母要分解成 \((1-\rho x)\) 而不是 \((x-\rho)\)，所以我们求分母多项式反转后的零点。
::::info[解释：]
简单的说就是：
::::
解方程 \(x^2-x-1=0\)，得到 \(x_1=\phi,x_2=\hat\phi\)。
我们设 \(\dfrac{x}{1-x-x^2}=\dfrac{a}{1-\phi x}+\dfrac{b}{1-\hat\phi x}\)，解得 \(a=\dfrac{1}{\sqrt 5},b=-\dfrac{1}{\sqrt 5}\)。
所以得到 \(F(x)\) 的 \(x^k\) 系数为 \(\dfrac{1}{\sqrt 5}(\phi^k-\hat\phi^k)\)，即 \(f_k=\dfrac{1}{\sqrt 5}\left(\left(\dfrac{1+\sqrt 5}{2}\right)^k-\left(\dfrac{1-\sqrt 5}{2}\right)^k\right)\)。
一般的递推数列的通项公式，步骤如下：
- 求出数列的普通生成函数 \(F(x)=\dfrac{P(x)}{Q(x)}\)。
- 求出 \(Q^R(x)\) 的零点。
- 待定系数求出 \(F(x)\) 的分式分解。
- 计算答案，化简答案。
普通生成函数的组合意义：
我们观察多项式乘法的过程，其实就是在两个括号中分别选取其中一项，然后系数相乘后再相加，指数是相加。
这种过程很像一种决策，例如二项式定理，系数 \({n\choose k}\) 是组合数，其实就是在 \(n\) 个括号中选取 \(k\) 个 \(x\)，和 \(n-k\) 个 \(1\)。
例题 AT_fps_24_a お菓子：
看一道例题，AT_fps_24_a お菓子。
每天可以选择 +1/3/4/6，这就像一种决策。
列出一天的决策的生成函数 \(x+x^3+x^4+x^6\)。
那么 \(D\) 天的生成函数就是 \((x+x^3+x^4+x^6)^D\)。
其中 \(x^N\) 的系数就是答案。
现在就是需要求出 \([x^N](x+x^3+x^4+x^6)^D\)。
注意到 \(x+x^3+x^4+x^6=x(1+x^2)(1+x^3)\)，暴力卷积就可以了。
例题 AT_fps_24_b 整数の組：
在看看另一道，AT_fps_24_b 整数の組。
将 \(a,b,c,d\) 的取值，理解成四种决策。
写出 \(a,b,c,d\) 的生成函数。
将这四个式子相乘，得到答案为 \([x^N]\dfrac{1}{(1-x)^2}=N+1\)。
习题 1：
AT_fps_24_i スコア，需要用到下面讲的多项式乘法。
P2000 拯救世界，这题需要快速高精度乘法，需要用到下面讲的多项式乘法。
指数生成函数 / EGF：
序列 \(b_k=\dfrac{a_k}{k!}\) 的生成函数被称为 \(\{a_k\}_{k=0}^{\infty}\) 的指数生成函数。
指数生成函数的组合意义：
我们考虑 \(\{a_k\}_{k=0}^{\infty},\{b_k\}_{k=0}^{\infty},\{c_k\}_{k=0}^{\infty}\) 的指数生成函数 \(A(x),B(x),C(x)\)，若 \(C(x)=A(x)B(x)\)，就有如下式子。
如果 \(a_i\) 表示 \(\text A\) 出现的方案数，\(b_i\) 表示 \(\text B\) 出现的方案数，那么 \(c_i\) 表示 \(\text A,\text B\) 共同出现 \(i\) 次的方案数。（乘以组合数是因为 \(\text A\) 和 \(\text B\) 不相同）
例题 AT_fps_24_f 色紙：
看一下这一道题，AT_fps_24_f 色紙。
红、黄、蓝三种颜色不同，且数量总和为 \(N\)，所以考虑指数生成函数的组合意义。
写出三种颜色的指数生成函数：
相乘，得到答案的指数生成函数为 \(\dfrac{e^{3x}-e^{-x}}4\)。
所以答案为 \(\dfrac{3^N-(-1)^N}4\)。
习题 2：
指数生成函数中 exp 的组合意义：
如果 \([\frac{x^N}{N!}]F(x)\) 表示恰有 \(N\) 个数的方案数，那么 \([\frac{x^N}{N!}]F^k(x)\) 的组合意义是将 \(N\) 个不同的数划分成 \(k\) 个集合的方案数，除以 \(k!\) 就相当于这 \(k\) 个集合相同。
对所有 \(k\) 求和，就可以得到 \(\exp F(x)\) 的组合意义是 \(N\) 个有标号元素构成的集合的生成集族的方案数。
这里讲的比较浅，后面实战运用的时候再讲细一点。
Lagrange 反演：
复合逆的定义和 Lagrange 反演公式：
现在有一个形式幂级数 \(F(x)\)，我们定义若形式幂级数 \(G(x)\)，满足 \(F\) 与 \(G\) 的复合 \(F(G(x))=x\) 或者 \(G(F(x))=x\)，那么称 \(F(x)\) 的复合逆为 \(G(x)\)。
一个形式幂级数存在复合逆，当且仅当其常数项为 \(0\)，且其一次项不为 \(0\)。证明直接考虑带入就行了。
Lagrange 反演公式是这样的。
如果形式幂级数 \(F(x)\) 存在复合逆且其复合逆为 \(G(x)\)。
那么：
当 \(k=1\) 时，更为常见：
Lagrange 反演的运用：
求解卡特兰数的通项公式：
卡特兰数的递推公式是 \(h_0=1,h_n=\sum_{i=0}^{n-1}h_ih_{n-1-i}(n>0)\)。
我们先考虑写出卡特兰数的生成函数 \(H(x)\)。
这是一个二次方程，直接用求根公式求出 \(H(x)\)，然后用广义二项式定理展开即可。
但是这里我们用 Lagrange 反演解决。
我们要求出 \([x^n]H(x)\)，可以先求出 \(H(x)\) 的复合逆，然后用 Lagrange 反演公式计算。
如果你真的直接开始计算 \(H(x)\)，你会发现最后的形式还有一个卷积，你根本不好算，甚至根本算不出来，问题出在哪里呢？
细心的你可能发现了，\([x^0]H(x)=1\) 啊，\(H(x)\) 根本就不存在复合逆！那咋办。
其实很简单，你令 \(F(x)=H(x)-1\)，计算 \(F(x)\) 的复合逆就可以了。将方程改写成如下形式：
设 \(G(x)=\dfrac{x}{(x+1)^2}\)，那么 \(G(F(x))=x\)，即 \(G(x)\) 为 \(F(x)\) 的复合逆。
那么：
所以当 \(n>0\) 时 \(h_n=[x^n]F(x)=\frac 1{n+1}{2n\choose n}\)，当 \(n=0\) 时，\(h_0=1=\frac 1{0+1}{0\choose 0}\)，所以可以不需要分类讨论，都是同一个通项公式。
例题 P2767 树的数量：
求大小为 \(n\) 的 \(m\) 叉树个数。
这其实和求卡特兰数的通项公式是一样的，当 \(m=2\) 时，就是卡特兰数的一种组合意义，答案就是 \(h_n\)。
我们再拿这道题练练手，熟悉一下 Lagrange 反演。
设 \([x^n]H(x)\) 为大小为 \(n\) 的 \(m\) 叉树个数。
那么 \(H(x)=1+xH^m(x)\)。接着去掉 \(H(x)\) 的常数项，设 \(F(x)=H(x)-1\)。
则 \(F(x)=x(F(x)+1)^m\)，即 \(\frac{F(x)}{(F(x)+1)^m}=x\)。
再设 \(G(x)=\frac{x}{(x+1)^m}\)，\(G(x)\) 为 \(F(x)\) 的复合逆。
那么：
所以当 \(n=0\) 时答案为 \(1\)，否则答案为 \(\frac 1n{mn\choose n-1}\)。
例题 P3978 [TJOI2015] 概率论：
这个题目会难一点，P3978 [TJOI2015] 概率论。
首先将期望转化为求叶子节点数量的总和，最后再除以卡特兰数即可。
设 \(a_{i,j}\) 表示大小为 \(i\) 的树且叶子节点数量为 \(j\) 的方案数。
考虑二元生成函数 \(T(x,y)=\sum_{i,j}a_{i,j}x^iy^j\)。
那么很容易列出方程 \(T=1+x(y+T^2-1)\)，最后 -1 是因为如果两棵子树都是空的，那么肯定是叶子节点，不能算入未新增叶子节点的方案中。
那么答案就为：
我们对 \(x\) 做 Lagrange 反演，先去除常数项，设 \(F(x,y)=T(x,y)-1\)。得到：
那么 \(F(x,y)\) 的复合逆为 \(G(x,y)=\frac{x}{y+x^2+2x}\)。
套用 Lagrange 反演公式：
由于 \(n\ge 1\)，所以不需要管常数项。
所以最终答案为 \(\frac{n+1}{2n\choose n}{2n-2\choose n-1}=\frac{n(n+1)}{4n-2}\)。
习题 3：
AT_abc222_h [ABC222H] Beautiful Binary Tree
多项式算法：
接下来讲讲 OI 中可能会用到的多项式算法。
多项式乘法：
普通的多项式乘法是每个元素挨个相乘在相加，这样很慢，时间复杂度是 \(O(n^2)\) 的。下面介绍一个更为快速的做法。
离散傅里叶变换 / DFT：
考虑换一种方式思考，我们知道 \(n+1\) 个横坐标不同的点确定一个 \(n\) 次多项式，那知道 \(n+1\) 个横坐标不同的点可以确定一个 \(n\) 次域上多项式吗？这是可以的，证明如下：
待定域上多项式 \(f(x)=a_0x^0+a_1x^1+a_2x^2+\cdots+a_nx^n\)，即证在域上的运算满足：
::::success[证明：]
第一个 \(n\times n\) 的矩阵是范德蒙德矩阵，其行列式为 \(\displaystyle\prod_{0\le j<i\le n}(x_i-x_j)\)，其显然不为 \(0\)，因为 \(x\) 互不相同，由线性方程的求解理论证毕。\(\Box\)
::::
后面我们会简称域上多项式为多项式。
考虑两个多项式 \(f(x),g(x)\in F[x]\)，设 \(f(x)\) 为 \(n\) 次多项式，\(g(x)\) 为 \(m\) 次多项式，则 \(f(x)g(x)\) 为 \(n+m\) 次多项式，由上述引理，我们只要知道 \(f(x_0)g(x_0),f(x_1)g(x_1),\cdots,f(x_{n+m})g(x_{n+m})\) 的值就可以确定 \(f(x)g(x)\) 这个多项式。
现在有两个困难。一：如何快速的通过 \(f(x)\) 的系数求解出 \(f(x_0),f(x_1),\cdots f(x_{n+m})\) 的值，即从 \(a\) 变到 \(y\)。二：如何快速的通过 \(f(x_0)g(x_0),f(x_1)g(x_1),\cdots,f(x_{n+m})g(x_{n+m})\) 求出 \(f(x)g(x)\) 的系数，即从 \(y\) 变到 \(a\)，其实也就是“一”的逆变换。
注意到，\(x_0,x_1,\cdots,x_{n+m}\) 我们可以任取，我们肯定想要 \(x\) 构成的范德蒙德矩阵的逆矩阵好求，这让我们想到了本原单位根！
我们设 \(N=n+m\)，\(\omega\) 为 \(N\) 次本原单位根，
则：
这样我们就不需要用 \(O(N^3)\) 的时间高斯消元求逆矩阵了，这非常的棒！
我们设列向量 \(a=(a_0,\cdots,a_{N-1})^{\top}\)，记 \(a\) 的离散傅里叶变换（DFT）为 \(\mathcal{F}a=P\times a\)，\(a\) 的离散傅里叶逆变换（IDFT）为 \(\mathcal{F}^{-1}a=P^{-1}\times a\)，我们有 \(\mathcal{F}^{-1}(\mathcal{F}a)=a\)。
可以发现离散傅里叶变换就是通过系数求出了函数值，离散傅里叶逆变换就是通过函数值求出了系数，且我们通过定义可以发现其满足一个性质：\(\mathcal{F}f\cdot\mathcal{F}g=\mathcal{F}(f\ast g)\)（这个式子特别重要，许多的卷积都需要构造出一种变换满足这个性质），其中 \(\cdot\) 是对应位置相乘，\(\ast\) 是卷积（在这里就是多项式乘法）。所以我们可以通过如下方式求出 \(f\ast g\)：\(f\ast g=\mathcal{F}^{-1}(\mathcal{F}(f\ast g))=\mathcal{F}^{-1}(\mathcal{F}f\cdot\mathcal{F}g)\)，\(\cdot\) 运算可以 \(O(N)\) 快速求，接下来的困难就是如何快速进行离散傅里叶变换和离散傅里叶逆变换。
虽然这个想法很天才，但是计算的时间复杂度并没有降低，进行矩阵乘法的时间复杂度还是 \(O(N^2)\) 的。
快速傅里叶变换 / FFT：
分治写法：
FFT 是 DFT 的一种基于分治思想的实现方式，可以将计算离散傅里叶变换的复杂度降低到 \(O(n\log n)\)（这里的 \(n\) 是两个多项式的次数之和）。
为了方便分治，不妨设 \(n=2^k\)，因为在多项式后面补 \(0\) 不会有影响。考虑 \(b=\mathcal{F}a\) 的第 \(k\) 个分量，\(b_k=\displaystyle\sum_{i=0}^{n-1}\omega^{ki}a_i=\sum_{i=0}^{\frac n2-1}(\omega^2)^{ki}a_{2i}+\omega^k\sum_{i=0}^{\frac n2-1}(\omega^2)^{ki}a_{2i+1}\)。注意到这个式子其实就是 \(b_k=(\mathcal{F}a')_{k\bmod \frac n2}+\omega^k(\mathcal{F}a'')_{k\bmod \frac n2}\)，其中 \(a'=(a_0,a_2,\cdots,a_{n-2})^{\top},a''=(a_1,a_3,\cdots,a_{n-1})^{\top}\)。
注意这里离散傅里叶变换的本原单位根并不相同。
\(k\) 需要模 \(\frac n2\) 的原因是因为当 \(k\ge\frac n2\) 时，\((\mathcal{F}a')_k\) 和 \((\mathcal{F}a'')_k\) 没有意义，但是 \((\omega^2)^{\frac n2}=1\)（\(\omega^2\) 是 \(\frac n2\) 次本原单位根），\((\omega^2)^k=(\omega^2)^{k\bmod\frac n2}\)。
将上面的式子变一下：\(b_k=(\mathcal{F}a')_k+\omega^k(\mathcal{F}a'')_k,b_{k+\frac n2}=(\mathcal{F}a')_k-\omega^k(\mathcal{F}a'')_k,0\le k<\frac n2\)。（因为 \(\omega^{\frac n2}=-1\)）
递归到最后一层就是直接返回，因为 \(1\) 次单位根就是 \(1\)，所以此时 \(\mathcal{F}a=a\)。
代码实现：
void FFT(int n, C *a, C g) {
if (n == 1) return;
static C tmp[kMaxN], t, x, y;
for (int i = 0; i < n / 2; i++) {
tmp[i] = a[i * 2], tmp[i + n / 2] = a[i * 2 + 1];
}
copy(tmp, tmp + n, a), FFT(n / 2, a, g * g), FFT(n / 2, a + n / 2, g * g), t = 1;
for (int i = 0; i < n / 2; i++, t = t * g) {
x = a[i], y = t * a[i + n / 2];
a[i] = x + y, a[i + n / 2] = x - y;
}
}
逆变换也很简单，逆变换其实就是将本原单位根 \(\omega\) 变成 \(\omega^{-1}\) 再跑一遍 FFT。
void Idft(int n, C *a, C g) {
FFT(N, a, {g.a, -g.b});
for (int i = 0; i < N; i++) {
a[i] = a[i] / N;
}
}
最终的代码就很容易写出来了。
#include <bits/stdc++.h>
using namespace std;
const int kMaxN = 4e6 + 5;
const double pi = acos(-1);
int n, m, N;
struct C { // 复数
double a, b;
C(double x = 0, double y = 0) { a = x, b = y; }
C operator+(const C &y) { return {a + y.a, b + y.b}; }
C operator-(const C &y) { return {a - y.a, b - y.b}; }
C operator*(const C &y) { return {a * y.a - b * y.b, a * y.b + b * y.a}; }
C operator*(const double &y) { return {a * y, b * y}; }
C operator/(const C &y) { return {(a * y.a + b * y.b) / (y.a * y.a + y.b * y.b), (y.a * b - a * y.b) / (y.a * y.a + y.b * y.b)}; }
C operator/(const double &y) { return {a / y, b / y}; }
} g, a[kMaxN], b[kMaxN];
void FFT(int n, C *a, C g) {
if (n == 1) return;
static C tmp[kMaxN], t, x, y;
for (int i = 0; i < n / 2; i++) {
tmp[i] = a[i * 2], tmp[i + n / 2] = a[i * 2 + 1];
}
copy(tmp, tmp + n, a), FFT(n / 2, a, g * g), FFT(n / 2, a + n / 2, g * g), t = 1;
for (int i = 0; i < n / 2; i++, t = t * g) {
x = a[i], y = t * a[i + n / 2];
a[i] = x + y, a[i + n / 2] = x - y;
}
}
void Idft(int n, C *a, C g) {
FFT(N, a, {g.a, -g.b});
for (int i = 0; i < N; i++) {
a[i] = a[i] / N;
}
}
int main() {
ios::sync_with_stdio(0), cin.tie(0);
cin >> n >> m, n++, m++;
for (int i = 0; i < n; i++) {
cin >> a[i].a;
}
for (int i = 0; i < m; i++) {
cin >> b[i].a;
}
N = pow(2, ceil(log2(n + m))), g = C(cosl(2 * pi / N), sinl(2 * pi / N));
FFT(N, a, g), FFT(N, b, g);
for (int i = 0; i < N; i++) {
a[i] = a[i] * b[i];
}
Idft(N, a, g);
for (int i = 0; i < n + m - 1; i++) {
cout << (int)(a[i].a + 0.5) << ' '; // 这里四舍五入是为了防止精度误差
}
return 0;
}
这份代码可以通过 P3803 【模板】多项式乘法（FFT）。
倍增写法：
倍增的写法比分治写法常数更小。
我们考虑分治中 \(a\) 是怎么变的：
- \((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)\)
我们设 \(r_i\) 为最后一次递归结束后 \(i\) 位置的元素是谁，通过观察我们可以发现 \(r_i\) 满足如下递推式：
其中 \(\lg=\log_2n\)。
有了 \(r\)，我们就可以不用递归，快速计算最后一层的 \(a\)，然后像递归一样向上归并即可。
void Init(int lg) {
for (int i = 0; i < 1 << lg; i++) {
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (lg - 1));
}
}
void FFT(int n, C *a, int type = 1) {
for (int i = 0; i < n; i++) {
if (i < r[i]) swap(a[i], a[r[i]]);
}
static C w, t, x, y;
for (int d = 2; d <= n; d *= 2) {
w = C(cosl(2 * pi / d), type * sinl(2 * pi / d));
for (int i = 0; i < n; i += d) {
t = C(1, 0);
for (int j = i; j < i + d / 2; j++, t = t * w) {
x = a[j], y = t * a[j + d / 2];
a[j] = x + y, a[j + d / 2] = x - y;
}
}
}
}
三次变两次优化：
由于 \((a+bi)^2=a^2-b^2+2abi\)，所以我们可以将 \(f\) 放在实部，\(g\) 放在虚部，那么平方后的多项式虚部除以二就是卷积后的结果。但是这样的精度误差比较大。
快速数论变换 / NTT：
由于 FFT 需要进行大量的复数运算，常数很大，并且需要用到浮点类型，精度也不高，所以我们引入快速数论变换 / NTT。NTT 其实就是在 \(\mathbb{F}_p\) 下的 FFT（\(p\) 是质数），但是复数如何取模呢？只需考虑本原单位根如何取模。想一想前面用到了本原单位根的什么性质？设 \(\omega_n\) 为 \(n\) 次单位根（\(n\) 为偶数）。
- \(\omega_n^n=1\)；
- \(1,\omega_n,\omega_n^2,\cdots,\omega_n^{n-1}\) 互不相同；
- \(\omega_n^{\frac n2}=-1\);
- \(\omega_{\frac n2}^{\frac k2}=\omega_n^k\)，其中 \(k\) 是偶数。
这让我们想到原根，设 \(g\) 是质数 \(p\) 的原根，我们令 \(g^{\frac{p-1}{n}}\) 为 \(n\) 次单位根，需要证明满足上述性质。
::::success[证明：]
第一点：\((g^{\frac{p-1}{n}})^n=g^{p-1}\equiv 1\pmod p\)。
第二点：设 \(0\le k,l<n\) 且 \(k\neq l\)，假设存在 \(k,l\) 使 \(g^{\frac{k(p-1)}{n}}\equiv g^{\frac{l(p-1)}{n}}\pmod p\)。
则由原根的性质得：\(k\frac{p-1}{n}\equiv l\frac{p-1}{n}\pmod {p-1}\)。
然后就有 \(k\equiv l\pmod{\dfrac{p-1}{(p-1,\dfrac{p-1}{n})}}\)，即 \(k\equiv l\pmod n\)，又有 \(0\le k,l<n\)，则 \(k=l\)，矛盾！
第三点：\(((g^{\frac{p-1}{n}})^{\frac n2})^2=g^{p-1}\equiv 1\pmod p\)，由二次剩余的知识知 \((g^{\frac{p-1}{n}})^{\frac n2}=\pm 1\)，但若 \((g^{\frac{p-1}{n}})^{\frac n2}=1\)，则 \(g^{\frac{p-1}{2}}=1\)，与 \(g\) 是原根矛盾！所以 \((g^{\frac{p-1}{n}})^{\frac n2}=-1\)。
第四点：\((g^{\frac{p-1}{\frac n2}})^{\frac k2}=(g^{\frac{p-1}{n}})^k\)。\(\Box\)
::::
这就说明了我们可以用原根替代本原单位根，我们可以取 \(g^{\frac{p-1}{n}}\) 为 \(n\) 次单位根，这需要保证 \(n\mid (p-1)\)，其中 \(n\) 是 \(2\) 的幂，我们可以让 \(p\) 取 \(998244353=119\times 2^{23}+1\)，就可以满足上述条件，而且记住 \(998244353\) 的原根为 \(3\)。这样我们就可以做 NTT 了，代码如下：
#include <bits/stdc++.h>
using namespace std;
const int kMaxN = 2.2e6 + 5, kM = 998244353;
int n, m, a[kMaxN], b[kMaxN], r[kMaxN];
int Fpow(int a, int b, int p, int r = 1) {
for (; b; a = 1LL * a * a % p, b >>= 1) {
if (b & 1) r = 1LL * r * a % p;
}
return r % p;
}
void Init(int lg) {
for (int i = 0; i < 1 << lg; i++) {
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (lg - 1));
}
}
void NTT(int n, int *a, int g) {
for (int i = 0; i < n; i++) {
if (i < r[i]) swap(a[i], a[r[i]]);
}
for (int d = 2, w; d <= n; d *= 2) {
w = Fpow(g, n / d, kM);
for (int i = 0; i < n; i += d) {
for (int j = i, t = 1, x, y; j < i + d / 2; j++) {
x = a[j], y = 1LL * t * a[j + d / 2] % kM;
a[j] = (x + y) % kM, a[j + d / 2] = (x - y + kM) % kM;
t = 1LL * t * w % kM;
}
}
}
}
void Idft(int n, int *a, int g, int invn = 0) {
NTT(n, a, Fpow(g, kM - 2, kM)), invn = Fpow(n, kM - 2, kM);
for (int i = 0; i < n; i++) {
a[i] = 1LL * a[i] * invn % kM;
}
}
int W(int n) { return powl(2, ceill(log2l(n))); }
int main() {
ios::sync_with_stdio(0), cin.tie(0);
cin >> n >> m, n++, m++;
for (int i = 0; i < n; i++) {
cin >> a[i];
}
for (int i = 0; i < m; i++) {
cin >> b[i];
}
int N = W(n + m), g;
Init(__lg(N)), g = Fpow(3, (kM - 1) / N, kM);
NTT(N, a, g), NTT(N, b, g);
for (int i = 0; i < N; i++) {
a[i] = 1LL * a[i] * b[i] % kM;
}
Idft(N, a, g);
for (int i = 0; i < n + m - 1; i++) {
cout << a[i] << ' ';
}
return 0;
}
这是 NTT 中可能会用到的素数和其原根，有需要的同学可以看一看：https://blog.miskcoo.com/2014/07/fft-prime-table。
FFT 和 NTT 都需要理解透彻，这样可以帮助你卡常。
多项式牛顿迭代：
多项式牛顿迭代过程：
多项式牛顿迭代可以用来处理下面的问题。
已知函数 \(G(x)\)，求一个多项式 \(F(x)\bmod x^n\)，满足方程 \(G(F(x))\equiv 0\pmod {x^n}\)。
当 \(n=1\)，需要单独求出满足 \(G(F(x))\equiv 0\pmod x\) 的 \(F(x)\)。
接下来假设我们已经求出了 \(F_0(x)\) 使得 \(G(F_0(x))\equiv 0\pmod {x^{\lceil\frac n2\rceil}}\)，考虑 \(G(F(x))\) 在 \(F_0(x)\) 上的泰勒展开：
由于 \(F_0(x)\equiv F(x)\pmod {x^{\lceil\frac n2\rceil}}\)，所以有 \((F(x)-F_0(x))^2\equiv 0\pmod {x^{2\lceil\frac n2\rceil}}\)，也就有 \((F(x)-F_0(x))^2\equiv 0\pmod {x^n}\)，更高次同理。
所以就有 \(G(F(x))\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))\pmod {x^n}\)，由于 \(G(F(x))\equiv 0\pmod {x^n}\)，转化得：\(F(x)\equiv F_0(x)-\dfrac{G(F_0(x))}{G'(F_0(x))}\pmod {x^n}\)。
有了这个式子，这道题就做完了。
下面我们来看几道例题：
多项式求逆：
给定一个多项式 \(F(x)\)，请求出一个多项式 \(G(x)\)，满足 \(F(x)\ast G(x)\equiv 1\pmod {x^n}\)。
设 \(H(G(x))=\dfrac{1}{G(x)}-F(x)\)，我们需要求出其零点多项式。
当 \(n=1\) 时，\(G(x)\equiv \dfrac{1}{F(x)}\pmod x\)，也即 \(F(x)\) 常数项的逆元。
假设已经求出满足 \(H(G_0(x))\equiv 0\pmod {x^{\lceil\frac n2\rceil}}\) 的 \(G_0(x)\)。因为 \(H'(G(x))=-\dfrac{1}{G(x)^2}\)（可以将 \(F(x)\) 看作常数项），所以 \(G(x)\equiv G_0(x)-\dfrac{\dfrac{1}{G_0(x)}-F(x)}{-\dfrac{1}{G_0(x)^2}}\pmod {x^n}\)，也即 \(G(x)\equiv G_0(x)(2-G_0(x)F(x))\pmod {x^n}\)，在每一层用多项式乘法即可。
其时间复杂度为 \(T(n)=T(\dfrac n2)+O(n\log n)=O(n\log n)\)。
代码：https://www.luogu.com.cn/record/198798739。
多项式 exp：
给出 \(n-1\) 次多项式 \(A(x)\)，求一个 \(\bmod{\:x^n}\) 下的多项式 \(B(x)\)，满足 \(B(x) \equiv \text e^{A(x)}\)。这里保证 \([x^0]A(x)=0\)。
设 \(F(B(x))=\ln B(x)-A(x)\)，则 \(F'(B(x))=\dfrac{1}{B(x)}\)，我们需要求出其零点多项式。
当 \(n=1\) 时，\(B(x)=1\)。
假设已经求出满足 \(F(B_0(x))\equiv 0\pmod {x^{\lceil\frac n2\rceil}}\) 的 \(B_0(x)\)。那么 \(B(x)\equiv B_0(x)-\dfrac{\ln B_0(x)-A(x)}{\dfrac{1}{B_0(x)}}\pmod {x^n}\)，即 \(B(x)\equiv B_0(x)(1-\ln B_0(x)+A(x))\pmod {x^n}\)，在每一层用多项式乘法和多项式 ln 即可。
其时间复杂度为 \(T(n)=T(\dfrac n2)+O(n\log n)=O(n\log n)\)。
这里简单说一下为什么要保证 \([x^0]A(x)=0\)，否则 \([x^0]e^{A(x)}=e^{A(0)}=e^{[x^0]A(x)}\) 是个超越数，无法在模 \(p\) 意义下表达。
这也说明了我们需要注意多项式多项式变换的常数项（后面的求导再积分也是）。
代码：https://www.luogu.com.cn/record/198806655。
多项式牛顿迭代还可以求多项式开根，原理都是一样的，这里就不再赘述。
多项式初等函数：
多项式求逆：
使用多项式牛顿迭代可解决。
多项式开方：
使用多项式牛顿迭代可解决。
多项式 ln：
由 \(\ln\) 的定义知，\([x^0]f(x)=1\)，否则没有定义，所以 \([x^0]\ln f(x)=0\)。
考虑对 \(\ln f(x)\) 求导后再积分：
多项式求导和积分很好处理，再用多项式求逆可解决。
注意 \(\ln f(x)\) 的常数项为 \(0\)。
多项式 exp：
使用多项式牛顿迭代可解决。
多项式三角函数：
需要保证 \([x^0]f(x)=0\)。
考虑使用欧拉公式 \(e^{ix}=\cos x+i\sin x\)。
那么就有：
可以直接使用多项式 exp 求解（\(i\) 其实就是 \(\omega_4\)）。
而其他的三角函数可以直接用 \(\sin,\cos\) 表示，都可以用多项式求逆得到。
多项式反三角函数：
对于 \(\arcsin\)，需要保证 \([x^0]f(x)=0\)。
对于 \(\arccos\)，需要保证 \([x^0]f(x)=1\)。
对于 \(\arctan\)，需要保证 \([x^0]f(x)=0\)。
将 \(\arcsin,\arccos,\arctan\) 求导再积分可得：
多项式快速幂：
\(f(x)^k=\exp(k\ln f(x))\)，但是这里需要保证 \([x^0]f(x)=1\)。
如果 \([x^0]f(x)\) 不等于 \(1\) 且不等于 \(0\)，我们可以将整个多项式除以常数项最后再乘上。
如果 \([x^0]f(x)\) 等于 \(0\)，我们可以将 \(x^k\) 提出最后再乘上，其中 \(k\) 是满足 \([x^k]f(x)\neq 0\) 的最小的 \(k\)。
时间复杂度：\( O(n\log n)\)（可是它有时常数太大，比 \(O(n\log n\log k)\) 的普通快速幂还慢）
多项式带余除法：
其实就是初中学的大除法。我们回忆一下这是怎么做的，将多项式次数从高到低写下来，然后类似整数除法列出竖式，每次将被除多项式的最高次幂解决，知道被除多项式的次数小于 \(m\) 即为余数多项式。
类似与整数的除法，如果不考虑余数，将这个竖式一直列下去就会得到这个分数的具体值。类比到多项式，那么多项式反转后求逆在截断就会得到商多项式的反转。
即若 \(F(x)=Q(x)G(x)+R(x)\)，则：
得出 \(Q(x)\) 后就很容易得出 \(R(x)\) 了。
下面具体证明一下这个结论。
::::success[证明：]
首先 \(F^R(x)=x^{\deg F}F\left(\dfrac 1x\right)\)，然后直接硬推。
\(\Box\)
::::
分治 FFT / NTT：
分治 FFT / NTT 可以处理一类半在线卷积的问题。
例题 P4721【模板】分治 FFT：
给定序列 \(g_{1\dots n - 1}\)，求序列 \(f_{0\dots n - 1}\)。
其中 \(f_i=\sum_{j=1}^if_{i-j}g_j\)，边界为 \(f_0=1\)。
多项式求逆做法：
注意到 \(F=FG+1\)，所以 \(F=(1-G)^{-1}\)，直接多项式求逆即可。
时间复杂度 \(O(n\log n)\)。
分治 FFT / NTT 做法：
我们考虑 cdq 分治，现在已经计算出 \([1,mid]\) 的 \(f\) 值，也计算出了 \([1,l-1]\) 对 \([mid+1,r]\) 的贡献，需要处理 \([l,mid]\) 对 \([mid+1,r]\) 的贡献。
直接将 \(f\) 的 \([l,mid]\) 这一段与 \(g\) 的 \([0,r-l]\) 这一段进行卷积，就可以算出贡献了，然而 \([l,mid]\) 已经计算出来了，所以我们将半在线卷积转换成了 NTT。
这两种做法都可以做，也说明了多项式求逆可以用分治 FFT / NTT 实现。
多项式连续点值平移：
给定次数为 \(n\) 的多项式 \(F(x)=\sum_{k=0}^nf_kx^k\)，你需要求出 \(G(x)=F(x+c)\) 的每一项系数。
仍旧暴力推式子：
系数是差卷积的形式，直接使用 \(O(n\log n)\) 的多项式乘法即可。
Chirp Z 变换：
给定一个 \(n\) 项多项式 \(F(x)=\sum_{k=0}^nf_kx^k\) 以及 \(c,m\)，请计算 \(F(c^0),F(c^1),\cdots,F(c^{m-1})\)。
多项式多点求值的特殊情况。
先带入求得 \(F(c^j)=\sum_{k=0}^nf_kc^{jk}\)。
\(jk\) 的形式不好直接卷积，注意到 \(jk={j+k\choose 2}-{j\choose 2}-{k\choose 2}\)，将 \(jk\) 变成了和的形式，那么就好卷积了。
这是也是差卷积的形式，直接 \(O((n+m)\log (n+m))\) 卷积即可。
常系数齐次线性递推：
给出初值 \(a_0,a_1,\cdots,a_{k-1}\) 和转移系数 \(f_1,f_2,\cdots,f_k\)，满足递推式 \(a_n=\sum_{i=1}^kf_ia_{n-i}\)。你需要快速求出 \(a_n\)。
多项式取模做法：
我们先考虑这样的暴力算法，我们要求 \(a_n\)，用递推式展开得 \(a_n=f_1a_{n-1}+f_2a_{n-2}+\cdots+f_ka_{n-k}\)。
一直对 \(a_i(i\ge k)\) 用递推式展开，就可以得到 \(a_n\) 等于 \(f_1,f_2,\cdots,f_k\) 的一个线性组合，求出这个系数就做出了这道题。
我们考虑将 \(a_n\) 写成 \(x^n\)，那么递推式就表示 \(x^n\) 变成 \(f_1x^{n-1}+f_2x^{n-2}+\cdots+f_kx^{n-k}\)。
想想这个东西像什么？
我们考虑 \(f\) 的特征多项式 \(F(x)=x^k-f_1x^{k-1}-f_2x^{k-2}-\cdots-f_k\)。那么每次操作就相当于将 \(x^n\gets x^n-x^{n-k}F(x)\)，最终会变成一个次数小于 \(k\) 的多项式。我们惊人的发现，这就是多项式取模，最后得到的线性组合的系数就是 \(x^n\bmod F(x)\)。
由于 \(n\) 很大，对它做快速幂即可。
时间复杂度 \(O(k\log k\log n)\)，但是因为需要多项式取模，常数巨大。
Bostan-Mori 算法：
Bostan-Mori 算法可以用来求解分式远处系数。
首先，根据之前的做法，常系数齐次线性递推数列的生成函数可以写作一个有理函数。
设 \(a_n\) 的生成函数为 \(A(x)\)，我们现在要构造 \(\frac{P(x)}{Q(x)}=A(x)\) 即 \(P(x)=Q(x)A(x)\)，我们考虑 \(n>\deg P\) 时 \([x^n]P(x)=0\)。
即：
直接令 \(q_0=1,q_i=-f_i(1\le i\le k)\) 即可满足上述条件，因为当 \(n>\deg P\) 时都要满足这个式子，我们尽量让 \(\deg P\) 小，所以取 \(P(x)=Q(x)A(x)\bmod x^k\) 即可。
所以题目便转化为求 \([x^n]\frac{P(x)}{Q(x)}\)，其中 \(n\) 很大。
将分式上下同乘 \(Q(-x)\)，那么分母中就会出现平方差公式，得到的多项式就没有了奇次项，我们设 \(V(x^2)=Q(x)Q(-x)\)。
由于分母没有奇次项，那么其逆元也没有。
所以我们将分子 \(P(x)Q(-x)\) 写成 \(S_1(x^2)+xS_2(x^2)\) 分别表示偶次项和奇次项。
那么答案就可以写成 \([x^n]\frac{S_1(x^2)}{V(x^2)}+[x^{n-1}]\frac{S_2(x^2)}{V(x^2)}\)，其中肯定有一个为 \(0\)，我们根据 \(n\) 的奇偶选择递归计算哪一边。
显然复杂度就是 \(O(k\log k\log n)\)，由于一轮只需要用到 5 次 DFT，常数其实很小。
多项式多点求值：
分治求解：
注意到 \(f(a)=f(x)\bmod (x-a)\)。
我们先利用多项式取模求出 \(g(x)=f(x)\bmod\prod_{i=0}^m(x-a_i)\)，接着递归计算 \(h(x)=f(x)\bmod\prod_{i=0}^{\lfloor\frac m2\rfloor}(x-a_i)=g(x)\bmod\prod_{i=0}^{\lfloor\frac m2\rfloor}(x-a_i)\)。
区间右半部分同理。
那么时间复杂度为 \(T(m)=2T(\frac m2)+O(m\log m)=O(m\log^2 m)\)。
转置原理：
这个作者还不会……
多项式快速插值：
这个作者还不会……
浙公网安备 33010602011771号