快速傅里叶变换浅析
问题描述
对于给定的\(n\)次多项式\(F_1(x)\)和\(m\)次多项式\(F_2(x)\)，求\(F_1(x)\)与\(F_2(x)\)的卷积。
关于多项式
卷积
通过两个函数\(F_1\),\(F_2\)生成第三个函数的数学算子，可以简单理解为多项式的积。
多项式的表示方法
系数表示法：通过数列\(a=\{a_1,a_2\ldots a_n\}\)表示n次多项式的系数，即\(F(x)=∑_{(i=0)}^na_i x^i\)。系数表示法下的n次多项式\(F(x)\)和\(m\)次多项式\(G(x)\)用乘法分配率进行运算的时间复杂度为\(O(nm)\)。
点值表示法：由于对于任意多项式\(F(x)=∑_{i=0}^na_i x^i\)都可以在平面上以曲线形式表示，而一个\(n\)次曲线可以被曲线上\((n+1)\)个点\(\{(x_i,y_i)|1≤i≤n+1\}\)唯一地表示，因此这\((n+1)\)个点的坐标可以唯一地表示\(F(x)\)。
复数运算
定义
虚数：定义与数轴上点相对应的数为实数；相应地，定义在数轴上没有对应点的数为虚数。负数的偶次方根在数学上被定义为纯虚数，记作\(z=bi(b∈R且b≠0,i=\sqrt{-1})\)，其中\(i\)为虚数单位，\(i=\sqrt{-1}\)。
复数：在数学上，定义形如\(z=a+bi(a,b∈R,i=\sqrt{-1})\)的数为复数，其中\(a\)称为实部，\(b\)称为虚部。当\(b=0\)时，\(z\)为实数；当\(b≠0\)时，\(z\)为纯虚数。
复平面：由于所有复数\(z=a+bi\)都可以表示成类似平面向量的形式，因此我们将复数表示在一个以实轴与虚轴构成的平面直角笛卡尔坐标系内，就可用平面向量运算法则进行复数运算。这样以实轴与虚轴构成的平面直角笛卡尔坐标系我们称之为复平面，其中横轴为实轴，竖轴为虚轴。
平面向量运算
平面向量：我们定义二维平面内既有方向又有大小的量为平面向量，向量在横纵轴上分别有其分量（即有向线段在两轴上的正投影）。
弧度制：指用\(\frac{弧长}{半径}\)度量对应圆心角度数的方式，记作\(rad\)（但通常省略）。\(1rad= \frac{180°}{π}\)，因此弧度制中通常使用\(π\)表示弧度。
平面向量加法运算法则：平面向量加法运算遵循平行四边形定则（三角形定则），即用两向量为边作平行四边形（三角形），向量和即为平行四边形的对角线（三角形的第三边）。在代数运算上，两向量的分量分别作代数加法形成和的分量。运算过程为
\(A=\left[\begin{matrix}x_1\\y_1\end{matrix}\right],B=\left[\begin{matrix}x_2\\y_2\end{matrix}\right],A+B=\left[\begin{matrix}x_1+x_2\\y_1+y_2\end{matrix}\right]\)。
平面向量乘法运算运算法则：平面向量作乘法运算时，将两向量模长（有向线段长度）作代数乘法后与两向量夹角的正弦值作代数乘法得到向量积的模长。运算过程为\(|A×B|=|A||B|sin<A,B>\)。向量积表示为一个空间向量，无法在平面中表示出来，但其模长与两向量形成平行四边形面积数值相等。
复数运算
复数加法：复数加法运算完全遵循平面向量的加法运算，即分量相加形成和的分量。运算过程为\(z_1=a_1+b_1 i,z_2=a_2+b_2 i,z_1+z_2=(a_1+a_2 )+(b_1+b_2)i\)。
复数乘法：复数乘法运算法则与平面向量乘法运算稍有差异，遵循模长相乘，辐角（有向线段与实轴的夹角）相加的原则。计算过程如下：
单位根
定义
在复平面上，以原点为圆心，作\(r=1\)的单位圆。以圆心为起点，圆的\(n(n=2^k,k∈Z)\)等分点为终点作\(n\)个向量。设其中辐角为正且最小的向量所对应的复数为\(n\)次单位根，记作\(ω_n\)。单位根辐角为\(\frac{2π}{n}\)。在代数中，若\(z^n=1\)，就称\(z\)为\(n\)次单位根。
根据复数乘法运算法则，其余\((n-1)\)个复数分别为\(ω_n^1,ω_n^2\ldots ω_n^n\)。其中\(ω_n^0=ω_n^n=1\)，对应复平面上以实轴为正方向的向量。
对于其余\((n-1)\)个复数的运算，可以用欧拉公式解决：
性质
对于任意\(ω_n^k\)，均有\(ω_n^k=\cosk\frac{2π}{n}+i\sink \frac{2π}{n}\)。
\(ω_2n^2k=ω_n^k\)。证明：\(ω_2n^2k=\cos2k\frac{2π}{2n}+i\sin2k\frac{2π}{2n}=ω_n^k\)。
\(ω_n^{k+\frac{n}{2}}=-ω_n^k\)。证明：\(ω_n^{\frac{n}{2}}=\cos\frac{n}{2}•\frac{2π}{n}+i\sin\frac{n}{2}=\cosπ+i\sinπ=-1\)。
\(ω_n^0=ω_n^n=1\)。由零次幂定义可知。
快速傅里叶变换
关于FFT的高效性
系数表示法下的多项式乘法：由于根据乘法分配律，需要将一个多项式中的每一项与另一多项式中所有项分别相乘，所以系数表示法下的多项式时间复杂度为\(O(nm)\)，运算过程如下：
点值表示法下的多项式乘法：对于两个化为点值表示法的多项式，只需要将其各点值对应相乘（运算法则参考平面向量乘法），所以时间复杂度为\(O(n)\)。
FFT具体原理
前文中我们提到了多项式的点值表示法，但我们需要知道求哪些点的点值，而单位根的概念就解决了这一问题。对于一个\(n\)次多项式，我们用\(ω_n^k(k∈Z且0≤k≤n)\)来确定这\((n+1)\)个点值。但上述做法的时间复杂度仍然是\(O(n^2)\)，因此需要引入分治思想优化。
观察\(n\)次（此处默认\(n\)为偶数）多项式\(F(x)\)的表示形式：
若直接利用单位根求其一个点值，时间复杂度是\(O(n)\)，总体时间复杂度仍然为\(O(n^2)\)。但若将其按指数奇偶性分类，可以得到如下形式：
再将两部分拆分为两个多项式，即可得：
再将两部分拆分为两个多项式，即可得：
根据函数奇偶性质可知\(G_1\)为奇函数（\(G_1(x)=-G_1(-x)\)），而\(G_2\)为偶函数（\(G_2(x)=G_2(-x)\)）。那么对于求出的一个点值，只需要利用奇偶性直接得出与其对称的另一个点值，由此可将计算1个点值的时间复杂度将为\(O(\frac{n}{2})\)。继续观察，不难发现若将\(G_1\)中的\(x\)提取出来，就可以得到：
将\(x^2\)作为整体代入，则得到：
不难发现，我们可以把\({G'}_1,G_2\)继续进行与上述步骤相同的操作，就可以将计算\({G'}_1,G_2\)点值的时间复杂度降为原来的\(\frac{1}{2}\)。由此，我们定义出了一个形如二叉树的结构，即可通过分治算法处理点值计算的问题。如此一来，计算\(F(x)\)的时间复杂度就优化到了\(O(n\log_2n)\)。
快速傅里叶逆变换
在前文中，笔者讨论了基于多项式点值表示法的快速傅里叶变换原理。但是由于需要求出原多项式卷积，所以计算结果必须转化回系数表示法下的多项式，这里就涉及到快速傅里叶逆变换。
我们设\(\{y_0,y_1…y_n\}\)为\(F(x)=∑_{i=0}^na_ix^i\)的傅里叶变换。设有另一向量\(\{A_0,A_1…A_n\}\)满足\(A_k=∑_{i=0}^ny_i(ω_n^{-k})^i\)，即为\(G(x)=∑_{i=0}^ny_i x^i\)的点值表示形式。对\(∑_{i=0}^na_ix^i\)行展开变形得到如下：
设\(S(x)=∑_{i=0}^nx^i\)，将单位根\(ω_n^k\)代入得到：
由上述推导可得：
将其带入前一个推导出的式子\(A_k=∑_{i=0}^na_i [∑_{j=0}^n(ω_n^{j-k})^i]\)，可得：
故得：
由此即可得到点值表示法下的多项式向系数表示法的傅里叶逆变换。
算法总结与补充
关于时间复杂度
FFT的实现依赖了多项式的点值表示法与单位根，两者的时间复杂度都为\(O(n\log_2 n)\)，而中途的点值乘法运算时间复杂度为\(O(n)\)，所以最终的时间复杂度为\(O(n\log_2n)\)。
FFT的分治思想主要体现在其求点值是对多项式按次数奇偶性的划分，通过分治可以将原本为\(O(n^2)\)的点值运算时间复杂度优化为\(O(n\log_2n)\)。
关于FFT相对于系数表示法下多项式乘法的具体优化过程，笔者选择用《算法导论》第30章中关于FFT的具体计算流程图进行阐释。
关于算法的语言实现（C++）
FFT的实现必然涉及负数运算。C++的STL库中提供了#include<complex>头文件，其中有现成的complex类。但由于算法竞赛对时间复杂度与空间复杂度控制要求较高，笔者在此不建议使用库函数。
递归写法的FFT由于栈开销过大，容易被各种毒瘤题目卡时间复杂度，因此笔者建议在实现过程中采用迭代写法。
点击查看代码
#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<cstring>
#include<queue>
#include<set>
#include<ctime>
#define RI register int
#define RC register complex
using namespace std;
const int maxn=1e7+10;
const double Pi=acos(-1.0);/*定义弧度制的π*/
struct complex/*复数*/
{
double x,y;
inline complex(double tx=0,double ty=0){x=tx;y=ty;}
}A[maxn],B[maxn],R[maxn];
/*复数运算（等同于复平面上的向量运算）*/
inline complex operator +(complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
inline complex operator -(complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
inline complex operator *(complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
int n,m,limit=1;
int l,r[maxn];
inline void FFT(complex *X,int opt);
int main()
{
ios::sync_with_stdio(false);cin.tie(0);
cin >> n >> m;
for(RI i=0;i<=n;i++) cin >> A[i].x;
for(RI i=0;i<=m;i++) cin >> B[i].x;
for(;limit<=m+n;limit<<=1,l++);
for(RI i=0;i<limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
FFT(A,1);FFT(B,1);/*对原多项式分别进行快速傅里叶变换*/
for(RI i=0;i<=limit;i++) R[i]=A[i]*B[i];/*点值计算*/
FFT(R,-1);/*对结果进行逆傅里叶变换 转为系数表示法*/
for(RI i=0;i<=n+m;i++)
printf("%d ",(int)(R[i].x/limit+0.5));
return 0;
}
inline void FFT(complex *X,int opt)
{
for(RI i=0;i<limit;i++)
if(i<r[i]) swap(X[i],X[r[i]]);/*求出迭代序列*/
for(RI mid=1;mid<limit;mid<<=1)
{
complex Ur(cos(Pi/mid),opt*sin(Pi/mid));/*求单位根*/
for(int len=mid<<1,i=0;i<limit;i+=len)/*确定区间*/
{
complex temp(1,0);
for(RI k=0;k<mid;k++,temp=temp*Ur)/*枚举区间左半部分*/
{
/*蝴蝶操作*/
RC x=X[i+k],y=temp*X[i+mid+k];
X[i+k]=x+y;
X[i+mid+k]=x-y;
}
}
}
return;
}
