【多项式】FFT

【多项式】FFT

Preface

本文对所有 $\LaTeX$ 编译后生成的文本共有大约 $7000$ 字，其中前半部分为前置知识部分，介绍了多项式的有关概念、运算法则以及复数的概念、运算法则以及单位根有关内容，并证明了蝴蝶操作所用到的有关复数的两个重要引理公式。如果你对上述内容已经有了解，可以跳过 Pre-knowledge 部分。Pre-knowledge 部分大约有 $2000$ 字。

由于内容比较长，本文还没有经审阅人完全审阅完成，如果您发现了文本中的错误，请私信或评论我指出。

Pre-knowledge

多项式

Definition

称一个关于 $x$ 的式子

$$f(x) = \sum_{i = 0}^{n} a_i \times x^i$$

为一个 $n$ 次多项式，其中 $a_i$ 为常数。称 $n$ 为 $f(x)$ 的次数。显然，$f(x)$ 可以看做一个关于 $x$ 的 $n$ 次函数 $y = f(x)$。

回忆初中解析几何最后一个大题的第一问，正常情况下都是给定三个点的坐标，求一个关于 $x$ 的二次曲线解析式方程。而类似的如果求一条直线的解析式，则需要给出两个点的坐标。

类似的，对于如果想要确定一个 $n$ 次函数的解析式，则需要 $n + 1$ 个点的坐标。这是因为一个 $n$ 次函数共有 $(n + 1)$ 个系数，根据某我忘了名字的基本定理，$n$ 元非无解方程组有唯一解的必要条件是有 $n$ 个本质不同的方程，因此这里需要 $(n + 1)$ 个点构造出 $(n + 1)$ 个方程，才能求出这 $(n + 1)$ 个系数。

也就是说，只要给定了 $(n + 1)$ 个点的坐标，就可以唯一确定一个 $n$ 次函数的解析式，进而也就可以确定这个多项式。

那么，给出 $(n + 1)$ 个互不相同的点的坐标来确定一个多项式的形式，称为多项式的点值表示法。对应的，写成 $f(x) = \sum_{i = 0}^{n} a_i \times x^i$ 形式的多项式被称为多项式的系数表示法。

Operation

在进行多项式运算时，如果两个多项式的次数不一样，注意到关于多项式次数的定义中没有规定最高次项的系不能为 $0$，因此可以认为次数较低的多项式的次数也为较高的次数多项式的次数。

如果使用点值表示法进行多项式运算，则必须保证两个多项式所给点值的横坐标一一对应相等。

以下设 $n$ 次多项式 $A(x) = \sum_{i = 0}^n a_i \times x^i$，$B(x) = \sum_{i = 0}^{n} b_i \times x^i$

多项式加减法：

两个 $n$ 次多项式 $A(x)$ 和 $B(x)$ 相加减的结果为

$$f(x)~ = ~A(x) \pm B(x)~ = ~\sum_{i = 0}^{n} (a_i \pm b_i) x^i$$

多项式乘法：

两个 $n$ 次多项式 $A(x)$ 和 $B(x)$ 相乘的结果是一个 $2n$ 次多项式

$$f(x)~=~A(x) \times B(x)~=~\sum_{i = 0}^n \sum_{j = 0}^n a_i \times b_j \times x^{i + j}$$

同时，上式显然可以写成这种形式

$$f(x) = \sum_{i = 0}^{2n} \sum_{j = 0}^{\min(i, n)} a_j \times b_{i - j} \times x^i$$

对于点值表示法，由于 $A(x) \pm B(x)$ 或 $A(x) \times B(x)$ 本身表示的就是两个多项式的值做加减乘法，因此直接将对应横坐标位置的纵坐标相加/减/乘即可。需要特别注意的是，做乘法时，需要 $A(x)$ 和 $B(x)$ 各给出 $2n$ 组点值，而不是 $n$ 组。

复数

Definition

定义常数 $i$，满足

$$i^2 = -1$$

则所有形如

$$z = a + b \times i,~~~a, b \in R$$

的数字 $z$ 构成的集合称为复数集，记为 $C$，$C$ 中的每个元素都称作复数。

对于 $z = a + bi$，称 $a$ 是 $z$ 的实部，$b$ 是 $z$ 的虚部

Geometric Interpretation

复平面是一个笛卡尔平面，有两条坐标轴，纵轴为虚轴，横轴为实轴，两轴相互垂直。

对于一个复数 $z = a + bi$，它在复平面上对应一个从原点指向 $(a,~b)$ 的向量。也即实轴坐标为实部值，虚轴坐标为虚部值的点。

显然复平面上从原点出发的任意一个向量也对应唯一的一个复数，也即向量 $(a, b)$ 对应一个复数 $z = a + b_i$。

易证复数与复平面上从原点出发的向量是一一对应关系。

以实轴正方向为始边，$z$ 所对应的向量 $Z$ 为终边的角 $\theta$ 称为复数 $z$ 的幅角。

Operation

复数的模：

复数 $z= a + b_i$ 的模为其在复平面上对应向量的长度，记做$|z|$。

$$|z| = \sqrt{a^2 + b^2}$$

共轭复数：

复数 $z$ 在复平面上对应的向量关于实轴对称后对应的复数称为 $z$ 的共轭复数，记做 $\overline{z}$

$z$ 与 $\overline z$ 满足如下关系：

设 $z$ 的幅角为 $\theta_0$，$\overline z$ 的幅角为 $\theta_1$

$$\theta_0 + \theta_1 = \pi$$

$$|z| = |\overline z|$$

且两者的实部相同，虚部互为相反数。

复数加减法：

两个复数 $z_1 = a_1 + b_1i,~z_2 = a_2 + b_2 i$ 相加减的结果为

$$z_0~=~z_1 \pm z_2~=~(a_1 \pm a_2) + (b_1 \pm b_2)i$$

在复平面上，两个复数相加减的结果为他们所对应的向量按照平行四边形定则相加减

复数乘法：

两个复数 $z_1 = a_1 + b_1i,~z_2 = a_2 + b_2 i$ 相乘的结果为

$$z_0~=~z_1 \times z_2~=~(a_1 + b_1i) \times (a_2 + b_2i) = a_1a_2 + a_1b_2i + a_2b_1i + b_1b_2i^2$$

又因为

$$i^2 = -1$$

所以

$$z_0~=~(a_1a_2 - b_1b_2) + (a_1b_2 + a_2b_1) i$$

在复平面上，$z_1,~z_2,~z_0$ 所对应的幅角 $\theta_1,~\theta_2,~\theta_0$ 有如下关系：

$$\theta_0 = \theta_1 + \theta_2$$

他们的模有如下关系

$$|z_0| = |z_1| \times |z_2|$$

考虑 $z = a + b_i$ 和它的共轭复数 $\overline z = a - b_i$

有

$$z \times \overline z~=~(a + bi) \times (a - bi) = a^2 + b^2$$

因此两个互为共轭复数的数之积一定是一个实数。

复数除法：

对于两个复数 $z_1 = a_1 + b_1i,~z_2 = a_2 + b_2 i$，他们相除的结果为

$$z_0 = \frac{z_1}{z_2}$$

考虑分数上下同时乘 $\overline{z_2}$，有

$$z_0~=~\frac{z_1 \overline {z_2}}{a_2^2 + b_2^2}$$

分母是一个实数，可以直接将分子的实部虚部除以分母。

复数指数幂：

有欧拉公式

$$e^{i\theta} = \cos \theta + i \sin \theta$$

其中 $e$ 是自然对数的底数

当取 $\theta = \pi$ 时，有

$e^{i\pi} = \cos \pi + i \sin \pi$

又因为 $\cos \pi = 1,~sin \pi = 0$

所以 $e^{i\pi} = 1$

单位根

Definition

在复数域下，满足 $x^n = 1$ 的 $x$ 被称为 $n$ 次单位根

根据代数基本定理，$n$ 次单位根一共有 $n$ 个。

经过计算可得，将所有的 $n$ 次单位根按照幅角大小排列，第 $k$ $(0 \leq k < n)$ 个 $n$ 次单位根为

$$x_k~=~e^{i \frac{2 k\pi}{n}}$$

Proof

$$x_k^n~=~e^{i 2 k\pi}$$

根据欧拉公式

$$e^{i2k\pi} = \cos 2k\pi + i \sin 2k\pi$$

又因为

$$\cos 2k\pi = 1,~~\sin 2k\pi~=~0$$

所以

$$x_k^n = e^{i\frac{2k\pi}{n}}$$

是原方程的一个解。

显然这 $n$ 个解是互不相同的，又根据代数基本定理，该方程有且仅有 $n$ 个解。因此该方程解的与 $x_k$ 一一对应。证毕。

Property

因为 $\sin^2 \theta + \cos^2 \theta~=~1$，所以所有的 $n$ 次单位根的模都是 $1$。

$n$ 个单位根在复平面上平分单位圆。他们与单位圆的 $n$ 等分线所在线段分别重合。

本原单位根

Definition

$0$ 到 $(n - 1)$ 次方的值能生成所有 $n$ 次单位根的 $n$ 次单位根称为为 $n$ 次本原单位根。

显然 $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$ 次本原单位根可能不止一个，但是下文的“本原单位根”特指 $e^{i \frac{2\pi}{n}}$

Property

设 $n$ 是一个正偶数，且 $n = 2m$，有：

$$(\omega_{n}^k)^2~=~\omega_m^k$$

$$\omega_n^{m + k}~=~-\omega_{n}^k$$

Proof

一式：

考虑两个 $\omega_n^k$ 相乘，幅角相加，为之前幅角的两倍，而模为 $1 \times 1 = 1$。

所以相乘的结果为一个幅角为 $\omega_n^k$ 幅角两倍的单位向量。容易验证 $\omega_m^k$ 的幅角是 $\omega_n^k$ 幅角的两倍，且模为 $1$。由于向量和复数是一一对应的，所以一式成立。

二式：

显然 $\omega_n^m$ 的幅角为 $\pi$。等式左边可以写成 $w_{n}^k \times w_n^m$，即为 $w_n^k$ 绕原点旋转 $\pi$ 弧度。根据平面解析几何定理，旋转前后的两个向量横纵坐标分别互为相反数。证毕。

Pre-knowledge部分结束了。

---

Algorithm

考虑对两个多项式做乘法，如果运用系数表示法，显然需要 $O(n^2)$ 的时间复杂度，而如果已知两个多项式的点值表示法，则只需要 $O(n)$ 的时间复杂度。因为只需要将对应的点值纵坐标相乘就可以了。

但是我们将一个多项式从系数表示法改为点值表示法（称为求值）需要 $O(n^2)$ 的复杂度（因为每个横坐标都需要 $O(n)$ 的时间去计算），而将一个点值表示法改为系数表示法（称为插值）则需要 $O(n^3)$ 的复杂度来做高斯消元。但是只要我们将这两步都优化至低于 $O(n^2)$ 的复杂度，就可以得到一个比直接用系数表示法乘更优的做法。

而求出一个 $n$ 次多项式在每个 $n$ 次单位根下的点值的过程，被称为离散傅里叶变换(Discrete Fourier Transform，DFT)，而将这些点值重新插值成系数表示法的过程，叫做离散傅里叶逆变换(Inverse Discrete Fourier Transform,IDFT)。

以下设进行变换的多项式为 $(n - 1)$ 次多项式 $A(x) = \sum_{i = 0}^{n - 1} a_i \times x^i$。其中 $n$ 是 $2$ 的整数幂。如果不是的话，则向 $A(x)$ 的更高次数位 $n$ 补充 $a_n = 0$ ，令其成为 $n$ 次多项式，一直进行直到其次数+1的值是 $2$ 的整数幂，取 $n$ 等于其次数，$m = \frac{n}{2}$。

DFT

考虑求出一个长度为 $n$ 数列 $\{b_i\}$，这个数列的第 $k$ 项为 $A(x)$ 在 $n$ 次单位根的 $k$ 次幂处的点值。

因此有

$$b_k~=~\sum_{i = 0}^{n - 1} a_i \times \omega_{n}^i$$

注意上式中的 $i$ 为 $\Sigma$ 循环的循环变量，而不是 $-1$ 的二次方根。

这个过程是 $O(n^2)$ 的，我们考虑使用快速傅里叶变换(Fast Fourier Transform，FFT)来优化这个过程。

FFT

我们考虑对 $A(x)$ 按照系数角标的奇偶性分类，即

$$A(x)~=~\sum_{i = 0}^{n - 1} a_i x^i~=~\sum_{i = 0}^m a_{2i} \times x^{2i} + \sum_{i = 0}^m a_{2i + 1} \times x^{2^i + 1}$$

对于上式的后半部分，提出一个 $x$，得到

$$A(x)~=~\sum_{i = 0}^{m - 1} a_{2i} x^{i2} + 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}$$

设 $A_0(x)$ 是一个 $(m - 1)$ 次多项式，满足

$$A_0~=~\sum_{i = 0}^{m - 1} a_{2i} x^i$$

设 $A_1(x)$ 是一个 $(m - 1)$ 次多项式，满足

$$A_1~=~\sum_{i = 0}^{m - 1} a_{2i + 1}x^i$$

联立以上三式，可以得到

$$A(x)~=~A_0(x^2) + x \times A_1(x^2)$$

如果求出了 $A_0$ 和 $A_1$ 在各点的点值，由于上式可以 $O(1)$ 计算，所以我们可以 $O(n)$ 计算 $A(x)$ 在各个点的点值了。而求 $A_0$ 和 $A_1$ 的过程和求 $A$ 的过程完全一致，因此可以递归处理。

但是很遗憾，我们 $A_0$ 和 $A_1$ 各需要递归一次，根据主定理，这样的时间复杂度是 $O(n^2)$ 的。

但是考虑我们求的是在 $n$ 次单位根的各个幂次下的点值，根据 Pre-Knowledge 的最后一部分，我们得到了公式

$$(\omega_{n}^k)^2~=~\omega_m^k$$

$$\omega_n^{m + k}~=~-\omega_{n}^k$$

第二个式子使得我们考虑小于 $m$ 次的点值和大于 $m$ 次点值之间的关系。

对于 $0 \leq k < m$，我们有

$$A(\omega_n^k)~=~A_0((\omega_n^k)^2) + w_n^kA_1((\omega_n^k)^2)$$

根据上面的第一个公式，化简得到

$$A(\omega_n^k)~=~A_0(\omega_m^k) + w_n^kA_1(\omega_m^k)$$

这只是前半部分的次数，我们考虑后半部分次数：

$$A(\omega_n^{m + k})~=~A_0((\omega_n^{m + k})^2) + \omega_n^{m + k} A_1((\omega_n^{m + k})^2)$$

根据第二个公式，化简得到

$$A(\omega_n^{m + k})~=~A_0((w_n^k)^2)+-\omega_n^k A_1((\omega_n^k)^2)$$

再应用第一个公式得到

$$A(\omega_n^{m + k})~=~A_0(\omega_m^k) - w_n^kA_1(\omega_m^k)$$

我们惊喜的发现，大于 $m$ 次的点值可以由 $A_0$ 和 $A_1$ 在小于 $m$ 次的点值求出。只要求出了 $A_0$ 和 $A_1$ 在小于 $m$ 次的点值，就可以线性求出 $A$ 在整个 $n$ 次幂处的点值。而求 $A_0$ 和 $A_1$ 在小于 $m$ 次的点值也是可以递归求解的。

考虑时间复杂度：递推关系为 $T(n)~=~2T(n / 2) + O(n)$。因为 $O(n)~=~\Theta(n)$，所以 $T(n)~=~\Theta(n^{\log_2^2} \log n) = \Theta(n \log n)$。

由此，得到了快速计算 DFT 的办法，并证明了其时间复杂度为 $O(n \log n)$。这种方法被即为 FFT。

上面推导出 $A(\omega_n^k)$ 和 $A(\omega_n^{m + k})$ 的值的式子被称为蝴蝶操作(Butterfly Operation)。这个名字的由来是如果将 $A$ 在各处的点值画成一个长条形的数组，在数组下面一次依次是 $n$ 次本原单位根的 $0 \sim m - 1$ 次幂，则求值过程可以画成 $\omega_n^0$ 连一条边向 $A(\omega_n^0),~A(\omega_n^{\frac{n}{2} - 1})$，$\omega_n^1$ 连向 $A(\omega_n^1),~A(\omega_n^\frac{n}{2})$，以此类推。画出的图形如同蝴蝶的翅膀。

IDFT

至此，我们已经有了 $O(n \log n)$ 的算法来计算两个系数表示法的多项式相乘后的点值表示。接下来我们只需要用 $O(n \log n)$ 的时间复杂度完成插值的过程，就可以得到一个完整的 $O(n \log n)$ 的系数型多项式乘法算法了。

我们目前已知一个 $(n - 1)$ 次多项式 $A(x)~=~\sum_{i = 0}^{n - 1} a_i x^i$ 进行了离散傅里叶变换后的点值 $\{b_i\}$，即

$$b_k~=~\sum_{i = 0}^{n - 1} a_i \times \omega_n^{ik}$$

现在试图还原系数数列 $\{a_i\}$。

推导过程比较复杂，我们直接给出结论：

$$a_k~=~\frac{1}{n} \sum_{i = 0}^{n - 1} b_i \omega_n^{-ki}$$

下面证明上面这个式子是 DFT 式子(即上面计算 $b_k$ 的式子)的逆变换。

我们首先将 DFT 的式子带入 $\sum_{i = 0}^{n - 1} b_i \omega_n^{-ki}$

得到上式为

$$\begin{align} \sum_{i = 0}^{n - 1} b_i \omega_n^{-ki}~ & =~\sum_{i = 0}^{n - 1} \sum_{j = 0}^{n - 1} a_j \times \omega_n^{ji} \times \omega_n^{-ki}\\ & =~\sum_{i = 0}^{n - 1}\sum_{j = 0}^{n - 1} a_j \times \omega_n^{ij - ki}\\ & =~\sum_{j = 0}^{n - 1} a_j \times \sum_{i = 0}^{n - 1} \omega_n^{i(j - k)} \end{align}$$

考虑分类讨论。

1、当 $j = k$ 时：

则 $j - k = 0$，因此 $\omega_{n}^{i(j - k)} = \omega_n^0 = 1$。于是

$$\sum_{i = 0}^{n - 1} \omega_n^{i(j - k)}~=~n \times 1 = n$$

2、当 $j \neq k$ 时：

显然 $|j - k| < n$，因此 $\omega_n^{j - k} \neq 1$，则 $\sum_{i = 0}^{n - 1} \omega_n^{(j - k) i}$ 是一个公比不为 $1$ 的等比数列前缀和。根据等比数列求和公式得到

$$\sum_{i = 0}^{n - 1} \omega_n^{(j - k) i}~=~\frac{1 - \omega_n^{(n - 1)(j - k)} \times \omega_n^{j - k}}{1 - \omega_{n}^{j - k}}~=~\frac{1 - \omega_n^{n(j - k)}}{1 - \omega_{n}^{j - k}}~=~\frac{1 - (\omega_n^{(j - k)})^n}{1 - \omega_{n}^{j - k}}$$

根据复数的几何性质，易证对于所有的 $n$ 次单位根 $T$，$T^{x + n}~=~T^{x}$，其中 $x$ 为任意整数。

$\sum_{i = 0}^{n - 1} \omega_n^{(j - k) i}~=~\frac{1 - (\omega_n^{(j - k)})^0}{1 - \omega_{n}^{j - k}}~=~\frac{1 - 1}{1 - \omega_n^{j - k}}~=~0$

第二种情况的证明被称为消去引理。

综上讨论，当 $j \neq k$ 时，因为另一个因数是 $0$，前面的 $\sum a_j$ 不会对式子产生贡献，而 $j = k$ 时，会对答案产生 $n$ 倍的贡献。

原式可以写成

$$\sum_{i = 0}^{n - 1} b_i \omega_n^{-ki}~=~~\sum_{j = 0}^{n - 1} a_j \times \sum_{i = 0}^{n - 1} \omega_n^{i(j - k)}~=~\sum_{j = 0}^{n - 1} a_j \times [j = k] \times n~=~a_k \times n$$

将上式带入 IDFT 的式子 $a_k~=~\frac{1}{n} \sum_{i = 0}^{n -1} b_i \omega_n^{-ki}$ 的右边，得到

$$\text{右边}~=~\frac{1}{n} a_k \times n~=~a_k~=~\text{左边}$$

类似的，可以证明将 IDFT 的式子带入 DFT 也是可以使等式成立的，这里略去。

由此，我们证明了变换

$$a_k~=~\frac{1}{n} \sum_{i = 0}^{n - 1} b_i \omega_n^{-ki}$$

是 DFT 的逆变换，称为 IDFT。我们可以通过这个式子求出这个多项式的系数表示法。

IFFT

下面的问题是如何用较低的复杂度计算 $B(x)~=~\sum_{i = 0}^{n - 1} b_i \times x^i$ 在 $w_n^{-ki}$，其中 $0 \leq k < n$ 处的点值。

而 $w_n^{-k}$ 可以看做 $n$ 次本原单位根每次逆时针旋转本原单位根幅角的弧度，因此 $\omega_n^{-k}$ 和 $\omega_n^k$ 是一一对应的。具体的，$w_n^{-k} = w_n^{k + n}$。因此我们只需要使用 FFT 的方法，求出 $B(x)$ 在 $\omega_n$ 各个幂次下的值，然后数组反过来，即令 $a_k~=~\frac{1}{n} \sum_{i = 0}^n B(w_n^{n - k})$ 即可。

这一步快速计算插值的过程叫做快速傅里叶逆变换(Inverse Fast Fourier Transform，IFFT)。

至此，我们得到了一个时间复杂度为 $O(n \log n)$ 的多项式乘法计算方法。

Code

根据上面的推导，我们可以很轻松的写出 FFT 的递归形式：

void FFT(std::complex<double> *A, int N) {
  if (N == 1) {
    return;
  }
  int M = N >> 1;
  std::complex<double> A0[M], A1[M];
  for (int i = 0; i < M; ++i) {
    A0[i] = A[i << 1];
    A1[i] = A[(i << 1) | 1];
  }
  FFT(A0, M); FFT(A1, M);
  auto W = std::complex<double>(cos(1.0 * PI / M), sin(1.0 * PI / M)), w = std::complex<double>(1.0, 0.0);
  for (int i = 0; i < M; ++i) {
    A[i] = A0[i] + w * A1[i];
    A[i + M] = A0[i] - w * A1[i];
    w *= W;
  }
}

以及 IFFT：

void IFFT(std::complex<double> *A, int N) {
  FFT(A, n)
  std::reverse(A + 1, A + N);
}

需要注意的是，因为 $\omega_n^0~=~\omega_n^n$，所以是第 $1$ 个点值和第 $N - 1$ 个交换，而不是第 $0$ 个点值和第 $N - 1$ 个交换。

optimization

上面这个 FFT 交到 luogu 上以后只有 $77$ 分。他的复杂度显然是 $O(n \log n)$ 的，但是递归和动态开空间带来的巨大常数让他难以通过 $10^6$ 的数据。

我们考虑优化上面的代码。

首先我们注意到我们动态申请和删除了 $O(n \log n)$ 的空间，但是我们同一时刻只最多需要 $2N$ 的空间，而递归调用又是栈式的，即先申请的后删除。这意味着如果用一个数组做内存池，没有被分配的内存总是连续的，被分配的内存也总是连续的。并且分配一个数组可以 $O(1)$ 而不是 $O(size)$ 完成。

然而这个优化并没有什么卵用。

我们考虑将上面的递归改成迭代。

既然递归是自上而下的，那么我们的迭代就是一个自下而上的合并过程。当 $n = 8$ 时，我们考虑递归时调用原多项式系数下标的递归树：

step 1: 0 1 2 3 4 5 6 7
step 2: 0 2 4 6,  1 3 5 7
step 3: 0 4,  2 6,  1 3,  5 7
step 4: 0,  4,  2,  6,  1,  3,  5,  7

上表的 step 代表了递归的层数，冒号后的数字代表数字的下标。我们考察最后算的一层，也就是第 $4$ 层的二进制值：

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

看起来还是没有头绪，但是我们将二进制值反过来

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

我们发现上面 $8$ 个二进制的排列是单调递增的，在十进制下分别是

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

因此我们只需要将数列 $\{a\}$ 按照下标二进制翻转后的大小排序，得到序列 $\{a'\}$。递归倒数第二层对于 $\{a\}$ 的蝴蝶操作就变成了对于 $\{a'\}$ 相邻两个数进行合并。而上面几层同理。

因此问题变成了如何在 $O(n \log n)$ 的时间内将序列按照下标二进制翻转后排序得到的序列。

显然对于任何一个数我们都可以 $O(\log n)$ 的运用进制转换来确定其二进制逆序值，但是这样因为涉及到了大量的除法和取模，常数很大，我们考虑规避这个做法。

我们考虑基数排序的操作，我们可以将一个数的数位分为两半，先对后半部分数位构成的数进行排序，然后再对前半部分数位构成的数进行排序。两个数字的前半部分数位如果相同，那么它们的先后顺序即为后半部分数位的先后顺序。这样由于后半部分本身是有序的，就可以自低位向高位对数列进行排序。

同样的，我们自低位向高位对每个数的二进制逆序排序。

初始时，序列里只有两个数

000, 100

他们逆序后为

000, 001

分别是最小的和次小的翻转值。

然后我们考虑将这两个数的第 $2$ 位(最右侧为最低位，最低位为第 $1$ 位)都置为 $1$，得到了两个数

010, 110

将这两个数也加入序列中，得到

000, 100, 010, 110

翻转值即为

000, 001, 010, 011

注意到第 $i$ 次操作的时候，加入了大于 $i$ 位的位置都是 $0$，且第 $i$ 位是 $1$ 的所有数，他们显然比序列中原来存在的大于等于 $i$ 为的位置都是 $0$ 的数要大，而比大于 $i$ 位存在 $1$ 的数小，并且去掉 $1$ 以后即为序列中原有数字的顺序。这样我们就证明了这个方法的正确性。

而我们只需要处理 $O(\log n)$ 位，每位的处理都是 $O(n)$ 的，因此总时间复杂度 $O(n \log n)$，空间复杂度 $O(n)$，事实上对于每位的处理是跑不满 $O(n)$ 的，并且规避了取模和除法操作，常数极小。

代码如下

void MakeRev(const int N) {
  int d = N >> 1, p = 0;
  tax[p++] = 0;
  tax[p++] = d;
  for (int w = 2; w <= N; w <<= 1) {
    d >>= 1;
    for (int p0 = 0; p0 < w; ++p0) {
      tax[p++] = tax[p0] | d;
    }
  }
}

然后考虑通过我们得到的 tax 数组来对原数组 $A$ 进行排序。

记 $rev(i)$ 为 $i$ 二进制逆序后的值，显然 $rev(rev(i)) = i$。因此对于一个下标 $p$，设它排序后在新序列中的下标为 $q$，那么原序列下标为 $q$ 的位置在新序列的下标一定是 $p$。所以我们只需要对于所有 $rev(i) > i$ 的位置 $i$，交换 $i$ 和 $rev(i)$ 位置的值即可。

这部分操作的名字叫做位逆序置换。

for (int i = 1; i < N; ++i) if (tax[i] > i) {
  std::swap(A[i], A[tax[i]]);
}

最后考虑迭代的蝴蝶操作过程。

对于求$A(x)$ 在 $n$ 次单位根的各幂次的点值时，$m = \frac{n}{2}$ 次单位根的各幂次在 $A_0$ 或 $A_1$ 处的点值已经被计算并存储在了 $A$ 数组中。我们令 $A_0(\omega_m^k)~=~A[k]$，$A_1(\omega_m^k)~=~A[k + m]$。那么 $A(\omega_n^k)~=~A[k] + \omega_n^k \times A[k + m]$，$A(\omega_n^{k + m})~=~A[k]~-\omega_n^k \times A[k + m]$，并将上面两个值分别存入 $A[k]$ 和 $A[k + m]$ 中即可。其中 $A(x)$ 代表多项式，$A[k]$ 代表数组。那么在更上面一层递归调用时，根据 $A_0$ 和 $A_1$的定义，当前的 $A[k]$ 即是上面一层所计算的多项式 $A'(x)$ 的导出多项式 $A'_0(x)$ 在 $\omega_n^k$ 处的点值，$A[k + m]$ 即是 $A'_1(x)$ 在 $\omega_n^k$ 处的点值。

至此，我们得到了一个迭代完成 FFT 的算法。

Code

#include <cmath>
#include <cstdio>
#include <cstring>
#include <complex>
#include <iostream>
#include <algorithm>

const int maxn = 6000006;
const double PI = acos(-1);

int n, m, k, sm;
int tax[maxn];
std::complex<double> F[maxn], G[maxn];

void MakeRev(const int N);
void FFT(std::complex<double> *A, int N);

int main() {
  freopen("1.in", "r", stdin);
  qr(n); qr(m);
  for (int i = 0, x; i <= n; ++i) {
    x = 0; qr(x); F[i] = x;
  }
  for (int i = 0, x; i <= m; ++i) {
    x = 0; qr(x); G[i] = x;
  }
  sm = n + m; k = 1;
  while (k <= sm) k <<= 1;
  MakeRev(k);
  FFT(F, k); 
  FFT(G, k);
  for (int i = 0; i < k; ++i) {
    F[i] *= G[i];
  }
  FFT(F, k);
  std::reverse(F + 1, F + k);
  for (int i = 0; i < sm; ++i) {
    qw(int((F[i].real()) / k + 0.5), ' ', true);
  }
  qw(int(F[sm].real() / k + 0.5), '\n', true);
  return 0;
}

void FFT(std::complex<double> *A, int N) {
  for (int i = 1; i < N; ++i) if (tax[i] > i) {
    std::swap(A[i], A[tax[i]]);
  }
  for (int len = 2, M = 1; len <= N; M = len, len <<= 1) {
    std::complex<double> W(cos(PI / M), sin(PI / M)), w(1.0, 0.0);
    for (auto L = 0, R = len - 1; R <= N; L += len, R += len) {
      auto w0 = w;
      for (auto p = L, lim = L + M; p < lim; ++p) {
        auto x = A[p] + w0 * A[p + M], y = A[p] - w0 * A[p + M];
        A[p] = x; A[p + M] = y;
        w0 *= W;
      }
    }
  }
}

void MakeRev(const int N) {
  int d = N >> 1, p = 0;
  tax[p++] = 0;
  tax[p++] = d;
  for (int w = 2; w <= N; w <<= 1) {
    d >>= 1;
    for (int p0 = 0; p0 < w; ++p0) {
      tax[p++] = tax[p0] | d;
    }
  }
}

Appreciation

十分感谢 @Dusker 抽出时间审阅这篇长文。