快速傅里叶变换（FFT）

本文主要介绍快速傅里叶变换在快速计算多项式乘法上的应用。

概述

快速傅里叶变换（fast Fourier transform，简称 FFT），即利用计算机计算离散傅里叶变换（DFT）的高效、快速计算方法的统称。比如，它可以以 $O(n \log n)$ 的复杂度计算多项式乘法。

为了方便，下面定义一些符号：

• 用 $\mathrm{F}(x)$ 表示一个多项式，叫做“多项式 $\mathrm{F}$”。
• 设 $n$ 次多项式 $\mathrm{F}$ 的系数为数列 $F$，则 $\mathrm{F}(x)=\sum\limits_{i=0}^nF_ix^i$。
• 多项式的乘法：设 $n$ 次多项式 $\mathrm{A}$，$m$ 次多项式 $\mathrm{B}$，$n+m$ 次多项式 $\mathrm{C}$，$\mathrm{C}=\mathrm{A}\mathrm{B}$。则 $\forall 0 \le i \le n+m,C_i=\sum\limits_{j+k=i}A_jB_k$。

模板题：Luogu P3803。

如果直接计算卷积，时间复杂度为 $O(nm)$，而通过 FFT 可以优化计算过程。

DFT 与 IDFT 的思想

FFT 的算法流程：把系数表达转换为点值表达 -> 点值表达相乘 -> 把点值表达转换为系数表达。其中“把系数表达转换为点值表达”的算法叫做 DFT，“把点值表达转换为系数表达”的算法叫做 IDFT。

单位根 $\omega$ 及其性质

前置知识：复数及其运算和几何意义。

单位根（unit root）：$n(n \in \mathbb{N}^+)$ 次单位根是 $n$ 次幂为 $1$ 的复数。换句话说，$n$ 次单位根就是方程 $x^n=1$ 的复数解。

如上图，在复平面内有一个单位圆。在单位圆上的所有复数的模长均为 $1$。

首先考虑复数的模长：考虑一个复数 $x$ 使得 $x^n=1$，则 $|x|^n=|1|=1$，所以 $|x|=1$，即 $x$ 在单位圆上。于是我们的探究范围缩小到了单位圆上。

然后考虑复数的幅角：令单位根 $x$ 的幅角为 $\mathrm{ang}(x)$，则 $\mathrm{ang}(x^n)=n\mathrm{ang}(x)=0$，即 $n$ 次单位根 $n$ 等分单位圆。我们将单位根沿着单位圆逆时针标号，记 $\omega_n^k(k \in \mathbb{Z},0 \le k<n)$ 为第 $k-1$ 个 $n$ 次单位根。

比如当 $n=4$ 时：

备注：当 $k \notin [0,n)$ 时，$\omega_n^k=\omega_n^{k \bmod n}$

单位根的一些性质：

1. $\omega_n^0=1$
2. $(\omega_n^1)^k=\omega_n^k$（每乘上一个 $\omega_n^1$，相当于把幅角增加 $\frac{1}{n}$ 个圆周。乘上 $k$ 次，相当于总共把幅角增加 $\frac{k}{n}$ 个圆周，即 $\omega_n^k$。）
3. $\omega_n^i\omega_n^j=\omega_n^{i+j}$（$\omega_n^i\omega_n^j=(\omega_n^1)^i(\omega_n^1)^j=(\omega_n^1)^{i+j}=\omega_n^{i+j}$）
4. $\omega_{2n}^{2k}=\omega_n^k$，同理 $\forall p \in \mathbb{N}^+,\omega_{pn}^{pk}=\omega_n^k$（$\frac{pk}{pn}=\frac{k}{n}$）
5. $\forall n \in \mathbb{N}^+,n \equiv 0 \pmod{2},\omega_n^{k+\frac{n}{2}}=-\omega_n^k$（转动 $\frac{\frac{n}{2}}{n}=\frac{1}{2}$ 个圆周，等同于关于原点对称。）

DFT 的理论

考虑 $n-1(n=2^k,k \in \mathbb{N}^+)$ 次（$n$ 项）多项式 $\mathrm{F}$，$\mathrm{F}(x)=\sum\limits_{i=0}^{n-1}F_ix^i=\sum\limits_{i=0}^{\frac{n}{2}-1}F_{2i}x^{2i}+\sum\limits_{i=0}^{\frac{n}{2}-1}F_{2i+1}x^{2i+1}$。

又考虑两个 $\frac{n}{2}-1$ 次（$\frac{n}{2}$ 项）多项式 $\mathrm{L},\mathrm{R}$，$\mathrm{L}(x)=\sum\limits_{i=0}^{\frac{n}{2}-1}F_{2i}x^i$，$\mathrm{R}(x)=\sum\limits_{i=0}^{\frac{n}{2}-1}F_{2i+1}x^i$，则 $\mathrm{F}(x)=\mathrm{L}(x^2)+x\mathrm{R}(x^2)$。

将单位根 $\omega_n^k,\omega_n^{\frac{n}{2}+k}(k<\frac{n}{2})$ 带入 $\mathrm{F}(x)$，得到：

$$\begin{aligned} \mathrm{F}(\omega_n^k) &=\mathrm{L}((\omega_n^k)^2)+\omega_n^k\mathrm{R}((\omega_n^k)^2) \\ &= \mathrm{L}(\omega_{\frac{n}{2}}^k)+\omega_n^k\mathrm{R}(\omega_{\frac{n}{2}}^k) \end{aligned}$$

$$\begin{aligned} \mathrm{F}(\omega_n^{\frac{n}{2}+k}) &= \mathrm{L}((\omega_n^{\frac{n}{2}+k})^2)+\omega_n^{\frac{n}{2}+k}\mathrm{R}((\omega_n^{\frac{n}{2}+k})^2) \\ &= \mathrm{L}(\omega_n^{n+2k})-\omega_{n}^k\mathrm{R}(\omega_n^{n+2k}) \\ &= \mathrm{L}(\omega_n^{2k})-\omega_{n}^k\mathrm{R}(\omega_n^{2k}) \\ &= \mathrm{L}(\omega_{\frac{n}{2}}^k)-\omega_{n}^k\mathrm{R}(\omega_{\frac{n}{2}}^k) \end{aligned}$$

比对后发现，$\mathrm{F}(\omega_n^k)$ 和 $\mathrm{F}(\omega_n^{\frac{n}{2}+k})$ 只有一个加减号的区别。

$\mathrm{L}$ 与 $\mathrm{R}$ 均为 $\mathrm{F}$ 的子问题，性质完全相同，可以一直分治下去，每一层的时间复杂度为 $O(n)$，共分治 $\log_2 n$ 层，总时间复杂度为 $O(n \log n)$。

DFT 的实现

复数运算

C++ 自带的 complex 常数大，一般手写复数运算（此处会用到加法、减法和乘法）。

struct complex{
	double x,y;
	complex(double a=0,double b=0){x=a,y=b;}
};
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);}

单位根的计算

tG 表示 $\omega_{len}^1$，buf 表示 $\omega_{len}^0$。

C++Copy
complex tG(cos(2*pi/len),sin(2*pi/len)),buf(1,0);

DFT 的实现

C++Copy
void dft(complex*f,int len){
	if(len==1)return;
	complex*fl=f,*fr=fl+(len>>1);
	for(int k=0;k<len;k++)h[k]=f[k];
	for(int k=0;k<len>>1;k++)fl[k]=h[k<<1],fr[k]=h[k<<1|1];
	dft(fl,len>>1);
	dft(fr,len>>1);
	complex tG(cos(2*pi/len),sin(2*pi/len)),buf(1,0);
	for(int k=0;k<len>>1;k++){
		h[k]=fl[k]+buf*fr[k];
		h[k+(len>>1)]=fl[k]-buf*fr[k];
		buf=buf*tG;
	}
	for(int k=0;k<len;k++)f[k]=h[k];
}

IDFT 的理论

设经过 DFT 处理后得到的点集 $G=\mathrm{DFT}(\mathrm{F}(x))$，则 $G_k=\sum\limits_{i=0}^{n-1}(F_i(\omega_n^k)^i)$。

结论为 $nF_k=\sum\limits_{i=0}^{n-1}(G_i(\omega_n^{-k})^i)$，证明如下：

$$\begin{aligned} \sum\limits_{i=0}^{n-1}(G_i(\omega_n^{-k})^i) &= \sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}(F_j(\omega_n^i)^j)(\omega_n^{-k})^i) \\ &= \sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}(F_j\omega_n^{ij}\omega_n^{-ik})) \\ &= \sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}(F_j\omega_n^{i(j-k)})) \\ &= \sum\limits_{j=0}^{n-1}(\sum\limits_{i=0}^{n-1}(F_j\omega_n^{i(j-k)})) \end{aligned}$$

分类讨论：

当 $j=k$ 时，$\sum\limits_{i=0}^{n-1}(F_j\omega_n^{i(j-k)})=\sum\limits_{i=0}^{n-1}F_j=nF_k$。

当 $j \ne k$ 时，设 $p=j-k$，则 $\sum\limits_{i=0}^{n-1}(F_j\omega_n^{i(j-k)})=\sum\limits_{i=0}^{n-1}(F_j\omega_n^{ip})=\sum\limits_{i=0}^{n-1}(\omega_n^p)^iF_{k+p}=\dfrac{(\omega_n^p)^n-1}{\omega_n^p-1}F_{k+p}=0F_{k+p}=0$。

证毕。

FFT 的初步实现

提交记录 备份

FFT 的代码优化

将 DFT 和 IDFT 合并成一个函数

提交记录 备份

预处理 $\sin$ 值和 $\cos$ 值

提交记录 备份

二进制反转优化

可以发现，每次将序列按照下标的奇偶性分成两部分，下标二进制表示下末尾为 $0$ 的排在前，末尾为 $1$ 的排在后。说白了，就是将数的二进制的最低 $k$ 位进行一个翻转，再排序。这个可以通过递推 $O(n)$ 求出。

迭代实现

提交记录 备份

例题