快速傅里叶变换(FFT)

OI-算法

---

部分内容参考文献
本文作者$\text{dxbdly}$，引用或转载请注明原文链接

1. 前言

$\text{FFT}$ 是一种用来实现 $\text{DFT}$（离散傅里叶变换）的高效算法。

它一般用于加速多项式乘法，可使其达到 $O(n\log n)$ 的速度，有时也用于加速大整数乘法。

注：$FFT$ 算法涉及较多前置数学知识，本文仅会介绍重要的前置知识，其他不清楚的请读者自行百度。

---

2. 前置知识1：多项式表示法

$1^o$ 系数表示法

通常，我们表示多项式的方式是系数表示法，即使用一个系数数列来表示多项式。

也就是形如：

$$f(x) = a_0 + a_1x + a_2x^2+\ldots+a_nx^n$$
即：
$$f(x) = \{a_0,a_1,a_2\ldots,a_n\}$$

但如果依照这种表示方式进行多项式乘法，乘积的系数数列会十分复杂，难以研究。

所以我们考虑使用其他方式表示多项式。

$2^o$ 点值表示法

我们把多项式看成一个函数，由于我们可以使用 $n + 1$ 个此函数上的点，来唯一表示此函数。

也就是说，这 $n + 1$ 个点同样可以唯一表示此多项式。

也就是当：

$$f(x_0) = a_0 + a_1x_0 + a_2x_0^2 + \ldots + a_nx_0^n\\ f(x_1) = a_0 + a_1x_1 + a_2x_1^2 + \ldots + a_nx_1^n\\ f(x_2) = a_0 + a_1x_2 + a_2x_2^2 + \ldots + a_nx_2^n\\ \vdots\\ f(x_n) = a_0 + a_1x_n + a_2x_n^2 + \ldots + a_nx_n^n$$

就有：

$$f(x) = \{(x_0, f(x_0)), (x_1, f(x_1)), (x_2, f(x_2)), \ldots , (x_n, f(x_n))\}$$

这样表示有何好处呢？

我们会发现，在点值表示法下，多项式的乘积将很好表示。

设：有两多项式
$f(x) = \{(x_0, f(x_0)), (x_1, f(x_1)), \ldots ,(x_n, f(x_n))\}$ 与
$g(x) = \{(x_0, g(x_0)), (x_1, g(x_1)), \ldots ,(x_n, g(x_n))\}$

那么乘积 $F(x)$ 可表示为：

$$F(x) = \{(x_0, f(x_0) \times g(x_0)), (x_1, f(x_1) \times g(x_1)), \ldots , (x_n, f(x_n)\times g(x_n))\}$$

也就是说，如果我们能将 $f(x)$ 与 $g(x)$ 表示为点值表示法，那么就能用 $O(n)$ 的时间求出 $F(x)$ 的点值表示。

那么留给我们的就是两个问题：

如何将 $f(x)$ 与 $g(x)$ 从系数表示变为点值表示
如何将 $F(x)$ 从点值表示变回系数表示

---

3. 前置知识2：单位复根

我们发现，在解决这两个问题之前，我们需要先找出 $n + 1$ 个点

并且它们适于我们进行转换，换句话说，我们希望这 $n + 1$ 个点满足某种规律，且对所有多项式都适用。

在这里，我们需要引入一个概念—— 单位复根

$1^o$ 单位复根的定义

PS：不知道复数是什么的读者请看。

下面附上手写的复数类型，包括复数的三种运算：

//The code is from dxbdly
struct Complex {
    double x, y;
    Complex(double x0 = 0, double y0 = 0) { x = x0, y = y0; }
}g1[maxn], g2[maxn], sum[maxn];
Complex operator + (Complex a, Complex b) { return Complex(a.x + b.x, a.y + b.y); }
Complex operator - (Complex a, Complex b) { return Complex(a.x - b.x, a.y - b.y); }
Complex operator * (Complex a, Complex b) { return Complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); }

我们称 $w^n = 1$ 在复数意义下的解为 $n$ 次复根。

显然，这样的复根有 $n$ 个，我们记它们为 $\{w_n^k | k = 0,1,2,\ldots,n - 1\}$

我们考虑这个方程的几何意义：

在复平面上作一个单位圆，则此单位圆上的复数向量模长都为 $1$
所以我们可以将 $w$ 理解成复平面上的一个单位向量，它的俯角的 $n$ 倍正好是 $360^o$
$w_n$ 就是 单位复根，且这 $n$ 次复根将这个圆周 $n$ 等分。

所以，我们可以利用俯角求出 $w_n$ ：

$$w_n = \cos(\frac{2\pi}{n}) + i \times \sin(\frac{2\pi}{n})$$

$2^o$ 单位复根的性质

单位复根有三个重要的性质：

对于任意正整数 $n$ 与整数 $k$，有：

• $w_n^n = 1$
  解释：相当于转一周回到原点。
• $w_n^k = w_{2n}^{2k}$
  解释：在 $n$ 份中取 $k$ 份与 在 $2n$ 中取 $2k$ 份显然相同
• $w_{2n}^{k + n} = -w_{2n}^{k}$
  解释：相当于转半圈，坐标取反
• $w_{n}^{a} + w_{n}^{b} = w_{n}^{a+b}$
  解释：分别取 $a$ 份与 $b$ 份显然就是共取 $a+b$ 份

后续的 $\text{DFT}$ 推导过程中，需要反复运用这三个式子，请好好理解。

---

4. $\text{DFT}$

确定了要带入的 $n + 1$ 个点后，我们先来解决系数表示变点值表示的问题。

这个过程称为 $\text{DFT}$ 。

如果直接带入点，复杂度仍是 $O(n ^ 2)$ 的，不能接受。

由于我们取的点是存在规律的，我们考虑能否一次性将所有点的函数值求出。

$1^o$ 分治 $\text{DFT}$

很自然的想到了分治。

这里积累一种套路：如果要对多项式分治，可以考虑尝试将奇数项与偶数项分开，构造子问题。

举个例子：

设：

$$F(x) = a_0 + a_1x + a_2x^2 + \ldots + a_7x^7$$

奇偶分组：

$$F(x) = (a_0 + a_2x^2 + a_4x^4 + a_6x^6) + (a_1x + a_3x^3 + a_5x^5 + a_7x^7)\\ \quad\quad=(a_0 + a_2x^2 + a_4x^4 + a_6x^6) + x(a_1 + a_3x^2 + a_5x^4 + a_7x^6)$$

可以发现两部分都是 $F(x)$ 的子结构。

设：

$$G(x) = a_0 + a_2x + a_4x^2 + a_6x^3\\ H(x) = a_1 + a_3x + a_5x^2 + a_7x^3$$

用 $G(x)$ 与 $H(x)$ 表示 $F(x)$ ：

$$F(x) = G(x^2) + xH(x^2)$$

接下来，考虑带入点 $w_n^k$ ：

$$F(w_n^k) = G((w_n^k)^2) + w_n^kH((w_n^k)^2)\\ \quad= G(w_n^{2k}) + w_n^kH(w_n^{2k})\\ \quad\quad= G(w_{n/2}^k) + w_n^kH(w_{n/2}^k)$$

由于要在 $O(n\log n)$ 的时间内，所以每次需要将求解的范围缩小。

考虑单位复根性质 $3$ ，则可尝试表示 $F(w_n^{k+2/n})$

$$F(w_n^{k+n/2}) = G(w_n^{2k+n}) + w_n^{k+n/2}H(w_n^{2k+n})\\ \quad\quad\;\;= G(w_n^{2k}) + w_n^{k+n/2}H(w_n^{2k})\\ \quad= G(w_n^{2k}) - w_n^kH(w_n^{2k})\\ \quad\quad\quad= G(w_{n/2}^{k}) - w_n^kH(w_{n/2}^{k})$$

这样，我们就可以用 $G(w_{n/2}^k)$ 与 $H(w_{n/2}^k)$ 同时算出两个值，保证了 $O(n\log n)$ 的复杂度。

$2^o$ 位逆序置换(蝴蝶变换)

上面的分治算法虽以达到 $O(n \log n)$ 的复杂度。

可实际上由于要使用递归，并且需传输大量的数组，指针，实现起来常数巨大。

我们考虑将其用迭代实现，这样，我们就可以从最后的状态倍增的向上合并。

也就是说，我们首先要想办法求出最后分组的状态，也就是确定合并的顺序。

举个例子：

原数列：$\{x_0,x_1,x_2,\ldots,x_8\}$

最终数列：$\{x_0\}\{x_4\}\{x_2\}\{x_6\}\{x_1\}\{x_5\}\{x_3\}\{x_7\}$

而这个序列是有规律的：

对于原来的序列，将每个数用二进制表示，然后将二进制翻转，就是对应的最终位置的下标
我们称这种翻转为位逆序置换，又称蝴蝶变换

根据定义，很容易想到用 $O(n \log n)$ 的方式求解每个数变换后的结果。

但是这还不够，实际上，我们可以通过递推求出 $1$ 到 $n$ 的位逆序置换值。

设：$R(x)$ 表示 $x$ 在 $1$ 到 $n$ 中的位逆序置换值

边界： $R(0) = 0$

考虑转移，将 $x$ 的翻转分为其个位的翻转与个位前的翻转，即：$x \bmod 2$ 的翻转与 $\lfloor\frac{x}{2}\rfloor$ 的翻转

则有：

$$R(x) = \lfloor\frac{R(\lfloor\frac{x}{2}\rfloor)}{2}\rfloor + (x \bmod 2) \times \frac{n}{2}$$

注意：

$\lfloor\frac{R(\lfloor\frac{x}{2}\rfloor)}{2}\rfloor$ 是因为去掉个位反转后，首位补的 $0$ 会翻转到个位，要去掉。
$(x \bmod 2) \times \frac{n}{2}$ 则是如果个位为 $1$ ，那么首位为 $1$ ，否则补 $0$ 。

位逆序置换部分代码：

//The code is from dxbdly
inline void Change(Complex f[], int len) {
    for(register int i = 0; i < len; ++i) {
        rev[i] = rev[i >> 1] >> 1;
        if(i & 1) rev[i] |= (len >> 1);
    }
    for(register int i = 0; i < len; ++i)
        if(i < rev[i]) swap(f[i], f[rev[i]]);
}

这样就可以用 $O(n)$ 的时间预处理，然后按照 $\text{DFT}$ 的过程将其合并就好，具体实现后面会讲。

---

5. $\text{IDFT}$

$\text{IDFT}$，即傅里叶反变换算法，它可以解决将点值表示变回系数表示的问题。

它与 $\text{DFT}$ 算法极为类似，实现时一般也写在一起（后面细讲）。

$1^o$ 线性代数意义下的 $\text{DFT}$ 与 $\text{IDFT}$

可以发现，$DFT$ 本质是一个线性变换，我们尝试将这个过程用矩阵表示出来：

$$\left[ \begin{matrix} y_0\\ y_1\\ y_2\\ \vdots\\ y_n \end{matrix} \right] = \left[ \begin{matrix} 1 & 1 & 1 & \ldots & 1\\ 1 & w_n^1 & w_n^2 & \ldots & w_n^{n - 1} \\ 1 & w_n^2 & w_n^4 & \ldots & w_n^{2(n - 1)} \\ \vdots & \vdots & \vdots &\ddots & \vdots \\ 1 & w_n^{n - 1} & w_n^{2(n - 1)} & \ldots & w_n^{(n - 1) ^ 2} \end{matrix} \right] \left[ \begin{matrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n - 1} \end{matrix} \right]$$

现在我们相当于知道 $y$ 点值矩阵 与 $w_n^i$ 点矩阵，要求 $a$ 系数矩阵。

同时乘上 $w_n^i$ 的逆矩阵，也就是这样：

$$\left[ \begin{matrix} y_0\\ y_1\\ y_2\\ \vdots\\ y_n \end{matrix} \right] \left[ \begin{matrix} \frac{1}{n} & \frac{1}{n} & \frac{1}{n} & \ldots & \frac{1}{n}\\ \frac{1}{n} & \frac{1}{w_n^1} & \frac{1}{w_n^2} & \ldots & \frac{1}{w_n^{n - 1}} \\ \frac{1}{n} & \frac{1}{w_n^2} & \frac{1}{w_n^4} & \ldots & \frac{1}{w_n^{2(n - 1)}} \\ \vdots & \vdots & \vdots &\ddots & \vdots \\ \frac{1}{n} & \frac{1}{w_n^{n - 1}} & \frac{1}{w_n^{2(n - 1)}} & \ldots & \frac{1}{w_n^{(n - 1) ^ 2}} \end{matrix} \right] = \left[ \begin{matrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n - 1} \end{matrix} \right]$$

这就是 $\text{IDFT}$ 的矩阵表示，可以发现，这是个类似 $\text{DFT}$ 的矩阵。

所以，我们只需要将所有 $w_n^k$ 改为 $\frac{1}{w_n^k}$ 就好了。

而根据定义：

$$\frac{1}{w_n^k} = \cos(\frac{2\pi}{k}) + i\times\sin(-\frac{2\pi}{k})$$

所以，我们只要在 $\text{DFT}$ 中带一个标记，就可以实现 $\text{IDFT}$ 的功能。

$2^o$ $\text{DFT}$ 与 $\text{IDFT}$ 的实现

自己写的时候容易晕，注解一下各个循环在枚举什么，先上代码：

//The code is from dxbdly
inline void DFT(Complex f[], int len, int op) {
    Change(f, len);
    for(register int h = 2; h <= len; h <<= 1) {
        Complex step(cos(2 * Pi / h), op * sin(2 * Pi / h));
        for(register int i = 0; i < len; i += h) {
            Complex w(1, 0);
            for(register int k = i; k < i + h / 2; ++k) {
                Complex G = f[k], H = w * f[k + h / 2];
                f[k] = G + H, f[k + h / 2] = G - H;
                w = w * step;
            }
        }
    }
    if(op == -1) {
        for(register int i = 0; i < len; ++i) f[i].x /= len;
    }
}

注意：
$op$ 为标志变量，用于区分 $\text{DFT}$ 与 $\text{IDFT}$。
第一层循环， $h$ 枚举的是当前合并层每一组的数量，其中 $step$ 为单位复数。
第二层循环， $i$ 枚举的是每一组的编号起点
第三层循环， $k$ 枚举的是当前要求解的下标（另一个是 $k + h / 2$）
注意这里用了滚动数组，并没有新记录 $G$ 函数与 $H$ 函数。

6. 例题

一道模板题：

模板题

代码为了图方便封装过了的。

代码：

//The code is from dxbdly
#include<bits/stdc++.h>
using namespace std;
inline int read() {
    int x = 0;
    char c = getchar();
    bool f = 0;
    while(!isdigit(c)) f |= (c == '-'), c = getchar();
    while(isdigit(c)) x = (x * 10) + (c ^ 48), c = getchar();
    return f ? -x : x;
}
namespace FFT {
    const double Pi = acos(-1);
    const int maxn = 3e6 + 5;
    int rev[maxn];
    struct Complex {
        double x, y;
        Complex(double x0 = 0, double y0 = 0) { x = x0, y = y0; }
    }g1[maxn], g2[maxn], sum[maxn];
    Complex operator + (Complex a, Complex b) { return Complex(a.x + b.x, a.y + b.y); }
    Complex operator - (Complex a, Complex b) { return Complex(a.x - b.x, a.y - b.y); }
    Complex operator * (Complex a, Complex b) { return Complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); }
    inline void Change(Complex f[], int len) {
        for(register int i = 0; i < len; ++i) {
            rev[i] = rev[i >> 1] >> 1;
            if(i & 1) rev[i] |= (len >> 1);
        }
        for(register int i = 0; i < len; ++i)
            if(i < rev[i]) swap(f[i], f[rev[i]]);
    }
    inline void DFT(Complex f[], int len, int op) {
        Change(f, len);
        for(register int h = 2; h <= len; h <<= 1) {
            Complex step(cos(2 * Pi / h), op * sin(2 * Pi / h));
            for(register int i = 0; i < len; i += h) {
                Complex w(1, 0);
                for(register int k = i; k < i + h / 2; ++k) {
                    Complex G = f[k], H = w * f[k + h / 2];
                    f[k] = G + H, f[k + h / 2] = G - H;
                    w = w * step;
                }
            }
        }
        if(op == -1) {
            for(register int i = 0; i < len; ++i) f[i].x /= len;
        }
    }
    inline void Work(int n, int m, int num1[], int num2[]) {
        int len = 1;
        while(len < n * 2 || len < m * 2) len <<= 1;
        for(register int i = 0; i < n; ++i) g1[i] = Complex(num1[i], 0);
        for(register int i = 0; i < m; ++i) g2[i] = Complex(num2[i], 0);
        for(register int i = n; i < len; ++i) g1[i] = Complex(0, 0);
        for(register int i = m; i < len; ++i) g2[i] = Complex(0, 0);
        DFT(g1, len, 1), DFT(g2, len, 1);
        for(register int i = 0; i < len; ++i) sum[i] = g1[i] * g2[i];
        DFT(sum, len, -1);
        for(register int i = 0; i < n + m - 1; ++i) printf("%d ", (int)(sum[i].x + 0.5));
    }
}
const int maxn = 1e6 + 5;
int n, m, num1[maxn], num2[maxn];
int main() {
    n = read() + 1, m = read() + 1;
    for(register int i = 0; i < n; ++i) num1[i] = read();
    for(register int i = 0; i < m; ++i) num2[i] = read();
    FFT::Work(n, m, num1, num2);
    return 0;
}

几个要注意的点：
$\text{FFT}$ 几乎所有的操作都要求在 $len$ 为 $2$ 的整数幂下进行，所以需要对多项式先补 $0$
$\text{STL}$ 的复数类型很慢，这里直接自己写的，三种运算，应该算基础吧。
最后关于精度，四舍五入输出，大部分时候 $\text{double}$ 够用。

7. $\text{The End}$

刚学的 $\text{FFT}$ ，没做什么题，文章有不对的地方请私信博主。

最后，$\text{FFT}$ 是一种很重要的多项式算法，之后的许多多项式有关的算法都可以用 $\text{FFT}$ 优化。

其中有许多算法思想更是值得我们学习。

解决 $\text{DFT}$ 问题的算法还有很多，例如 快速数论变换(NTT)，$\text{MTT}$（任意模数 $\text{NTT}$），$\text{FWT}$（快速沃尔什变换）

但这些算法都是以 $\text{FFT}$ 为基础衍生而来，所以读者需要将 $\text{FFT}$ 理解清楚后进行深入的学习

希望本文对大家有所帮助，感谢观看。

部分内容参考文献
本文作者$\text{dxbdly}$，引用或转载请注明原文链接

$$\text{Q.E.D}$$