这篇文章写得很烂，而且大部分来自于算导，如有感兴趣者可以直接阅读算导或参考资料QAQ！本文省略大量证明，请各位谨慎食用QAQ！

概述

快速傅里叶变换（FFT）是一种快速计算序列的离散傅里叶变换（DFT）及其逆变换（IDFT）的算法，其复杂度为O(n \log n)，比朴素的DFT计算更优。OI中常用它来计算多项式乘法等需要应用卷积的地方。
除了本文章中介绍的数学知识，你还需要知道：

• 数学笔记：复数 | KSkun’s Blog
• 数学笔记：矩阵、行列式 | KSkun’s Blog

多项式的表达法

系数表达

一个n次多项式可以写成A(x) = \sum_{i=0}^{n} a_ix^i的形式，每一项对应有系数，我们把这些系数组成的向量a = (a_0, a_1, \ldots, a_n)叫做这个多项式的系数表达。

点值表达

一个n次多项式A(x)的点值表达是一个由n+1个点值对组成的集合\{ (x_0, y_0), (x_1, y_1), \ldots, (x_n, y_n) \}，其中，x值互不相同，且y_i = A(x_i)。
对于n+1个点值能表达唯一的n次多项式的证明，参见[学习笔记] 多项式与快速傅里叶变换(FFT)基础 – rvalue – 博客园。
有了点值表达，我们来想想多项式相乘的结果：A(x) \cdot B(x)应该怎么计算。假设A多项式的次数是n，B的次数是m，那么得到的结果一定是一个n+m次多项式，因此我们需要n+m+1个点值才能确定它。因此我们需要把A和B的点值表示扩展为n+m+1个点值，然后对它们的点值表示中的y值对应相乘，结果就是新多项式的点值表示。

单位根的使用与多项式乘法的步骤

我们知道，朴素进行多项式乘法需要O(n^2)时间。如果我们能快速进行系数转点值（DFT）及其逆运算，通过点值表达做乘法只需要O(n)的时间，就能够实现加速这一过程的目的。我们可以通过选取精妙的点值来实现这一加速。
还记得复数中的单位根吗？它有一个折半引理，也就证明了 (\omega_n^k)^2 = \omega_{\frac{n}{2}}^k 。这个结论会在FFT中用到。这决定了我们可以分治实现O(n \log n)的系数转点值及其逆运算。这种加速后的算法被称为FFT。
那么我们计算多项式乘法步骤就是：

1. 系数转点值，一次FFT，O(n \log n)
2. 计算乘积的点值，O(n)
3. 点值转系数，一次FFT的逆运算，O(n \log n)

离散傅里叶变换（DFT）

离散傅里叶变换是一种完成系数表达转换为点值表达的变换。这一变换是对多项式A在n个n次单位根处求点值，也就是说，是在完成下面这样一个事情：
多项式的系数表达给出为a = (a_0, a_1, \ldots, a_n)，求y = (y_0, y_1, \ldots, y_n)，其中y_i = A(\omega_{n+1}^i)。
但是这样看起来其实也是一个O(n^2)的算法啊，教练加速去哪里了？往下看↓

快速傅里叶变换（FFT）

下面的叙述中，n代表多项式的次数界，次数界是指严格大于一个多项式的次数的正整数。为了方便，我们默认n是2的某次幂。
上面我们提到了要把这个问题分治解决，那么就要构造出适合分治解决的形式。首先是怎么构造分治的分支。我们考虑把系数按奇偶性分开，构造下面这样两个多项式
\begin{aligned} A_1(x) &= a_0 + a_2x + a_4x^2 + \cdots \\ A_2(x) &= a_1 + a_3x + a_5x^2 + \cdots \end{aligned}
那么就有A(x) = A_1(x^2) + A_2(x^2)x。
那么我们考虑那个折半引理的式子 (\omega_n^k)^2 = \omega_{\frac{n}{2}}^k ，要求的一堆平方实际上变成了\frac{n}{2}个数字，那么我们就分别计算这\frac{n}{2}个数处的点值，然后合并到当前点值上就可以了。这样我们成功实现了缩小问题规模的目的。
具体而言，我们需要对\frac{n-1}{2}之前和之后的点值分别讨论
y_k = A(\omega_n^k) = A_1(\omega_n^{2k}) + \omega_n^kA_2(\omega_n^{2k}) = y^{[1]}_k + \omega_n^ky^{[2]}_k \ (k = 0, 1, \ldots, \frac{n}{2} - 1)
y_{k+\frac{n}{2}} = A(\omega_n^{k+\frac{n}{2}}) = A_1(\omega_n^{2k+n}) + \omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n}) = A_1(\omega_n^{2k}) - \omega_n^kA_2(\omega_n^{2k}) = y^{[1]}_k - \omega_n^ky^{[2]}_k \ (k = 0, 1, \ldots, \frac{n}{2} - 1)
对于递归边界，因为\omega_1^0 = 1，所以y_0 = a_0\omega_1^0 = a_0，直接不做处理即可。

逆变换：点值转系数（IDFT）

我们发现我们可以用下面这张矩阵来表示DFT的过程
\begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_{n-1} \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n & \omega_n^2 & \cdots & \omega_n^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)^2} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix}
其实中间那个大块头矩阵就是一个范德蒙矩阵V_n，实际上我们还有一个定理来求这个矩阵的逆矩阵： (V_n^{-1})_{j, k} = \frac{\omega_n^{-kj}}{n} 。那么我们就可以把逆矩阵往点值上一乘，点值就乖乖地变回系数啦~具体而言，变化是这样的a_j = \frac{1}{n} \sum_{k=0}^{n-1} y_k\omega_n^{-kj}。我们发现这个玩意和FFT在做的事情很相似，差就差在一个除以n还有单位根指数上面的负号。负号我们可以通过函数参数传进去，除以n可以在函数外部做，那么直接用FFT的过程不就好了。

FFT的优化方式

蝴蝶操作

我们发现上面的代码中\omega_n^ky_k^{[2]}的计算似乎是重复的，我们可以临时存起来而不是计算两次。

迭代实现

我们来观察FFT的调用树：

（图片来自《算法导论（第3版）》，参见参考资料）
我们发现我们可以从最下层逐对合并结果，而需要把最下层的顺序处理出来。通过观察我们发现神奇的事：最下层第i个位置的下标刚好是它的二进制倒过来。那么交换一个位置和它二进制倒过来位置的元素即可。二进制倒过来的数可以预处理好。

不要用STL

STL比手写complex慢一些。

代码

本代码可以通过【P3803】【模板】多项式乘法（FFT） – 洛谷一题。

// Code by KSkun, 2018/3
#include <cstdio>
#include <cmath>

#include <algorithm>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    char c = fgc();
    while(c < '0' || c > '9') {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(c >= '0' && c <= '9') {
        res = res * 10 + c - '0';
        c = fgc();
    }
    return res * neg;
}

const int MAXN = 1 << 22;
const double PI = acos(-1);

struct Complex {
    double real, imag;
    Complex(double real = 0, double imag = 0): real(real), imag(imag) {}
    inline Complex operator+(const Complex &rhs) const {
        return Complex(real + rhs.real, imag + rhs.imag);
    }
    inline Complex operator-(const Complex &rhs) const {
        return Complex(real - rhs.real, imag - rhs.imag);
    }
    inline Complex operator*(const Complex &rhs) const {
        return Complex(real * rhs.real - imag * rhs.imag, real * rhs.imag + imag * rhs.real);
    }
    inline Complex& operator*=(const Complex &rhs) {
        Complex old = *this;
        real = old.real * rhs.real - old.imag * rhs.imag;
        imag = old.real * rhs.imag + old.imag * rhs.real;
        return *this;
    }
};

int n, m, len, rev[MAXN];
Complex a[MAXN], b[MAXN];

inline void fft(Complex *arr, int f) {
    for(int i = 0; i < n; i++) {
        if(i < rev[i]) std::swap(arr[i], arr[rev[i]]);
    }
    for(int i = 1; i < n; i <<= 1) {
        Complex wn(cos(PI / i), f * sin(PI / i));
        for(int j = 0; j < n; j += i << 1) {
            Complex w(1, 0);
            for(int k = 0; k < i; k++) {
                Complex x = arr[j + k], y = w * arr[j + k + i];
                arr[j + k] = x + y;
                arr[j + k + i] = x - y;
                w *= wn;
            }
        }
    }
}

int main() {
    n = readint(); m = readint();
    for(int i = 0; i <= n; i++) {
        a[i].real = readint();
    }
    for(int i = 0; i <= m; i++) {
        b[i].real = readint();
    }
    m += n;
    for(n = 1; n <= m; n <<= 1) len++;
    for(int i = 0; i < n; i++) {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
    }
    fft(a, 1);
    fft(b, 1);
    for(int i = 0; i <= n; i++) {
        a[i] *= b[i];
    }
    fft(a, -1);
    for(int i = 0; i <= m; i++) {
        printf("%d ", int(a[i].real / n + 0.5));
    }
    return 0;
}

参考资料

• [学习笔记] 多项式与快速傅里叶变换(FFT)基础 – rvalue – 博客园
• 《算法导论（第3版）》，机械工业出版社，9787111407010