多项式（Ⅰ）：基础工业

对 Alex_wei 博客的抄写。

复数与单位根

复数

跳出实数域 \(\mathbb{R}\)，定义 \(i^2 = -1\)，即 \(i=\sqrt{-1}\)，并在此基础上定义复数 \(a+bi\)，其中将 \(b\not= 0\) 的称为虚数。复数域记为 \(\mathbb{C}\)。

我们根据上面的定义一下复数的四则运算。

• 加法：\((a+bi)+(c+di) = (a+c)+(b+d)i\)

• 减法：\((a+bi)-(c+di) = (a-c)+(b-d)i\)

• 乘法：\((a+bi)(c+di) = (ac-bd)+(ad+bc)i\)

• 除法：\(\displaystyle \frac{a+bi}{c+di}= \frac{(a+bi)(c-di)}{(c+di)(c-di)}=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i\)

复数的除法其实是一个分母有理化的过程，我们乘上一个分母的共轭复数，就是实部相等，虚部为相反数的复数，来使得分母有理化。

复平面与棣莫弗定理

描述一个复数 \(a+bi\) 需要两个值 \(a\) 和 \(b\)，其实 \(a\) 表示实部，\(b\) 表示虚部。这表明我们可以把它放到平面直角坐标系上进行描述，称为复平面。其在复平面上的坐标为 \((a,b)\)，实部 \(a\) 为横坐标，虚部 \(b\) 为纵坐标。

一个复数可以唯一对应一个复平面上的向量。我们将向量起点平移至远点，那么它的重点就指向与其对应的复数。我们把平面向量的一些性质应用到复平面上，就能得到一些定义：

• 定义复数 \(z=a+bi\) 的模为 \(|z|=\sqrt{a^2+b^2}\)。

• 定义复数 \(z=a+bi\) 的辐角为 \(\text{Arg} \ z=\theta\)，其中 \(\tan \theta=\frac{a}{b}\)。满足 \(-\pi < \theta \leq \pi\) 的 \(\theta\) 称为辐角主值，记作 \(\arg z\)，即 \(\arg z =\arctan\frac{a}{b}\)。说这么多其实就是这个向量与 \(x\) 轴的夹角。

• 辐角确定了 \(z\) 所在的直线，模确定了 \(z\) 在直线上的长度。有了这两个，我们可以用另一种方法来描述复数。

根据 \(z = a+bi\) 的模 \(r\) 和辐角 \(\theta\)，可知 \(z\) 的实部 \(a=r\cos\theta\)，虚部 \(b=r\sin\theta\)，据此我们定义复数的三角形式 \(z = r(\cos\theta + i\sin\theta)\)。

利用 \(\sin,\cos\) 的和角公式可得 \(z_1z_2 = r_1r_2(\cos(\theta_1+\theta_2)+i\sin(\theta_1+\theta_2))\)。该等式称为棣莫弗定理，它说明复数相乘，模长相乘，辐角相加。

根据棣莫弗定理，我们得到对于虚数单位 \(i\) 的一种直观理解：将一个复数 \(z\) 乘以 \(i\) 相当于将其逆时针旋转 \(\frac{\pi}{2}\) 弧度。我们再看 \(i\) 在复平面上是在 \((0,1)\) 的位置，它其实就是让 \(1\) 逆时针旋转了 \(\frac{\pi}{2}\) 弧度。所以说，\(i\) 就代表了旋转。

单位根

当 \(r=1\) 时，\(z=\cos\theta + i\sin\theta\) 在单位圆上。此时根据棣莫弗定理有 \(z^n = \cos(n\theta)+i\sin(n\theta)\)，这个式子很美妙，它在复数旋转和复数乘法之间构建了桥梁：\(z\) 的 \(n\) 次幂相当于从 \((1,0)\) 开始，以 \(z\) 辐角的角度在单位圆上旋转 \(n\) 次得到的结果，称为将 \(z\) 旋转 \(n\) 次。

引用一下 wls 的图片。

在这里，我们将单位圆 \(n\) 等分，取任意 \(n\) 等分点 \(P_k(0\leq k < n)\)，将其旋转 \(n\) 次均得到 \(1\)，也就是 \(P_k^n = 1\)，我们探究一下为什么。

因为 \(P_k\) 相当于从 \(1\) 开始在单位圆上旋转 \(\frac{2k\pi}{n}\) 弧度，因此 \(P_k = \cos(\frac{2k\pi}{n})+i\sin(\frac{2k\pi}{n})\)，你会发现旋转 \(n\) 次后，\(P_k^n = \cos(2k\pi)+i\sin(2k\pi) = 1\)。我们称所有 \(P_k\) 为 \(n\) 次单位根，将 \(P_1\) 记作 \(\omega_n\)，则有 \(P_k = P_1^k = \omega_n^k\)。

根据 \(P_k^n=1\) 我们可知道任意 \(n\) 次单位根 \(\omega_n^k\) 均为 \(x^n = 1\) 的一个解。一般来说，我们一般将 \(n\) 次单位根直接代指 \(\omega_n\)，即从 \(1\) 开始逆时针方向的第一个单位根。

下面我们讨论一下单位根的一些性质，它是我们实现 FFT 的关键：

• 循环性：由 \(\omega_n^n =1\) 可知，\(\omega_n^k\) 的循环节是 \(n\)，这也就是你在复平面上一次转 \(\frac{2\pi}{n}\) 弧度，转 \(n\) 次就会回到原点。所以说，\(\omega_n^k = \omega_n^{k+tn}(t\in \mathbb{Z})\)。

• \(\omega_n^k = \omega_{cn}^{ck}(c > 0)\)，这个放到复平面上更容易理解，你把整个圆分成 \(cn\) 份，它的第 \(ck\) 个单位根就是你把圆分成 \(n\) 份，它的第 \(k\) 个单位根。

• 当 \(n\) 为偶数时，将 \(\omega_n^k\) 取反相当于将其逆时针（或顺时针）转半圈，所以 \(\omega_n^k = -\omega_n^{k\pm\frac{n}{2}}(2\mid n)\)。

• 单位根的对称性：因为 \(n\) 次单位根将单位圆 \(n\) 等分，均匀分布在圆周，所以它们的中心就在原点，即 \(\displaystyle \sum_{i=0}^{n-1}\omega_n^i = 0\)。

• 若 \(\gcd(k,n) = 1\)，则 \(\omega_n^k\) 称为本原单位根。所有本原单位根的 \(0\sim n-1\) 次幂互不相同，你会发现，它似乎和原根有着极大的相似性。

单位根与原根

阐释了单位根的性质以后，我们发现，单位根和原根之间似乎有着奇妙的关联，可以说，原根就是模 \(p\) 意义下整数域的单位根。

设 \(n=\varphi(p)\) 且 \(p\) 存在原根 \(g\)，则 \(g^0,g^1,\dots,g^{n-1},g^n = g^0,g^{n+1}=g^1\) 这样的循环和 \(n\) 次单位根的循环一模一样。这使得在模 \(p\) 意义下涉及 \(n\) 次单位根的运算时，可直接用原根 \(g\) 代替。也就是说，对于 \(d\mid n,g^{\frac{n}{d}}\) 可以直接代替模 \(p\) 意义下的 \(d\) 次单位根。

结合群论来看，单位根和原根都是对应域上一个大小为 \(n=\varphi(p)\) 的循环群的生成元，它们均满足 \(n\) 次幂时对应域的单位元，且它们的 \(0\sim n-1\) 次幂互不相同。换言之，它们同构。　

这就是快速傅里叶变换 FFT 和快速数论变化 NTT 之间的内在联系。

欧拉公式

欧拉公式的内容：

\[e^{ix} = \cos x + i\sin x \]

即单位圆上从 \((1,0)\) 开始旋转 \(x\) 弧度得到的复数，也即大小为 \(x\) 弧度的角的终边与单位圆的交点。

推荐观看 在 3.14 分钟内理解 \(e^{i\pi}\)-3Blue1Brown。

有了这个公式，我们可以更简介的描述一个复数：\(re^{i\theta}\) 表示一个模长为 \(r\)，辐角为 \(\theta\) 的复数，使记号变得更简洁。

将该表示法应用于单位根，可得 \(\omega_n=e^{\frac{2\pi}{n}i}\)。

代入 \(t=\pi\)，得到著名等式

\[e^{i\pi} = -1 \]

代入 \(t = 2\pi = \tau\)，得

\[e^{i\tau} = 1 \]

这说明对于任意 \(k\in\mathbb{Z},(e^{2\pi i})^{k+\frac{t}{n}}\) 相等恰对应 \(\omega_n^{kn+t}\) 相等。

多项式

基本概念

形如 \(\displaystyle \sum_{i=0}^n a_ix^i\) 的有限和式称为多项式。记作 \(\displaystyle f(x)= \sum_{i=0}^na_ix^i\)。

其中，\(a_i\) 称为 \(i\) 次项前的系数，也称 \(x^i\) 前的系数，记作 \([x^i]f(x)\)。超过最高次数 \(n\) 的系数 \(a_i(i>n)\) 视为 \(0\)。

当项数无限时，形如 \(\displaystyle \sum_{i=0}^{\infty}a_ix^i\) 的和式，称为形式幂级数，它在生成函数的部分出现过，在这里我们不予讨论。

• 多项式系数非零的最高次项的次数称为该多项式的度，或者叫次数，记作 \(\deg f\)。

• 使得 \(f(x) = 0\) 的所有 \(x\) 称为多项式的根。

• 若 \(a_i\) 均为实数，则称 \(f\) 为实系数多项式。若 \(a_i\) 可以均为复数，则称 \(f\) 为复系数多项式。

• 代数基本定理：任何非零一元 \(n\) 次复系数多项式恰有 \(n\) 个复数根。这些复数根可能重合。证明？略/cy。

系数表示法和点值表示法

\(\displaystyle f(x)=\sum_{i=0}^na_ix^i\) 这样的式子给出了所有 \(i\) 次项前的系数，这种描述多项式的方法称为系数表示法，这也是我们最常见的表示方法。

将 \(x=x_i\) 代入，得到 \(y_i=f(x_i)\)，称 \((x_i,y_i)\) 为 \(f\) 在 \(x_i\) 处的点值。用若干点值 \((x_i,y_i)\) 描述多项式的方法称为点值表示法。

我们可以探究一下，给出 \(n\) 个 \(x_i\) 互不相同的点值 \((x_0,y_0),\dots,(x_{n-1},y_{n-1})\)，它所唯一确定的多项式的最高次数。

由两点确定一条直线入手，我们猜测是 \(n-1\)，事实上，的确如此。因为 \(n-1\) 次多项式需要 \(n\) 个系数描述，而 \(n\) 个点值恰好就能提供 \(n\) 个信息。

证明就考虑我们把点值代入到 \(n-1\) 次多项式中，得到 \(n\) 元线性方程组，我们把它代入以后，该线性方程组的系数矩阵就是一个 \(n\) 行 \(n\) 列的范德蒙德矩阵，因为 \(x_i\) 互不相同，所以行列式非零，方程组的解唯一。

所以我们得到一个结论：\(n\) 个点值唯一确定的多项式的最高次数为 \(n-1\)。

• 从系数表示法转为点值表示法称为求值。

• 从点值表示法转为系数表示法称为插值。

多项式的运算

多项式运算

设 \(\displaystyle f(x)=\sum_{i=0}^na_ix^i,g(x) = \sum_{i=0}^mb_ix^i\)

系数表示法

• 设 \(h = f+g\)，则

\[h(x)=f(x)+g(x)=\sum_{i=0}^{\max(n,m)}(a_i+b_i)x^i \]

也就是两多项式相加，对应系数相加，且 \(\deg (f+g) = \max(\deg f,\deg g)\)。

• 设 \(h = f * g\)，则

\[h(x)=f(x)g(x) = \sum_{i=0}^{n+m}(\sum_{j=0}^ia_jb_{i-j})x^i \]

也就是两多项式相乘，每两个系数相乘贡献至次数之和的系数，\(\deg(f * g) = \deg f + \deg g\)。

可以看出在系数表示法下，多项式相加的复杂度为 \(O(\max(n,m))\)，多项式相乘的复杂度为 \(O(nm)\)。

点值表示法

• 根据 \((f+g)(x) = f(x)+g(x)\)，可知两个多项式相加时，对应点值相加。

• 根据 \((fg)(x) = f(x)g(x)\)，可知两个多项式相乘时，对应点值相乘。

因此，在点值表示法下，计算两个多项式相加需要 \(\max(n,m)+1\) 个点值，计算两个多项式相乘也需要同样多的点值，复杂度均为 \(O(n+m)\)

注：我们常用 \(f * g\) 和 \(fg\) 表示多项式相乘，即进行加法卷积；用 \(f \cdot g\) 表示多项式点乘，即对应系数相乘。

实际上，FFT 等快速多项式乘法就是运用的这个性质来快速相乘。先转化成点值表示法，相乘后在转回系数表示法，时间复杂度 \(\mathcal{O}((n+m)\log(n+m))\)，我们在后面部分会讲述这个算法。

拉格朗日插值

刚才我们得到一个结论：\(n\) 个点值唯一确定的多项式的最高次数为 \(n-1\)。接下来我们思考如何在点值表示法和系数表示法之间转化。

从系数表示法转为点值表示法，最简单也是最常见的是 \(O(n^2)\) 直接代入 \(n\) 个点每次 \(O(n)\) 求值。\(O(n\log^2n)\) 的多项式多点求值是高级科技，以后我学了会补充上，现在先不讨论。

从点值表示法转为系数表示法，首先因为可以先设出系数表示法的形式，然后把 \(x_i\) 代入进去，这就变成了 \(n\) 个 \(n\) 元一次方程组的形式，显然通过高斯消元就可以 \(O(n^3)\) 的求解此问题了。

还有一种方法，就是我们现在要提到的拉格朗日插值法，它能在 \(O(n^2)\) 的时间复杂度内求解此问题。

算法内容

拉格朗日插值的构造方法其实和中国剩余定理的构造方式有着异曲同工之处。

简单来说，就是我们需要构造一个 \(n\) 次多项式，使得它满足在已知的点 \(x_i\) 处取值为 \(y_i\)。

构造 \(f_i(x)\)，使得其仅在 \(x_i\) 处取值为 \(1\)，其他 \(x_j,j\not= i\) 处取值为 \(0\)。

那么也就是说 \(y_if_i(x)\)，仅在 \(x_i\) 处取值为 \(y_i\)，其他 \(x_j,j\not= i\) 处取值为 \(0\)。

接下来我们考虑如何构造 \(f_i(x)\)：

• 如果要在 \(x_j,j\not=i\) 处取值为 \(0\)，可以让 \((x-x_j),j\not=i\) 乘起来，即 \(\displaystyle \prod_{j\not=i}(x-x_j)\)，我们把它记为 \(g(x)\)。有了这个式子，对于不是 \(x_i\) 的地方，取值就都是 \(0\) 了。

• 如果要在 \(x_i\) 处取值为 \(1\)，那么就让 \(f_i(x) = \frac{g(x)}{g(x_i)}\)，这样当 \(x=x_i\) 时，\(f_i(x)\) 的取值就为 \(1\) 了，因为当 \(i\not=j\) 时，\(x_i\not=x_j\)，所以不会出现除以 \(0\) 的情况。

综上，采取一个类似于 CRT 的构造方法，我们就构造出了唯一的符合要求的 \(n\) 次多项式：

\[\sum_{i=1}^ny_i\prod_{j\not=i}\frac{x-x_j}{x_i-x_j} \]

为得到 \(f\) 的各项系数，我们先 \(O(n^2)\) 的算出 \(\displaystyle F(x) = \prod_{i=0}^{n-1}(x-x_i)\)，对每个 \(i\) 暴力 \(O(n)\) 的除掉一次式 \((x-x_i)\)，算出 \(\frac{F(x)}{x-x_i}\) 的各项系数。这一步的目的就是求分式上面那一部分，然后我们要做的就是在乘以 \(\frac{y_i}{\prod_{j\not=i}x_i-x_j}\) 得到 \(f_i\)，这个分式的下面我们可以 \(O(n^2)\) 预处理出来，最后的答案就是 \(\displaystyle f = \sum_{i=0}^{n-1}f_i\)，最后的时间复杂度就是 \(O(n^2)\)。

通常情况下，以模板题为例，它会让你求出 \(f(x)\) 在给定某个 \(x\) 处的取值，此时我们不把 \(x\) 看成一个未知量，而是直接代入 \(x\)，时间复杂度仍为 \(O(n^2)\)。

signed main()
{
	n = read() , k = read();
	for(re int i=1;i<=n;i++) x[i] = read() , y[i] = read();
	for(re int i=1;i<=n;i++)
	{
		nume = deno = 1;
		for(re int j=1;j<=n;j++)
		{
			if(i == j) continue;
			nume = nume * ((k - x[j] + mod) % mod) % mod;
			deno = deno * ((x[i] - x[j] + mod) % mod) % mod;
		}
		ans = (ans + y[i] * nume % mod * ksm(deno,mod-2) % mod) % mod;
	}
	cout << ans;
	return 0;
}

多项式快速插值在 \(O(n\log^2n)\) 的时间内将点值表示法转化为系数表示法，这个也是以后再讨论。

连续取值插值

对于给定点值横坐标为连续整数时，我们有 \(O(n)\) 插值的方法。

以 \(0\sim n-1\) ，即 \(x_i =i\) 为例：

\[f(x)=\sum_{i=0}^{n-1}y_i\prod_{j\not=i}\frac{x-j}{i-j} \]

• 分子是 \(\prod(x-i)\) 少了一项 \((x-i)\)，我们同时维护一个前缀后缀积即可。设 \(\displaystyle p_i = \prod_{j=0}^{i-1}(x-j),s_i = \prod_{j=i+1}^{n-1}(x-j)\)

• 分母对于 \(\displaystyle i>j,\prod_{j=0}^{i-1}(i-j) =i!\)；对于 \(\displaystyle i<j,\prod_{j=i+1}^{n-1}(i-j) =(-1)(-2)\dots(i-n+1)\)，将所有负号提出来，得到 \((-1)^{n-i-1}(n-i-1)!\)，因此，最终结果就是

\[f(x) = \sum_{i=0}^{n-1}y_i\frac{p_is_i}{(-1)^{n-i-1}i!(n-i-1)!} \]

预处理阶乘逆元，时间复杂度 \(O(n)\)。

快速傅里叶变换

快速傅里叶变换（Fast Fourier Transfrom，FFT）是一切多项式算法的根基。要想真正搞懂它，必须先真正理解单位根和对线性代数有一定的了解。因为我比较菜，所以有些本质性的东西可能会略去，有些地方感性理解即可。

求值的本质

设 \(f(x) = \displaystyle \sum_{i=0}^{n-1}a_ix^i\)，将 \(x_0\) 代入，得 \(f(x_0) = \displaystyle \sum_{i=0}^{n-1}a_ix_0^i\)。将其写成点积的形式：

\[f(x_0) = \sum_{i = 0} ^ {n - 1} a_i x_0 ^ i = \begin{bmatrix} x_0 ^ 0 & x_0 ^ 1 & \cdots & x_0 ^ {n - 1} \end{bmatrix} \times \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n - 1} \end{bmatrix} \]

这样，如果有 \(x_0,x_1,\dots,x_{m-1}\) 需要求值，整个过程就可以写成 \(m\times n\) 维矩阵乘以 \(n\) 维列向量的形式：

\[\begin{bmatrix} x_0 ^ 0 & x_0 ^ 1 & \cdots & x_0 ^ {n - 1} \\ x_1 ^ 0 & x_1 ^ 1 & \cdots & x_1 ^ {n - 1} \\ \vdots & \vdots & \ddots & \vdots \\ x_{m - 1} ^ 0 & x_{m - 1} ^ 1 & \cdots & x_{m - 1} ^ {n - 1} \\ \end{bmatrix} \times \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n - 1} \end{bmatrix} = \begin{bmatrix} f(x_0) \\ f(x_1) \\ \vdots \\ f(x_{m - 1}) \end{bmatrix} \]

左侧矩阵就是范德蒙德矩阵。

当 \(n=m\) 时为范德蒙德方阵，如果 \(x_i\) 互不相同，那么矩阵逆存在，这帮助我们快速从点值表示法转回系数表示法。

朴素计算求值的复杂度为 \(O(nm)\)。快速傅里叶变换即在离散傅里叶变换的基础上通过选取合适的 \(x_i\)，使得可以快速求值。

离散傅里叶变换和离散傅里叶逆变换

在介绍 FFT 之前，我们先介绍一下离散傅里叶变换（Discrete Fourier Transfrom，DFT）和离散傅里叶逆变换（Inverse Discrete Fourier Transfrom，IDFT）

在 OI 中，离散傅里叶变换可以看成是多项式 \(f(x)\) 在 \(n\) 个单位根处求值，也就是把系数表达转换为点值表达。

而离散傅里叶逆变换，便是从点值表达转化为系数表达。

而我们所要学习的 FFT，便是用来实并加速 DFT 的。

快速傅里叶变换

先放一张图。

快速傅里叶变化目的就是让我们在 \(O(n\log n)\) 的时间复杂度内实现系数表达转化为点值表达，在 \(O(n)\) 的点值乘法之后，应用 IFFT 将点值表达转化为系数表达之后，我们便可以把 \(O(n^2)\) 的多项式乘法优化到 \(O(n\log n)\)。

接下来我们的目标就明确了，如何在 \(O(n\log n)\) 的时间内，求出 \(n-1\) 次多项式 \(f(x) = \displaystyle\sum_{i=0}^{n-1}a_ix^i\) 的至少 \(n\) 个点值。

简化情况

我们试着把问题问题简化，由浅入深。

我们都知道函数的性质是有限的：奇偶性，单调性，周期性等。一般函数没有单调性和周期性，但是一般函数总能被表示成一个偶函数和一个奇函数之和，这给我们提供了思路。

我们将任意的多项式拆成偶函数 \(f_e(x)\) 和奇函数 \(f_o(x)\) 之和，则有

\[\begin{cases} f(x) = f_e(x) + f_o(x) \\ f(-x) = f_e(x) - f_o(x) \end{cases} \]

选择 \(\lceil\frac n 2\rceil\) 对两两互为相反数的值 \((x_i,-x_i)\)，求出所有 \(x_i\) 在 \(f_e(x)\) 和 \(f_o(x)\) 处的取值。

我们设 \(n\) 是偶数，不难发现，\(f_e\) 是 \(n-2\) 次多项式，\(f_o\) 是 \(n-1\) 次多项式，本质上还是求了 \(n-1\) 次多项式的 \(n\) 个点值，对时间复杂度没有优化。

但是 \(f_e\) 和 \(f_o\) 的项数减半，我们尝试利用这个性质。

因为 \(f_e = a_0x^0 + a_2x^2 + \dots\)，我们可以自然的想到换元 \(u = x^2\)，这样 \(f(u) = a_0u^0 + a_2u^2 + \dots\)。我们设 \(f_e'(x)=a_0x^0+a_2x^1+\dots\)，则 \(f_e(x) = f_e'(x^2)\)。

同理，我们从 \(f_o\) 中提出一个 \(x\)，设 \(f_o'(x)=a_1x^0+a_3x^1+\dots\)，则\(f_o(x) = xf_o'(x^2)\)，因此我们有

\[\begin{cases} f(x) = f'_e(x ^ 2) + xf'_o(x ^ 2) \\ f(-x) = f'_e(x ^ 2) - xf'_o(x ^ 2) \end{cases} \]

这才是真正意义上的规模减半。这个问题的复杂度就是 \(T(n) = 2T(\frac{n}{2})+O(n) = O(n\log n)\)。

单位根的引入

问题来了，如何保证递归后我们取得的 \(x\) 值也满足两两互为相反数呢？

考虑一开始的 \((x_i,-x_i)\)，这说明一开始存在 \(i'\) 使得 \(x_i^2 = -x_{i'}^2\)，它们互不相同但是它们的 \(4\) 次方相等。

再进一步考虑，因为问题会递归 \(w = \lceil\log_2n\rceil\) 层，所以我们需要找到 \(k=2^w\) 个互不相等的 \(x\)，但它们的 \(k\) 次幂相等，它们可以等于任何数，方便起见，设为 \(1\)，即 \(x^k=1\)，那么 \(x\) 就是所有的 \(k\) 次单位根，单位根将在这里发挥巨大的作用。

递归求解

得到大致框架后，我们具体地描述整个算法流程：

首先将 \(n\) 补齐到不小于 \(n\) 的最小的 \(2\) 的整数幂，即 \(2^{\lceil \log_2n\rceil}\)，这一步的目的就是让我们在递归的时候不会出错，因为每次会划分成两个子问题，每个子问题含这个问题的一半。

为了方便起见，我将 \(f_e\) 设为 \(A_1\)，\(f_o\) 设为 \(A_2\)，\(f\) 设为 \(A\)。

那么就有 \(A(x) = A_1(x^2) + xA_2(x^2)\)

我们将 \(\omega_n^k(k<\frac{n}{2})\) 代入得：

\[A(\omega_n^k)=A_1(\omega_n^{2k})+A_2(\omega_n^{2k})\omega_n^k \]

\[=A_1(\omega_{\frac{n}{2}}^{k})+A_2(\omega_{\frac{n}{2}}^{k})\omega_n^k \]

将 \(\omega_n^{k+\frac{n}{2}}(k<\frac{n}{2})\) 代入得：

\[A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+A_2(\omega_n^{2k+n})\omega_n^{k+\frac{n}{2}} \]

\[=A_1(\omega_{\frac{n}{2}}^{k})-A_2(\omega_{\frac{n}{2}}^{k})\omega_n^k \]

这也就是说，如果我们知道了两个多项式 \(A_1,A_2\) 分别在 \(\omega_{\frac{n}{2}}^0,\omega_{\frac{n}{2}}^1,\dots,\omega_{\frac{n}{2}}^{\frac{n}{2}-1}\) 处的点值表示。

套用第一个公式，我们就能求出 \(A\) 在 \(\omega_{n}^0,\omega_{n}^1,\dots,\omega_{n}^{\frac{n}{2}-1}\) 处的点值表示；套用第二个公式，我们就能知道 \(A\) 在 \(\omega_{n}^{\frac{n}{2}},\omega_{n}^{\frac{n}{2}+1},\dots,\omega_{n}^{n}\) 处的点值表示。

我们又发现 \(A_1\) 和 \(A_2\) 两个多项式是原规模一半的子问题，它们的性质完全相同。

因此就可以一直分治下去，直到多项式只剩下一项为止。

代码实现

我们手写一个复数类，然后应用刚才提到的这两个式子：

\[A(\omega_n^k) =A_1(\omega_{\frac{n}{2}}^{k})+A_2(\omega_{\frac{n}{2}}^{k})\omega_n^k \]

\[A(\omega_n^{k+\frac{n}{2}}) =A_1(\omega_{\frac{n}{2}}^{k})-A_2(\omega_{\frac{n}{2}}^{k})\omega_n^k \]

由于我们已经知道了 \(A_1,A_2\)，所以我们只需要代入前 \(\frac{n}{2}\) 个单位根便能得到整个 \(A\) 的点值。

complex operator + (complex a,complex b) { return complex{a.x+b.x,a.y+b.y}; }

complex operator - (complex a,complex b) { return complex{a.x-b.x,a.y-b.y}; }

complex operator * (complex a,complex b) { return complex{a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y}; }

il void FFT(complex A[],int n)
{
	if(n == 1) return ;
	complex A1[n] , A2[n];
	for(re int i=0;i<n/2;i++) A1[i] = A[i<<1] , A2[i] = A[(i<<1)+1];
	FFT(A1,n/2) , FFT(A2,n/2);
	complex Wn = {cos(2*Pi/n),sin(2*Pi/n)};
	complex Wk = {1,0};
	for(re int i=0;i<n/2;i++)
	{
		A[i] = A1[i] + A2[i] * Wk;
		A[i+(n>>1)] = A1[i] - A2[i] * Wk;
		Wk = Wk * Wn;
	}
}

这里面，\(A\) 数组的初值实部存系数，虚部存 \(0\)。最终返回多项式的点值。

其实这就类似于一个合并的过程，到最底层 \(n=1\) 时 \(A(\omega_n^0) = A(1)\)，也就是该多项式的系数就是此时的点值，然后不断按照式子合并成一个大多项式的点值，最后就变成了我们想要求的多项式的点值。

快速傅里叶逆变换

我们求得点值后，因为我们最终想要的还是系数表示法，所以我们还要把点值转化为系数，也就是我们的 IFFT。

由于我的知识面浅薄，线性代数的角度还不能理解，所以我就只说这是怎么构造的，证明就略去了。

从另一个角度看 DFT, 我们实际上是把一个矩阵 \(W=[\omega^{ij}]\) 左乘在了系数向量上, 那么实现 IDFT 实际上我们需要找到这个单位根矩阵的逆。通过计算，可以得出 \(W'=[\omega^{-ij}]\)。乘起来是一个快速傅里叶变换矩阵的常数倍。

不加证明的给出最终的矩阵：

\[\mathcal F ^ {-1} = \frac 1 n \begin{bmatrix} (\omega ^ {-0}) ^ 0 & (\omega ^ {-0}) ^ 1 & \cdots & (\omega ^ {-0}) ^ {n - 1} \\ (\omega ^ {-1}) ^ 0 & (\omega ^ {-1}) ^ 1 & \cdots & (\omega ^ {-1}) ^ {n - 1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega ^ {-(n - 1)}) ^ 0 & (\omega ^ {-(n - 1)}) ^ 1 & \cdots & (\omega ^ {-(n - 1)}) ^ {n - 1} \\ \end{bmatrix} \]

\[W^{-1}=\displaystyle \frac{W'}{n} \]

\[\frac{1}{\omega_n^k}=\frac{1}{\cos\theta + i\sin\theta} = \cos(-\theta)+i\sin(-\theta)=\omega_n^{-k} \]

因此，对一个序列做 IFFT，只需将 FFT 递归公式里的 \(\omega_n^k\) 换成 \(\omega_n^{-k}\)，并在最后除以一个 \(n\) 即可。

上面从线性代数角度来说的，如果学过的话可能理解起来比较容易，下面是从代数角度构造这个解，可能更适合没学过线性代数的同学，但其实这个也是根据线代得出的。

设多项式 \(A(x) = a_0 + a_1x + a_2x^2 + \dots + a_{n-1}x^{n-1}\)

我们通过 FFT，已经得到了一组点值 \((\omega_n^0,y_0),(\omega_n^1,y_1),\dots,(\omega_n^{n-1},y_{n-1})\)

其中 \(y_i=\displaystyle \sum_{j=0}^{n-1}a_j(\omega_n^i)^j\)

我们构造多项式 \(B(x) = y_0 + y_1x+\dots+y_{n-1}x^{n-1}\)

把 \(n\) 个单位根的倒数 \(\omega_n^0,\omega_n^{-1},\dots,\omega_n^{-(n-1)}\) 代入 \(B(x)\) 得到 \(n\) 个新点值，设为 \((z_0,z_1,\dots,z_{n-1})\)，\(z_k\) 和 \(a_k\) 具有一定的联系，我们探讨一下：

\[z_k = \sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i =\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^i)^j(\omega_n^{-k})^i \]

\[= \sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i \]

当 \(j=k\) 时，内层的和式等于 \(n\)

当 \(j\not=k\) 时，内层是个等比数列，套入公式得到 \(\displaystyle \frac{(\omega_n^{j-k})^n-1}{\omega_n^{j-k}-1}=\frac{(\omega_n^n)^{j-k}-1}{\omega_n^{j-k}-1}= 0\)

所以 \(z_k = na_k\)，即 \(a_k = \displaystyle \frac{z_k}{n}\)！

最后的结果是一样的。

代码实现

il void IFFT(complex A[],int n)
{
	if(n == 1) return ;
	complex A1[n] , A2[n];
	for(re int i=0;i<n/2;i++) A1[i] = A[i<<1] , A2[i] = A[(i<<1)+1];
	FFT(A1,n/2,type) , FFT(A2,n/2,type);
	complex Wn = {cos(2*Pi/n),-sin(2*Pi/n)};//变化的只有这里的负号
	complex Wk = {1,0};
	for(re int i=0;i<n/2;i++)
	{
		A[i] = A1[i] + A2[i] * Wk;
		A[i+(n>>1)] = A1[i] - A2[i] * Wk;
		Wk = Wk * Wn;
	}
}

你会发现，变化的只有一个负号，所以我们其实可以合并 FFT 和 IFFT 为一个函数，用一个标记判断是进行的 FFT 还是 IFFT 即可。

总代码实现

模板题代码

#include<bits/stdc++.h>
//#define int long long
#define ll long long
#define complex comple
#define re register
#define il inline
const int N = 4e6 + 5;
const double Pi = acos(-1.0);
using namespace std;
int max(int x,int y){return x > y ? x : y;}
int min(int x,int y){return x < y ? x : y;}

int n,m;
struct complex{
	double x,y;
	complex (double xx = 0 , double yy = 0) { x = xx; y = yy; }
}a[N],b[N];

il int read()
{
	int f=0,s=0;
	char ch=getchar();
	for(;!isdigit(ch);ch=getchar()) f |= (ch=='-');
	for(; isdigit(ch);ch=getchar()) s = (s<<1) + (s<<3) + (ch^48);
	return f ? -s : s;
}

complex operator + (complex a,complex b) { return complex{a.x+b.x,a.y+b.y}; }

complex operator - (complex a,complex b) { return complex{a.x-b.x,a.y-b.y}; }

complex operator * (complex a,complex b) { return complex{a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y}; }

il void FFT(complex A[],int n,int type)
{
	if(n == 1) return ;
	complex A1[n] , A2[n];
	for(re int i=0;i<n/2;i++) A1[i] = A[i<<1] , A2[i] = A[(i<<1)+1];
	FFT(A1,n/2,type) , FFT(A2,n/2,type);
	complex Wn = {cos(2*Pi/n),1.0 * type * sin(2*Pi/n)};
	complex Wk = {1,0};
	for(re int i=0;i<n/2;i++)
	{
		A[i] = A1[i] + A2[i] * Wk;
		A[i+(n>>1)] = A1[i] - A2[i] * Wk;
		Wk = Wk * Wn;
	}
}

signed main()
{
	n = read() , m = read();
	for(re int i=0;i<=n;i++) scanf("%lf",&a[i].x);
	for(re int i=0;i<=m;i++) scanf("%lf",&b[i].x);
	m = n + m , n = 1;
	while(n <= m) n <<= 1;
	FFT(a,n,1) , FFT(b,n,1);
	for(re int i=0;i<n;i++) a[i] = a[i] * b[i];//将求得的点值相乘
	FFT(a,n,-1);
	for(re int i=0;i<=m;i++) printf("%d ",(int)(a[i].x/n+0.5));//注意精度，四舍五入
	return 0;
}

迭代实现

蝴蝶变换

我们观察一下递归过程中的下标变化

我们发现

我们发现一个规律：每个位置的元素下标都做了二进制翻转，这个变换称为位逆序变换（也称蝴蝶变换）。

我们应用位逆序变换得到后序列，然后从下向上合并就行了，递归就变成了迭代，无疑是一个很大的优化。

蝴蝶变换的实现

蝴蝶变换可以 \(O(n)\) 递推实现。 首先还是一样的令 \(n=2^{\lceil \log_2^n\rceil}\)，然后设 \(R_x\) 表示 \(x\) 二进制翻转后的数，我们要求出 \(R_0\sim R_{n-1}\)。

由于我们是递推求解，所以在求 \(R_x\) 时，\(R_{\frac{x}{2}}\) 的值已经求出，所以我们先把 \(x\) 右移一位变成 \(\frac{x}{2}\)，然后翻转 \(\frac{x}{2}\)，再右移一位，这样做的目的就是先把 \(x\) 除了最后一位翻转出来，假设 \(x\) 二进制位有 \(5\) 位，翻转过后再右移一位是因为 \(R_{\frac{x}{2}}\) 是把这 \(5\) 位翻转，而我们想要让它翻转的其实只是后 \(4\) 位，所以把它右移。

最后再判断 \(x\) 的最后一位是否是 \(1\)，如果是的话，要再在翻转后的最高位上加上 \(1\)，即 \(\frac{n}{2}\)。这样递推式就出来了。

\[R_x = \frac{R_{\frac{x}{2}}}{2}+[x\& 1] \times \frac{n}{2} \]

常数优化

三次变两次优化

按照正常做法，我们会进行两次 FFT，分别求出点值后相乘，然后再做一次 IDFT。来回共调用了三次函数。

我们可以把它优化成两次，设两个多项式分别是 \(f(x),g(x)\)

具体操作是，我们可以把 \(g(x)\) 放到 \(f(x)\) 的虚部上去，FFT 求出点值之后，求出 \(f(x)^2\)，把虚部取出来除以 \(2\) 就是点值相乘后的结果了，正确性证明也很简单：

\[(a+bi)^2 = (a^2-b^2)+(2ab)i \]

这样我们就能把常数优化到原来的 \(\displaystyle \frac{2}{3}\) 了。

	FFT(a,n,1);
	for(re int i=0;i<n;i++) a[i] = a[i] * a[i];//将求得的点值相乘
	FFT(a,n,-1);
	for(re int i=0;i<=m;i++) printf("%d ",(int)(a[i].y/2/n+0.5));//注意精度，四舍五入

快速数论变换

快速数论变换（Number-Theoretic Transform, NTT），是 DFT 在数论基础上实现的，它和 FTT 一样，都是为了加速 DFT 而出现的。

前置知识：原根。

有了 FFT 的基础后，其实 NTT 就比较容易理解了。

原根的引入

根据阶和原根的定义，我们可以发现原根一个很优美的性质：那就是 \(\{g,g^2,\dots,g^{\varphi(p)}\}\) 在模 \(p\) 意义下互不相同，也就是说，\(g\) 就是这个循环群的生成元，它在模 \(p\) 意义下会出现循环，这与单位根十分相似。

根据原根存在定理，只有 \(2,4,p^\alpha,2p^\alpha\) 才有原根，所以我们的模数应为质数，此时 \(\varphi(p) = p-1\)。若 \(p-1\) 的分解中含有较多的 \(2\) 这个因子，我们就有足够的余地选取对称的 \(n(n=2^k)\) 个原根来进行递归了。

为了多次二分，模数 \(p\) 应该选取形如 \(q\times 2^k + 1\) 的质数，其中 \(k\) 为整数。为了保证可以二分，多项式相乘后的长度不能超过 \(2^k\)。

常见的模数有：

• \(998244353(119\times 2^{23}+1)\)，它的最大长度为 \(2^{23}\)，它的一个原根 \(g=3\)

• \(469762049(7\times 2^{26}+1)\)，它的最大长度为 \(2^{26}\)，它的一个原根 \(g=3\)

• \(2281701377(17\times 2^{27}+1)\)，它的最大长度为 \(2^{27}\)，它的一个原根 \(g=3\)

NTT 和 INTT

为了实现 \(n\) 次原根乘 \(n\) 次方会出现循环的性质，并且保证对称，我们直接用 \(g^{\frac{\varphi(p)}{n}}\)（\(p\) 为质数时就是 \(g^{\frac{p-1}{n}}\)） 来替换 \(n\) 次原根。这样的话，假设 \(p\) 是质数，那么设 \(g_n^0 = 1,g_n^1 = g^{\frac{p-1}{n}},\dots,g_n^k = g^{\frac{p-1}{n}k}\)。

我们来验证一下我们选取的原根具有的性质（以下地方省略 \(\bmod \ p\)）：

• 循环性：\(g_n^{k+n} = g_n^k\)。

  证：\(g_n^{k+n} = g^{p-1}g_n^k = g_n^k\)

• 可加性：\(g_n^kg_n^m = g_n^{k+m}\)

  证：\(g_n^kg_n^m = g^{\frac{p-1}{n}k}g^{\frac{p-1}{n}m} =g_n^{k+m}\)

• \(g_n^k = g_{cn}^{ck}(c > 0)\)，这个感性理解以下就行，你把整个圆分成 \(cn\) 份，它的第 \(ck\) 个原根就是你把圆分成 \(n\) 份，它的第 \(k\) 个原根。

  证也很好证：\(g_n^k = g^{\frac{p-1}{n}k}=g^{\frac{p-1}{cn}ck}=g_{cn}^{ck}\)

• 当 \(n\) 为偶数时，将 \(g_n^k\) 取反相当于将其逆时针（或顺时针）转半圈，所以 \(g_n^k = -g_n^{k\pm\frac{n}{2}}(2\mid n)\)。

  证：因为 \((g_n^{\frac{n}{2}})^2 = g_n^n = 1\)，而根据原根的定义 \(g_n^{\frac{n}{2}}\) 是不能等于 \(1\) 的，所以 \(g_n^{\frac{n}{2}}=-1\)。

• 对称性：因为 \(n\) 次原根将圆 \(n\) 等分，均匀分布在圆周，所以它们的中心就在原点，即 \(\displaystyle \sum_{i=0}^{n-1}g_n^i = 0\)。

放个图，你会发现，原根和单位根的相似之处。

可以证明，\(n\) 个原根的倒数也是与单位根相似的 \(g_n^0,g_n^{-1}\)，因此我们在 INTT 的时候取 \(g_n^1\) 的逆元即可。

代码实现

这里还是以 FFT 模板题为例，因为题目里保证了每项的系数不会超过 \(9\)，所以多项式相乘最终得到的系数不会超过模数 \(p\)，我们就可以近似的看作这个乘法是在模 \(p\) 意义下实现的。

递归版本

#include<bits/stdc++.h>
//#define int long long
#define ll long long
#define next nxt
#define re register
#define il inline
const int N = 4e6 + 5;
const ll mod = 998244353;
const ll g = 3;
using namespace std;
int max(int x,int y){return x > y ? x : y;}
int min(int x,int y){return x < y ? x : y;}

int n,m;
ll invg,invn;
ll a[N],b[N];

il int read()
{
	int f=0,s=0;
	char ch=getchar();
	for(;!isdigit(ch);ch=getchar()) f |= (ch=='-');
	for(; isdigit(ch);ch=getchar()) s = (s<<1) + (s<<3) + (ch^48);
	return f ? -s : s;
}

ll ksm(ll a,ll b)
{
	ll res = 1;
	while(b)
	{
		if(b&1) res = res * a % mod;
		a = a * a % mod;
		b >>= 1;
	}
	return res;
}
ll inv(ll x) { return ksm(x,mod-2); }

il void NTT(ll A[],int n,int type)
{
	if(n == 1) return ;
	ll A1[n/2] , A2[n/2];
	for(re int i=0;i<n/2;i++) A1[i] = A[i<<1] , A2[i] = A[(i<<1)+1];
	NTT(A1,n/2,type) , NTT(A2,n/2,type);
	ll g1 = ksm(type==1 ? g : invg ,(mod-1)/n);
	ll gk = 1;
	for(re int i=0;i<n/2;i++)
	{
		A[i] = (A1[i] + A2[i]*gk) % mod;
		A[i+(n>>1)] = ((A1[i] - A2[i] * gk) % mod + mod) % mod;
		gk = gk * g1 % mod;
	}
}

signed main()
{
	n = read() , m = read();
	for(re int i=0;i<=n;i++) a[i] = read();
	for(re int i=0;i<=m;i++) b[i] = read();
	m = n + m , n = 1;
	while(n <= m) n <<= 1;
	invg = inv(g) , invn = inv(n);
	NTT(a,n,1) , NTT(b,n,1);
	for(re int i=0;i<n;i++) a[i] = a[i] * b[i] % mod;
	NTT(a,n,-1);
	for(re int i=0;i<=m;i++) cout << a[i] * invn % mod << " ";
	return 0;
}

迭代版本

#include<bits/stdc++.h>
//#define int long long
#define ll long long
#define next nxt
#define re register
#define il inline
const int N = 4e6 + 5;
const ll mod = 998244353;
const ll g = 3;
using namespace std;
int max(int x,int y){return x > y ? x : y;}
int min(int x,int y){return x < y ? x : y;}

int n,m;
ll invg,invn;
ll a[N],b[N],R[N];

il ll read()
{
	ll f=0,s=0;
	char ch=getchar();
	for(;!isdigit(ch);ch=getchar()) f |= (ch=='-');
	for(; isdigit(ch);ch=getchar()) s = (s<<1) + (s<<3) + (ch^48);
	return f ? -s : s;
}

ll ksm(ll a,ll b)
{
	ll res = 1;
	while(b)
	{
		if(b&1) res = res * a % mod;
		a = a * a % mod;
		b >>= 1;
	}
	return res;
}
ll inv(ll x) { return ksm(x,mod-2); }

il void NTT(ll A[],int n,int type)
{
	for(re int i=0;i<n;i++) if(i < R[i]) swap(A[i],A[R[i]]);
	for(re int blo=2;blo<=n;blo<<=1)
	{
		ll g1 = ksm(type==1 ? g : invg , (mod-1)/blo);
		for(re int i=0;i<n;i+=blo)
		{
			ll gk = 1;
			for(re int j=0;j<blo/2;j++)
			{
				ll x = A[i+j] , y = A[i+j+(blo>>1)] * gk; 
				A[i+j] = (x + y) % mod;
				A[i+j+(blo>>1)] = ((x-y) % mod + mod) % mod;
				gk = gk * g1 % mod;
			}
		}
	}
}

signed main()
{
	n = read() , m = read();
	for(re int i=0;i<=n;i++) a[i] = read();
	for(re int i=0;i<=m;i++) b[i] = read();
	m = n + m , n = 1;
	while(n <= m) n <<= 1;
	invg = inv(g) , invn = inv(n);
	for(re int i=0;i<n;i++) R[i] = (R[i>>1]>>1) + ((i&1) ? n/2 : 0);
	NTT(a,n,1) , NTT(b,n,1);
	for(re int i=0;i<n;i++) a[i] = a[i] * b[i] % mod;
	NTT(a,n,-1);
	for(re int i=0;i<=m;i++) cout << a[i] * invn % mod << " ";
	return 0;
}

任意模数多项式乘法

MTT，是由毛啸在 2016 年国家集训队论文《再探快速傅里叶变换》 中提出的，取其首字母，因而叫 MTT。

参考资料：Oxide-[学习笔记]拆系数 FFT

如果现在要求两个整数多项式 \(F(x),G(x)\) 的卷积，序列长度 \(n \leq 10^5\)，多项式系数 \(A_i,B_i\leq 10^9\)，答案对 \(p\leq 10^9\) 取模。

你会发现这题的数据范围是 \(10^9 \times 10^9 \times 10^5 = 10^{23}\)(因为有 \(a_{x+y}\leftarrow a_xa_y\)，所以这样的 \(x,y\) 的组合是 \(O(n)\) 的) 级别的，long double 也存不下。使用 FFT 会掉精度，而 NTT 会因为模数的性质而失去作用，这时候我们就需要用到 MTT 了。

MTT 有两种方法，一种是三模数 NTT，另一种是拆系数 FFT。

其中，三模 NTT 精度优秀，但常数较大，会调用 \(9\) 次 FFT 函数；而拆系数 FFT 则相反，朴素的拆系数 FFT 要调用 \(7\) 次 FFT 函数，优化后只需要 \(4\) 次。

下面，我们分别介绍这两种算法。

拆系数 FFT

在我的观点看来，拆系数 FFT 相较于三模 NTT 更容易理解。

算法流程

我们设 \(\text{base}\) 为一个 \(\sqrt p\) 级别的数，类似除法的，我们可以将 \(F(x)\) 分解成 \(\text{base} \cdot A(x) + B(x)\)，这样拆出来的系数就是 \(<\text{base}\) 的，因此最终答案就是 \(10^{14}\) 次方级别的，可以接受。

故

\[F * G = (\text{base} \cdot A + B)(\text{base} \cdot C +D) \]

\[=\text{base}^2\cdot A * C + \text{base} \cdot(A * D + B * C) + B * D \]

这样以后显然就有一个四次 DFT 和三次 IDFT 的朴素做法：对 \(A,B,C,D\) 进行 DFT，\(\hat{A} * \hat{C},\hat{A} * \hat{D} + \hat{B} * \hat{C},\hat{C} * \hat{D}\) 进行 IDFT，共 \(7\) 次。这里 \(\hat{A}\) 表示的是点值表示法后的多项式 \(A\)。

类似于 FFT 三次转两次优化，我们可以利用虚部来存储。（下午将 \(f(\omega_n^k)\) 简写为 \(f(k)\))我们令

\[f(k) = A(k) + i \times B(k) \]

\[g(k) = A(k) - i\times B(k) \]

\[h(k) = C(k) + i \times D(k) \]

推一下式子后，你会发现：\(f(k)\) 和 \(g(n-k)\) 是共轭的，这里 \(n\) 还是最小的比多项式相乘后次数高的 \(2\) 的整数幂。下面的式子推导来源于 Oxide，我就直接贴上了，简单看一下就行

\[\begin{aligned} f(k)&=A(k)+i\times B(k) \\&=\sum_{j=0}^{n-1}a_j(\omega_n^k)^j+i\times \sum_{j=0}^{n-1}b_j(\omega_n^k)^j\\&=\sum_{j=0}^{n-1}(a_j+i\times b_j)(\omega_n^k)^j\\&=\sum_{j=0}^{n-1}(a_j+i\times b_j)\left(\cos\left(\frac{2\pi kj}{n}\right)+i\sin \left(\frac{2\pi kj}{n}\right)\right)\\&=\sum_{j=0}^{n-1}\left(a_j\times \cos\left(\frac{2\pi kj}{n}\right)-b_j\times \sin \left(\frac{2\pi kj}{n}\right)\right)+i\left(b_j\times \cos\left(\frac{2\pi kj}{n}\right)+a_j\times \sin \left(\frac{2\pi kj}{n}\right)\right) \end{aligned} \]

\[\begin{aligned} g(n-k)&=A(n-k)-i\times B(n-k)\\&=\sum_{j=0}^{n-1}a_j(\omega_n^{-k})^j-i\times \sum_{j=0}^{n-1}b_j(\omega_n^{-k})^j\\&=\sum_{j=0}^{n-1}(a_j-i\times b_j)(\omega_n^{-k})^j\\&=\sum_{j=0}^{n-1}(a_j-i\times b_j)\left(\cos\left(\frac{2\pi kj}{n}\right)-i\sin \left(\frac{2\pi kj}{n}\right)\right)\\&=\sum_{j=0}^{n-1}\left(a_j\times \cos\left(\frac{2\pi kj}{n}\right)-b_j\times \sin \left(\frac{2\pi kj}{n}\right)\right)-i\left(b_j\times \cos\left(\frac{2\pi kj}{n}\right)+a_j\times \sin \left(\frac{2\pi kj}{n}\right)\right) \end{aligned} \]

需要注意一个特判：\(f(k)\) 和 \(g(n-k)\) 共轭也就意味着 \(f(n)\) 和 \(g(0)\) 配对，其实就是 \(f(0)\) 和 \(g(0)\) 配对，所以 \(g(0)\) 需要特殊计算一下。

综上，我们只需要计算 \(f,h\) 的点值，两次 FFT，然后利用 \(f\) 来 \(O(n)\) 求解 \(g\)，便可求出三个函数的点值。

令 \(p = f * h , q = g * h\)，那么就有

\[p = A * C - B * D + i(A * D + B * C) \]

\[q = A * C + B * D + i(A * D - B * C) \]

把 \(p,q\) 转会系数多项式需要 \(2\) 次 FFT，将 \(p,q\) 对应项相加就可得到 \(A * C, A * D\)，进而我们也能知道 \(B * C,B * D\)，从而得到最终的系数。

综上所述，我们只需要 \(4\) 次 FFT。

代码实现

在实现代码的时候，为了不掉精度，建议把 double 换成 long double，并且要注意时刻取模！

\(7\) 次 FFT

#include<bits/stdc++.h>
#define int long long
#define ll long long
#define complex comple
#define double long double
#define re register
#define il inline
const int N = 4e5 + 5;
const double Pi = acos(-1.0);
const ll bas = (1<<15) , bas2 = bas * bas;
using namespace std;
int max(int x,int y){return x > y ? x : y;}
int min(int x,int y){return x < y ? x : y;}

int n,m,mod;
int R[N];
ll x,x1,x2,x3;

il int read()
{
	int f=0,s=0;
	char ch=getchar();
	for(;!isdigit(ch);ch=getchar()) f |= (ch=='-');
	for(; isdigit(ch);ch=getchar()) s = (s<<1) + (s<<3) + (ch^48);
	return f ? -s : s;
}

namespace Poly
{
	struct complex
	{
		double x,y;
		complex (double xx = 0 , double yy = 0) { x = xx; y = yy; }
	}f0[N],f1[N],g0[N],g1[N],A[N],B[N],C[N];
	complex operator + (complex a,complex b) { return complex{a.x+b.x,a.y+b.y}; }
	complex operator - (complex a,complex b) { return complex{a.x-b.x,a.y-b.y}; }
	complex operator * (complex a,complex b) { return complex{a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y}; }
	void FFT(complex A[],int n,int type)
	{
		for(int i=0;i<n;i++) if(i < R[i]) swap(A[i],A[R[i]]);
		for(int blo=2;blo<=n;blo<<=1)
		{
			complex Wn = {cos(2*Pi/blo),1.0 * type * sin(2*Pi/blo)};
			for(int i=0;i<n;i+=blo)
			{
				complex Wk = {1,0};
				for(int j=0;j<blo/2;j++)
				{
					complex x = A[i+j] , y = A[i+j+(blo>>1)] * Wk;
					A[i+j] = x + y , A[i+j+(blo>>1)] = x - y;
					Wk = Wk * Wn;
				}
			}
		}
	}
	void init()
	{
		m = n + m , n = 1;
		while(n <= m) n <<= 1;
		for(re int i=0;i<n;i++) R[i] = (R[i>>1]>>1) + ((i&1) ? n/2 : 0);
		FFT(f0,n,1) , FFT(f1,n,1) , FFT(g0,n,1) , FFT(g1,n,1);
		for(re int i=0;i<n;i++) A[i] = f0[i] * g0[i] , B[i] = f0[i]*g1[i] + f1[i]*g0[i] , C[i] = f1[i] * g1[i];//将求得的点值相乘
		FFT(A,n,-1) , FFT(B,n,-1) , FFT(C,n,-1);
		for(re int i=0;i<=m;i++)
		{
			x1 = (ll)(A[i].x/n+0.5) % mod , x2 = (ll)(B[i].x/n+0.5) % mod, x3 = (ll)(C[i].x/n+0.5) % mod;//时刻取模
			x1 = x1 * bas2 % mod , x2 = x2 * bas % mod;
			cout << (x1 + x2 + x3) % mod << " ";
		}
	}
}
using namespace Poly;

signed main()
{
	n = read() , m = read() , mod = read();
	for(re int i=0;i<=n;i++)
	{
		x = read() % mod;
		f0[i].x = x / bas , f1[i].x = x % bas;
	}
	for(re int i=0;i<=m;i++)
	{
		x = read() % mod;
		g0[i].x = x / bas , g1[i].x = x % bas;
	}
	init();
	return 0;
}

\(4\) 次 FFT

#include<bits/stdc++.h>
#define int long long
#define ll long long
#define complex comple
#define double long double
#define re register
#define il inline
const int N = 4e5 + 5;
const double Pi = acos(-1.0);
const ll bas = (1<<15) , bas2 = bas * bas;
using namespace std;
int max(int x,int y){return x > y ? x : y;}
int min(int x,int y){return x < y ? x : y;}

int n,m,mod;
int R[N];
ll x,AC,AD,BC,BD;

il int read()
{
	int f=0,s=0;
	char ch=getchar();
	for(;!isdigit(ch);ch=getchar()) f |= (ch=='-');
	for(; isdigit(ch);ch=getchar()) s = (s<<1) + (s<<3) + (ch^48);
	return f ? -s : s;
}

namespace Poly
{
	struct complex
	{
		double x,y;
		complex (double xx = 0 , double yy = 0) { x = xx; y = yy; }
	}f[N],g[N],h[N];
	complex operator + (complex a,complex b) { return complex{a.x+b.x,a.y+b.y}; }
	complex operator - (complex a,complex b) { return complex{a.x-b.x,a.y-b.y}; }
	complex operator * (complex a,complex b) { return complex{a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y}; }
	void FFT(complex A[],int n,int type)
	{
		for(int i=0;i<n;i++) if(i < R[i]) swap(A[i],A[R[i]]);
		for(int blo=2;blo<=n;blo<<=1)
		{
			complex Wn = {cos(2*Pi/blo),1.0 * type * sin(2*Pi/blo)};
			for(int i=0;i<n;i+=blo)
			{
				complex Wk = {1,0};
				for(int j=0;j<blo/2;j++)
				{
					complex x = A[i+j] , y = A[i+j+(blo>>1)] * Wk;
					A[i+j] = x + y , A[i+j+(blo>>1)] = x - y;
					Wk = Wk * Wn;
				}
			}
		}
	}
	void init()
	{
		m = n + m , n = 1;
		while(n <= m) n <<= 1;
		for(re int i=0;i<n;i++) R[i] = (R[i>>1]>>1) + ((i&1) ? n/2 : 0);
		FFT(f,n,1) , FFT(h,n,1);
		g[0] = {f[0].x,-f[0].y};
		for(re int i=1;i<n;i++) g[i] = {f[n-i].x,-f[n-i].y};
		for(re int i=0;i<n;i++) h[i].x /= n , h[i].y /= n , f[i] = f[i] * h[i] , g[i] = g[i] * h[i];
		FFT(f,n,-1) , FFT(g,n,-1);
		for(re int i=0;i<=m;i++)
		{
			AC = ((ll)((f[i].x+g[i].x)/2+0.5)) % mod;
			AD = ((ll)((f[i].y+g[i].y)/2+0.5)) % mod;
			BD = ((ll)(round(g[i].x-AC)) % mod + mod) % mod;
			BC = ((ll)(round(f[i].y-AD)) % mod + mod) % mod;
			AC = bas2 % mod * AC % mod , AD = bas % mod * (AD+BC) % mod;
			cout << (AC + AD + BD) % mod << " ";
		}
	}
}
using namespace Poly;

signed main()
{
	n = read() , m = read() , mod = read();
	for(re int i=0;i<=n;i++)
	{
		x = read() % mod;
		f[i].x = x / bas , f[i].y = x % bas;
	}
	for(re int i=0;i<=m;i++)
	{
		x = read() % mod;
		h[i].x = x / bas , h[i].y = x % bas;
	}
	init();
	return 0;
}

三模 NTT

前置知识：exCRT。

算法流程

选取三个常用的 NTT 模数，分别算出 \(F * G\) 在这些模数下的结果，然后使用中国剩余定理合并即可。

一般来说，我们选择 \(p_1 =998244353,p_2=1004525809,p_3=469762049\)，设结果分别为 \(x_1,x_2,x_3\)，则有

\[\begin{cases} & x \equiv x_1(\bmod \ p_1)\\ & x \equiv x_2(\bmod \ p_2)\\ & x \equiv x_3(\bmod \ p_3)\\ \end{cases} \]

注：后面 \(p_1,p_2,p_3\) 用 \(A,B,C\) 表示。

如果要用 CRT 合并的话，\(p_1p_2p_3\) 的积超过了 long long 的范围，所以需要 __int128。

而使用 exCRT 合并则不需要 __int128：

• 先把前两个合并：

  \[ \begin{aligned} &x_1 + k_1A = x_2 + k_2B\\ &x_1 + k_1A \equiv x_2(\bmod \ B)\\ &k_1 \equiv \frac{x_2-x_1}{A}(\bmod \ B) \end{aligned} \]

• 于是我们就求出了 \(k_1\)，那就求出了 \(x\equiv x_1+k_1A(\bmod \ AB)\)，记 \(x_4=x_1 + k_1A\)

  \[ \begin{aligned} &x_4 + k_4AB = x_3 + k_3C\\ &x_4 + k_4AB \equiv x_3(\bmod \ C)\\ &k_4 \equiv \frac{x_3-x_4}{AB}(\bmod \ C) \end{aligned} \]

求出了 \(k_4,x\equiv x_4+k_4AB(\bmod \ ABC)\)，因为 \(x_4+k_4AB < ABC\)，所以 \(x=x_4+k_4AB\)，求出答案后，最后再对题目给出的模数取模即可。

从上述过程中可以看出，三模 NTT 要调用 \(9\) 次 NTT 函数（\(6\) 次 DFT，\(3\) 次 IDFT），所以说效率上要比拆系数 FFT 低，但是我的拆系数 FFT 好像写丑了，模板题跑了 3.5s，而三模 NTT 只跑了 1.3s，令人感叹。

代码实现

#include<bits/stdc++.h>
//#define int long long
#define ll long long
#define next nxt
#define re register
#define il inline
const int N = 4e6 + 5;
const ll mod1 = 998244353;
const ll mod2 = 1004535809;
const ll mod3 = 469762049;
const ll g = 3;
using namespace std;
int max(int x,int y){return x > y ? x : y;}
int min(int x,int y){return x < y ? x : y;}

int n,m;
ll invg,inv1,inv2,invmod1,invmod2,invmod3,x1,x2,x3,x,k,mod;
ll a1[N],b1[N],a2[N],b2[N],a3[N],b3[N],R[N];

il ll read()
{
	ll f=0,s=0;
	char ch=getchar();
	for(;!isdigit(ch);ch=getchar()) f |= (ch=='-');
	for(; isdigit(ch);ch=getchar()) s = (s<<1) + (s<<3) + (ch^48);
	return f ? -s : s;
}

namespace Poly
{
	ll ksm(ll a,ll b,ll mod)
	{
		ll res = 1;
		while(b)
		{
			if(b&1) res = res * a % mod;
			a = a * a % mod;
			b >>= 1;
		}
		return res;
	}
	ll inv(ll x,ll mod) { return ksm(x,mod-2,mod); }
	il void NTT(ll A[],int n,int type,ll mod)
	{
		invg = inv(g,mod);
		for(re int i=0;i<n;i++) if(i < R[i]) swap(A[i],A[R[i]]);
		for(re int blo=2;blo<=n;blo<<=1)
		{
			ll g1 = ksm(type==1 ? g : invg , (mod-1)/blo , mod);
			for(re int i=0;i<n;i+=blo)
			{
				ll gk = 1;
				for(re int j=0;j<blo/2;j++)
				{
					ll x = A[i+j] , y = A[i+j+(blo>>1)] * gk; 
					A[i+j] = (x + y) % mod;
					A[i+j+(blo>>1)] = ((x-y) % mod + mod) % mod;
					gk = gk * g1 % mod;
				}
			}
		}
	}
	il void init()
	{
		m = n + m , n = 1;
		while(n <= m) n <<= 1;
		for(re int i=0;i<n;i++) R[i] = (R[i>>1]>>1) + ((i&1) ? n/2 : 0);
		NTT(a1,n,1,mod1),NTT(b1,n,1,mod1),NTT(a2,n,1,mod2),NTT(b2,n,1,mod2),NTT(a3,n,1,mod3) , NTT(b3,n,1,mod3);
		for(re int i=0;i<n;i++) a1[i] = a1[i] * b1[i] % mod1 , a2[i] = a2[i] * b2[i] % mod2 , a3[i] = a3[i] * b3[i] % mod3;
		NTT(a1,n,-1,mod1),NTT(a2,n,-1,mod2),NTT(a3,n,-1,mod3);
		invmod1 = inv(n,mod1) , invmod2 = inv(n,mod2) , invmod3 = inv(n,mod3);//最后别忘了再不同的模意义下除以n
		for(re int i=0;i<n;i++) a1[i] = a1[i] * invmod1 % mod1 , a2[i] = a2[i] * invmod2 % mod2 , a3[i]= a3[i] * invmod3 % mod3;
		inv1 = inv(mod1,mod2) , inv2 = inv(1ll * mod1 * mod2 % mod3,mod3);
		for(re int i=0;i<=m;i++)
		{
			k = ((a2[i]-a1[i]) % mod2 + mod2) % mod2 * inv1 % mod2;
			x = k * mod1 + a1[i];//第一步合并
			k = ((a3[i]-x) % mod3 + mod3) % mod3 * inv2 % mod3;
			x = (x + k * mod1 % mod * mod2 % mod) % mod;//第二步合并
			cout << x << " ";
		}
	}
}
using namespace Poly;

signed main()
{
	n = read() , m = read() , mod = read();
	for(re int i=0;i<=n;i++) a1[i] = read() % mod , a2[i] = a1[i] , a3[i] = a1[i];
	for(re int i=0;i<=m;i++) b1[i] = read() % mod , b2[i] = b1[i] , b3[i] = b1[i];
	Poly::init();
	return 0;
}

应用

学了 FFT/NTT 有一段时间了，但是直到现在才知道它最基本的用途。

一般来说，FFT/NTT 是用来优化一些卷积形式的式子的。

我们都知道多项式 \(f\) 和 \(g\) 相乘后，\(x^i\) 的系数是

\[\sum_{j+k=i}a_j\times b_k \]

其中 \(a_j\) 是 \(f\) 里 \(x^j\) 的系数，\(b_k\) 是 \(g\) 里 \(x^k\) 的系数。这个也很好理解。而 FFT/NTT 的用途便是去加速这种形式的式子。

举个例子，比方说我们想求出 \(E_i\) 在 \(0\sim n\) 的取值，其中

\[E_i = \sum_{j+k=i}f_j \times g_k \]

\(f,g\) 的值已经给出。

你会发现，求 \(E_i\) 的这个式子就很像多项式的卷积形式，这时候，我们就可以采用 FFT/NTT 来加速。

具体的，虽然 \(f,g\) 并不是一个多项式，但是我们可以采取类似于形式幂级数的思想，把它想成是一个多项式，也就是 \(f(x) = f_0 + f_1x+f_2x^2+\dots,g(x)=g_0 + g_1x + g_2x^2\)，我们不关心 \(x\) 的取值，令这两个多项式相乘，\(x^i\) 前面的系数，便是我们想要的 \(E_i\) 的答案，这就是 FFT/NTT 的实际应用。在做题中，我们往往就需要对题目中给出的式子进行一系列的推导，最终化成我们想要的卷积形式。

例题

P3338 [ZJOI2014]力

给出 \(n\) 个数 \(q_1,q_2, \dots q_n\)，定义

\[F_j~=~\sum_{i = 1}^{j - 1} \frac{q_i \times q_j}{(i - j)^2}~-~\sum_{i = j + 1}^{n} \frac{q_i \times q_j}{(i - j)^2} \]

\[E_i~=~\frac{F_i}{q_i} \]

对 \(1 \leq i \leq n\)，求 \(E_i\) 的值(\(n \leq 10^5,0 < F_i < 10^9\))。

既然题目让我们求 \(E_i\)，我们就直接代进去，发现 \(q_j\) 没了

\[F_j=\sum_{i = 1}^{j - 1} \frac{q_i}{(i - j)^2}-\sum_{i = j + 1}^{n} \frac{q_i}{(i - j)^2} \]

我们在两个 \(\sum\) 里都加上 \(i=j\) 这一项，我们假定它是有意义的，因为这一项左右两边相等会被消掉，所以是可以加上的。

\[F_j=\sum_{i = 1}^j \frac{q_i}{(i - j)^2}-\sum_{i = j}^{n} \frac{q_i}{(i - j)^2} \]

这时候我们就可以抽象出来一个函数了，设 \(f_i=q_i,g_i = \displaystyle \frac{1}{i^2}\)，那么原式就变为了

\[F_j = \sum_{i=1}^j f_i \times g_{j-i} - \sum_{i = j}^{n} f_i\times g_{i-j} \]

假设 \(f_0 = 0,g_0 = 0\)，这样我们就能消去我们加上的那一项的影响了，并且左边的 \(\sum\) 范围也可以再做一下变换

\[F_j = \sum_{i=0}^j f_i \times g_{j-i} - \sum_{i = j}^{n} f_i\times g_{i-j} \]

你会发现，左边的 \(\sum\) 已然是一个卷积 的形式，这意味着我们可以用 FFT 来求出，接下来我们看右边的怎么化简。

我们将右边式子展开你会发现 \(f_j \times g_0,f_{j+1} \times g_1,\dots,f_{j+(n-j)} \times g_{n-j}\)，于是式子可以化成

\[\sum_{i=0}^{n-j}f_{j+i}\times g_i \]

下面有一个套路，就是翻转的套路：设 \(f_{n-i}' = f_i\)，那么原式就变为

\[\sum_{i=0}^{n-j}f_{n-j-i}' \times g_i \]

自此你会发现，这部分也被我们转化成卷积的式子了。那么这个题的思路就出来了：

• 维护 \(f,g,f'\) 三个数组，运用 FFT 求得 \(f * g,f' * g\)。最后 \(E_i\) 就是 \([x^i](f * g) - [x^{n-i}](f' * g)\)。

代码。

• P5488 差分与前缀和

可能这个题是我理解生成函数和多项式的关系的第一道题。

给定一个长为 \(n\) 的序列 \(a\)，求出其 \(k\) 阶差分或前缀和。
结果的每一项都需要对 \(1004535809\) 取模。若 \(t=0\) 表示求前缀和，\(t=1\) 表示求差分。

在一些多项式题中，我们往往会把一个序列抽象成一个生成函数，然后按照题目的要求去乘上另一个生成函数，而将生成函数转化成封闭形式往往就是解题关键。虽然是生成函数有无穷项，可根据题目的要求，我们大多数情况下只需要对前 \(n\) 项进行多项式卷积即可，这样就能得到我们最终的答案。

首先给个公式：\(\displaystyle \binom{n}{m} = \frac{\frac{n!}{(n-m)!}}{m!}\)，也很简单。

在这个题中，我们将序列 \(a\) 看成是一个 OGF：

\[F(x) = \sum_{i=0}^{\infty}a_ix^i \]

我们先来看一下前缀和的式子：

\[s_i = \sum_{j=1}^i a_j \]

我们把它和卷积的式子比较一下：

\[s_i = \sum_{j+k=i}a_j\times b_k \]

你会发现对 \(a\) 做一次前缀和，相当于和一个系数全是 \(1\) 的多项式做一次卷积。我们把这个多项式也想成是一个 OGF，那么就是

\[G(x) = 1 + x + x^2 + \dots \]

而它可能是我们最常用的生成函数了，易知它的封闭形式就是 \(\displaystyle \frac{1}{1-x}\)。

所以说做一次前缀和后，答案就是

\[F(x) \times G(x) = F(x) \times \frac{1}{1-x} \]

而二阶前缀和，就是在此基础上再乘上一个 \(\displaystyle \frac{1}{1-x}\)，以此类推， \(k\) 阶前缀和的式子就是

\[F(x) \times (\frac{1}{1-x})^k \]

然后我们需要用广义二项式定理把它拆开

\[\begin{aligned} &\ \ \ \ \ \ (\frac{1}{1-x})^k\\ &=(1-x)^{-k}\\ &=\sum_{i=0}^{\infty}\binom{-k}{i}(-x)^i\\ &=\sum_{i=0}^{\infty}(-1)^i\frac{\frac{(-k)!}{(-k-i)!}}{i!}x^i \end{aligned} \]

\[=\sum_{i=0}^{\infty}(-1)^i\frac{(-k)(-k-1)\dots(-k-i+1)}{i!}x^i \]

然后你会发现，无论 \(i\) 是奇数还是偶数，负号都能消掉，于是式子变成：：

\[=\sum_{i=0}^{\infty}\frac{k(k+1)\dots(k+i-1)}{i!}x^i \]

\[=\sum_{i=0}^{\infty}\frac{(k+i-1)!}{(k-1)!i!}x^i \]

\[=\sum_{i=0}^{\infty}\binom{k+i-1}{i}x^i \]

到此，你就求得了 \(k\) 阶前缀和前面的系数，因为 \(k\) 过大，可以考虑递推求解，递推公式也很显然：

\[\binom{k+i-1}{i} = \binom{k+i-2}{i-1}\times \frac{k+i-1}{i} \]

有了多项式的系数，和 \(F(x)\) 卷一下答案就出来了。接下来我们考虑差分。

因为差分是前缀和的逆运算，所以差分的 OGF 就是 \((1-x)\)，之后的推导就和上面差不多了：

\[F(x) \times (1-x)^k = F(x) \times \sum_{i=0}^{\infty}\binom{k}{i}(-x)^i \]

\[= F(x) \times \sum_{i=0}^{\infty}(-1)^i\binom{k}{i}x^i \]

这就是最简形式，然后稍微推一推推出来递推式：

\[\binom{k}{i} = \binom{k}{i-1}\times \frac{(k-i+1)}{i} \]

求出前 \(n\) 项即可。

于是，这个题就做完了，
代码。

通过这个题，我们可以了解生成函数和多项式结合的初步应用，这需要我们对生成函数的封闭形式以及一些组合数学的知识十分熟悉。感觉现在才刚刚入了多项式的门，不得不说，数学真是奇妙。