【算法笔记】快速 Fourier 变换（FFT）

合集 - OI算法笔记(13)

1.【算法笔记】K.M.P.字符串匹配2022-01-182.【算法笔记】Tarjan 算法 · 上2022-01-213.【算法笔记】Manacher 算法2022-01-314.【算法笔记】二维凸包2022-01-305.【算法笔记】扩展 Euclid 算法2022-01-276.【算法笔记】Tarjan 算法 · 下2022-01-247.【算法笔记】Lagrange 插值法2022-04-148.【算法笔记】中国剩余定理（CRT）与扩展中国剩余定理（exCRT）2022-04-129.【算法笔记】Baby Step Giant Step（BSGS）及其扩展2022-04-0910.【算法笔记】Gauss 消元 · 行列式2022-04-0311.【算法笔记】Euler-Fermat 定理2022-03-19

12.【算法笔记】快速 Fourier 变换（FFT）2022-03-13

13.【算法笔记】乘法逆元2022-02-01

收起

• 本文总计约 20000 字，阅读大约需要 60 分钟。
• 警告！警告！警告！本文有大量的 $\LaTeX$ 公式渲染，可能会导致加载异常缓慢！

前言

我终于学会 FFT 了 QwQ！虽然我非常的菜，但是还是辛苦地把这些令人费解的知识给理解下来了，太难了~

或许初学 FFT 的人都会有和我一样的感觉：这在说的是啥？什么叫 DFT？蝴蝶变换又是什么奇葩的东西？Vandermonde 矩阵又是什么？

但是当我们逐渐平静下来，再来仔细地研究这个算法时，就会发现，FFT 真的是很美的一个算法。当我们看见单位根带入了奇偶分组的多项式时，仿佛真的看见了一个由系数组成的蝴蝶在翩翩起舞；当一个多项式被从中间一分为二的时候，我们更能感受到分治地美感之所在。

前人的智慧，总是会让我们叹服叫绝。但我们又应该如何去应用呢？如何把这些智慧内化于心，并且改进出更好的理论呢？我们又不能不感到任重而道远。

正好最近在研读《论语》，用其中的一言以蔽：这不正是做学问做应该有的，『仰之弥高，钻之弥坚。瞻之在前，忽焉在后』的灵魂吗？

大概，科学真正的内涵，就寓于这十六个字里吧。

题目引入

• $\text{Take 1}:$ 给定两个正整数 $a, b$，求 $a\times b$ 的值。其中 $1\le a,b\le 10^9$。
  『CaO 你这不是糊弄人吗，A*B Problem 还用得着你来讲吗？』

#include <cstdio>
using namespace std;

long long a, b;
int main(void) {
	scanf("%lld%lld", &a, &b);
	printf("%lld", a * b);
	return 0;
}
//by CaO

• $\text{Take 2}:$ 给定两个正整数 $a, b$，求 $a\times b$ 的值。其中 $1\le a,b\le 10^{1000}$。
  『嗯……看样子不能直接 A*B 了，写个高精罢。』
  『然而我会用 Python（逃）』

#include <iostream>
#include <cstring>
using namespace std;
const int MAXN = 10001;

char A[MAXN], B[MAXN];
int n, m, high;
int a[MAXN], b[MAXN], ans[MAXN];

int main(void) {
    cin >> A >> B;

    n = strlen(A), m = strlen(B);

    for(int i = 0; A[i]; ++i) a[i] = A[n - 1 - i] - '0';
    for(int i = 0; B[i]; ++i) b[i] = B[m - 1 - i] - '0';

    for(int i = 0; i < n; ++i) {
        for(int j = 0; j < m; ++j) {
            ans[i + j] += a[i] * b[j];
            ans[i + j + 1] += ans[i + j] / 10;
            ans[i + j] %= 10;
        }
    }

    high = 10000;
    while(!ans[high] && high) --high;
    while(high >= 0) cout << ans[high--];
    return 0;
}
//by CaO

• $\text{Take 3}:$ 给定两个正整数 $a, b$，求 $a\times b$ 的值。其中 $1\le a,b\le 10^{10000000}$。
  『……』

如果像解决 $\text{Take 2}$ 那样写简单的高精，那么时间复杂度为 $\Theta(n^2)$，其中 $n$ 代表 $a,b$ 位数的几何平均值。在 $a,b$ 都有 $10^7$ 位的时候，再用 $\Theta(n^2)$ 的乘法就会 TLE 了。

那么，有没有更高效的，像 $\Theta(n \log n)$ 时间复杂度的乘法算法呢？

前置知识

学习 FFT 的前置知识非常~非常~得多 QwQ，主要包含下面的这几条：

• 多项式的系数表示
• 多项式的点值表示
• 复数的乘法与单位根

如果您已经保证了上面的这些知识都已经掌握了，那么请您跳过这一段吧 QwQ~

多项式

我们在初中的时候就已经接触过了多项式相关的问题，但是我还是应该给出多项式的定义，形如：$F(x)=\sum_{i=0}^n{a_ix^i}=a_0+a_1x+a_2x^2+\cdots+a_nx^n$ 的式子，我们称为多项式。其中的 $a_i(i=1,2,\cdots,n)$ 为 $x$ 的 $i$ 次项系数，$a_0$ 为多项式的常数。

我们定义两个多项式 $F(X)$ 与 $G(X)$ 的卷积为 $(F* G)(x)$，如果有 $F(x)=\sum_{i=0}^{n}{a_ix^i}$，$G(x)=\sum_{i=0}^m{b_ix^i}$，那么我们有 $(F* G)(x)=\sum_{i=0}^{n+m}{c_ix^i}$。其中，$c_i=\sum_{j=0}^n{a_jb_{i-j}}$。

实际上，两个多项式的卷积就是它们的乘积。例如 $F(x)=x+1,G(x)=x^2-3x+1$，那么就有 $(F* G)(x)=F(x)\cdot G(x)=x^3-2x^2-2x+1$。

点值表示法

我们平时接触的多项式，都是用上面的，类似于 $F(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n$ 的形式表示的。像这样，用 $F(x)$ 的 $(n-1)$ 个系数来表示 $F(x)$ 的表示方法，叫做系数表示法。

假如我们向这个多项式中代入 $(n+1)$ 个值，记做 $x_0, x_1, \cdots, x_n$，我们就可以得到 $(n+1)$ 个值，记作 $F(x_0), F(x_1), F(x_2), \cdots, F(x_n)$。

由 Lagrange 插值定理，这 $(n+1)$ 个点 $(x_0, F(x_0)), (x_1, F(x_1)), (x_2, F(x_2)), \cdots, (x_n, F(x_n))$，可以唯一确定一个多项式 $F(x)=\sum_{i=0}^na_ix^i$。

那么，用这 $(n+1)$ 个点来表示一个多项式的表示方法，叫做点值表示法。

如果用点值表示法来表示两个最高次数均为 $n$ 的多项式 $F(x), G(x)$ 的话，我们分别代入 $(2n+1)$ 个值 $x_0, x_2, x_2, \cdots, x_{2n}$（为什么要代入 $(2n+1)$ 个值呢 QwQ？留给读者思考吧），那么$F(x)$ 和 $G(x)$ 都可以被表示为：

$$\begin{cases} F(x)=(x_0, F(x_0)), (x_1, F(x_1)), \cdots, (x_{2n}, F(x_{2n})); \\ G(x)=(x_0, G(x_0)), (x_1, G(x_1)), \cdots, (x_{2n}, G(x_{2n})). \end{cases}$$

如何求这两个多项式的卷积呢？我们可以直接拿这 $(2n+1)$ 个点的数值一一对应相乘，像这样：

$$(F* G)(x)=(x_0, F(x_0)\cdot G(x_0)),(x_1, F(x_1)\cdot G(x_1)),\cdots,(x_{2n}, F(x_{2n})\cdot G(x_{2n})).$$

这样求卷积的时间复杂度为 $\Theta(n)$。

『哇，那么就是说求多项式的卷积的时间复杂度被优化到线性的了吗？』

错！因为我们对于每个多项式 $F(x)$，至少要代入 $(2n+1)$ 个点才能将其转化为可以用于求卷积的点值表示法。而对每个点来说，我们可以应用秦九韶算法在 $\Theta(n)$ 内求单个点的点值，这个过程，也就是『系数转点值』总复杂度达到了 $\Theta(n^2)$。

而且，就算我们真的求出了 $(F* G)(x)$ 的结果，也是用点值表示法表示的，但我们平时使用的都是系数表示法。想象一下，如果我问你 $F(x)=(x+1)(x+2)$ 的结果，你会回答结果是 $((-1,0),(0,2),(1,6))$ 吗？

我们当然可以用 Gauss 消元法把点值在转化为系数表示，但它的时间复杂度为 $\mathcal{O}(n^3)$，还不如直接求卷积呢！

『呃……』

复数

虽然直接用点值表示法来计算多项式卷积看上去效率并不高，但是实际上，它为我们提供了该选择代入哪些值的自由。也就是说，如果我们选得点足够好，可能会让点值法求卷积的效率发生质的飞跃。

我们考虑引入复数的概念。

定义虚数单位 $\text i^2=-1$，那么形如 $z=a+b\text i$ 的数，我们称之为复数，其中的 $a$ 称作复数 $z$ 的实部，$b$ 称作复数 $z$ 的虚部。

由所有复数组成的集合记为 $\mathbb{C}$，即 $\mathbb{C}=\{a+b\text i|a,b\in\mathbb R\}$。

复数的四则运算

与实数的加减法类似，我们可以定义复数的四则运算：对于任意的两个复数 $z_1=a+b\text i,z_2=c+d\text i$，我们定义：

$$z_1\pm z_2=(a+b\text i)\pm(c+d\text i)=(a\pm c)+(b\pm d)\text i,\\ \\ z_1z_2=(a+b\text i)(c+d\text i)=(ac-bd)+(bc+ad)\text i,\\ \\ \dfrac{z_1}{z_2}=\dfrac{a+b\text i}{c+d\text i}=\dfrac{(a+b\text i)(c-d\text i)}{(c+d\text i)(c-d\text i)}=\dfrac{(ac+bd)+(bc-ad)\text i}{c^2+d^2}.$$

复数的坐标表示

既然复数有实部和虚部两个参数，那么我们就可以用这两个参数作为一个点的横纵坐标，在平面直角坐标系中表示出任意一个复数了。

换言之，对于复数 $z=x+y\text i$，我们可以用平面中的一个点 $(x,y)$ 来代表复数 $z$。如下图，点 $A(3,2)$ 就代表了复数 $z=3+2\text i$。

我们记坐标原点为 $O$，有向线段 $|OA|$ 的长度记为 $r$，$OA$ 与坐标系中横轴正半轴的夹角记为 $\theta$，那么对于复数 $z=x+y\text i$，都有转换关系： $\begin{cases}r=\sqrt{x^2+y^2},\\ x=r\cos\theta,\\ y=r\sin\theta\end{cases}$。

我们定义上述描述中的 $r$ 为复数 $z$ 的模长，$\theta$ 为复数 $z$ 的辐角。显然，$r\ge0,\theta\in[0,2\pi)$。

复数的三角形式

对于复数 $z_1,z_2$，设 $z_1$ 的模长和辐角分别为 $r_1,\alpha$，$z_2$ 的模长和辐角为 $r_2,\beta$，那么就有

$$\begin{cases} z_1=r_1\cos\alpha+\text ir_1\sin\alpha=r_1(\cos\alpha+\text i\sin\alpha),\\ z_2=r_2\cos\beta +\text ir_2\sin\beta =r_2(\cos\beta +\text i\sin\beta). \end{cases}$$

像这样表示复数的方式，称为复数的三角形式。

对于上面的两个复数 $z_1,z_2$，我们可以计算出来：

$$\begin{aligned} z_1z_2 &= r_1(\cos\alpha+\text i\sin\alpha)\cdot r_2(\cos\beta+\text i\sin\beta)\\ &= r_1r_2(\cos\alpha+\text i\sin\alpha)(\cos\beta+\text i\sin\beta)\\ &= r_1r_2[(\cos\alpha\cos\beta-\sin\alpha\sin\beta)+\text i(\cos\alpha\sin\beta+\sin\alpha\cos\beta)]\\ &= r_1r_2(\cos(\alpha+\beta)+\text i\sin(\alpha+\beta))\\ &= r_1r_2\cos(\alpha+\beta) +\text ir_1r_2\sin(\alpha+\beta). \end{aligned}$$

从这里，我们也可以知道，对于任意的两个复数，将它们的乘的结果，即为两复数的辐角相加，模长相乘。

当然，伟大的数学家 $\mathcal{Euler}$ 也给出了大名鼎鼎的 Euler 定理，即：

$$\text e^{\text i\theta}=\cos\theta+\text i\sin\theta,\quad \text{which e} =\lim_{n\rightarrow+\infin}\left(1+\dfrac{1}{n}\right)^n.$$

这就让我们能够用更简便的方式来证明上面的结论：若 $z_1=r_1\text e^{\text i\alpha},z_2=r_2\text e^{\text i\beta}$，那么 $z_1z_2=r_1r_2\text e^{\text i(\alpha+\beta)}$ 当然这不重要啦 QwQ。

单位根

根据 Gauss 代数基本定理，$n$ 次方程 $x^n=1$ 在复数域内有且仅有 $n$ 个根，通过上面的结论，我们不妨设 $x=r\text e^{\text i\theta}$，那么一定有：$r^n=1,n\theta=2k\pi$，其中 $k=0,1,2,\cdots,n-1$，我们取 $k=1$ 时的解，记做 $\omega_n^1=\text e^\frac{2\pi\text i}{n}$，显然，这 $n$ 个根分别可以记做 $\omega_n^0,\omega_n^1\cdots,\omega_n^{n-1}$，这 $n$ 个复数，我们定义为 $n$ 次单位根。

让我们更直观地看一下，这是 $5$ 次单位根 $\omega_5^0,\omega_5^1,\omega_5^2,\omega_5^3,\omega_5^4$ 在平面直角坐标系上的表示：

它们均匀地分布在单位圆上，并且将单位圆 $5$ 等分。

显然，单位根有如下的两个性质（划重点！后面要用到）：

1. $\omega_{2n}^{n+k}=-\omega_{2n}^k$，这个性质又被称为消去引理；
2. $\omega_{2n}^{2k}=\omega_n^k$，这个性质又被称为折半引理。

这样的性质，使得单位根有成为了点值表示法的点值的优势。为什么这么说？请往下看。

快速 Fourier 变换（FFT）

从 DFT 变换谈起

对于多项式 $F(x)=\sum_{i=0}^{n-1}a_ix^i$，我们考虑分别代入 $n$ 个 $n$ 次单位根，即 $\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}$，那么我们就可以获得 $n$ 个点：

$$F(x)=((\omega_n^0,F(\omega_n^0)),(\omega_n^1,F(\omega_n^1)),\cdots,(\omega_n^{n-1},F(\omega_n^{n-1}))).$$

这个代入的名称叫做离散 Fourier 变换（Discrete Fourier Transform，DFT）。尽管这个过程并不是那么令人惊讶，当然，朴素 DFT 的时间复杂度依旧是 $\Theta(n^2)$。

但是，如果我们把 DFT 的过程做一个小小的变换，那么就可以得到很大的优化。

快速 Fourier 变换的思想

前方高能，系好安全带 QwQ！

我们假定 $n$ 是一个偶数，对于多项式 $F(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}$，那么我们可以按照次数的奇偶性，把系数分组，如下：

$$\begin{aligned} F(x) &= (a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+\cdots+a_{n-1}x^{n-1})\\ &= (a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots+a_{n-1}x^{n-2}). \end{aligned}$$

我们不妨设 $A_1(x)=a_0+a_2x+\cdots+a_{n-2}x^{\frac{n}{2}-1},A_2(x)=a_1+a_3x+\cdots+a_{n-1}x^{\frac{n}{2}-1}$，那么

$$F(x)=A_1(x^2)+xA_2(x^2).$$

假如我们代入 $\omega_n^k\left(k=0,1,\cdots,\dfrac{n}{2}-1\right)$，则有：$F(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})$；

假如我们代入 $\omega_n^{k+\frac{n}{2}}\left(k=0,1,\cdots,\dfrac{n}{2}-1\right)$，则有：$F(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n})$。

由消去引理，则有 $F(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k})-\omega_n^{k}A_2(\omega_n^{2k})$；
由折半引理，则有

$$\begin{cases} F(\omega_n^k)&=A_1(\omega_\frac{n}{2}^k)+\omega_n^kA_2(\omega_\frac{n}{2}^k),\\ F(\omega_n^{k+\frac{n}{2}})&=A_1(\omega_\frac{n}{2}^k)-\omega_n^{k}A_2(\omega_\frac{n}{2}^k). \end{cases}$$

我们注意到上面这两个式子中，只有 $A_2$ 的系数的符号有所不同。所以，我们带入 $\omega_n^0,\omega_n^1,\cdots,\omega_n^{\frac{n}{2}-1}$ 的值，也就是一半的点值，就能够得到所有的点值了。比起朴素的 DFT，问题的规模便减少了一半！

既然 $A_1(x),A_2(x)$ 都是多项式，我们就可以像上面那样分治地完成这样一个过程。这时，时间复杂度为 $T(n)=2T\left(\dfrac{n}{2}\right)+\text O(n)$，解得 $T(n)=\Theta(n\log n)$。这样，多项式的乘法就比暴力的 $\Theta(n^2)$ 的乘法高效了很多！因此，这个操作又被称为快速 Fourier 变换（Fast Fourier Transform，FFT）。

怎么转回系数表示法呢？

怎么转回系数表示？

我们不妨定义多项式的系数组成了一个列向量，如下：

$$A=\left[ \begin{matrix} a_0\\ a_1\\ \vdots\\ a_{n-1} \end{matrix} \right].$$

又有 Vandermonde 矩阵：

$$B=\left[ \begin{matrix} \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\\ \omega_n^0 & \omega_n^1 & \cdots & \omega_n^{n-1}\\ \vdots & \vdots & \ddots & \vdots \\ \omega_n^0 & \omega_n^{n-1} & \cdots & \omega_n^{(n-1)^2} \end{matrix} \right].$$

可知：

$$A\times B=\left[ \begin{matrix} f(\omega_n^0) \\ f(\omega_n^1) \\ \vdots \\ f(\omega_n^{n-1}) \end{matrix} \right].$$

经计算（虽然我太蒻了 QwQ，不会算）得：

$$B^{-1}=\dfrac{1}{n}\left[ \begin{matrix} \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\\ \omega_n^0 & \omega_n^{-1} & \cdots & \omega_n^{-(n-1)}\\ \vdots & \vdots & \ddots & \vdots \\ \omega_n^0 & \omega_n^{-(n-1)} & \cdots & \omega_n^{-(n-1)^2} \end{matrix} \right].$$

那么，就可以得知：

$$A=\left[ \begin{matrix} f(\omega_n^0) \\ f(\omega_n^1) \\ \vdots \\ f(\omega_n^{n-1}) \end{matrix} \right]\cdot B^{-1}.$$

当然，如果上面的东西你看不懂，那就直接记结论吧：把 $n$ 个点值分别记作 $n$ 新多项式系数，代入 $\omega_n^{0},\omega_n^{-1},\cdots,\omega_n^{-(n-1)}$ 再做一次 FFT，得到新的 $n$ 个点值，把这些点值都除以 $n$，就得到了原多项式的各项系数。

上面的这步操作，叫做快速 Fourier 逆变换（Inverse Fast Fourier Transform，IFFT）。时间复杂度嘛，自然与 FFT 是一样的，也是 $\Theta(n\log n)$。

所以，我们是可以做到在 $\Theta(n\log n)$ 的复杂度内完成求多项式的卷积的。

代码

C++ 自带了 <complex> 库，但是封装的 complex 类常数大，在实战的时候容易被卡常，不如自己手动重载的好。

这个是递归的 FFT 代码。

#include <cstdio>
#include <cmath>
#define reg register  //register，卡常技巧来一个 QwQ
using namespace std;
const int MAXN = 5000001;
const double PI = acos(-1.0);

inline long long read(void) {  //为了防卡常，打一个快读 QwQ
    long long s = 1, w = 0; char ch = getchar();

    while(ch < '0' || ch > '9') {if(ch == '-') s = -1; ch = getchar();}
    while(ch >= '0' && ch <= '9') {w = w * 10 + ch - '0'; ch = getchar();}

    return s * w;
}

int n, m, siz = 1;
struct Complex {  //手动重载 Complex 类
    double Re, Im;
    Complex(double R = 0.0, double I = 0.0) {Re = R, Im = I;}
} f[MAXN], g[MAXN];

Complex operator+ (Complex a, Complex b) {return Complex(a.Re + b.Re, a.Im + b.Im);}
Complex operator- (Complex a, Complex b) {return Complex(a.Re - b.Re, a.Im - b.Im);}
Complex operator* (Complex a, Complex b) {return Complex(a.Re * b.Re - a.Im * b.Im, a.Re * b.Im + a.Im * b.Re);}

void fft(Complex* arr, int siz, int tp) {  //tp 代表变换的类型，如果 tp = 1 代表 FFT，tp = -1 代表 IFFT
    if(siz == 1) return ;  //只有一个常数项时，递归达到出口

    Complex arr1[siz >> 1], arr2[siz >> 1];

    for(reg int i(0); i < siz; i += 2) {arr1[i >> 1] = arr[i], arr2[i >> 1] = arr[i + 1];}  //奇偶分组
    fft(arr1, siz >> 1, tp); fft(arr2, siz >> 1, tp);  //分治做 FFT

    Complex Wn = Complex(cos(2.0 * PI / siz), sin(2.0 * PI / siz) * tp), w = Complex(1, 0);  //用 Wn 代表单位根
    for(reg int i(0); i < (siz >> 1); ++i) {
        arr[i] = arr1[i] + w * arr2[i]; arr[i + (siz >> 1)] = arr1[i] - w * arr2[i];  //一步内转移出两个点值
        w = w * Wn;
    }
}

int main(void) {
    n = read(), m = read();

    while(siz <= n + m) siz <<= 1;  //一定要让多项式的最高次数可以表示为 2^{k}-1 的形式

    for(reg int i(0); i <= n; ++i) f[i].Re = read();
    for(reg int j(0); j <= m; ++j) g[j].Re = read();

    fft(f, siz, 1); fft(g, siz, 1);  //分别对 F(x) 和 G(x) 做一遍 FFT

    for(reg int i(0); i <= siz; ++i) f[i] = f[i] * g[i];  //卷起来 QwQ
    fft(f, siz, -1);  //IFFT 还原回去

    for(reg int i(0); i <= n + m; ++i) printf("%d ", int(f[i].Re / siz + 0.5));  //不要忘了除以 n

    return 0;
}
//by CaO

蝴蝶变换

上面的递归版 FFT 或许非常快，但还不够快。

这是我在洛谷上提交 P3803【模板】多项式乘法（FFT） 的提交记录，我们可以发现，我们的代码在开 $\text O2$ 的情况下，最慢也已经跑到了 $1000\text{ms}$ 以上，虽然本题的时限开到了 $2000\text{ms}$，还是可以 AC 的，但是如果真的在实战中，很有可能会被卡常。

因为在递归的过程中，会有大量的函数入栈出栈的过程。众所周知，函数信息的入栈和出栈会造成大量的时间浪费，所以，把 FFT 函数设法变得非递归，可以产生非常有效的优化。

怎么样才能做到呢？

我们观察一下被转化为点值表示的多项式，这里以 $n=8$ 时为例：$f(x)=a_0+a_1x+a_2x^2+\cdots+a_7x^7$。

把它的系数写成一排，为：$[a_0\quad a_1\quad a_2\quad a_3\quad a_4\quad a_5\quad a_6\quad a_7]$；

第一次奇偶分组后，变为：$[a_0\quad a_2\quad a_4\quad a_6]\ \ [a_1\quad a_3\quad a_5\quad a_7]$；

第二次奇偶分组后，变为：$[a_0\quad a_4]\ \ [a_2\quad a_6]\ \ [a_1\quad a_5]\ \ [a_3\quad a_7]$；

第三次奇偶分组后，变为：$[a_0]\ \ [a_4]\ \ [a_2]\ \ [a_6]\ \ [a_1]\ \ [a_5]\ \ [a_3]\ \ [a_7]$。

我们再观察一下分组之后，每个系数的下标的规律：

分组	第 $0$ 项下标	第 $1$ 项下标	第 $2$ 项下标	第 $3$ 项下标	第 $4$ 项下标	第 $5$ 项下标	第 $6$ 项下标	第 $7$ 项下标
分组前	$(000)_2$	$(001)_2$	$(010)_2$	$(011)_2$	$(100)_2$	$(101)_2$	$(110)_2$	$(111)_2$
分组后	$(000)_2$	$(100)_2$	$(010)_2$	$(110)_2$	$(001)_2$	$(101)_2$	$(011)_2$	$(111)_2$

发现什么规律没？

没错！分组后，每个系数的下标均为分组前系数下标的二进制反转！这是更直观的图解。

『好漂亮！就跟蝴蝶的翅膀一样呢！』

所以，这样的变换，名字叫做蝴蝶变换。

我们可以事先预处理出所有下标的二进制反转，然后套一层循环就可以了。这样，就没必要进行递归了。

代码如下：

#include <bits/stdc++.h>
#define reg register  //继续卡常 QwQ
using namespace std;
const int MAXN = 5000001;
const double PI = acos(-1.0);

inline long long read(void) {  //卡常 QwQ
    long long s = 1, w = 0; char ch = getchar();

    while(ch < '0' || ch > '9') {if(ch == '-') s = -1; ch = getchar();}
    while(ch >= '0' && ch <= '9') {w = w * 10 + ch - '0'; ch = getchar();}

    return s * w;
}

int n, m, siz = 1;
int bit, rev[MAXN];

struct Complex {  //手动重载 Complex 类
    double Re, Im;
    Complex(double R = 0.0, double I = 0.0) {Re = R, Im = I;}
} f[MAXN], g[MAXN];

Complex operator+ (Complex a, Complex b) {return Complex(a.Re + b.Re, a.Im + b.Im);}
Complex operator- (Complex a, Complex b) {return Complex(a.Re - b.Re, a.Im - b.Im);}
Complex operator* (Complex a, Complex b) {return Complex(a.Re * b.Re - a.Im * b.Im, a.Re * b.Im + a.Im * b.Re);}

void fft(Complex* arr, int tp) {  //tp 代表变换的类型，如果 tp = 1 代表 FFT，tp = -1 代表 IFFT
    for(reg int i(0); i < siz; ++i) {
        if(i < rev[i]) swap(arr[i], arr[rev[i]]);  //蝴蝶变换
    }

    for(reg int mid(1); mid < siz; mid <<= 1) {  //枚举每次合并区间的半长度
        Complex Wn = Complex(cos(PI / mid), tp * sin(PI / mid));

        for(reg int i(0); i < siz; i += (mid << 1)) {  //枚举每个区间的左端点
            Complex w = Complex(1, 0);
            for(reg int j(0); j < mid; ++j) {  //FFT 代值
                Complex A1 = arr[i + j], A2 = w * arr[i + j + mid];
                arr[i + j] = A1 + A2, arr[i + j + mid] = A1 - A2;
                w = w * Wn;
            }
        }
    }
}

int main(void) {
    n = read(), m = read();

    while(siz <= n + m) siz <<= 1, ++bit;  //一定要让多项式的最高次数可以表示为 2^{k}-1 的形式
    for(reg int i(0); i < siz; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));  //预处理二进制反转，思考一下为什么这样能够行得通 QwQ

    for(reg int i(0); i <= n; ++i) f[i].Re = read();
    for(reg int j(0); j <= m; ++j) g[j].Re = read();

    fft(f, 1); fft(g, 1);  //分别对 F(x) 和 G(x) 做一遍 FFT

    for(reg int i(0); i <= siz; ++i) f[i] = f[i] * g[i];  //卷起来 QwQ
    fft(f, -1);  //IFFT 还原回去

    for(reg int i(0); i <= n + m; ++i) printf("%d ", int(f[i].Re / siz + 0.5));  //不要忘了除以 n

    return 0;
}
//by CaO

这份代码在洛谷的模板题上能够达到最慢 $500\text{ms}$ 的好成绩。

应用

FFT 的用途可多了去了 QwQ，几乎所有的学习多项式的 OIer，都是从 FFT 开始学习的。

当然，FFT 可能会因为卡精度而 WA 掉，所以从 FFT 又衍生出了快速数论变换（Number Theory Transform，NTT）。然后又衍生出了类如拆系数 FFT，任意模数 NTT（MTT）等黑科技，虽然我都不会 QwQ。

当然也可以用来加速高精度乘法 QwQ。

除此以外，FFT 在真正的计算机工程中也有大量的应用，从神经网络，到信号处理，应用非常广泛。

因此，学习 FFT 是一件很有用的事！

例题

本题目列表会持续更新。

• [P3803] 【模板】多项式乘法（FFT）
• [P1919] 【模板】A*B Problem 升级版（FFT 快速傅里叶变换）