多项式 I：拉格朗日插值与快速傅里叶变换

1. 复数和单位根

前置知识：弧度制，三角函数。

1.1 复数的引入

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

像这种从 $a$ 变成 $a+bx$ 的 扩域 操作并不少见，例如初中学习 “平方根” 时，经常用 $a+b\sqrt{x}(x>0)$ 表示一个数。这类数的加减乘都是容易的，除法即考虑平方差公式 $(c+d\sqrt{x})(c-d\sqrt{x})=c^2-d^{2}x$，因此 $\frac{a+b\sqrt{x}}{c+d\sqrt{x}}=\frac{(ac-bdx)+(bc-ad)\sqrt{x}}{c^2-d^{2}x}$。

将 $x$ 替换为 $-1$，复数四则运算可由实数四则运算结合 $i$ 的定义直接推广得到。

• 加法：$(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$。
• 除法：$\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$。

1.2 复平面与棣莫弗定理

描述一个复数 $a+bi$ 需要两个值 $a$ 和 $b$，其中 $a$ 表示 实部，$b$ 表示 虚部。这启发我们将其放在平面直角坐标系上描述，称为 复平面。其在复平面上的坐标为 $(a,b)$，实部 $a$ 为横坐标，虚部 $b$ 为纵坐标。

一个复数唯一对应一个平面向量，因为它们都可以用有序实数对描述。将向量起点平移至原点，则其终点指向与其对应的复数。考虑平面向量的一些特征，如其长度与所在直线斜率。将这些概念应用在复数上，我们得到如下定义：

• 定义复数 $z=a+bi$ 的 模 为 $|z|=\sqrt{a^2+b^2}$。
• 定义复数 $z=a+bi$ 的 辐角 为 $\mathrm{Arg}z=\theta$，其中 $\tan \theta =\frac{b}{a}$。满足 $-\pi <\theta \le \pi$ 的 $\theta$ 称为 辐角主值，记作 $\arg z$，即 $\arg z=\arctan \frac{b}{a}$。
• 辐角确定了 $z$ 所在直线，模确定了 $z$ 在直线上的长度。对比实部和虚部，模和辐角主值以另一种有序实数对的形式描述复数。

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

利用 $\cos$ 和 $\sin$ 的和角公式可得 $z_1{z}_2=r_1{r}_2(\cos ({\theta }_1+{\theta }_2)+i\sin ({\theta }_1+{\theta }_2))$。该等式称为 棣莫弗定理，它说明 复数相乘，模长相乘，辐角相加。

• 根据棣莫弗定理，我们得到对虚数单位 $i$ 的一种直观理解：将一个复数 $z$ 乘以 $i$ 相当于将其 逆时针旋转 $\frac{\pi }{2}$ 弧度。实际上，考虑虚数单位 $i$ 本身在复平面上的位置，发现就是 $1$ 逆时针旋转 $\frac{\pi }{2}$ 度。一般地，有旋转的地方就有 $i$ 存在，$i$ 即旋转。推荐观看：虚数的来源 - Veritasium。

1.3 单位根的定义

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

考虑将单位圆 $n$ 等分（如下图），取任意 $n$ 等分点 $P_k(0\le k<n)$，将其旋转 $n$ 次均得到 $1$，因为跨过的 $\frac{1}{n}$ 单位圆弧数为 $n$ 的倍数。这说明 $P_k^n=1$。

探究 $P_k$ 的表达式：因为它相当于从 $1$ 开始在单位圆上旋转 $\frac{2k\pi }{n}$ 弧度，因此 $P_k=\cos \left(\frac{2k\pi }{n}\right)+\sin \left(\frac{2k\pi }{n}\right)$。我们称所有 $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=0$ 的解。除特殊说明外，一般 $n$ 次单位根直接代指 ${\omega }_n$，即从 $1$ 开始逆时针方向的第一个单位根。

可以观看 3b1b 的视频 从 6:00 到 7:30 的部分加深对上述内容的理解。

• 单位根的循环性：由 ${\omega }_n^n=1$ 单位根的指数可对 $n$ 取模。换言之，考虑从 $1$ 开始不断乘以 ${\omega }_n$，我们将得到 $1,{\omega }_n,{\omega }_n^2,\cdots ,{\omega }_n^{n-1},{\omega }_n^n=1,{\omega }_n^{n+1}={\omega }_n,\cdots$，循环节为 $n$。想象从 $1$ 开始不断旋转 $\frac{2\pi }{n}$ 弧度，旋转 $n$ 次后我们将落入循环。换言之，${\omega }_n^k={\omega }_n^{k+tn}(t\in \mathbb{Z})$。
• ${\omega }_n^k={\omega }_{cn}^{ck}(c>0)$。对单位根有可视化认知后容易理解这一点。
• 当 $n$ 为偶数时，将 ${\omega }_n^k$ 取反相当于将其逆时针（或顺时针）转半圈，所以 $-{\omega }_n^k={\omega }_n^{k\pm \frac{n}{2}}(2\mid n)$。
• 单位根的对称性：因为 $n$ 次单位根将单位圆 $n$ 等分，均匀分布在圆周，所以它们的重心就在原点，即 $\sum _{i=0}^{n-1}{\omega }_n^i=0$。
• 若 $gcd(k,n)=1$，则 ${\omega }_n^k$ 称为 本原单位根。所有本原单位根的 $0\sim n-1$ 次幂互不相同。

1.4 单位根与原根

前置知识：原根，详见 初等数论学习笔记 I：同余相关。

单位根和原根有极大的相似性，可以说原根就是模 $P$ 下的整数域的单位根。

设 $n=\phi (P)$ 且 $P$ 存在原根 $g$，则 $g^0,g^1,\cdots ,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=\phi (P)$ 的 循环群 的 生成元。它们均满足 $n$ 次幂是对应域的单位元，且它们的 $0\sim n-1$ 次幂互不相同。换言之，它们 同构。

快速傅里叶变换 FFT 和快速数论变换 NTT 的联系恰在于此。

1.5 欧拉公式

前置知识：自然对数 $\ln$ 的底数 $e$ 及其相关性质。

这部分为接下来可能用到的符号做一些补充。

欧拉公式指出

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

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

欧拉公式的严谨证明超出了讨论范围，相关资料可以参考百度百科。我们给出一个对欧拉公式的感性理解方式，以加深读者对欧拉公式的直观印象与理解，来自 在 3.14 分钟内理解 $e^{i\pi }$ - 3Blue1Brown。这个视频相当有启发性。

根据 $(e^t{)}^{\prime }=e^t$，考虑 $e^t$ 随着 $t$ 增大在数轴上的取值。$e^t$ 以每秒钟 $t$ 均匀增加 $1$ 的速度向数轴正方向移动，则 $e^t$ 的速度就是它本身的位置。它的起始点为 $e^0=1$。

根据 $(e^{kt}{)}^{\prime }=k{e}^t$，类似可知 $e^{kt}$ 的速度等于 $k$ 倍它本身的位置。

当 $k=i$ 时，速度相当于将本身的位置逆时针旋转 $\frac{\pi }{2}$ 弧度，与位置垂直。这样，根据基础的高中物理知识，我们知道 $e^{it}$ 随着 $t$ 增大，在单位圆上做匀速圆周运动，且每秒移动 $1$ 弧度。

这样，$e^{it}$ 就等于模长为 $1$，辐角为 $t$ 的复数，即 $\cos t+i\sin t$。这使得我们可以用 $r{e}^{i\theta }$ 描述模长为 $r$，辐角为 $\theta$ 的复数，记号变得更简洁。

将该表示法应用至单位根，可知 ${\omega }_n=e^{\frac{2\pi i}{n}}$。

读者应当在 $r{e}^{i\theta }$，$r(\cos \theta +i\sin \theta )$ 及其在复平面上的位置建立直观联系，有助于理解接下来的内容。

带入 $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}$ 相等。

2. 多项式

2.1 基本概念

形如 $\sum _{i=0}^n{a}_i{x}^i$ 的 有限和式 称为 多项式，记作 $f(x)=\sum _{i=0}^n{a}_i{x}^i$。其中，$a_i$ 称为 $i$ 次项的 系数，也称 $x^i$ 前的系数，记作 $[x^i]f(x)$。超出最高次数 $n$ 的系数 $a_i(i>n)$ 视作 $0$。

当项数无限时，和式形如 $\sum _{i=0}^{\mathrm{\infty }}{a}_i{x}^i$，称为 形式幂级数，它暂时超出了我们的讨论范围。

• 多项式 系数非零 的最高次项的次数称为该多项式的 度，也称次数，记作 $\mathrm{deg}f$。如 $f(x)=\sum _{i=0}^n{a}_i{x}^i$ 其中 $a_n\ne 0$，则 $f$ 为 $n$ 次多项式，记作 $\mathrm{deg}f=n$。
• 使得 $f(x)=0$ 的所有 $x$ 称为多项式的 根。
• 若 $a_i$ 均为实数，则称 $f$ 为实系数多项式。若 $a_i$ 可以均为复数，则称 $f$ 为复系数多项式。
• 代数基本定理：任何非零一元 $n$ 次复系数多项式恰有 $n$ 个复数根。这些复数根可能重合。证明略。

2.2 系数表示法和点值表示法

$f(x)=\sum _{i=0}^n{a}_i{x}^i$ 给出了所有 $i$ 次项前的系数，这种描述多项式的方法称为 系数表示法。

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

考虑这两种表示法之间的联系。我们尝试探究在给定 $n$ 个点值 $(x_0,y_0),(x_1,y_1),\cdots ,(x_{n-1},y_{n-1})$ 其中 $x_i$ 互不相同时，所唯一确定的多项式的最高次数。

一个自然的猜测是 $n-1$，因为 $n-1$ 次多项式需要 $n$ 个系数描述，具有 $n$ 单位信息，而 $n$ 个点值同样具有 $n$ 单位信息。

证明即考虑直接带入，得到 $n$ 元线性方程组：对于任意 $0\le j<n$，$\sum _{i=0}^{n-1}{a}_i{x}_j^i=y_j$。该线性方程组的系数矩阵为

$$\left[\begin{array}{ccccc}1& x_0& x_0^1& \cdots & x_0^{n-1}\\ 1& x_1& x_1^1& \cdots & x_1^{n-1}\\ 1& x_2& x_2^1& \cdots & x_2^{n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_{n-1}& x_{n-1}^2& \cdots & x_{n-1}^{n-1}\end{array}\right]$$

因 $x_i$ 互不相同，所以该范德蒙德矩阵的行列式非零，方程组有唯一解。类似的，从线性方程组的角度出发，容易证明 $n$ 个点值不可能唯一确定 $n$ 次或更高次的多项式。

综上，我们得到重要结论：$n$ 个点值唯一确定的多项式最高次数为 $n-1$。

• 从系数表示法转为点值表示法称为 求值（Evaluation）。
• 从点值表示法转为系数表示法称为 插值（Interpolation）。

2.3 多项式运算

2.3.1 基本运算

设 $f(x)=\sum _{i=0}^n{a}_i{x}^i$，$g(x)=\sum _{i=0}^m{b}_i{x}^i$。

• 设 $h=f+g$，则 $h(x)=(\sum _{i=0}^n{a}_i{x}^i)+(\sum _{j=0}^m{b}_j{x}^j)=\sum _{i=0}^{max(n,m)}(a_i+b_i)x^i$，可知两多项式相加，对应系数相加，$\mathrm{deg}(f+g)=max(\mathrm{deg}f,\mathrm{deg}g)$。
• 设 $h=f\ast g$，则 $h(x)=(\sum _{i=0}^n{a}_i{x}^i)(\sum _{j=0}^m{b}_j{x}^j)=\sum _{i=0}^{n+m}(\sum _{j=0}^i{a}_j{b}_{i-j})x^i$，可知两多项式相乘，每两个系数相乘贡献至次数之和的系数，$\mathrm{deg}(f\ast g)=\mathrm{deg}f+\mathrm{deg}g$。

因此，在系数表示法下，计算两个多项式相加的复杂度为 $\mathcal{O}(n+m)$，计算两个多项式相乘的复杂度为 $\mathcal{O}(nm)$。

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

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

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

• 用 $f\ast g$ 和 $fg$ 表示多项式相乘，即进行加法卷积；用 $f\cdot g$ 表示多项式 点乘，即 对应系数相乘。

在系数表示法下计算两个多项式相乘，我们首先将其转化为点值表示法，相乘后再转回系数表示法。时间复杂度 $\mathcal{O}((n+m)\log (n+m))$。相关算法将在第四章介绍。

3. 拉格朗日插值

在 2.2 小节我们得到结论：$n$ 个点值唯一确定的多项式最高次数为 $n-1$。那么，如何在点值表示法和系数表示法之间快速转换呢？

从系数表示法转为点值表示法，最常用的方法是直接带入，时间复杂度 $\mathcal{O}(n^2)$。$\mathcal{O}(n{\log }^{2}n)$ 的多项式多点求值则需要高级多项式技巧，此处不作介绍。

从点值表示法转为系数表示法，最直接的方法是高斯消元，时间复杂度 $\mathcal{O}(n^3)$。接下来我们将介绍常用的拉格朗日插值法。

3.1 算法详解

拉格朗日插值的核心思想在于利用点值的可加性，每次只考虑一个点值，其它点值均视为 $0$，由此构造 $n$ 个多项式 $f_i(x)$，使得它们在 $x_i$ 处取值为 $y_i$，$x_j(j\ne i)$ 处取值为 $0$，则 $f=\sum _{i=0}^{n-1}{f}_i$。中国剩余定理 使用了类似的思想。

为了让 $f_i(x_j)=0(i\ne j)$，$f_i$ 应包含 $x-x_j$ 作为因式，因此设 $f_i(x)=\prod _{i\ne j}(x-x_j)$。但是此时 $f_i(x_i)$ 不一定等于 $y_i$，我们发现可以调整 $f_i$ 的系数，给 $f_i$ 乘上 $\frac{y_i}{f_i(x_i)}$。综上，我们得到 $f_i$ 的表达式

$$f_i(x)=y_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}$$

对 $f_i$ 求和得 $f$，我们得到拉格朗日插值公式

$$f(x)=\sum _{i=0}^{n-1}{y}_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}$$

为得到 $f$ 的各项系数，只需 $\mathcal{O}(n^2)$ 求出 $F(x)=\prod _{i=0}^{n-1}(x-x_i)$，对每个 $i$ 暴力 $\mathcal{O}(n)$ 将 $F$ 除掉一次式 $x-x_i$ 算出 $\frac{F(x)}{x-x_i}$ 的各项系数，再乘以 $\frac{y_i}{\prod _{j\ne i}{x}_i-x_j}$ 得到 $f_i$，则 $f=\sum _{i=0}^{n-1}{f}_i$。时间复杂度 $\mathcal{O}(n^2)$。

通常情况下，题目要求我们求出 $F(x)$ 在给定某个 $x$ 处的取值，此时我们不把 $x$ 看做函数的元，而是直接将 $x$ 带入上式。时间复杂度仍为 $\mathcal{O}(n^2)$。

多项式快速插值在 $\mathcal{O}(n{\log }^{2}n)$ 的时间内将点值表示法转化为系数表示法，这超出了我们的讨论范围。

3.2 连续取值插值

当给定点值横坐标为连续整数时，我们有快速单点插值的方法。

以 $0,1,\cdots ,n-1$ 即 $x_i=i$ 为例：

$$\begin{array}{r}f(x)=\sum _{i=0}^{n-1}{y}_i\prod _{j\ne i}\frac{x-j}{i-j}\end{array}$$

分子是 $\prod (x-i)$ 挖掉一个项，经典维护前缀后缀积。设 $p_i=\prod _{j=0}^{i-1}x-j$，$s_i=\prod _{j=i+1}^{n-1}x-j$。

分母对于 $i>j$，$\prod _{j=0}^{i-1}(i-j)=i!$。对于 $i<j$，$\prod _{j=i+1}^{n-1}(i-j)=(-1)(-2)\cdots (i-n+1)$，将所有负号提出来，得 $(-1{)}^{n-i+1}(i-n+1)!$。

因此

$$f(x)=\sum _{i=0}^{n-1}{y}_i\frac{p_i{s}_i}{(-1{)}^{n-i+1}i!(i-n+1)!}$$

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

3.3 应用

• 计算范德蒙德方阵的逆矩阵，详见 4.4.1 小节。

3.4 例题

P4781 【模板】拉格朗日插值

#include <bits/stdc++.h>
using namespace std;
constexpr int N = 2e3 + 5;
constexpr int mod = 998244353;
int ksm(int a, int b) {
  int s = 1;
  while(b) {
    if(b & 1) s = 1ll * s * a % mod;
    a = 1ll * a * a % mod, b >>= 1;
  }
  return s;
}
int n, k, ans, x[N], y[N];
int main() {
  cin >> n >> k;
  for(int i = 1; i <= n; i++) cin >> x[i] >> y[i];
  for(int i = 1; i <= n; i++) {
    int deno = 1, nume = 1;
    for(int j = 1; j <= n; j++) {
      if(i == j) continue;
      deno = 1ll * deno * (x[i] + mod - x[j]) % mod;
      nume = 1ll * nume * (k + mod - x[j]) % mod;
    }
    ans = (ans + 1ll * y[i] * nume % mod * ksm(deno, mod - 2)) % mod;
  }
  cout << ans << "\n";
  return 0;
}

4. 快速傅里叶变换

快速傅里叶变换（Fast Fourier Transform，FFT）是多项式算法的根基。想要真正理解它，必须先真正理解单位根，还需要对线性代数有基本了解。

推荐观看：

• 这个算法改变了世界 - Veritasium。
• The Fast Fourier Transform (FFT): Most Ingenious Algorithm Ever? - Reducible & B 站搬运。

4.1 求值的本质

设 $f(x)=\sum _{i=0}^{n-1}{a}_i{x}^i$，将 $x_0$ 带入，得 $f(x_0)=\sum _{i=0}^{n-1}{a}_i{x}_0^i$。考虑将其写成向量内积（点积）的形式：

$$f(x_0)=\sum _{i=0}^{n-1}{a}_i{x}_0^i=\left[\begin{array}{cccc}{x}_0^0& x_0^1& \cdots & x_0^{n-1}\end{array}\right]\times \left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\end{array}\right]$$

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

$$\left[\begin{array}{cccc}{x}_0^0& x_0^1& \cdots & x_0^{n-1}\\ x_1^0& x_1^1& \cdots & x_1^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ x_{m-1}^0& x_{m-1}^1& \cdots & x_{m-1}^{n-1}\end{array}\right]\times \left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}f(x_0)\\ f(x_1)\\ ⋮\\ f(x_{m-1})\end{array}\right]$$

左侧矩阵就是著名的 范德蒙德矩阵。

当 $m=n$ 时为范德蒙德方阵，$x_i$ 互不相同时其逆存在，帮助我们快速从点值表示法转回系数表示法。$m>n$ 时任取 $x_i$ 互不相同的 $n+1$ 行可以求逆。$m<n$ 时无法还原系数。这体现出 $n+1$ 个点值唯一确定最高次数不超过 $n$ 的多项式。

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

4.2 离散傅里叶变换

在介绍 FFT 之前，我们先给出离散傅里叶变换（Discrete Fourier Transform，DFT）的概念。

DFT 在工程中是将离散信号从时域转为频域的过程。碰巧的是，其表达式刚好可以用来对多项式进行多点求值，只不过这些点值是固定的 单位根 处的点值，但对于求值做多项式乘法已经足够了。

DFT 将一个长度为 $N$ 的复数序列 $x_0,x_1,\cdots ,x_{N-1}$ 通过如下公式转化为另一个长为 $N$ 的复数序列 $X_0,X_1,\cdots ,X_{N-1}$：

$$X_k=\sum _{n=0}^{N-1}{x}_n{e}^{-\frac{2\pi i}{N}kn}$$

也即

$$X_k=\sum _{n=0}^{N-1}{x}_n{\omega }_N^{-kn}$$

设多项式 $f(x)=\sum _{n=0}^{N-1}{x}_n{x}^i$，易知

$$X_k=\sum _{n=0}^{N-1}{x}_n({\omega }_N^{-k}{)}^n=f({\omega }_N^{-k})$$

根据上式，离散傅里叶变换可以看成多项式 $f(x)=\sum _{n=0}^{N-1}{x}_n{x}^i$ 在 $N$ 个单位根处求值。

没看懂？没关系。接下来我们将从另一角度出发，得到一个差别微小但更容易理解的算法 —— FFT。

4.3 算法详解

首先我们得搞清楚，DFT 是一个变换，而 FFT 是用于实现 DFT 的算法。在 FFT 的推导过程中，其所实现的变换和 DFT 变换有细微的差别，因此笔者也用 FFT 表示该算法实现的变换。

按理说学习一个算法时需要搞清楚它的用途，但如果直接从 DFT 角度入手尝试简化流程，那么为了让动机看起来自然，我们还要先学习 DFT 的实际含义，这涉及到大量前置知识。

但是，从计算机科学界尤为重要且为各位 OIer 熟知的多项式乘法的优化出发，我们通过自然推理得到一个算法，而该算法恰好可以快速实现 DFT（有一些细小的差别，详见 4.3.5）。

我们明确接下来的目标：给定次数 $n-1$ 的多项式 $f(x)=\sum _{i=0}^{n-1}{a}_i{x}^i(a_{n-1}\ne 0)$，快速求出它的至少 $n$ 个点值。

下文中，$f(x)$ 也被视为关于 $x$ 的 $n-1$ 次函数。

4.3.1 简化情况

解决一个普遍性的问题，最重要的思想就是 从简化情况入手，分析问题的性质。

函数的性质无非就那几个：奇偶性，单调性，周期性等。一般函数没有单调性和周期性，但总可以表示为一个偶函数和一个奇函数之和，这一定是突破点。

• 对于偶函数 $f(x)$，即所有奇数次项系数为 $0$，带入两个相反数 $w$ 和 $-w$ 时，它们的值相等。

• 对于奇函数 $f(x)$，即所有偶数次项系数为 $0$，带入两个相反数 $w$ 和 $-w$ 时，它们的值互为相反数。

因此，将任意多项式 $f(x)$ 拆成偶函数 $f_e(x)$ 和奇函数 $f_o(x)$ 之和，则

$$\{\begin{array}{l}f(x)=f_e(x)+f_o(x)\\ f(-x)=f_e(x)-f_o(x)\end{array}$$

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

不妨设 $n$ 为偶数，则 $f_e$ 是 $n-2$ 次多项式，$f_e$ 是 $n-1$ 次多项式，本质上依然相当于求出 $n-1$ 次多项式的 $n$ 个点值，对时间复杂度没有优化。

但是 $f_e$ 和 $f_o$ 的项数减半，尝试利用该性质。

因为 $f_e$ 只有偶数次项 $a_0{x}^0+a_2{x}^2+\cdots$，故考虑换元 $u=x^2$，则 $f_e(u)=a_0{u}^0+a_2{u}^1+\cdots$。换言之，我们设 $f_e^{\prime }(x)=a_0{x}^0+a_2{x}^1+\cdots$，则 $f_e(x)=f_e^{\prime }(x^2)$。

同理，从 $f_o$ 中提出一个 $x$，设 $f_o^{\prime }(x)=a_1{x}^0+a_3{x}^1+\cdots$，则 $f_o(x)=x{f}_o^{\prime }(x^2)$。

因此，

$$\{\begin{array}{l}f(x)=f_e^{\prime }(x^2)+x{f}_o^{\prime }(x^2)\\ f(-x)=f_e^{\prime }(x^2)-x{f}_o^{\prime }(x^2)\end{array}$$

这样才是真正意义上的规模减半。若问题可递归，则时间复杂度为 $T(n)=2T\left(\frac{n}{2}\right)+\mathcal{O}(n)$，解得 $T(n)=\mathcal{O}(n\log n)$。

4.3.2 引入单位根

问题来了，如何保证递归得到的问题也满足两两互为相反数呢？

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

进一步推演，因为问题会递归 $w=\lceil {\log }_{2}n\rceil$ 层，所以我们需要找到 $k=2^w$ 个互不相等的 $x$，但它们的 $k$ 次幂相等。这个相等的 $k$ 次幂可以任意选择，方便起见设为 $1$，则 $x^k=1$，$x$ 为所有 $k$ 次单位根。

逆推得到结果后，我们再顺着检查一遍：初始令 $x$ 为所有 $k$ 次单位根，每层递归会将这些数平方，得到所有 $\frac{k}{2},\frac{k}{4}\cdots$ 次单位根。因为 $\frac{k}{2},\frac{k}{4},\cdots$ 都是偶数，所以它们可以两两正负配对，直到递归 $w$ 层的平凡情况：$\frac{k}{2^w}=1$ 次单位根即 $x=1$。

由此可得通常使用（补齐到 $2$ 的幂）的快速傅里叶变换的范德蒙德方阵形式：

$$\left[\begin{array}{cccc}({\omega }_n^0{)}^0& ({\omega }_n^0{)}^1& \cdots & ({\omega }_n^0{)}^{n-1}\\ ({\omega }_n^1{)}^0& ({\omega }_n^1{)}^1& \cdots & ({\omega }_n^1{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{n-1}{)}^0& ({\omega }_n^{n-1}{)}^1& \cdots & ({\omega }_n^{n-1}{)}^{n-1}\end{array}\right]\times \left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}f({\omega }_n^0)\\ f({\omega }_n^1)\\ ⋮\\ f({\omega }_n^{n-1})\end{array}\right]$$

4.3.3 递归公式

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

首先将 $n$ 补齐到不小于 $n$ 的最小的 $2$ 的幂，即 $2^{\lceil {\log }_{2}n\rceil }$。

设当前需要求值的多项式 $f$ 长度为 $L(L=2^w,w\in \mathbb{N})$，若 $L=1$ 则直接返回 $a_0$。否则我们需求出 $f(x)$ 在所有 $L$ 次单位根 ${\omega }_L^k(0\le k<L)$ 处的取值。

将 $f(x)$ 写成 $f_e(x^2)+x{f}_o(x^2)$，则对于 $0\le k<\frac{L}{2}$，有

$$\begin{array}{rl}f({\omega }_L^k)& =f_e({\omega }_L^{2k})+{\omega }_L^k{f}_o({\omega }_L^{2k})\\ f({\omega }_L^{k+\frac{L}{2}})& =f_e({\omega }_L^{2k+L})+{\omega }_L^{k+\frac{L}{2}}{f}_o({\omega }_L^{2k+L})\end{array}$$

根据单位根的性质（在单位根部分有介绍）：

• ${\omega }_n^k={\omega }_{2n}^{2k}$。
• ${\omega }_n^k={\omega }_n^{k+tn}(t\in \mathbb{Z})$。
• ${\omega }_n^k=-{\omega }_n^{k+\frac{n}{2}}(2\mid n)$。

得

$$\begin{array}{rl}f({\omega }_L^k)& =f_e({\omega }_{\frac{L}{2}}^k)+{\omega }_L^k{f}_o({\omega }_{\frac{L}{2}}^k)\\ f({\omega }_L^{k+\frac{L}{2}})& =f_e({\omega }_{\frac{L}{2}}^k)-{\omega }_L^k{f}_o({\omega }_{\frac{L}{2}}^k)\end{array}$$

这样，我们只需求出 $\frac{L}{2}$ 次多项式 $f_e$ 和 $f_o$ 在所有 $\frac{L}{2}$ 次单位根处的取值。

4.3.4 蝴蝶变换

递归处理比较慢，我们希望像位运算卷积一样通过递推实现整个过程。

因为在边界条件下，每个位置的取值与其对应的系数相关，所以递归转递推的关键在于考察系数的变化。

考虑 $n=8$ 的情况，初始为 $(a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7)$。

进行第一层递归时，将 $f_e$ 的系数写在左半部分，$f_o$ 的系数写在右半部分，得 $(a_0,a_2,a_4,a_6),(a_1,a_3,a_5,a_7)$。

进行第二层递归时，类似地将每个子问题的 $f_e$ 和 $f_o$ 的系数分别写在左右两侧，得 $(a_0,a_4),(a_2,a_6),(a_1,a_5),(a_3,a_7)$。

进行第三层递归时，共有 $4$ 个大小为 $2$ 的子问题，且进行上述操作后系数的位置不发生改变。

我们看到，如果一个系数在规模为 $2n$ 的问题中的位置为 $p$，若 $p$ 是偶数，则递归至左半部分，若 $p$ 是奇数，则递归至右半部分，且在规模为 $n$ 的问题中的位置为 $\lfloor \frac{p}{2}\rfloor$。

进一步地，一个系数在第 $i$ 次递归时的方向决定了它最终下标在二进制下第 $w-i(n=2^w)$ 位的取值。若向左递归则为 $0$，向右递归则为 $1$。而它递归的方向又受到 $\lfloor \frac{p}{2^{i-1}}\rfloor$ 的奇偶性的影响，即 $p$ 在二进制下第 $i$ 位的取值，若为 $0$ 则向左递归，为 $1$ 则向右递归。

这就说明，$p$ 在二进制下第 $i$ 位的取值，等于它对应的系数的最终下标在二进制下第 $w-i$ 位的取值。

因此我们断言，系数 $a_p$ 在 $w$ 次递归后的下标等于 $p$ 二进制翻转后的值，设为 $r_p$。这里翻转指翻转第 $0\sim w-1$ 位的值。

$r_p$ 可递推求得：$r_0=0$。对于 $r_i(i>0)$，先右移一位得到 $i^{\prime }=\lfloor \frac{i}{2}\rfloor$，则 $r_i$ 的低 $w-1$ 位（第 $0\sim w-2$ 位）的值可由 $r_{i^{\prime }}$ 右移一位得到。第 $w-1$ 位的值只需检查 $i$ 的奇偶性。因此，$r_i=\lfloor \frac{r_{i^{\prime }}}{2}\rfloor +\frac{n}{2}(imod2)$，其中 $i^{\prime }=\lfloor \frac{i}{2}\rfloor$。

因此，在算法一开始先对系数数组 $a$ 执行 蝴蝶变换，即同时令 $a_i\to a_{r_i}$，然后类似 FWT，枚举问题规模 $2d(1\le d<n)$，枚举每个子问题 $2di(0\le 2di<n)$，再枚举当前子问题的所有对应位置 $(x=2di+k,y=2di+k+d)(0\le k<d)$ 执行变换，即设 $x$ 和 $y$ 对应位置的当前值为 $f_x$ 和 $f_y$，则 $f_x^{\prime }=f_x+{\omega }_{2d}^k{f}_y$，$f_y^{\prime }=f_x-{\omega }_{2d}^k{f}_y$。

进一步地，根据 $r$ 的定义，我们有 $r_{r_i}=i$。因此，执行蝴蝶变换只需对所有 $i<r_i$ 执行 $\mathrm{swap}(a_i,a_{r_i})$。

这就是 FFT，整个过程称为 对多项式 $f$ 做长度为 $n$ 的快速傅里叶变换，时间复杂度 $\mathcal{O}(n\log n)$。代码在 4.5 小节。

4.3.5 DFT 和 FFT

对比 DFT 和 FFT 矩阵：

$$\begin{gathered} {\mathcal{F}}_{\mathcal{D}}=\left[\begin{array}{cccc}({\omega }_N^0{)}^0& ({\omega }_N^0{)}^1& \cdots & ({\omega }_N^0{)}^{N-1}\\ ({\omega }_N^{-1}{)}^0& ({\omega }_N^{-1}{)}^1& \cdots & ({\omega }_N^{-1}{)}^{N-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{-(N-1)}{)}^0& ({\omega }_n^{-(N-1)}{)}^1& \cdots & ({\omega }_N^{-(N-1)}{)}^{N-1}\end{array}\right] \\ {\mathcal{F}}_{\mathcal{F}}=\left[\begin{array}{cccc}({\omega }_n^0{)}^0& ({\omega }_n^0{)}^1& \cdots & ({\omega }_n^0{)}^{n-1}\\ ({\omega }_n^1{)}^0& ({\omega }_n^1{)}^1& \cdots & ({\omega }_n^1{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_n^{n-1}{)}^0& ({\omega }_n^{n-1}{)}^1& \cdots & ({\omega }_n^{n-1}{)}^{n-1}\end{array}\right] \end{gathered}$$

可以发现 DFT 和 FFT 基本一致，它们的差别在于：

• 朴素 FFT 要求 $n$ 是 $2$ 的幂，但 DFT 序列长度可以是任意正整数。
• DFT 和 FFT 带入单位根的顺序是相反的。

注意这些细微差别，不要把 DFT 和 FFT 搞混了。

4.3.6 循环卷积

4.4 离散傅里叶逆变换

离散傅里叶逆变换（Inverse Discrete Fourier Transform，IDFT）可以视为单位根处插值的过程。即给出 $n=2^w$ 个在所有 $n$ 次单位根处的点值 $P_k=({\omega }_n^k,f({\omega }_n^k))(0\le k<n)$，要求还原 $f$ 的各项系数，其中 $f$ 的次数不大于 $n-1$。

类似地，IDFT 和 IFFT 之间也存在一些差异。想必各位读者根据之前的内容已经可以猜出这种差异是什么了。

4.4.1 范德蒙德方阵逆

考虑范德蒙德方阵

$$\mathcal{A}=\left[\begin{array}{cccc}{x}_0^0& x_0^1& \cdots & x_0^{n-1}\\ x_1^0& x_1^1& \cdots & x_1^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ x_{n-1}^0& x_{n-1}^1& \cdots & x_{n-1}^{n-1}\end{array}\right]$$

这玩意怎么求逆？我们考虑最暴力的方法：拉格朗日插值！

因为范德蒙德方阵可以看成多项式多点求值，根据

$$\left[\begin{array}{cccc}{x}_0^0& x_0^1& \cdots & x_0^{n-1}\\ x_1^0& x_1^1& \cdots & x_1^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ x_{n-1}^0& x_{n-1}^1& \cdots & x_{n-1}^{n-1}\end{array}\right]\times \left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}f(x_0)\\ f(x_1)\\ ⋮\\ f(x_{n-1})\end{array}\right]$$

再结合拉格朗日插值公式

$$f(x)=\sum _{i=0}^{n-1}f(x_i)\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}$$

可知从 $f(x_j)$ 贡献到 $a_i$ 的系数为 $({\mathcal{A}}^{-1}{)}_{ij}=[x^i]\prod _{k\ne i}\frac{x-x_k}{x_j-x_k}$。

这就是范德蒙德方阵逆当中每个元素的表达式。

4.4.2 算法介绍

很显然，$f$ 在经过快速傅里叶变换后，再进行快速傅里叶逆变换，仍得到 $f$。

因此，只需对快速傅里叶变换的矩阵求逆，即可得到快速傅里叶逆变换的矩阵。

设 $\omega$ 表示 ${\omega }_n$，则

$$\mathcal{F}=\left[\begin{array}{cccc}({\omega }^0{)}^0& ({\omega }^0{)}^1& \cdots & ({\omega }^0{)}^{n-1}\\ ({\omega }^1{)}^0& ({\omega }^1{)}^1& \cdots & ({\omega }^1{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }^{n-1}{)}^0& ({\omega }^{n-1}{)}^1& \cdots & ({\omega }^{n-1}{)}^{n-1}\end{array}\right]$$

则 $({\mathcal{F}}^{-1}{)}_{ij}=[x^i]\prod _{k\ne j}\frac{x-{\omega }^k}{{\omega }^j-{\omega }^k}$。

分子和分母均形如 $g(x)=\frac{\prod _{0\le k<n}(x-{\omega }^k)}{x-{\omega }^j}$，我们研究该函数的性质。

首先，因为 ${\omega }^k(0\le k<n)$ 为 $x^n-1=0$ 的 $n$ 个互不相同的根，根据因式定理，$\prod _{0\le k<n}(x-{\omega }^k)$ 就等于 $x^n-1$。感性理解：将 $\prod _{0\le k<n}(x-{\omega }^k)$ 展开，再应用单位根的 对称性。

模拟短除法 $\frac{x^n-1}{x-{\omega }^j}$，可知

$$g(x)=x^{n-1}+{\omega }^j{x}^{n-2}+({\omega }^j{)}^2{x}^{n-3}+\cdots +({\omega }^j{)}^{n-1}{x}^0$$

因此

$$({\mathcal{F}}^{-1}{)}_{ij}=\frac{[x^i]g(x)}{g({\omega }^j)}=\frac{({\omega }^j{)}^{n-1-i}}{n({\omega }^j{)}^{n-1}}=\frac{({\omega }^{-j}{)}^i{\omega }^{-j}}{n{\omega }^{-j}}=\frac{{\omega }^{-ij}}{n}$$

这就说明

$${\mathcal{F}}^{-1}=\frac{1}{n}\left[\begin{array}{cccc}({\omega }^{-0}{)}^0& ({\omega }^{-0}{)}^1& \cdots & ({\omega }^{-0}{)}^{n-1}\\ ({\omega }^{-1}{)}^0& ({\omega }^{-1}{)}^1& \cdots & ({\omega }^{-1}{)}^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }^{-(n-1)}{)}^0& ({\omega }^{-(n-1)}{)}^1& \cdots & ({\omega }^{-(n-1)}{)}^{n-1}\end{array}\right]$$

因此，对一个序列做 IFFT，只需将 FFT 递归公式里面的 ${\omega }_L^k$ 换成 ${\omega }_L^{-k}$，并在最后将所有数除以 $n$ 即可。

这就引出了 IDFT 公式及其对应矩阵：

$$\begin{gathered} x_n=\frac{1}{N}\sum _{k=0}^{N-1}{X}_k{e}^{\frac{2\pi i}{N}kn}=\sum _{k=0}^{N-1}{X}_k{\omega }_N^{kn} \\ {{\mathcal{F}}_{\mathcal{D}}}^{-1}=\frac{1}{N}\left[\begin{array}{cccc}({\omega }_N^0{)}^0& ({\omega }_N^0{)}^1& \cdots & ({\omega }_N^0{)}^{N-1}\\ ({\omega }_N^1{)}^0& ({\omega }_N^1{)}^1& \cdots & ({\omega }_N^1{)}^{N-1}\\ ⋮& ⋮& \ddots & ⋮\\ ({\omega }_N^{N-1}{)}^0& ({\omega }_N^{N-1}{)}^1& \cdots & ({\omega }_N^{N-1}{)}^{N-1}\end{array}\right] \end{gathered}$$

4.5 快速多项式乘法

通过 FFT 和 IFFT，我们可以在 $\mathcal{O}(n\log n)$ 的时间内实现 $n-1$ 次多项式在系数表示法和点值表示法之间的转换。

计算两个次数分别为 $n-1$ 和 $m-1$ 的多项式 $a,b$ 相乘，设结果为 $c$，则 $c$ 是 $n+m-2$ 次多项式，我们需要 $n+m-1$ 个点值才能确定它。根据点值的性质，首先找到不小于 $n+m-1$ 的最小的 $2$ 的幂 $L$，对系数表示法的 $a,b$ 分别做长度为 $L$ 的 FFT，将对应点值相乘得到 $\hat{c}$，再对 $\hat{c}$ 做 IFFT 还原 $c$。

---

首先我们需要实现一个复数类，它支持复数的加减乘运算。也可以使用 C++ 自带复数类 std::complex<T>，用法见 CPP reference。

FFT 和 IFFT 大体上一致，只有一些细小差别。我们可以类似实现 FWT 和 IFWT 那样，用同一个函数实现它们，并用一个参数区别。

等式 ${\omega }_n=\cos (\frac{2\pi }{n})+i\sin (\frac{2\pi }{n})$ 帮助我们求出 $n$ 次单位根。

• 注意浮点数的精度。当多项式系数的绝对值较大时，需使用 long double 甚至 5.2 小节的任意模数卷积。

模板题 P3803 代码：

#include <bits/stdc++.h>
using namespace std;

constexpr int N = 1 << 21;
constexpr double pi = acos(-1);

struct comp {
  double a, b; // a + bi
  comp operator + (const comp &x) const {return {a + x.a, b + x.b};}
  comp operator - (const comp &x) const {return {a - x.a, b - x.b};}
  comp operator * (const comp &x) const {return {a * x.a - b * x.b, a * x.b + b * x.a};}
} f[N], g[N], h[N];
int n, m, r[N];

void FFT(int L, comp *f, bool type) { // L 表示长度, type = 1 表示 FFT, type = 0 表示 IFFT
  for(int i = 1; i < L; i++) {
    r[i] = (r[i >> 1] >> 1) + (i & 1 ? L >> 1 : 0);
    if(i < r[i]) swap(f[i], f[r[i]]);
  }
  for(int d = 1; d < L; d <<= 1) {
    comp wd = {cos(pi / d), sin(pi / d)}; // 2d 次单位根
    if(!type) wd.b = -wd.b; // IFFT 单位根取倒数, 相当于沿 x 轴对称
    static comp w[N];
    w[0] = {1, 0};
    for(int j = 1; j < d; j++) w[j] = w[j - 1] * wd; // 用数组记录 0 ~ d-1 次单位根, 减少复数乘法次数
    for(int i = 0; i < L; i += d << 1) {
      for(int j = 0; j < d; j++) {
        comp x = f[i | j], y = w[j] * f[i | j | d];
        f[i | j] = x + y, f[i | j | d] = x - y;
      }
    }
  }
  if(!type) for(int i = 0; i < L; i++) f[i].a /= L; // IFFT 时所有数要除以长度 L, 我们只用到了实部所以只需将实部除以 L
}

int main() {
  cin >> n >> m;
  for(int i = 0; i <= n; i++) cin >> f[i].a;
  for(int i = 0; i <= m; i++) cin >> g[i].a;
  int L = 1;
  while(L <= n + m) L <<= 1;
  FFT(L, f, 1), FFT(L, g, 1);
  for(int i = 0; i < L; i++) h[i] = f[i] * g[i];
  FFT(L, h, 0);
  for(int i = 0; i <= n + m; i++) cout << (int) (h[i].a + 0.5) << (i < n + m ? ' ' : '\n');
  return 0;
}

4.6 快速数论变换

前置知识：原根和阶。

我们将 DFT 的数值范围从复数域推广至任意存在 $n$ 次单位根 $\alpha$ 的环 $R$，其中 $\alpha$ 满足 ${\alpha }^k(0\le k<n)$ 互不相同，对 $x_0,x_1,\cdots ,x_{n-1}$ DFT 得 $X_0,X_1,\cdots ,X_{n-1}$，则

$$X_i=\sum _{j=0}^{n-1}{x}_j{\alpha }^{ij}$$

也可以视作 $X_i=f({\alpha }^i)$，其中 $f(t)=\sum _{j=0}^{n-1}{x}_j{t}^j$。

类似可知 IDFT

$$x_j=\frac{1}{n}\sum _{i=0}^{n-1}{X}_i{\alpha }^{-ij}$$

即 DFT 和 IDFT 对序列进行的变换的本质抽象。

4.6.1 算法介绍

快速数论变换即在模 $p$ 意义下的整数域 $F=\mathbb{Z}/p$ 进行的快速傅里叶变换。

设变换长度为 $n$，则需存在 $n$ 次单位根 $a$ 满足 ${\delta }_p(a)=n$。大部分题目的模数 $p$ 满足：

• $p$ 为质数。
• $2^k\mid p-1$，其中 $2^k$ 不小于最大的可能的 NTT 长度。

第一条性质保证 $p$ 存在原根 $g$ 且 $n$ 可逆，第二条性质保证存在 $n$ 次单位根。注意它不是充要条件，只是充分条件。

根据原根的性质，${\delta }_p(g)=p-1$，即 $g$ 的 $0\sim p-2$ 次幂互不相同，则 $g$ 在 $F$ 上的性质和复数域上 $p-1$ 次单位根的性质完全一致：$g^k(0\le k<p-1)$ 和 ${\omega }_{p-1}^k(0\le k<p-1)$ 形成的域是同构的，$g^k$ 等价于 ${\omega }_{p-1}^k$。

因此，$q=g^{\frac{p-1}{n}}$ 等价于 ${\omega }_{p-1}^{\frac{p-1}{n}}$ 即 ${\omega }_n$，它的 $0\sim n-1$ 次幂互不相同，即 ${\delta }_p(q)=n$，也可以通过阶的性质 ${\delta }_p(g^k)=\frac{{\delta }_p(g)}{gcd({\delta }_p(g),k)}$ 说明。

常见 NTT 模数有：

• $998244353=7\times 17\times 2^{23}+1$，有原根 $3$。它可以用来做长度不超过 $2^{23}$ 的 NTT，也是最常见的模数。
• $1004535809=479\times 2^{21}+1$，有原根 $3$。
• $469762049=7\times 2^{26}+1$，有原根 $3$。

如果不是常见模数，我们需要检验 $p$ 是否为形如 $q{2}^k+1$ 的质数，$2^k$ 是否足够大，再运用原根相关的知识枚举并判定找到任意一个原根，就可以做 NTT 了。

注意以上只是模数 $p$ 可 NTT 的充分条件，它的更弱的条件是存在 ${\delta }_p(a)=n$ 和 $n^{-1}$。如果 NTT 是正解的一部分，那么一个合格的出题人应将 $p$ 设为常见模数，因为模数不是考察重点。

4.6.2 代码实现

朴素 NTT 已经比 FFT 快了不少，因为复数运算很耗时。

接下来加入一些常数优化：

• 预处理 $2d$ 次单位根的 $0\sim d-1$ 次幂，这样对每个 $i$ 枚举 $j$ 的时候就不需要重复计算单位根的幂。
• 用 unsigned long long 和 $16$ 次模一次的技巧，减少取模次数。类似技巧也用于优化矩阵乘法。

综上，写出如下代码（依然是 模板题）：尽管题目没有要求取模，但可视为在模 $p$ 意义下进行多项式乘法，其中 $p$ 大于答案的每一项。

#include <bits/stdc++.h>
using namespace std;

using ull = unsigned long long;
constexpr int N = 1 << 21;
constexpr int mod = 998244353;
constexpr int ivg = (mod + 1) / 3;

int ksm(int a, int b) {
  int s = 1;
  while(b) {
    if(b & 1) s = 1ll * s * a % mod;
    a = 1ll * a * a % mod, b >>= 1;
  }
  return s;
}

int n, m, r[N], f[N], g[N], h[N];
void FFT(int L, int *a, bool type) {
  static ull f[N], w[N];
  for(int i = 0; i < L; i++) f[i] = a[r[i] = (r[i >> 1] >> 1) | (i & 1 ? L >> 1 : 0)];
  for(int d = 1; d < L; d <<= 1) {
    int wd = ksm(type ? 3 : ivg, (mod - 1) / (d + d));
    for(int j = w[0] = 1; j < d; j++) w[j] = 1ll * w[j - 1] * wd % mod;
    for(int i = 0; i < L; i += d << 1) {
      for(int j = 0; j < d; j++) {
        int y = w[j] * f[i | j | d] % mod;
        f[i | j | d] = f[i | j] + mod - y, f[i | j] += y;
      }
    }
    if(d == (1 << 16)) for(int p = 0; p < L; p++) f[p] %= mod;
  }
  int inv = ksm(L, mod - 2);
  for(int i = 0; i < L; i++) a[i] = f[i] % mod * (type ? 1 : inv) % mod;
}
int main() {
  cin >> n >> m;
  for(int i = 0; i <= n; i++) cin >> f[i];
  for(int i = 0; i <= m; i++) cin >> g[i];
  int L = 1;
  while(L <= n + m) L <<= 1;
  FFT(L, f, 1), FFT(L, g, 1);
  for(int i = 0; i < L; i++) h[i] = 1ll * f[i] * g[i] % mod;
  FFT(L, h, 0);
  for(int i = 0; i <= n + m; i++) cout << h[i] << (i < n + m ? ' ' : '\n');
  return 0;
}

5. 应用与技巧

FFT 和 NTT 是所有快速多项式操作的基础。

设多项式 $f$ DFT 得到 $\hat{f}$，也记为 $\mathrm{DFT}(f)$。可以发现，DFT 是线性变换，因此它具有 线性性，这是它最重要且最常用的一个性质：

• $c\mathrm{DFT}(f)+d\mathrm{DFT}(g)=\mathrm{DFT}(cf+dg)$。

5.1 常数优化

在分析变换次数的时候，一般不区分 DFT 和 IDFT。

5.1.1 三次变两次优化

计算两 实系数 多项式 $f,g$ 相乘。

根据 $(a+bi{)}^2=(a^2-b^2)+2abi$，将 $f+gi$ 平方后取出虚部再除以 $2$ 即可。

这样，三次 DFT 变成了两次 DFT，对常数有显著优化。代码。

这个优化被接下来稍复杂的技巧完全包含，它们的核心思想都是利用实系数的性质。

5.1.2 合并两次实系数 DFT

同时计算两 实系数 多项式 $f,g$ 的 DFT。

类似地，我们设 $u=f+ig$，$v=f-ig$。先求出 $u$ 的 DFT $\hat{u}$，考虑能否利用 $\hat{u}$ 直接求出 $\hat{v}$。

在给出做法之前，我们还要引入一些复数相关的概念：

• 定义：对于两个复数 $z_1$ 和 $z_2$，若它们实部相等，虚部互为相反数，则称 $z_1,z_2$ 互为 共轭复数。$z_2$ 是 $z_1$ 的共轭，$z_1$ 是 $z_2$ 的共轭。
• 表示：复数 $z$ 的共轭记作 $\bar{z}$ 或 $\mathrm{conj}(z)$。
• 运算性质：两个共轭复数的辐角互为相反数，即 $\arg z_1=-\arg z_2$。根据棣莫弗定理，可知 复数积的共轭，等于它们共轭的积。这样理解：求共轭相当于把整个复平面沿 $x$ 轴翻转，求积再翻转等价于翻转再求积。
• 易知复数和的共轭等于它们共轭的和，且一个复数的共轭的共轭等于它本身。

首先，$v({\omega }_n^k)=\sum _{i=0}^{n-1}{v}_i({\omega }_n^k{)}^i$。为了让它和 $u$ 产生联系，因为 $f,g$ 为实系数多项式，所以 $u$ 和 $v$ 的各项系数互为共轭，我们对整个结果求两次共轭，并将一次共轭根据共轭的性质下放至 $v_i$ 和 ${\omega }_n^k$。

$$\begin{array}{rl}{\hat{v}}_k& =v({\omega }_n^k)\\ & =\sum _{i=0}^{n-1}{v}_i({\omega }_n^k{)}^i\\ & =\mathrm{conj}\left(\mathrm{conj}(\sum _{i=0}^{n-1}{v}_i({\omega }_n^k{)}^i)\right)\\ & =\mathrm{conj}(\sum _{i=0}^{n-1}\mathrm{conj}(v_i({\omega }_n^k{)}^i))\\ & =\mathrm{conj}(\sum _{i=0}^{n-1}\mathrm{conj}(v_i)\mathrm{conj}({\omega }_n^k{)}^i)\\ & =\mathrm{conj}(\sum _{i=0}^{n-1}\mathrm{conj}(v_i)\mathrm{conj}({\omega }_n^k{)}^i)\\ & =\mathrm{conj}(\sum _{i=0}^{n-1}{u}_i({\omega }_n^{n-k}{)}^i)\\ & =\{\begin{array}{ll}\mathrm{conj}({\hat{u}}^{n-k})& (1\le k<n)\\ \mathrm{conj}({\hat{u}}_0)& (k=0)\end{array}\end{array}$$

求得 $\hat{v}$ 之后，因为 $\hat{u}=\hat{f}+i\hat{g}$，$\hat{v}=\hat{f}-i\hat{g}$，所以 $\hat{f}=\frac{\hat{u}+\hat{v}}{2}$，$\hat{g}=\frac{\hat{u}-\hat{v}}{2i}=\frac{(\hat{v}-\hat{u})i}{2}$。

5.1.3 合并两次实系数 IDFT

这里实系数指 IDFT 后的各项系数为实数。如果是 IDFT 前的各项系数为实数，直接使用上一小节的技巧即可。

设需要 IDFT 的两个多项式为 $\hat{f}$ 和 $\hat{g}$。计算 $\mathrm{IDFT}(\hat{f})$，各项系数均为实数，虚部信息被浪费了。考虑计算 $\mathrm{IDFT}(\hat{f}+i\hat{g})$，则各项系数的实部即 $f$ 的系数，虚部即 $g$ 的系数。

5.2 任意模数卷积

给定两多项式 $f,g$，在系数模 $p$ 意义下求出它们的卷积 $h=f\ast g$。

若 $p$ 不是 NTT 模数，我们不能朴素地使用 FFT 求解该问题，因为 $h$ 的系数可达 $n{p}^2$。取 $n={10}^6$，$p={10}^9$，则 $n{p}^2={10}^{24}$，long double 也无法接受。

接下来介绍两种常见的解决该问题的方法。它们也被称为 MTT（非正式），得名于毛啸 2016 年的集训队论文。

5.2.1 拆系数 FFT

设 $B=\sqrt{p}$，将所有系数表示为 $kB+r(0\le r<B)$ 的形式，得到四个多项式 $f=f_{0}B+f_1$，$g=g_{0}B+g_1$，计算它们两两相乘的结果，则 $fg=(f_0{g}_0)B^2+(f_0{g}_1+f_1{g}_0)B+f_1{g}_1$。

显然有一个四次 DFT 和三次 IDFT 的朴素做法：对 $f_0,f_1,g_0,g_1$ 进行 DFT，${\hat{f}}_0\cdot {\hat{g}}_0,{\hat{f}}_0\cdot {\hat{g}}_1+{\hat{f}}_1\cdot {\hat{g}}_0,{\hat{f}}_1\cdot {\hat{g}}_1$ 进行 IDFT。

使用 6.1.2 和 6.1.3 的技巧，可以做到两次 DFT 和两次 IDFT。系数值域 $n{B}^2=np={10}^{15}$，勉强可以接受。

模板题 P4245 任意模数多项式乘法 的 代码。注意用 long double，double 会被卡精度，或者换一种精度较高的 FFT 写法。

5.2.2 三模数 NTT

前置知识：（扩展）中国剩余定理。

选取三个常用 NTT 模数，分别算出 $f\ast g$ 在这些模数下的结果，再使用中国剩余定理合并即可。

我们选择 $p_1=998244353$，$p_2=1004535809$ 和 $p_3=469762049$，设结果分别为 $h_1,h_2,h_3$。

如果使用 CRT，则需要 __int128，因为 $h$ 的真实值为

$$(h_1{p}_2{p}_3\times \mathrm{inv}(p_2{p}_3,p_1)+h_2{p}_1{p}_3\times \mathrm{inv}(p_1{p}_3,p_2)+h_3{p}_1{p}_2\times \mathrm{inv}(p_1{p}_2,p_3))mod(p_1{p}_2{p}_3)$$

使用 exCET 就不需要 __int128：

• 合并 $h\equiv h_1(\mathrm{mod}{p}_1)$ 和 $h\equiv h_2(\mathrm{mod}{p}_2)$。设 $h=h_1+y{p}_1$，则 $h_1+y{p}_1\equiv h_2(\mathrm{mod}{p}_2)$，解得 $y\in [0,p_2)$ 之后得到 $h\equiv h_1+y{p}_1(\mathrm{mod}{p}_1{p}_2)$，记作 $h\equiv x(\mathrm{mod}{p}_1{p}_2)$。
• 合并 $h\equiv x(\mathrm{mod}{p}_1{p}_2)$ 和 $h\equiv h_3(\mathrm{mod}{p}_3)$。设 $h=x+y{p}_1{p}_2$，类似解得 $y\in [0,p_3)$，得到 $h\equiv x+y{p}_1{p}_2(\mathrm{mod}{p}_1{p}_2{p}_3)$。因为 $x+y{p}_1{p}_2<p_1{p}_2{p}_3$，所以它就是真实值，可以直接取模。

效率比拆系数 FFT 低不少，因为进行了九次 DFT。代码。

5.3 分治 NTT

前置知识：CDQ 分治。

分治 NTT 在信息竞赛界应用广泛，这归功于他所解决的问题：形如 $f_i=\sum _{j=0}^{i-1}{f}_j{g}_{i-j}$ 的 半在线卷积。

5.3.1 算法介绍

形如给定 $g_1\sim g_{n-1}$ 和 $f_0$，求 $f_1\sim f_{n-1}$ 满足 $f_i=\sum _{j=0}^{i-1}{f}_j{g}_{i-j}$ 的问题被称为 半在线卷积。因为 $f_i$ 很大，一般会对 NTT 模数 $p$ 取模。

我们不能通过简单的 NTT 解决这个问题，因为 $f$ 的每一项均和之前项有关，这要求我们在尝试计算 $f_i$ 时必须已经求出 $f_0\sim f_{i-1}$，而 NTT 做不到这一点。或者说，它们解决的问题形式不同。

这是一个在线的问题，考虑使用 CDQ 分治 转化为静态问题。

• 设当前分治区间为 $[l,r]$，其中 $f_0\sim f_{l-1}$ 对 $f_l\sim f_r$ 的贡献已经计算完毕。
• 若 $l=r$，说明已经求出 $f_l$，返回。
• 否则，令 $m=\lfloor \frac{l+r}{2}\rfloor$，先向左递归 $[l,m]$ 求解 $f_l\sim f_m$。
• 为了求解 $f_{m+1}\sim f_r$，我们需要计算 $f_l\sim f_m$ 对它们的贡献：$f_i=\sum _{j=l}^m{f}_j{g}_{i-j}$。因为 $f_l\sim f_m$ 已知，所以总的贡献相当于两个已知的多项式相乘的结果。具体地，令 $f_j^{\prime }=f_{j+l}(0\le j\le m-l)$，计算 $h=f^{\prime }\ast g$，则 $f_i$ 受到 $h_{i-l}$ 的贡献。
• 向右递归 $[m+1,r]$ 求解 $f_{m+1}\sim f_r$。

因为 $f,g$ 的长度均不超过区间长度，所以时间复杂度 $\mathcal{O}(n{\log }^{2}n)$。

模板题 P4721 分治 FFT 代码：

#include <bits/stdc++.h>
using namespace std;

using ll = long long;
using ull = unsigned long long;

constexpr int N = 1 << 17;
constexpr int mod = 998244353;
void add(int &x, int y) {x += y, x >= mod && (x -= mod);}
int ksm(int a, int b) {
  int s = 1;
  while(b) {
    if(b & 1) s = 1ll * s * a % mod;
    a = 1ll * a * a % mod, b >>= 1;
  }
  return s;
}

int n, f[N], g[N];
void solve(int l, int r) {
  if(l == r) return;
  int m = l + r >> 1, L = 1;
  solve(l, m);
  while(L < r - l + 1) L <<= 1;
  static int a[N], b[N], c[N];
  memset(a, 0, L << 2);
  memcpy(a, f + l, m - l + 1 << 2);
  memcpy(b, g, L << 2);
  NTT(L, a, 1), NTT(L, b, 1);
  for(int i = 0; i < L; i++) c[i] = 1ll * a[i] * b[i] % mod;
  NTT(L, c, 0);
  for(int i = m + 1; i <= r; i++) add(f[i], c[i - l]);
  solve(m + 1, r);
}

int main() {
  cin >> n;
  for(int i = 1; i < n; i++) scanf("%d", &g[i]);
  f[0] = 1, solve(0, n - 1);
  for(int i = 0; i < n; i++) printf("%d%c", f[i], i + 1 < n ? ' ' : '\n');
  return 0;
}

5.3.2 阶梯网格路径计数

阶梯网格路径计数是一类可以使用分治 NTT 解决的经典问题。

给定 $n+1$ 列 $m+1$ 行的网格图，左下角和右上角分别记为 $(0,0)$ 和 $(n,m)$。从 $(0,0)$ 出发，每次只能向右或向上走，求走到 $(n,m)$ 的方案数。让我们回忆：方案数为 $(\genfrac{}{}{0ex}{}{n+m}{n})$。

当然，问题没有这么简单。我们限制第 $i$ 列不能走到行数大于 $c_i$ 的点，其中 $c_i$ 单调不降，且 $c_n=m$。换言之，求出在一个阶梯网格从左下角走到右上角的方案数。

显然有 $\mathcal{O}(nm)$ 的动态规划：设 $f_{i,j}$ 表示走到 $(i,j)$ 的方案数，则 $f_{i,j}$ 转移至 $f_{i+1,j}$ 和 $f_{i,j+1}$。

考虑一段 $c_i$ 相同的极长区间 $[l,r](l<r)$ 从 $f_{l,j_1}$ 转移到 $f_{r,j_2}$ 的贡献系数：为防止重复计数，从 $(l,j_1)$ 出发的第一步应当向右，因此有系数 $(\genfrac{}{}{0ex}{}{(r-l-1)+(j_2-j_1)}{r-l-1})$。设 $g_i$ 表示 $(\genfrac{}{}{0ex}{}{r-l-1+i}{i})$，则 $f_{l,j_1}{g}_i\to f_{r,j_1+i}$，为卷积形式。

在不受 $c_i$ 影响的时候，我们确实可以这样做。但是对于每个 $i$，它内层的所有 $j$ 之间也有转移，这让我们不方便借助分治和卷积加快整个过程。考虑进行一些等价转换便于分层。

稍作思考，设 $f_{k,j}$ 表示走到 $(i,j)$ 其中 $i+j=k$ 的方案数，则 $f_{k,j}$ 转移至 $f_{k+1,j}$ 和 $f_{k+1,j+1}$。

进一步地，为了避免在 $k>n$ 时受到 $j\ge k-n$ 的限制，考虑从 $(n,0)$ 开始沿左上方向把整个阶梯网格劈成两半，分别计算从 $(0,0)$ 和 $(n,m)$ 到斜对角线上的点（$i+j=n$）的方案数，而这两个问题是对称的。

综合上述分析，我们将问题转化为：从 $(0,0)$ 出发，每次可以向右或右上走一步，求走到最右侧一列每个点的方案数。其中，在第 $i$ 列不能走到行数大于 $c_i$ 的点。这个 $c_i$ 可通过原问题的 $c_i$ 进行简单转化得到，且容易证明其仍不降且满足很好的性质：$c_{i+1}-c_i\le 1$。

因此，从 $f_l$ 推到 $f_r$ 的时候，我们会发现对于 $j\le c_r-(r-l)$，设不等号右侧的数为 $x$，$f_{l,j}$ 对 $f_r$ 的贡献不会受到 $c$ 的影响：因为 $c_{i+1}-c_i\le 1$，所以从 $(l,j)$ 开始，就算每一步都往右上走，也不会突破限制。这样，$f_{l,0}\sim f_{l,x}$ 对 $f_r$ 的贡献可直接卷积计算。对于 $f_{l,x+1}\sim f_{l,c_l}$，分治下放至左右子区间 $(l,m]$ 和 $(m,r]$ 进行递归。

有了大致思路，设计算法就很简单了：设 $\mathrm{solve}(l,r,\mathrm{\Delta },F)$ 表示当前区间为 $(l,r]$，传入多项式 $F$ 的第 $i$ 项表示 $f_{l,i+\mathrm{\Delta }}$，返回多项式 $G$ 的第 $i$ 项表示 $f_{r,i+\mathrm{\Delta }}$。也可视作暂时将 $c_l\sim c_r$ 全部减去 $\mathrm{\Delta }$，传入 $f_l=F$，返回 $G=f_r$，接下来就使用这种视角。

对于 $j\le c_r-(r-l)$，设不等号右侧的数为 $x$，$F_0\sim F_x$ 和 $H_0\sim H_{r-l}$ 卷积求出 $F_0\sim F_x$ 转移得到的 $G_0\sim G_{c_r}$，其中 $H_i$ 即转移系数 $(\genfrac{}{}{0ex}{}{r-l}{i})$。

对于 $j>x$，分治下放：设 $F^{\prime }=F_{x+1}\sim F_{c_l}$，$mid=\lfloor \frac{l+r}{2}\rfloor$，则先递归左半部分 $F^{\prime }\leftarrow \mathrm{solve}(l,mid,\mathrm{\Delta }+x+1,F^{\prime })$，再递归右半部分 $F^{\prime }\leftarrow \mathrm{solve}(mid,r,\mathrm{\Delta }+x+1,F^{\prime })$，则此时 $F^{\prime }$ 表示从 $F_{x+1}\sim F_{c_l}$ 转移得到的 $G_{x+1}\sim G_{c_r}$。

两部分相加即得 $G$，返回即可。

边界 $r-l=1$ 的处理是平凡的。

两次下传 $F^{\prime }$ 时，$F^{\prime }$ 的长度显然不超过 $c_r-x=r-l$，因此在处理区间 $(l,r]$ 时涉及到的所有多项式的长度均不超过其父区间长度。设 $n,m$ 同级，则时间复杂度为 $\mathcal{O}(n{\log }^{2}n)$。

接下来对问题进行一些扩展：

• 如果没有 $c_{i+1}-c_i\le 1$ 的性质，但 $c_i$ 单调不降，还能做吗？答案是可以，只需将 $x$ 的定义改为 $c_l-(r-l)$，此时 $(l,r]$ 涉及到的多项式长度不超过父区间的长度加上端点处 $c$ 的差值，可以接受。

5.4 例题

CF1770G Koxia and Bracket

相比于 F，这道题就略显套路了。

将左括号视为 $1$，右括号视为 $-1$，找到最后一个使得前缀和小于 $0$ 的位置 $lst$，则 $s_1\sim s_{lst}$ 的每个左括号都有用，必然不会删去。同理，$s_{lst+1}\sim s_n$ 的每个右括号也都有用。

对于 $s_1\sim s_{lst}$ 的每个右括号 $s_i$，如果它对应的前缀和为 $c_i$，则为保证前缀和非负，在 $s_1\sim s_i$ 中至少需要删掉 $-c_i$ 个右括号。我们只关心 $-c_i$ 的前缀最大值，因为若这些位置满足条件，则所有位置满足条件，而每个前缀最大值一定会比先前的前缀最大值大 $1$，所以我们将问题转化为：给定长为 $m$ 的序列，其中有 $c$ 个位置被打上了标记，求出选择 $c$ 个位置的方案数，使得对于每个前缀，选择的位置数不小于被打上标记的位置数。

设 $f_{i,j}$ 表示考虑到第 $i$ 个位置，当前选择位置数减去打标记位置数的数量为 $j$。对于没有被打标记的位置 $p$，$f_{p-1,j}$ 转移到 $f_{p,j/j+1}$，否则转移到 $f_{p,j-1/j}$。

非常明显的阶梯格路计数问题，稍有变形，不过解法大差不差，核心思想是一致的：考虑 $f_l$ 转移到 $f_r$，用卷积描述一部分转移，剩下来 $\mathcal{O}(r-l)$ 个位置分别递归两侧处理。

对于本题，就是 $f_{l,j}(j\ge cnt)$ 对 $f_r$ 的贡献用卷积描述，其中 $cnt$ 表示 $l+1\sim r$ 的被打上标记的位置。剩余部分递归，长度只有 $cnt$，而且因为截取的是低位，我们甚至不需要记录偏移量 $\mathrm{\Delta }$。

$r-l=1$ 的边界情况是平凡的。

$s_{lst+1}\sim s_n$ 的左括号类似处理即可。

每个分治区间涉及到的多项式长度不超过其父区间长度，因此时间复杂度为 $\mathcal{O}(n{\log }^{2}n)$。代码。

参考资料

第一章：

• 复数 - OI Wiki。
• 复数 - 百度百科。
• [Video] 虚数的来源 - Veritasium。
• [Video] 【官方双语】微分方程概论 - 第五章：在 3.14 分钟内理解 $e^{i\pi }$ - 3Blue1Brown。

第二章：

• 多项式与生成函数简介 - OI Wiki。
• 多项式 - 百度百科。
• 代数基本定理 - OI Wiki。

第四章：

• 快速傅里叶变换 - OI Wiki。
• FFT 理性瞎扯 - xtx1092515503。
• Vandermonde matrix - Wikipedia。
• Discrete Fourier transform - Wikipedia。
• 浅谈范德蒙德 (Vandermonde) 方阵的逆矩阵的求法以及快速傅里叶变换 (FFT) 中 IDFT 的原理 - Deadecho。
• 范德蒙德矩阵的逆矩阵怎么计算? - FFjet 的回答 - 知乎。
• Discrete Fourier transform over a ring - Wikipedia。
• NTT 与多项式全家桶 - command_block。
• [Video] 这个算法改变了世界 - Veritasium。
• [Video] The Fast Fourier Transform (FFT): Most Ingenious Algorithm Ever? - Reducible & B 站搬运。

第五章：

• P3803 题解 - NaCly_Fish。
• 2016 集训队论文《再探快速傅里叶变换》—— 毛啸。