快速傅里叶变换FFT学习小记

Posted 2020-08-05 新博客↓

FFT学得还是有点模糊，原理那些基本还是算有所理解了吧，不过自己推这个推不动。

看的资料主要有这两个：

http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform

https://www.zybuluo.com/397915842/note/37965

这儿简单做做笔记。

多项式点值表示

首先$FFT$用来快速计算两个多项式的乘积。

一个$n$次多项式（最高次为$n$），可以用系数表示法表示和点值表示法表示。

• 系数表示法：$A(x)=\sum_{i=0}^{n}c_ix^i$，可以知道用系数表示法进行多项式乘法时间复杂度是$\Theta (n^2)$。
• 点值表示法：用$n+1$个点$(x_k,y_k)，k=0\dots n$，其中$x_k$互不相同，$y_k=A(x_k)$。

这$n+1$个点是可以唯一表示一个$n$次多项式的，因为：$$\begin{bmatrix} 1 & x_0 & x_0^2 & ... & x_0^n \\ 1 & x_1 & x_1^2 & ... & x_1^n \\\vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & ... & x_n^n \\\end{bmatrix}\begin{bmatrix} c_0 \\ c_1 \\ \vdots \\c_n \\\end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\y_n \\\end{bmatrix}$$

前面那个方阵是范德蒙德矩阵，其行列式可以知道是非0，所以该方阵存在逆矩阵，方程有解且有唯一解。

 

现在考虑两个n次多项式$A(x)$和$B(x)$，其乘积的多项式结果是$C(x)$，那么$C(x)$的次数为$2n$。

不妨给$A$和$B$取$2n+1$个点，这样求得两个多项式的点值表示为$(x_k,A(x_k))、(x_k,B(y_k))，k=0\dots 2n$，而$C(x)$的点值表示在$\Theta (n)$就能直接求出，即为$(x_k,A(x_k)B(y_k))$。

 

这样进行多项式乘法运算，原本用系数表示多项式需要$\Theta (n^2)$的时间复杂度，而点值则只需要$\Theta (n)$。

 

不过将一个多项式转化成点值的形式，需要枚举$2n$个点，然后对于每一个点$\Theta (n)$求得多项式的值，这样时间复杂度是$\Theta (n^2)$的；另外，求得结果后也是一个点值的表示，需要还原回系数形式，直接解上面矩阵表示的方程，这要立方的时间复杂度。

这时就是$FFT$做的事了。前面转化为点值形式称为$DFT$，后面逆转化回系数形式称为$IDFT$。

 

DFT

对于多项式$A(x)$，比如我们要取$n$个点$(x_k,y_k)$，$x_k$选取的是$n$次单位根（满足$z^n=1$的复数解$z$），即：$$e^{\frac{2\pi ki}{n}}, k = 0\cdots n - 1, i=\sqrt{-1}$$

另外，计算这个可以用欧拉公式：$$e^{\theta i}=\cos\theta + i\sin\theta$$

记$\omega_n=e^{\frac{2\pi i}{n}}$，那么点值表示即为$$(\omega_n^0, A(\omega_n^0)), (\omega_n^1, A(\omega_n^1)),\cdots ,(\omega_n^{n-1}, A(\omega_n^{n-1}))$$

 

现在如何求$A(\omega_n^0), A(\omega_n^1), \cdots, A(\omega_n^{n-1})$？

看下面：$$A(x) = c_0 + c_1x + c_2x^2 + c_3x^3 + \cdots + c_nx^{n-1}\\A^{[0]}(x) = c_0 + c_2x + c_4x^2 + \cdots + c_{n - 2}x^{n/2 - 1} \\A^{[1]}(x) = c_1 + c_3x + c_5x^2 + \cdots + c_{n - 1}x^{n/2 - 1}$$

然后，可以发现$$A(x) = A^{[0]}(x^2) + xA^{[1]}(x^2)$$

而这边是将$\omega_n^k，k=0\cdots n-1$代入$x$，其中$x^2$也就是$(\omega_n^k)^2$可以化简：$$(\omega_n^k)^2=\left(e^{\frac{2\pi ki}{n}}\right)^2=e^{\frac{2\pi ki}{n/2}}=\omega_{\frac{n}{2}}^k$$而$\omega_{\frac{n}{2}}^k$复数根有$n/2$个。

那么可以看到，把$A(x)$分成两个子问题，而这两个子问题选取的点也缩减了一半。这样时间复杂度就是$\Theta (nlogn)$了。于是可以用递归去实现这个分治的算法。另外一开始$n$拓展成$2$的次幂，多出来的补$0$即可。

 

不过，事实上可以改用迭代实现，而不用递归。画出递归的每一次划分。。$此处没图$。。

然后对比最初和最后一层各个下标的二进制，可以发现相当于是每一个下标的二进制反转了。

那么可以一开始直接反转，从最后一层开始，一层层一个个迭代上去去求得结果。

另外，有个现成的$Rader$雷德算法可以直接求得反转结果。

 

IDFT

上面就在$\Theta (nlogn)$下将系数表示转化成点值表示；而现在需要逆过来，把点值表示还原成所需要的结果，系数表示。$$ \begin{bmatrix} (\omega_n^0)^0 & (\omega_n^0)^1 & \cdots &(\omega_n^0)^{n-1} \\ (\omega_n^1)^0 & (\omega_n^1)^1 & \cdots & (\omega_n^1)^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{n-1})^0 & (\omega_n^{n-1})^1 & \cdots & (\omega_n^{n-1})^{n-1} \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ \vdots \\ c_{n-1} \end{bmatrix} = \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ \vdots \\ A(\omega_n^{n-1}) \end{bmatrix}$$

其实就是求$c_k$。

而那个范德蒙德矩阵的逆也有结论的，就是：$$\frac{1}{n}\begin{bmatrix} (\omega_n^{-0})^0 & (\omega_n^{-0})^1 & \cdots & (\omega_n^{-0})^{n-1} \\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 & \cdots & (\omega_n^{-1})^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 & \cdots & (\omega_n^{-(n-1)})^{n-1} \end{bmatrix}$$

证明过程。。两个矩阵相乘。。我就不搬运了= =。。

这样的话等式两边同时乘上方阵的逆：$$\begin{bmatrix} c_0 \\ c_1 \\ \vdots \\ c_{n-1} \end{bmatrix} = \frac{1}{n} \begin{bmatrix} (\omega_n^{-0})^0 & (\omega_n^{-0})^1 & \cdots & (\omega_n^{-0})^{n-1} \\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 & \cdots & (\omega_n^{-1})^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 & \cdots & (\omega_n^{-(n-1)})^{n-1} \end{bmatrix} \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ \vdots \\ A(\omega_n^{n-1}) \end{bmatrix}$$

那么可以发现，这可以用$DFT$一样的过程去求了，不同的是选的点是$\omega_n^{-k}$，最后结果除以$n$。

 

总之

总的过程大概就是：

1. 进行一些预处理，把次数拓展成2的次幂，高位多出来的设置成0；
2. 用$DFT$在$\Theta (nlogn)$分别将两个多项式$A(x)$和$B(x)$转化成点值的形式；
3. 利用点值表示在$\Theta (n)$直接就能求得$A(x)\times B(x)$的点值表示；
4. 用$IDFT$在$\Theta (nlogn)$将乘积结果的点值表示还原成系数表示。

 

代码的实现照着别人的，稍微改了下点，写了一遍：

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
#define MAXN 133333
const double PI=acos(-1.0);

struct Complex{
	double real,imag;
	Complex(){}
	Complex(double _real,double _imag):real(_real),imag(_imag){}
	Complex operator-(const Complex &cp) const{
		return Complex(this->real-cp.real,this->imag-cp.imag);
	}
	Complex operator+(const Complex &cp) const{
		return Complex(this->real+cp.real,this->imag+cp.imag);
	}
	Complex operator*(const Complex &cp) const{
		return Complex(this->real*cp.real-this->imag*cp.imag,this->real*cp.imag+this->imag*cp.real);
	}
	void setValue(double _real=0,double _imag=0){
		real=_real; imag=_imag;
	}
};

int len;
Complex wn[MAXN],wn_anti[MAXN];

void FFT(Complex y[],int op){
	// Rader, 位逆序置换 
	for(int i=1,j=len>>1,k; i<len-1; ++i){
		if(i<j) swap(y[i],y[j]);
		k=len>>1;
		while(j>=k){
			j-=k;
			k>>=1;
		}
		if(j<k) j+=k;
	}
	// h=1，Wn^0=1 
	for(int h=2; h<=len; h<<=1){
		// Wn = e^(2*PI*i/n)，如果是插值Wn = e^(-2*PI*i/n)，i虚数单位 
		Complex Wn = (op==1 ? wn[h] : wn_anti[h]);
		for(int i=0; i<len; i+=h){
			Complex W(1,0);
			for(int j=i; j<i+(h>>1); ++j){
				Complex u=y[j],t=W*y[j+(h>>1)];
				y[j]=u+t;
				y[j+(h>>1)]=u-t;
				W=W*Wn;
			}
		}
	}
	if(op==-1){
		for(int i=0; i<len; ++i) y[i].real/=len;
	}
}
void Convolution(Complex A[],Complex B[],int n){
	// 初始化 
	for(len=1; len<(n<<1); len<<=1);
	for(int i=n; i<len; ++i){
		A[i].setValue();
		B[i].setValue();
	}
	// e^(θi) = cosθ+isinθ -> Wn = cos(2*PI/n)+isin(2*PI/n)
	for(int i=0; i<=len; ++i){ // 预处理 
		wn[i].setValue(cos(2.0*PI/i),sin(2.0*PI/i)); 
		wn_anti[i].setValue(wn[i].real,-wn[i].imag);
	}
	FFT(A,1); FFT(B,1);
	for(int i=0; i<len; ++i){
		A[i]=A[i]*B[i];
	}
	FFT(A,-1);
}