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

• 本文总计约 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 快速傅里叶变换）