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$。
总之
总的过程大概就是:
- 进行一些预处理,把次数拓展成2的次幂,高位多出来的设置成0;
- 用$DFT$在$\Theta (nlogn)$分别将两个多项式$A(x)$和$B(x)$转化成点值的形式;
- 利用点值表示在$\Theta (n)$直接就能求得$A(x)\times B(x)$的点值表示;
- 用$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); }