FFT快速傅里叶变换

快速傅里叶变换——FFTyL2013-9现在是1页\一共有51页\编辑于星期日我们关心的问题快速解决多项式乘法问题衍生问题——高精度乘法现在是2页\一共有51页\编辑于星期日问题的描述记一个多项式次数界为n的多项式A(x)则其中a为每一项的系数注意最高次系数为n-1现在是3页\一共有51页\编辑于星期日问题的描述两个多项式相乘我们记一般形式为C的次数界为A与B次数界的和普通的时间复杂度为O(n^2)现在是4页\一共有51页\编辑于星期日PART1——中心思想1点值与插值2N次单位复数根3FFT4演示现在是5页\一共有51页\编辑于星期日转换思路普通的相乘方法提出概念：点值，插值1点值与插值2N次单位复数根3FFT4演示现在是6页\一共有51页\编辑于星期日点值一个次数界为n的多项式A(x)的点值表达就是一个由n个点值对所组成的集合：其中每一个x都不相同，且E.G.多项式的一个合法点值表达是1点值与插值2N次单位复数根3FFT4演示现在是7页\一共有51页\编辑于星期日插值插值运算是点值运算的逆运算假设我们得到了一个有n个点值对的点值表达那我们可以确定唯一的一个次数界为n的多项式1点值与插值2N次单位复数根3FFT4演示现在是8页\一共有51页\编辑于星期日多项式乘法我们来探究一下如何用点值与插值来完成多项式乘法我们确定一组x，求得A与B的点值表达那我们可以得知C的点值表达通过插值运算，我们可以知道多项式C的系数1点值与插值2N次单位复数根3FFT4演示现在是9页\一共有51页\编辑于星期日注意的地方设A与B的次数界均为n则C的次数界为2n则我们要找出2n个x来求点值表达否则不可以进行插值运算1点值与插值2N次单位复数根3FFT4演示现在是10页\一共有51页\编辑于星期日算法流程对于次数界均为n的多项式A与B1点值运算构造长度为2n的点值表达2逐点相乘计算出C的点值表达3插值运算通过C的点值表达求出多项式C的每项系数1点值与插值2N次单位复数根3FFT4演示现在是11页\一共有51页\编辑于星期日时间复杂度可以证明，若选取n个x计算点值与插值的时间复杂度均为O(N^2)本质上没有解决时间的问题但我们可以巧妙的选择这些数来优化时间复杂度。1点值与插值2N次单位复数根3FFT4演示现在是12页\一共有51页\编辑于星期日PART2——N次单位复数根及其性质1点值与插值2N次单位复数根3FFT4演示现在是13页\一共有51页\编辑于星期日N次单位复数根n次单位复数根是满足的复数w。n次单位复数根根恰好有n个，对于k=0,1,...,n-1，这些根是。为了解释这个表达式，我们利用复数的指数形式的定义：下一页图说明n个单位复数根均匀地分布在以复平面的原点为圆心的单位半径的圆周上。1点值与插值2N次单位复数根3FFT4演示现在是14页\一共有51页\编辑于星期日N次单位复数根1点值与插值2N次单位复数根3FFT4演示现在是15页\一共有51页\编辑于星期日性质我们需要N次单位复数根我们首先来探究这些根的性质1点值与插值2N次单位复数根3FFT4演示现在是16页\一共有51页\编辑于星期日性质1——主n次单位根我们称为主n次单位根同时注意到，n次单位复数根都是经过旋转而得到的每次旋转都是一定角度n次单位复数根可视为公比为主n次单位根的等比数列1点值与插值2N次单位复数根3FFT4演示现在是17页\一共有51页\编辑于星期日性质2——群的性质因为所以推论1点值与插值2N次单位复数根3FFT4演示现在是18页\一共有51页\编辑于星期日性质3——消去引理&折半引理....消去引理：推论：折半引理：如果n>0为偶数，那么n次单位复数根的平方的集合就是n/2次单位复数根的集合。证明：可以知道的平方相同。1点值与插值2N次单位复数根3FFT4演示现在是19页\一共有51页\编辑于星期日性质4——求和引理求和引理：对于任意整数n≥1和不能被n整除的非负整数k，有等比数列求和所以注意k不能是n的倍数，否则分母为01点值与插值2N次单位复数根3FFT4演示现在是20页\一共有51页\编辑于星期日PART3——FFT及其关键算法1点值与插值2N次单位复数根3FFT4演示现在是21页\一共有51页\编辑于星期日DFT——离散傅里叶变换我们希望计算次数界为n的多项式在n次单位复数根处的值（共n个）接下来定义结果yy即为a的离散傅里叶变换（DFT）我们也可记为1点值与插值2N次单位复数根3FFT4演示现在是22页\一共有51页\编辑于星期日FFT——快速傅里叶变换大前提：n为2的整数幂（方便计算）利用复数单位根复数根的特殊性质我们可以在时间内计算出FFT利用了分治策略1点值与插值2N次单位复数根3FFT4演示现在是23页\一共有51页\编辑于星期日PART3.1——点值运算1点值与插值2N次单位复数根3FFT4演示现在是24页\一共有51页\编辑于星期日分治策略如何求出单个数x的函数值A(x)？我们定义两个新多项式观察两个多项式的特点1分别拥有奇数下标的系数与偶数下标的系数2次数界变为n/2（缩小了一半）1点值与插值2N次单位复数根3FFT4演示现在是25页\一共有51页\编辑于星期日分治策略对于一个数x，求A(x)则根据上两个多项式1点值与插值2N次单位复数根3FFT4演示现在是26页\一共有51页\编辑于星期日分治策略至此我们成功的转换了问题原问题：求一个多项式A(x)在的函数值。现问题：求两个多项式在的函数值。1点值与插值2N次单位复数根3FFT4演示现在是27页\一共有51页\编辑于星期日分治策略现问题：求两个多项式在的函数值。根据折半引理，并不是n个不同的值，而是由n/2个n/2次单位复数根组成，而每个根恰好出现2次于是，我们递归地对n/2的多项式A[0](x)与A[1](x)在n/2个n/2次单位复数根进行求值1点值与插值2N次单位复数根3FFT4演示现在是28页\一共有51页\编辑于星期日程序实现1点值与插值2N次单位复数根3FFT4演示现在是29页\一共有51页\编辑于星期日我们关心的程序实现问题1/2：规定程序递归出口3/4/12：定义主n次单位根，更新w值5/6/7/8：定义两个多项式并递归求解13：返回DFT9/10/11：递归结束后的处理工作1点值与插值2N次单位复数根3FFT4演示现在是30页\一共有51页\编辑于星期日递归结束后的处理工作10：11：1点值与插值2N次单位复数根3FFT4演示现在是31页\一共有51页\编辑于星期日FFT时间复杂度对于运行时间有以下的递归式所以采用FFT，我们可以在O(nlgn)时间内实现点值运算（求出次数界为n的多项式在n次单位复数根处的值）。1点值与插值2N次单位复数根3FFT4演示现在是32页\一共有51页\编辑于星期日PART3.2——插值运算1点值与插值2N次单位复数根3FFT4演示现在是33页\一共有51页\编辑于星期日矩阵乘积我们可以把点值运算表示成一个矩阵方程我们把DFT写成矩阵乘积y=Vna1点值与插值2N次单位复数根3FFT4演示现在是34页\一共有51页\编辑于星期日逆矩阵到此为止我们把插值运算改写成y与的逆矩阵的乘积1点值与插值2N次单位复数根3FFT4演示现在是35页\一共有51页\编辑于星期日逆矩阵定理：证明比较冗长。我们证明，其中In为n*n的单位矩阵。考虑中(j,j')处的元素：如果j=j'，则此和为1否则根据求和引理，此和为01点值与插值2N次单位复数根3FFT4演示现在是36页\一共有51页\编辑于星期日逆DFT给定矩阵，可以推导出：通过比较DFT的运算我们可以看到，对FFT作以下修改就可以计算出逆DFT：把a与y互换，用，再把计算结果的每个元素除以n。于是我们也可以在O(nlgn)时间内算出逆DFT。1点值与插值2N次单位复数根3FFT4演示现在是37页\一共有51页\编辑于星期日PART4——标程演示1点值与插值2N次单位复数根3FFT4演示现在是38页\一共有51页\编辑于星期日程序实现优化因为像伪代码中，赋值数组是不切实际的但我们发现，a中的应用的系数是有规律的所以根据位移与起始点的不同来优化FFT的实现1点值与插值2N次单位复数根3FFT4演示现在是39页\一共有51页\编辑于星期日程序实现优化(a0,a1,a2,a3,a4,a5,a6,a7)(a0,a2,a4,a6)(a1,a3,a5,a7)(a0,a4)(a2,a6)(a1,a5)(a3,a7)(a0)(a4)(a2)(a6)(a1)(a5)(a3)(a7)1点值与插值2N次单位复数根3FFT4演示现在是40页\一共有51页\编辑于星期日type Node=record x,y:double; end;arr=array[0..200000]ofNode;operator+(a,b:Node)c:Node;beginc.x:=a.x+b.x;c.y:=a.y+b.y;end;operator-(a,b:Node)c:Node;beginc.x:=a.x-b.x;c.y:=a.y-b.y;end;operator*(a,b:Node)c:Node;beginc.x:=a.x*b.x-a.y*b.y;c.y:=a.x*b.y+a.y*b.x;end;//定义复数类型，复数运算procedureDft(vara:arr;s,t:longint);//a答案数组，s起始量，t位移量begin if(n>>t)=1thenexit; Dft(a,s,t+1);Dft(a,s+1<<t,t+1); fori:=0ton>>(t+1)-1do begin j:=s+i<<(t+1); wt:=w[i<<(t)]*a[j+1<<t]; tt[i]:=a[j]+wt; tt[i+n>>(t+1)]:=a[j]-wt; end; fori:=0ton>>t-1doa[s+i<<t]:=tt[i];end;1点值与插值2N次单位复数根3FFT4演示现在是41页\一共有51页\编辑于星期日begin n:=1; cin(a);cin(b); fori:=0ton-1dow[i].x:=cos(pi*i*2/n); fori:=0ton-1dow[i].y:=sin(pi*i*2/n); Dft(a,0,0);Dft(b,0,0);//DFT fori:=0ton-1doa[i]:=a[i]*b[i]; fori:=0ton-1dow[i].y:=-w[i].y; Dft(a,0,0);//逆DFT fori:=0ton-1do begin ans[i]:=ans[i]+round(a[i].x/n); ans[i+1]:=ans[i]div10; ans[i]:=ans[i]mod10; end; while(i>0)and(ans[i]=0)dodec(i); forj:=idownto0do write(ans[j]); writeln;end.//感谢wwt同学提供的程序1点值与插值2N次单位复数根3FFT4演示现在是42页\一共有51页\编辑于星期日#include<iostream>#include<cstdio>#include<string.h>#include<complex>#defineN50010//usingnamespacestd;//出处：/ysjjovo/article/details/6299359感谢之constdoublepi=acos(-1);constcomplex<double>I(0,1);constdoubleeps=1e-6;charaa[N],bb[N];intans[4*N];//charans[4*N];!!!!!!!complex<double>a[4*N],b[4*N];//an-1,an-2,……,a1,a0intn;voidinitial(){intlenA=strlen(aa),lenB=strlen(bb);n=max(lenA,lenB);doublet=log2(n);inttt=int(t);if(t-tt>0)tt++;n=1<<(tt+1);//doublelengthinti;for(i=0;i<lenA;i++)a[i]=aa[lenA-1-i]-'0';while(i<n)a[i++]=0;for(i=0;i<lenB;i++)b[i]=bb[lenB-1-i]-'0';while(i<n)b[i++]=0;}1点值与插值2N次单位复数根3FFT4演示现在是43页\一共有51页\编辑于星期日voidbitReverse(complex<double>*a){for(inti=1;i<n-1;i++){intj=0;for(intk=1,tmp=i;k<n;j=((j<<1)|(tmp&1)),k<<=1,tmp>>=1);if(j>i)swap(a[i],a[j]);}}voiditerative_FFT(complex<double>*a,intsig){bitReverse(a);for(intm=2;m<=n;m<<=1)//是把a[]按m的长度分组，{intmh=m>>1;for(inti=0;i<mh;i++){complex<double>wi=exp(i*sig*pi/mh*I);for(intj=i;j<n;j+=m){intk=j+mh;complex<double>u=a[j],t=wi*a[k];;a[j]=u+t;a[k]=u-t;}}}if(sig==-1)for(inti=0;i<n;i++)a[i]/=n;}1点值与插值2N次单位复数根3FFT4演示现在是44页\一共有51页\编辑于星期日voidconvolution(){iterative_FFT(a,1);iterative_FFT(b,1);for(inti=0;i<n;i++)a[i]*=b[i];//a*biterative_FFT(a,-1);}voidoutput(){intk=0;ans[0]=0;for(inti=0;i<n;i++)//{inttmp=a[i].real()+eps;//fourignorefiveincreaseans[i]+=tmp;if(ans[i])k=i,ans[i+1]=ans[i]/10,ans[i]%=10;elseans[