【数学】快速傅里叶变换（FFT）

快速傅里叶变换（FFT）

FFT 是之前学的，现在过了比较久的时间，终于打算在回顾的时候系统地整理一篇笔记，有写错的部分请指出来啊 qwq。

卷积

卷积、旋积或褶积（英语：Convolution）是通过两个函数 $f$ 和 $g$​​ 生成第三个函数的一种数学算子。

定义

设 $f,g$​ 在 $R1$​ 上可积，那么 $h(x) = \int_{-∞}^∞f(\tau)g(x-\tau)d\tau$ 称为 $f$ 与 $g$​ 的卷积。

对于整系数多项式域，$n-1$ 次多项式 $A,B$ 的相乘可以得到 $h(x) = A(x)B(x) = \sum_{i=0}^{2n-2} \sum_{j+k=i}a_j b_k = \sum_{i=0}^{2n-2} \sum_{j=0}^i a_j b_{i-j} x^i$，对应的卷积为 $C_i = \sum A_j B_{i-j}$。

系数表示法

即用多项式各项系数来刻画这个多项式，例如 $n-1$ 次多项式就可以写成这样：$A(x) = a_0 + a_1x+...+a_{n-1}x^{n-1}$

点值表示法

我们知道，$n$​​​​ 个不同的点可以确定一个 $n-1$​​​​ 次的多项式，所以我们可以使用 $n$​​​​ 个（不同）点来刻画一个 $n-1$​​​​ 次多项式。

这样做会有什么方便呢？

例如 $f(x) = (x_0, f(x_0)),...(x_n,f(x_n)),g(x) = (x_0, g(x_0)),...(x_n, g(x_n))$，那么它们的卷积 $h(x) = (x_0, f(x_0)g(x_0)),...(x_n, f(x_n)g(x_n))$。

这意味着在系数表示法中需要 $O(n^2)$ 次的乘法运算在点值表示法中只需要 $O(n)$ 次。

系数表示法转点值表示法（DFT）

下面考虑如何将 $n-1$ 次多项式从系数表示法转为点值表示法。

因为用普通的方法选取 $n$​​​​​ 个点然后将系数表示法转为点值表示法的复杂度为 $O(n^2)$​​​​​（因为需要选 $n$​​​​​ 个点，然后对于每个点 $x$ 需要计算共 $n$ 项的结果），我们考虑如何优化这一步。

注意到满足 $w^n=1$​​ 的单位根 $w$​ 有 $n$​ 个，故从这里入手。

我们记方程 $w^n = 1$ 的第 $k$ 个单位根为 $w_n^k$​。

方便起见，设 $n$​ 为 $2$​ 的幂（就算不是也可以看作是高次项的系数为 $0$​）。

将 $A(x)$ 按照次数的奇偶性分别分成两组 $F(x),G(x)$，并表示为 $A(x) = F(x^2) + xG(x^2)$

例如 $A(x) = a_0+a_1x+a_2x^2+a_3x^3$​，那么 $F(x) = a_0+a_2x,G(x)=a_1+a_3x$​。

将 $x=w_n^k$​ 代入 $A(x)$​，由复数的性质，

$A(w_n^k) = F(w_{\frac{n}{2}}^k) + w_n^k G(w_{\frac{n}{2}}^k)$​ ，类似地 $A(w_{n}^{k+\frac{n}{2}}) = F(w_{\frac{n}{2}}^k) - w_n^k G(w_{\frac{n}{2}}^k)$​。​

推导：

$A(w_{n}^{k+\frac{n}{2}}) = F(w_n^{2k+n}) + w_n^{k+\frac{n}{2}}G(w_n^{2k+n}) \\= F(w_{\frac{n}{2}}^k) + w_n^{k+\frac{n}{2}} G(w_{\frac{n}{2}}^k) \\= F(w_{\frac{n}{2}}^k) - w_n^k G(w_{\frac{n}{2}}^k)$​​​

可以发现对于两个相应的单位根 $w_n^k,w_n^{\frac{n}{2}+k}$​，可以用对应的 $F,G$​ 算出（可以递归地实现这个过程），而且计算的范围折半，所以一共需要计算 $O(logN)$​ 层，每一层执行 $O(n)$​ 次运算，所以复杂度为 $O(NlogN)$​。

点值表示法转系数表示法（IDFT）

下面考虑如何将 $n-1$​ 次多项式从点值表示法转为系数表示法。

因为对于每个点值 $y_i = \sum_{k=0}^{n-1}w_n^{ki}$，其中 $i\in[0, n-1]$​​​​，我们可以写出等式：

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

现在我们已经有向量 $y$​​ 了（就是右式），因此，如果要得到向量 $a$​​，只需要两边乘上 $w$​​​ 矩阵的逆即可。

这里的 $w$​​ 矩阵正是著名的范德蒙矩阵，它的逆正好是每一项都取倒数，然后除以 $n$​​。

因此有 $a_i = \frac{1}{n}\sum_{k=0}^{n-1}w_n^{-ki}$，其中 $i\in[0, n-1]$​。

有没有发现 $a_i,y_i$​ 的形式非常接近？据此，我们可以在实现的时候在同一个函数中写出逆变换和正变换，然后在得到的结果 res 中除以 $n$ 就可以了。（参照下面的代码）

至此，FFT 的基本原理讲述完毕，下面是优化。

位逆序置换

按照上文的讲述，如果不看下面的代码，那么编写出来的是递归版本，但是这个版本的常数太大了，因此运行起来的效果不好，故使用位逆序置换来降低常数。

我们看看递归过程是什么样的，以 $n=8$ 为例：

$$\{x_0, x_1, x_2,x_3,x_4,x_5,x_6,x_7\}\\ \{x_0, x_2, x_4,x_6\},\{x_1,x_3,x_5,x_7\}\\ \{x_0, x_4\}, \{x_2,x_6\},\{x_1,x_5\},\{x_3,x_7\}\\ \{x_0\},\{x_4\},\{x_2\},\{x_6\},\{x_1\},\{x_5\},\{x_3\},\{x_7\}\\ $$

这里就有一个非常神奇的规律：在最后一行中，原下标所对应的二进制数翻转正好是在最后一行的序数。例如 $x_6$ 的下标是 $6 = 110_{(2)}$，那么它的序数正好是 $011_{(2)} = 3$。

据此，可以处理出 rev 数组，它记录的正是最后一行所有元素对应的下标。

简单地说，递归形式是自上而下地做 FFT，而利用位逆序置换我们可以自下而上地做 FFT，它们在实际运行中有着常数上的区别。

模板题及代码

https://www.luogu.com.cn/problem/P3803

给定一个 $n$ 次多项式 $F(x)$，和一个 $m$ 次多项式 $G(x)$。

请求出 $F(x)$ 和 $G(x)$ 的卷积。

#include<bits/stdc++.h>
using namespace std;

const int N=3e5+5;
const double pi=acos(-1);

int n, m;

// 复数类
struct Complex{
	double x, y;
	Complex operator + (const Complex &o)const { return {x+o.x, y+o.y}; }
	Complex operator - (const Complex &o)const { return {x-o.x, y-o.y}; }
	Complex operator * (const Complex &o)const { return {x*o.x-y*o.y, x*o.y+y*o.x}; }
};

Complex a[N], b[N];
int res[N];

int rev[N], bit, tot;

void fft(Complex a[], int inv){ // inv 指示正变换、逆变换。
	for(int i=0; i<tot; i++) if(i<rev[i]) swap(a[i], a[rev[i]]);
	
	for(int mid=1; mid<tot; mid<<=1){
		auto w1=Complex({cos(pi/mid), inv*sin(pi/mid)});
		for(int i=0; i<tot; i+=mid*2){
			auto wk=Complex({1, 0});
			for(int j=0; j<mid; j++, wk=wk*w1){
				auto x=a[i+j], y=wk*a[i+j+mid];
				a[i+j]=x+y, a[i+j+mid]=x-y;
			}
		}
	}
}

int main(){
	cin>>n>>m;
	for(int i=0; i<=n; i++) cin>>a[i].x;
	for(int i=0; i<=m; i++) cin>>b[i].x;
	
	while((1<<bit)<n+m+1) bit++; // 结果次数分布在 [0, n+m] 内，一共有 n+m+1 位。
	
	tot=1<<bit; // 得到上文所说的 n
	
	for(int i=0; i<tot; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
	
	fft(a, 1), fft(b, 1); // 正变换 DFT
	
	for(int i=0; i<tot; i++) a[i]=a[i]*b[i];
	
	fft(a, -1); // 逆变换 IDFT
	
	for(int i=0; i<=n+m; i++) res[i]=(int)(a[i].x/tot+0.5), printf("%d ", res[i]);
	
	return 0;
}