基础知识：FFT - 简单入门

七海

算法数学

发布于：2021年3月29日

次浏览

一句话简介：FFT 即快速傅里叶变换，是一种可以在 $n\log n$ 的时间内完成多项式乘法的算法——的一部分。

前置知识

在了解 FFT 之前，你可能需要先知道的一些东西：

多项式

对于数域 $\mathbb{F}$，若有 $\forall i\in\{1,2,3,\cdots,n \}$，则： \[ f(x)=a_0+a_1x+a_2x^2+\cdots+a_nx=\sum_{i=0}^n a_ix^i \ (a_n \neq 0) \] 为数域 $\mathbb{F}$ 上的一个多项式；该多项式的次数即为其中最高次项的次数为 $n$，记作 $\deg f(x) = n$；

表示方法

通常有两种表示方法

系数表示法

即上文定义中的表示方法，是一个向量；可以映射为系数向量 $\vec{a}$： \[ f(x) = \sum_{i = 0}^n a_ix^i ⇔ \vec{a} = (a_0, a_1, \cdots ,a_n) \ (a_n \neq 0) \] 我们将向量 $\vec{a}$ 成为多项式 $f(x)$ 的系数表示。

点值表示法

设 $\deg f(x) = n$；众所周知，确定一个 $n$ 次函数的解析式，需要 $n + 1$ 个点的坐标；同理，这样也可以确定一个有 $n + 1$ 个系数的 $n$ 次多项式。

因此，对于 $\forall x_i \in \mathbb{F}, i \in [0, n]$，有 $y_i = f(x_i)$；如果对于 $\forall i \neq j$ 总是满足 $x_i \neq x_j$，那么我们可以在数域 $\mathbb{F}$ 中绘制 $n+1$ 个点，获得 $n+1$ 个不同的方程，从而确定这个多项式；

那么，这 $n + 1$ 个不同的点构成的点集，就是多项式的点值表示。

如果将这 $n + 1$ 个不同的点代入表达式中进行运算，可以得到下面的矩阵等式： \[ \begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n \\ 1 & x_1 & x_1^2 & \cdots & x_1^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_n \end{bmatrix} \] 我们将最左侧的 $n + 1$ 阶矩阵称为范德蒙矩阵。

那么如何转换多项式的两种表达方式呢？

系数表达式 => 点值表达式
这个还是非常的简单的；只需要选取 $n + 1$ 个不同的 $x_i$，分别带入解析式中求出 $y_i$ 即可。我们假定使用快速幂计算乘方，那么这样做的复杂度就是 $O(n^2\log n)$ 的。
点值表达式 => 系数表达式
这个过程又被称为多项式插值；朴素地，我们可以将多项式看作一个 $n$ 次函数，将点值表达中的 $n + 1$ 个点代入其中得到等数目的方程，然后求解这个方程组得到每个系数，组成系数表达式。
形式化地说，我们可以利用范德蒙矩阵进行高斯消元，就可以在 $O(n^3)$ 的时间内求解；但是这样实在是太慢了，时间上无法接受；因此我们也可以用拉格朗日插值法求解： \[ f(x) = \sum_{i = 1}^n y_i \prod_{i = j} \frac{x-x_j}{x_i-x_j} \] 这样，我们可以在 $O(n^2)$ 的时间复杂度内完成多项式插值。

虽然但是，我们可以看到想要在这两种形式中转化，一般会消耗大量的时间。

操作

令参与操作的两个多项式为 $f(x)$ 和 $g(x)$；若 $\deg f(x) \neq \deg g(x)$，我们可以为其中度数较小的一方增加系数为 $0$ 的高次项，从而使得两个多项式齐次；因此，下面的讨论都仅针对于两个多项式齐次的情况。

此外，如果两个多项式均仅提供点值表示，那么要求提供的点值的横坐标一一对应且相等。

现在假定 $\deg f(x) = \deg g(x) = n$，系数表示分别为 $\vec{a}、\vec{b}$；

加减法

对于系数表示，若 $h(x) = f(x) ± g(x)$，则有： \[ h(x) = f(x) ± g(x) = \sum_{i = 0}^n (a_i ± b_i)x^i \] 且 $h(x)$ 的系数表示记为 $\vec{c}$，那么 $c_i = a_i ± b_i$；

对于点值表示，只需要将对应的纵坐标相加/相减即可。

乘法（卷积）

同上，对于系数表示，可以得到多项式乘法的定义如下： \[ h(x) = (f \cdot g)(x) = f(x) \times g(x) = \sum_{i = 0}^n \sum_{j = 0}^n a_ib_jx^{i + j} \] 显然，两个 $n$ 次的多项式相乘会得到一个 $2n$ 次的多项式；若 $h(x)$ 的系数表示记为 $\vec{c}$，那么也有： \[ c_i= \sum_{j+k=i}a_jb_k = \sum_{j = i}^{n} a_jb_{i - j} \] 那么系数向量 $\vec{c}$ 为稀疏向量 $\vec{a}、\vec{b}$ 的卷积，记作：$\vec{c} = \vec{a} ∗ \vec{b}$；

对于点值表示，依然只需要将对应的纵坐标相乘即可；但是因为新得到的多项式次数更高，所以每个因子多项式都需要提供 $2n$ 个点参与运算；

综上所述，我们可以看到：因为多项式运算本质上还是多项式的值进行运算，所以对于点值表示法而言，这些运算都是可以 $O(n)$ 完成的（即仅将对应点值的纵坐标进行运算即可）；唯一需要注意的就是多项式乘法，点值表示法需要提供更多组点值进行运算。

复数

~~高等数学和复变函数的完全败北~~

首先定义虚数单位 $i^2 = -1$；则所有形如： \[ z = a + bi, \ \ a,b \in \mathbb{R} \] 的数字为复数，它们构成的集合成为复数集，记为 $\mathbb{C}$；上面的表达式中，$a$ 被称为实部，$b$ 被称为虚部。

复平面

复平面是一个笛卡尔平面，有两条坐标轴，纵轴为虚轴，横轴为实轴，两轴相互垂直；对于每一个复数 $z = a + bi$，它都可以在复平面上被表示为一个从原点指向点 $(a, b)$ 的向量；显然，复数和复平面上从原点出发的向量之间的映射关系是双射。

和一般向量一样，以实轴正方向为始，$\vec{z}$ 为终的角 $\theta$ 称为复数 $z$ 的幅角。

基本操作

模：即复数在复平面上所对应的向量的模，记作 $|z|$；$|z| = \sqrt {a^2 + b^2}$
共轭复数：复数在复平面上所对应的向量关于实轴对称后对应的复数，记为 $\bar{z}$；
- $z$ 和 $\bar{z}$ 的实部相同，虚部互为相反数
- 两者的幅角满足：$\theta_{z} + \theta_{\bar{z}} = \pi$
- 两者的模满足：$|z| = |\bar{z}|$

运算

令参与运算的两个复数为 $z_1 = a_1 + b_1i$ 和 $z_2 = a_2 + b_2i$；

加减法：$z = z_1 ± z_2$
- $z = (a_1 ± a_2) + (b_1 ± b_2)i$
- 对于向量：按照平行四边形定则
乘法：$z = z_1 \times z_2$
- $z = (a_1 + b_1i) (a_2 + b_2i) = $ $(a_1a_2 - b_1b_2) + (a_1b_2 + a_2b_1)i$
- 对于向量：它们的幅角满足 $\theta_z = \theta_{z_1} + \theta_{z_2}$
- 对于向量：它们的模满足 $|z| = |z_1| \times|z_2|$
- 若 $z_2 = \bar{z_1}$，那么 $z = a_1^2 + b_1^2$，显然是一个实数
除法：$z = \frac{z_1}{z_2}$
- 通分 $\bar{z_2}$，可以得到：$z = \frac{z_1\bar{z_2}}{a_2^2 + b_2^2}$
- 显然，分母是实数，所以可以直接除进分子的实部和虚部
指数幂：首先，有欧拉公式 $e^{i\theta} = \cos\theta + i\sin\theta$
- 在数学分析和复变函数中，欧拉公式是一个把复指数函数与三角函数联系起来的一个公式： \[ e^{i\theta} = \cos\theta + i\sin\theta \] 这个恒等式就是欧拉公式，它是数学里最令人着迷的一个公式——它将数学里最重要的几个数字联系到了一起：
  - 两个超越数：自然对数的底 $e$，圆周率 $\pi$；
  - 两个单位：虚数单位 $i$ 和自然数的单位 $1$ ；
  - 以及被称为人类伟大发现之一的 $0$。
  它将指数函数的定义域扩大到复数，建立了三角函数和指数函数的关系，被誉为“数学中的天桥”。因此数学家们评价它是“上帝创造的公式”。
- 任何复数都可以表示为欧拉公式的形式，即： \[ z = |z|\cos\theta_z + i|z|\sin\theta_z = |z|e^{i\theta_z} \]
- 特殊值：当 $\theta = \pi$ 时，有 $e^{i\pi} = -1$
- 因此，复数的幂只需要将它转化为欧拉公式的形式，然后取幂即可

和多项式一样，复数的运算也需要根据不同的运算在不同的表示形式中反复切换。

单位根

在 $\mathbb{C}$ 中，满足 $x^n = 1$ 的 $x$ 被称为 $n$ 次单位根；根据代数基本定理可知，$n$ 次单位根共有 $n$ 个，它们分别是： \[ x_k = e^{i\frac{2k\pi}{n}}, \ k \in [0, n) \] 显然，$x_k^n = 1$；可以用上一节讲到的欧拉公式来证明。

性质

因为 $\sin\theta + \cos\theta = 1$，所有的 $n$ 次单位根的模均为 $1$.

在复平面上，$n$ 个 $n$ 次单位根平分了单位圆，且和单位圆的 $n$-等分线重叠。

本原单位根

定义：$0$ 到 $(n - 1)$ 次幂的值能生成全部 $n$ 个 $n$ 次单位根的 $n$ 次单位根称为为 $n$ 次本原单位根。

显然，对于任何 $n > 1$，$x_1 = e^{i\frac{2\pi}{n}}$ 是一个 $n$ 次本原单位根。

特别地，我们记 $\omega_n = e^{i \frac{2\pi}{n}} = \cos \frac{2\pi}{n} + i \sin \frac{2\pi}{n}$ 为 $n$ 次本原单位根；虽然它未必唯一，但是接下来的所有讨论中的“$n$ 次本原单位根”都特指 $\omega_n$。

$n$ 次本原单位根 $\omega_n$ 有一些特殊的运算性质：

折半定理：$(\omega_n^k)^2 = \omega_{\frac{n}2}^k$
令 $n$ 为偶数：$\omega_n^{\frac{n}2 + k} = -\omega_n^k$
相消定理：$\omega_{dn}^{dk} = \omega_n^k$
求和定理：可以描述为如下公式：
\[ \sum_{i = 0}^{n - 1} (\omega_n^k)^i = \begin{cases} 0, \ \ & k \neq mn \ , \ m \in \mathbb{Z} \\ n, \ \ & k = mn \ , \ m \in \mathbb{Z} \end{cases} \] 在 $k \neq mn$ 的场合下，由于公比不为 $1$，我们可以用等比数列的求和公式推导： \[ \sum_{i = 0}^{n - 1} (\omega_n^k)^i = \frac{1 - (\omega_n^k)^n}{1 - \omega_n^k} = \frac{1 - 1}{1 - \omega_n^k} = 0 \] 否则，公比为 $\omega_n^k = 1$，上述推导并不成立；但由于 $\omega_n^n = \omega_n^0 = 1$，所以上计算结果为 $n$；

显然，第一个式子是成立的：每份角度扩大一倍，那么可以等分的份数就减为一半，这很合理；

显然，$\omega_n^{\frac{n}2}$ 的幅角是 $\pi$；由和差角公式可以知道幅角加上 $\pi$ 意味着正弦和余弦值都变为原来的相反数，所以第二个式子也非常的正确；证毕（

相消定理也非常的好理解：两边的表达式转化为欧拉公式的形式是一样的。

~~不会有人真的连高中数学学的三角函数相关变换的忘记了吧？不会这人就是我吧？~~

到现在为止，我们已经把必要的前置知识~~数学基础~~给介绍完了。

分析过程

根据上面的描述，我们已经知道了对于多项式的系数表示，计算其乘法需要 $O(n^2)$ 的时间；而对于点值表示法，这个运算只需要 $O(n)$ 的时间；但是如果我们想要通过两个多项式的系数表示得到它们乘积的系数表示，我们需要遵循以下的步骤：

求值：将待乘的多项式转化为点值表示，需要 $O(n^2)$
运算：点值表示多项式乘法，需要 $O(n)$
插值：将得到的结果转化为系数表示，朴素算法需要 $O(n^3)$，不低于 $O(n^2)$

可以看出，如果仅是将多项式转化为点值表示，并不能做到算法复杂度的降低——但是它为我们提供了思路：如果我们能将求值和插值步骤优化到低于 $O(n^2)$，就可以降低多项式乘法的整体复杂度。

离散傅里叶变换

特别地，现在对上面的步骤进行一些定义：

离散傅里叶变换（DFT）：指求出一个 $n - 1$ 次多项式在每个 $n$ 次单位根下的点值的过程
离散傅里叶逆变换（IDFT）：将上面求出的那些点值重新插值得到系数表示的过程

对于多项式 $f(x)$，若 $\deg f(x) = n - 1$，那么对它进行的 DFT 被称为长度为 $n$ 的离散傅里叶变换，记作 $\text{DFT}_n(f)$。

显然，这个操作就是朴素做法，具有 $O(n^2)$ 的时间复杂度；现在我们考虑使用 FFT 优化它。

快速傅里叶变换

现在，我们假设 $n$ 是偶数，且有 $m = \frac{n}2$；有多项式 $f(x) = \sum_{i = 0}^{n - 1} a_ix^i$，显然 $\deg f(x) = n - 1$；

分治

考虑分治：我们对 $f(x)$ 的系数角标按照奇偶性进行分类，可以得到： \[ \begin{align} f(x) &= \sum_{i = 0}^{m - 1} a_{2i}x^{2i} + \sum_{i = 0}^{m - 1} a_{2i+1}x^{2i+1} \\ &= \sum_{i = 0}^{m - 1} a_{2i}x^{2i} + x(\sum_{i = 0}^{m - 1} a_{2i+1}x^{2i}) \\ &= \sum_{i = 0}^{m - 1} a_{2i}(x^2)^i + x(\sum_{i = 0}^{m - 1} a_{2i+1}(x^2)^i) \end{align} \] 我们对于前半部分（偶数部分）记作 $f_0(x)$，后半部分（奇数部分）记为 $f_1(x)$： \[ \begin{align} f_0(x) &= \sum_{i = 0}^{m - 1} a_{2i}x^i \\ f_1(x) &= \sum_{i = 0}^{m - 1} a_{2i + 1}x^i \\ \end{align} \] 那么，$f(x)$ 可以被表示为如下形式： \[ f(x) = f_0(x^2) + x \cdot f_1(x^2) \] 如果知道了 $f_0(x)$ 和 $f_1(x)$ 在 $n$ 个不同的位置的值，那么我们只需要计算 $O(n)$ 次上式，就可以得到 $f(x)$ 的点值表示了；然后可以继续递归处理 $f_{0/1}$；

但是仅仅是这压根是不够的；根据主定理，因为 $f_{0/1}$ 各需要递归一次，所以复杂度还是 $O(n^2)$ 的。

使用单位根

因此，我们考虑代入单位根求 $\text{DFT}_n(f)$。根据单位根的性质，我们知道了 $\omega_n^{m + k} = -\omega_n^k$：这提示我们需要考虑小于 $m$ 次的点值和大于等于 $m$ 次的点值之间存在的关系；因此我们以此进行分类讨论：

对于 $k \in [0, m)$，我们可以利用单位根的折半定理进行简化，有： \[ \begin{align} f(\omega_n^k) &= f_0((\omega_n^k)^2) + \omega_n^k \cdot f_1((\omega_n^k)^2) \\ &= f_0(\omega_m^k) + \omega_n^k \cdot f_1(\omega_m^k) \end{align} \] 在考虑 $k + m \in [m, n)$，我们同样可以利用单位根的性质得到下面的化简： \[ \begin{align} f(\omega_n^{k + m}) &= f_0((\omega_n^{k + m})^2) + \omega_n^{k + m} \cdot f_1((\omega_n^{k + m})^2) \\ &= f_0((-\omega_n^k)^2) - \omega_n^k \cdot f_1((-\omega_n^k)^2) \\ &= f_0((\omega_n^k)^2) - \omega_n^k \cdot f_1((\omega_n^k)^2) \\ &= f_0(\omega_m^k) - \omega_n^k \cdot f_1(\omega_m^k) \end{align} \] 上面的两个式子又被称为蝴蝶操作，式子中的 $\omega_n^k$ 又被称为旋转因子；

至于叫“蝴蝶操作”这个名字的原因，现在我们可以考虑两层：一层是我们将要求的 $f(x)$，它的 DFT 长度为 $n$；另一层是我们已经求好的（或者将要递归去求的）$f_{0/1}(x)$，它们的 DFT 长度均为 $m$；对于 $f_{0/1}(\omega_m^k)$，它们能影响到的位置只有 $f(\omega_n^k)$ 和 $f(\omega_n^{k + m})$；所以可以画出下面的图：

图中实边贡献的权值为 +1，虚边的权值是 -1；结合推导的式子很容易理解，形状也许确实是蝴蝶（

综上所述，我们可以发现：$k \geq m$ 部分的值可以在求 $k < m$ 部分的值时一并求出；因此，只要知道了$f_0(x)$ 和 $f_1(x)$ 在小于 $m$ 次位置的值（显然，这是 $m = \frac{n}2$ 次），就可以在 $O(n)$ 的时间内求出 $f(x)$ 在整全部 $n$ 个位置的值；显然，和朴素做法一样，这个过程也是递归的。

分析算法的时间复杂度：$T(n) = T(\frac{n}2) + O(n)$，典型的 $O(n\log n)$；

综上所述，得到的快速计算 $\text{DFT}_n(f)$ 的方法，就被称为 FFT。

插值方法

在上一个部分，我们已经将朴素求值的过程用 FFT 取代，获得了更优的时间复杂度；但是朴素插值的过程仍然是 $O(n^2)$ 的；这个部分将说明朴素 IDFT 的过程，并尝试压缩它的时间复杂度。

离散傅里叶逆变换

问题
有 $\deg f(x) = n - 1$ 的多项式 $f(x)$ ；已知其点值表示 $\{(x_k, \ y_k) \ | \ k \in [0, n)\}$；其中： \[ \text{DFT:} \ \ y_k = \sum_{i = 0}^{n - 1} a_ix_i \ , \ x_i = \omega_n^{ik} \] 现在我们要从这个点集反向求出 $f(x)$ 的系数表示 $\vec{a}$。

至于这个问题的推导结果，~~因为我不会推导~~我们先给出它的结论： \[ \text{IDFT:} \ \ a_k = \frac1n \sum_{i = 0}^{n - 1}y_i\omega_n^{-ki} \] 然后我们进行反向推导，来证明它确实是 $\text{DFT}$ 式子的逆变换：

将 $\text{IDFT}$ 式子中的 $y_i$ 展开为 $\text{DFT}$ 式： \[ \begin{align} a_k = \frac1n \sum_{i = 0}^{n - 1}y_i\omega_n^{-ki} &= \frac1n \sum_{i = 0}^{n - 1}\sum_{j = 0}^{n - 1} a_j\omega_n^{ji}\omega_n^{-ki} \\ &= \frac1n \sum_{i = 0}^{n - 1}\sum_{j = 0}^{n - 1} a_j\omega_n^{i(j - k)} \\ &= \frac1n \sum_{j = 0}^{n - 1} a_j\sum_{i = 0}^{n - 1}(\omega_n^{j - k})^i \end{align} \]
现在我们使用单位根的求和定理来考虑等式右侧的 $\sum_{i = 0}^{n - 1}(\omega_n^{j - k})^i$：
- $j = k$ 时，$\omega_n^{i(j - k)} = \omega_n^0 = 1$，故： \[ \sum_{i = 0}^{n - 1}(\omega_n^{j - k})^i = 1 \times n = n \]
- $j \neq k$ 时，由求和定理可知： \[ \sum_{i = 0}^{n - 1}(\omega_n^{j - k})^i = 0 \] 既可以使用消去定理，也可以使用单位根的定义得到这个值为 0；
- 综上所述，仅当 $j = k$ 时，外侧的求和会对结果产生 $n$ 倍的贡献。
因此，我们可以继续推导第一步我们得到的式子： \[ \begin{align} a_k = \frac1n \sum_{i = 0}^{n - 1}y_i\omega_n^{-ki} &= \frac1n \sum_{j = 0}^{n - 1} a_j\sum_{i = 0}^{n - 1}(\omega_n^{j - k})^i \\ &= \frac1n a_k \times n \end{align} \] 显然，等式的两侧是相等的。

所以，我们证明了提供的结论确实是 $\text{DFA}$ 的逆变换。

如果愿意，也可以把式子写成矩阵的形式，然后对对应的矩阵求逆，也能得到这个结果。

分析这种做法的复杂度：对于每一个 $k \in [0, n)$ 我们都需要 $O(n)$ 的时间去求出等式右侧的求和；所以 $\text{IDFT}$ 公式的复杂度毫无疑问也是 $O(n^2)$ 的。

快速傅里叶逆变换

那么现在，我们就要考虑如何在更短的复杂度下求解 IDFT。显然 IDFT 公式也是一个多项式，我们定义： \[ a_k = F(x_k) = \frac1n \sum_{i = 0}^{n - 1}y_ix_k^i, \ \ x_k = \omega_n^{-k} \] 在回忆一下我们前面刚推到完的 FFT： \[ y_k = f(x_k) = \sum_{i = 0}^{n - 1}a_ix_k^i, \ \ x_k = \omega_n^k \] 我们就会发现它们出奇的相似！（好吧，我也是推到这一步才豁然开朗，，只能说我确实把复变函数的一些知识忘得精光，完全的不掌握了把）；$\omega_n^k$ 和 $\omega_n^{-k}$ 很显然是一对共轭复数——在复平面的角度来看它们只是幅角的旋转方向不同，实际上还是一一对应的；所以我们完全可以把 FFT 的方法拿到这里来，求出 $F(x)$ 在 $\omega_n$ 不同幂下的值之后，翻转结果就可以得到答案（$\omega_n^{-k} = \omega_n^{n-k}$）。

~~是不是突然有些明白了在蝴蝶操作中 $\omega_n^k$ 被称为旋转因子的原因了（~~

综上所述，FFT + 事后翻转 = $\text{IFFT}$；这就是快速傅里叶逆变换插值，时间复杂度显然是 $O(n\log n)$。

综上所述，我们使用 $\text{FFT}$ 和 $\text{IFFT}$ 取代求值和插值步骤朴素的 $\text{DFT}$ 和 $\text{IDFT}$ 算法，就可以把多项式的点值表达和系数表达之间的互相转化的复杂度从 $O(n^2)$ 降为 $O(n\log n)$，从而得到了一个总体复杂度为 $O(n\log n)$ 的多项式乘法的算法。

所以严格上来说，快速的多项式乘法算法并不是 FFT；FFT 只是其中非常重要的组成部分罢了。

算法实现

好了，你已经学会了 FFT 了！快上！（指上机写代码

朴素的递归实现

然后我们就写出了像上面说的那样的递归 FFT：

complex getW(uint n, uint k = 1)
{
    const auto angle = (k % n) * PI * 2 / n;
    return complex(cos(angle), sin(angle));
}

void subprocess(addr f, size_t n)
{
    if (n == 1) return;
    const auto m = n / 2;
    vector<complex> &sub = _alloc(n);
    auto f0 = sub.begin(), f1 = f0 + m;
    for (int i = 0; i < m; ++ i)
        f0[i] = f[i * 2],
        f1[i] = f[i * 2 + 1];
    subprocess(f0, m);
    subprocess(f1, m);
    auto unit = getW(n), w = getW(1, 0);
    for (int i = 0; i < m; ++ i)
        f[i] = f0[i] + w * f1[i],
        f[i + m] = f0[i] - w * f1[i],
        w *= unit;
    _free_mem_stack_frame();
}

上面的代码是最核心的递归过程；一些包装的函数在下面：

void start(vector<complex> &p, bool inv)
{
    subprocess(p.begin(), p.size());
    if (inv) reverse(++ p.begin(), p.end());
}

void fft(vector<complex> &p)
{start(p, false);}

void ifft(vector<complex> &p)
{
    start(p, true);
    const auto n = p.size();
    for (auto &i : p) _post_ifft(i, n);
}

如果是系数表达变点值表达，只需要开始递归即可；反之，在得到结果之后还需要进行后处理：首先将答案数组的 $[1, N - 1]$ 位置翻转（因为 $\omega_n^0 = \omega_n^n$，所以交换的是 $1$ 和 $n - 1$ 的位置），然后根据推导的表达式，我们还要对每个位置都除 $N$；需要注意的是这一步也会引入误差，需要 round。

然后我们再为它添加一些头和尾的代码，用来去做模板题检验正确性：【模板】多项式乘法（FFT）；增加的代码就像下面这样：

signed main()
{
    ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    std::cout.tie(nullptr);
    int n = scanner.nextInt(),
        m = scanner.nextInt();
    vector<FFT::complex> a, b;
    const auto siz = FFT::stretch(n + m);
    a.resize(siz), b.resize(siz);
    FFT::init_mem_pool(siz);
    for (int i = 0; i <= n; ++ i)
        a[i].real(scanner.nextInt());
    for (int i = 0; i <= m; ++ i)
        b[i].real(scanner.nextInt());
    FFT::fft(a), FFT::fft(b);
    for (int i = 0; i < siz; ++ i)
        a[i] *= b[i];
    FFT::ifft(a);
    for (int i = 0; i <= n + m; ++ i)
        cout << (int)a[i].real() << " \n"[n + m == i];
    return 0;
}

然后交上去就会发现…… A 了？和别人博客里说的不一样啊？

实际上是因为递归带来的较大的常数开销，很大一部分都来自于每层递归重复开辟内存空间用来做分治；因此，我写的时候就使用了一个内存池来处理递归中带来的内存分配和释放，就可以有效的降低它的常数。

当然，这种做法还是不够快的；所以我们还要对它继续优化。

优化迭代实现

因为递归是从上到下的过程，那么我们的迭代实现就是一个从下而上的一个过程。

我们把长度为 $8$ 的情况下，上述朴素的递归做法的调用树画出来，如下图：

Layer 1:		   0 1 2 3 4 5 6 7
				  /				  \
Layer 2:	   0 2 4 6		   1 3 5 7
			  /		  \		  /		  \
Layer 3:	 0 4	 2 6	 1 5	 3 7
			 / \ 	 / \	 / \	 / \
Layer 4:	0	4	2	6	1	5	3	7

当然，只看这个我们什么也看不出来；现在我们将焦点集中在最后一层上，并将它们的二进制和二进制反转以及二进制反转后得到的数字写出来，如下表：

Layer 4:	0	4	2	6	1	5	3	7
Binary:		000 100 010 110 001 101 011 111
Reverse:	000 001 010 011 100 101 110 111
Decimal:	0	1	2	3	4	5	6	7

可以看出最后一层经过一番操作之后就变成了最开始的数组。因此，我们只需要在最开始的数组上进行这样的一番操作，就可以得到最后一层的这个新的数组，同时第一次合并也只是对相邻的位置进行合并。那么这样上面那几层就可以这样一直合并上去得到最终的数组了。

那么我们要怎么样快速的进行这一波神秘操作呢？现在我们考虑反转后的二进制表示（$1, 2, \dots,2^k-1$），我们可以按照这样的顺序构造出这个序列：

1
2
3

round 1:												(1)001	(0)000	
round 2:								(3)011	(2)010	(1)001	(0)000	
round 3:	...	(6)110	(5)101	(4)100	(3)011	(2)010	(1)001	(0)000

第一轮，所有的位都是 0；第一个位置是 0 已经满足了要求；我们选择最低位，放置 1，得到了序列 0, 1；第二轮，已经有满足要求的序列长度 2，我们选择次低位，将这一位置 1，那么又可以填两个数字，这两个数字的相对位置就和最开始满足要求的序列那样…… 最终，重复 n 次就可以构造出 $0,\dots,2^k-1$ 的序列。

那么我们现在构造的是二进制反向的序列，那么只需要从高位开始，重复这个过程就可以了：

vector<uint> &reverse(uint length)
{
    static vector<uint> rev(length);
    uint msk = length >> 1u, i = 0;
    rev[i ++] = 0, rev[i ++] = msk;
    for (uint w = 2; w < length; w <<= 1)
    {
        msk >>= 1u;
        for (uint j = 0; j < w; ++ j)
            rev[i ++] = rev[j] | msk;
    }
    return rev;
}

得到的数组是关于原数组的一个排序。我们只需要利用这个排列进行位逆序置换——因为每一个 i 和 rev[i] 都是一一对应的，它们本是有序的，那么只要交换所有逆序的对，就可以得到适用迭代的数组。

计算在 $n$ 次单位根的各幂次点值的时候，令 $m=\frac{n}2$，那么 $f_0$ 和 $f_1$ 已经由下层（迭代意义上的上层）计算并储存在 $A[k]$ 和 $A[k+m]$ 处，只需取出来计算后更新对应的值就完成了一层的迭代。到此，我们已经完成了对迭代优化的实现的过程描述，核心代码实现如下：

void start(vector<complex> &p, bool inv)
    {
        const auto siz = p.size();
        const auto rev = reverse(siz);
        for (int i = 1; i < siz; ++ i)
            if (rev[i] > i) swap(p[rev[i]], p[i]);
        for (uint n = 2, m = 1; n <= siz; m = n, n <<= 1)
        {
            auto unit = get_w(n), w = get_w(1, 0);
            for (uint l = 0, r = n - 1; r <= siz; l += n, r += n)
            {
                complex w0 = w, d, tmp0, tmp1;
                for (uint i = l, lim = l + m; i < lim; ++ i)
                    d = w0 * p[i + m], tmp0 = p[i] + d, tmp1 = p[i] - d,
                    p[i] = tmp0, p[i + m] = tmp1, w0 *= unit;
            }
        }
        if (inv) reverse(++ p.begin(), p.end());
    }

这样也可以通过上面的那个模板题：

signed main()
{
    ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    std::cout.tie(nullptr);
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
//  freopen("out.txt", "w", stdout);
#endif
    int n = scanner.nextInt(),
        m = scanner.nextInt();
    vector<FFT2::complex> a, b;
    const auto siz = FFT2::stretch(n + m);
    a.resize(siz), b.resize(siz);
    for (int i = 0; i <= n; ++ i)
        a[i].real(scanner.nextInt());
    for (int i = 0; i <= m; ++ i)
        b[i].real(scanner.nextInt());
    FFT2::fft(a), FFT2::fft(b);
    for (int i = 0; i < siz; ++ i)
        a[i] *= b[i];
    FFT2::ifft(a);
    for (int i = 0; i <= n + m; ++ i)
        cout << (int)a[i].real() << " \n"[n + m == i];
    return 0;
}

这样，FFT 就简要的介绍完了。

参考资料

更新于：2021年3月29日

算法

ACM

数论

Codeforces Round#712 div.2 回顾

补题链接：https://codeforces.com/contest/1504 记录只能说进入四月一来一直颓的可以，，不仅没怎么训练没怎么准备考研，倒是莫名其妙地出去混了一天又一天，真的...

Codeforces Round#710 div.3 回顾

比赛链接：https://codeforces.ml/contest/1506 记录比赛的时候可以说还是有点捞的，，竟然没有时间写完最后一题，属是不是很行（但是精彩还是在比完之后： _...