下面的所有知识均属多项式范畴。

---

前言

本文略去多项式的定义，毕竟初中数学都讲过。

本文会出现一些可能没有见过的表达：

• 对于一个关于 $x$ 的多项式 $f(x)$，称其最高次项次数为 $f(x)$ 的度数，记做 $\mathrm{deg}f$。
• 多项式同样定义了四则运算，同时也衍生出了同余符号 $\equiv$、对数函数 $\ln (f(x))$，指数函数 $\exp (f(x))$ 等等。
• 下文一般采用 $f(x)\mathrm{m}\mathrm{o}\mathrm{d}{x}^n$ 来表示多项式 $f(x)$ 的所有次数小于 $n$ 的项之和。

注意：本文的知识环环相扣，请勿乱序阅读。

---

拉格朗日插值

插值，通俗地讲，就是指通过多项式在某些点的值求出整个多项式的表达式。

而拉格朗日插值就是计算插值的一种方法，它支持通过 $n+1$ 个点插出一个度数为 $n$ 的多项式。

按惯例，献上模板题。

首先对于多项式 $f(x)$，必然满足 $f(x)\equiv f(a)(\mathrm{m}\mathrm{o}\mathrm{d}(x-a))$。

证明：

转化为 $f(x)-f(a)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}(x-a))$。

设左侧的多项式为 $g(x)$，显然其常数项为 $0$。

又因为其第 $i$ 项为 $p_i(x^i-a^i)$（$p_i$ 为系数），并且 $x^i-a^i$ 必定可以分解为 $(x-a)\cdot q(x)$ 的形式（$q(x)$ 为另一个多项式），所以得证。

所以我们可以得到一个关于 $f(x)$ 线性同余方程组：
$$\{\begin{array}{c}f(x)\equiv y_1(\mathrm{m}\mathrm{o}\mathrm{d}(x-x_1))\\ f(x)\equiv y_2(\mathrm{m}\mathrm{o}\mathrm{d}(x-x_2))\\ f(x)\equiv y_3(\mathrm{m}\mathrm{o}\mathrm{d}(x-x_3))\\ ⋮\\ f(x)\equiv y_n(\mathrm{m}\mathrm{o}\mathrm{d}(x-x_n))\end{array}$$
这里可以使用中国剩余定理解方程组（这里就没有必要关注它们是否互质了）。

我们可以设 $M={\prod }_{i=1}^n(x-x_i)$，$M_i=\frac{M}{(x-x_i)}={\prod }_{j\ne i}(x-x_j)$。

则 $M_i$ 在模 $x-x_i$ 意义下的逆元为 ${\prod }_{j\ne i}\frac{1}{x_i-x_j}$。

（因为 $M_i{M}_i^{-1}={\prod }_{j\ne i}\frac{x-x_j}{x_i-x_j}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}(x-x_i))$，则对于每一个 $j\ne i$ 都存在 $ x - x_j \equiv x_i - x_j~(\bmod (x - x_i))$，移项可得 $x-x_i\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}(x-x_i))$，这显然是成立的。）

所以有 $f(x)={\sum }_{i=1}^n{y}_i{\prod }_{j\ne i}\frac{x-x_j}{x_i-x_j}(\mathrm{m}\mathrm{o}\mathrm{d}M)$，它在模意义下是唯一的。

（实际上有另外的证明方式证明了其在通常意义下也是唯一的，这里略去。）

本题中，我们只需要代入 $k$ 求解即可，时间复杂度为 $O(n^2)$。不过更多时候，我们需要求出其所有系数，也可以做到 $O(n^2)$，甚至还存在更低的复杂度，这里就不展开介绍了。

代码就不贴了。

---

多项式乘法

多项式乘法的定义，这里就不必多言了，初中数学都讲过。

代数形式：

设有多项式 $f(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}_n$，$g(x)=b_0+b_{1}x+b_2{x}^2+\cdots +b_m{x}_m$。

则它们的乘法积为 $t(x)=c_0+c_{1}x+c_2{x}^2+\cdots +c_{n+m}{x}^{n+m}$。

其中，对于 $t(x)$ 的系数 $c_k$ 有公式：$c_k={\sum }_{i=0}^k{a}_i{b}_{k-i}$。

可以发现这个公式的形式十分类似卷积。但请注意：它们不是同一个东西，不要弄混淆了。

事实上，多项式乘法就是用于求卷积的。如果您对卷积感兴趣，建议百度。

---

现在，给出两个多项式 $f(x)$ 和 $g(x)$ 的系数，求出它们的乘法积 $t(x)$。$f(x)$ 和 $g(x)$ 的次数小于 $10^5$。

根据上面的公式，我们可以得到一个 $O(n^2)$ 的优秀算法。

我们考虑换一种思想计算这个乘法积。

众所周知，有这样一条著名的定理：任意 $k$ 个点可以唯一确定一个至多为 $k-1$ 次的多项式。

所以我们得到了一个新的思路：先求出足够多的 $f(x)$ 和 $g(x)$ 的点值，将它们用 $O(n)$ 的复杂度乘起来（横坐标不变，纵坐标相乘），再将这些点转化为一个多项式（上面提到了这样的过程叫做插值）。

这种方法叫做离散傅里叶变换（即著名的 $\mathrm{DFT}$）。

问题是，如果钦定一些 $x$ 去求对应的 $f(x)$ 的值，时间复杂度还是 $O(n^2)$ 的。

所以，下面介绍快速傅里叶变换，同时会介绍快速数论变换（即著名的 $\mathrm{FFT}$ 和 $\mathrm{NTT}$）。

---

前置知识：单位根，原根。

还记得二次函数 $y=x^2$ 的性质吗？

它关于 $y$ 轴对称，所以我们如果代入 $x=k$，同样可以得到 $x=-k$ 时 $y$ 的值。

这启发我们，可以选取一些特殊的 $x$，使得我们可以在求出一些点值的同时，快速地求出另外的一些点值。

这个时候，单位根就派上用场了。

先讲一些使得单位根能够应用于此的性质：

1. 记 ${\omega }_n$ 为 $n$ 次单位根，即满足 ${\omega }_n^k=e^{\frac{2k\pi i}{n}}$，而 $e^{\frac{2k\pi i}{n}}=\cos \frac{2k\pi }{n}+i\sin \frac{2k\pi }{n}$，$i$ 为虚数单位。

（注意 $n$ 次单位根的 $n$ 一般取 $2$ 的正整数次幂，而下文皆如此。）

这个计算公式不仅可以让我们计算出单位根的具体数值，并且可以得出以下性质：

• ${\omega }_n^n=\cos 2\pi +i\sin 2\pi =1={\omega }_n^0$。

这条性质说明 $n$ 次单位根存在长度为 $n$ 的指数循环节；

• ${\omega }_n^{\frac{n}{2}}=\cos \pi +i\sin \pi =-1$。

这条性质间接地证明了 $n$ 次单位根的最小指数循环节长度为 $n$，同时为下文折半引理的证明做了铺垫。

2. 单位根消去引理：${\omega }_{2n}^{2k}={\omega }_n^k$。

根据公式有 ${\omega }_{2n}^{2k}=e^{\frac{4k\pi i}{2n}}=e^{\frac{2k\pi i}{n}}={\omega }_n^k$，得证。

这条性质可以做为沟通不同次单位根的桥梁，同时也证明了 $n$ 次单位根的值域是 $2n$ 次单位根值域的真子集。

3. 单位根折半引理：${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$。

根据上面的性质有 ${\omega }_n^{k+\frac{n}{2}}={\omega }_n^k\cdot {\omega }_n^{\frac{n}{2}}=-{\omega }_n^k$，得证。

这条性质可以做为沟通同次不同指数的单位根的桥梁，为 $\mathrm{DFT}$ 的分治优化创造了可能。

而我们该如何运用单位根的这些性质呢？

我们将 $f(x)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3$ 做一个小小的改动：
$$\begin{array}{rl}f(x)& =a_0+a_{1}x+a_2{x}^2+a_3{x}^3\\ & =a_0+a_2{x}^2+a_{1}x+a_3{x}^3\\ & =(a_0+a_2{x}^2)+x(a_1+a_3{x}^2)\\ & =f_1(x^2)+x{f}_2(x^2)\end{array}$$
可以发现，这样的变换，会使得原多项式被拆分为两个规模缩小了一半多项式。

所以，我们又有了一个新的算法模型：

我们考虑将 $f(x)$ 拆分为两个规模缩小了一半的多项式 $f_1(x),f_2(x)$，然后对它们进行分治，接着，我们可以用 $f_1(x),f_2(x)$ 在 $x=k^2$ 处的点值计算出 $f(x)$ 在 $x=k$ 处的点值。如果分治时每一层的复杂度都能够做到 $O(n)$，那么总的时间复杂度就是 $O(n{\log }_{2}n)$。

接下来，我们考虑利用单位根的性质，将每一层转化点值的操作复杂度变为 $O(n)$。

如果我们将单位根 ${\omega }_n^k$ 和 ${\omega }_n^{k+\frac{n}{2}}$ 分别代入多项式 $f(x)$：
$$\begin{gathered} f({\omega }_n^k)=f_1({\omega }_n^{2k})+{\omega }_n^k{f}_2({\omega }_n^{2k})=f_1({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k{f}_2({\omega }_{\frac{n}{2}}^k) \\ f({\omega }_n^{k+\frac{n}{2}})=f_1({\omega }_n^{2k+n})+{\omega }_n^{k+\frac{n}{2}}{f}_2({\omega }_n^{2k+n})=f_1({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k{f}_2({\omega }_{\frac{n}{2}}^k) \end{gathered}$$
可以发现，只要我们求出了 $f_1(x),f_2(x)$ 在 $x={\omega }_{\frac{n}{2}}^k$ 处的点值，我们就可以直接得出 $f(x)$ 在 $x={\omega }_n^k,x={\omega }_n^{k+\frac{n}{2}}$ 两处的点值！

具体来说，我们只要在所有长度为 $n$ 的多项式中求出在所有 $n$ 次单位根处的点值即可。

这样，做到 $O(n)$ 的复杂度已经易如反掌了。

我们再来理一遍这个算法吧：（同时加入了一些很多细节）

1. 我们找到一个 $p$，使得 $p$ 是 $2$ 的正整数次幂且 $p\ge n+m+1$。

这一步操作之前没有提到，但是很有必要，它为分治算法的进行提供了方便，这样我们就无需考虑多项式的度数；同时，这一步操作不会影响到点值的大小，也不会影响到最终插出的多项式的度数 （想一想，为什么）。

代码：（$q$ 的作用下文会提及）

int p = 1, q = -1;
while(p <= n + m) ++q, p <<= 1;

2. 接着，我们对两个多项式 $f(x)$ 和 $g(x)$ 分别跑 $\mathrm{FFT}$。在 $\mathrm{FFT}$ 运行过程中，如果当前多项式长度为 $1$，直接返回系数；否则，将当前的多项式 $f(x)$ 按照系数底数的奇偶性将项分类（这里请读者自己想一想为什么不可以直接对半分类），将其拆分为两个多项式，然后分治。

但是，如果这一步直接采用递归，时间似乎会有些吃不消。

我们不妨来探讨一下系数究竟是怎样分配的：（以 $n=8$ 为例）
{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\} \\ \downarrow \\ \{a_0,a_2,a_4,a_6\}~\{a_1,a_3,a_5,a_7\} \\ \downarrow \\ \{a_0, a_4\}~\{a_2, a_6\}~\{a_1, a_5\}~\{a_3, a_7\} \\ \downarrow \\ \{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​}
把每个系数的底数分离：

$$0,4,2,6,1,5,3,7$$

转换成二进制：

$$000,100,010,110,001,101,011,111$$

我们再把它们翻转过来：
$$000,001,010,011,100,101,110,111$$

这不就是从小到大排序的吗？

所以我们就有了如下代码：（设 $r_i$ 为分配完系数后第 $i$ 个位置上的系数）

for(register int i = 0; i < p; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << q);

解释一下：第 $i$ 个位置的 $r$，不就等于第 $\lfloor \frac{i}{2}\rfloor$ 个位置的 $r$ 右移一位吗？

如果 $i$ 是奇数，我们再在最高位补上一个 $1$ 即可（$q$ 的作用就在这里）。

以上过程称为 $\mathrm{Rader}$ 排序。

3. 用拆分出的两个多项式的点值，求出当前多项式的点值。

代码：（$\mathrm{Complex}$ 为复数类，定义了构造函数，重定义了乘法）

inline void FFT(Complex *f)
{for(register int i = 0; i < p; i++) if(i < r[i]) swap(f[i], f[r[i]]);for(register int len = 2; len <= p; len <<= 1){int k = len >> 1;Complex _w(cos(PI / k), sin(PI / k));for(register int i = 0; i < p; i += len){Complex w(1.0, 0.0);for(register int j = 0; j < k; j++){Complex u = f[i + j], v = w * f[i + j + k];f[i + j] = u + v;f[i + j + k] = u - v;w = w * _w;}}}return;
}

4. 接着，我们将 $f$ 和 $g$ 对应的点值相乘，时间复杂度为 $O(n)$。

你以为结束了？并没有，我们还有第五步：

5. 最后，我们将多项式插值出每一项的系数，算法结束。我们称这一步叫 $\mathrm{IFFT}$。

“算法结束” 四个字就在这里，近在咫尺，却又远在天边。

如果用拉格朗日插值计算，时间复杂度还是 $O(n^2)$，功亏一篑。

难道前面的分析过程都只能打水漂？

不不不，我们冷静分析一波。

看看下面这张图：

多项式乘法的过程实际上类似上图中的矩阵乘法。

那么我们想得到后面的系数矩阵，不就只要把中间的这个矩阵求逆就行了吗！

而这个逆矩阵非常特殊：

它等于原矩阵的每个值求倒数，再乘上一个 $\frac{1}{n}$。

（然而笔者线性代数太菜了，无法解释其原因。）

既然无法解释，那就考虑推式子证明吧：

设 $y_k$ 是 $f(x),g(x)$ 在傅里叶变换后，在 $x={\omega }_p^k$ 处点值相乘的结果。

我们设多项式 $s(x)=y_0+y_{1}x+y_2{x}^2+\cdots +y_{p-1}{x}^{p-1}$ 在 $x={\omega }_n^0,{\omega }_n^{-1},\cdots \cdots ,{\omega }_n^{-(p-1)}$ 处的点值分别为 $u_{0\sim p-1}$，即满足：
$$u_k=\sum _{i=0}^{p-1}{y}_i({\omega }_p^{-k}{)}^i$$
现在我们要用它证明等式 $c_k=\frac{u_k}{p}$ 成立，其中 $c_k$ 指最终 $f(x)$ 和 $g(x)$ 的乘法积的系数。

我们直接莽式子：
$$\begin{array}{rl}{u}_k& =\sum _{i=0}^{p-1}{y}_i({\omega }_p^{-k}{)}^i\\ & =\sum _{i=0}^{p-1}(\sum _{j=0}^{p-1}{c}_j({\omega }_p^i{)}^j)({\omega }_p^{-k}{)}^i\\ & =\sum _{i=0}^{p-1}\sum _{j=0}^{p-1}{c}_j({\omega }_p^{j-k}{)}^i\\ & =\sum _{j=0}^{p-1}{c}_j\sum _{i=0}^{p-1}({\omega }_p^{j-k}{)}^i\end{array}$$
然后对 ${\sum }_{i=0}^{p-1}({\omega }_p^{j-k}{)}^i$ 这一部分分讨，设 $v=j-k$，则该式等于 ${\sum }_{i=0}^{p-1}({\omega }_p^v{)}^i$。

在这个式子的意义下，这就是一个等比数列求和。

如果公比 ${\omega }_p^v=1$，即 $v=0$，则此时原式等于 $p$；

否则，套用等比数列公式，原式等于：
$$\begin{array}{rl}& \frac{({\omega }^v{)}^p-1}{{\omega }_p^v-1}\\ =& \frac{({\omega }^p{)}^v-1}{{\omega }_p^v-1}\\ =& \frac{1-1}{{\omega }_p^v-1}\\ =& 0\end{array}$$
故当且仅当 $v=j-k=0$，即 $j=k$ 时，后面的这一部分才会产生贡献，则有 $u_k=p{c}_k$，即 $c_k=\frac{1}{p}{u}_k$。

得证。

然后可以发现，$\mathrm{FFT}$ 和 $\mathrm{IFFT}$ 的代码可以合为一份：

inline void FFT(Complex *f, double t)
{for(register int i = 0; i < p; i++) if(i < r[i]) swap(f[i], f[r[i]]);for(register int len = 2; len <= p; len <<= 1){int k = len >> 1;Complex _w(cos(PI / k), t * sin(PI / k));for(register int i = 0; i < p; i += len){Complex w(1.0, 0.0);for(register int j = 0; j < k; j++){Complex u = f[i + j], v = w * f[i + j + k];f[i + j] = u + v;f[i + j + k] = u - v;w = w * _w;}}}for(register int i = 0; i < p; i++) f[i].r = (int)(f[i].r / p + 0.5);return;
}

当运行 $\mathrm{FFT}$ 时，我们令 $t=1$ 即可；当运行 $\mathrm{IFFT}$ 时，我们令 $t=-1$ 即可。

（注意：由于 $\cos$ 和 $\sin$ 的精度本身就十分地不稳定，加上复数乘法的运算次数较多，我们最后计算答案的时候需要处理这个误差。）

这样，算法真正大功告成啦！

---

优化：

可以发现，对于一个复数 $a+bi$，有 $(a+bi{)}^2=(a^2-b^2)+2abi$。

如果我们把 $g(x)$ 放到 $f(x)$ 的虚部上，对 $f(x)$ 跑一次 $\mathrm{FFT}$，然后将 $f(x)$ 的每个点值变为其平方。那么，虚部除以 $2$ 不就是我们想要的值吗？

这样，我们就可以只跑一遍 $\mathrm{FFT}$，一遍 $\mathrm{IFFT}$ 即可，常数优化至原来的 $\frac{2}{3}$。

---

然而，$\mathrm{FFT}$ 的精度仍然是一个令人头疼的问题。我们想考虑一种方式，使得整个运算建立在整数域上。

整数，想到数论，有一个东西恰好符合其所有要求——原根。

关于原根的定义和性质，可以参考我的博客整除和同余相关理论。

可以发现，使得单位根可用于 $\mathrm{FFT}$ 的性质，原根同样满足：

1. $n$ 次原根计算公式为 $g_n^k=g^{\frac{k(p-1)}{n}}\mathrm{m}\mathrm{o}\mathrm{d}p$，其中 $p$ 为模数，$g$ 特指模 $p$ 的最小原根。

这里的 $p$ 一般取 $998244353$，它的最小原根等于 $3$，并且有 $998244353=2^{23}\times 7\times 17+1$，其中 $2^{23}=8388608$，其大小完全可以满足我们的需求。

同样，我们也可以得出 $g_n^n=g_n^0\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 和 $g_n^{\frac{n}{2}}\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 以及它们的推论。

2. 原根同样存在消去引理和折半引理。

所以，在 $\mathrm{FFT}$ 中对多项式的变换，在 $\mathrm{NTT}$ 中同样适用。同时，也需要 $\mathrm{Rader}$ 排序。

同理，对于 $\mathrm{INTT}$，原根同样可以在原来的基础上求倒数再乘 $\frac{1}{n}$，只需要求逆元即可。

代码：

#define MOD 998244353
#define G 3
#define InvG 332748118inline int GetMOD(int k) { if(k >= MOD) k -= MOD; return k; }inline int Pow(int k, int p)
{int res = 1;while(p){if(p & 1) res = 1LL * res * k % MOD;k = 1LL * k * k % MOD;p >>= 1;}return res;
}inline void NTT(int *f, bool Inverse = false)
{for(register int i = 0; i < p; i++) if(i < r[i]) swap(f[i], f[r[i]]);for(register int len = 2; len <= p; len <<= 1){int k = len >> 1, g = Pow(Inverse ? InvG : G, (MOD - 1) / len);for(register int i = 0; i < p; i += len){int w = 1;for(register int j = 0; j < k; j++){int u = f[i + j], v = 1LL * w * f[i + j + k] % MOD;f[i + j] = GetMOD(u + v);f[i + j + k] = GetMOD(u - v + MOD);w = 1LL * w * g % MOD;}}}return;
}

---

后记：

强推神佬的视频：快速傅里叶变换($\mathrm{FFT}$)——有史以来最巧妙的算法？

强推神佬的博客：$\mathrm{FFT}$ $\mathrm{NTT}$

---

$ToBeContinued\cdots \cdots$