FFT快速傅里葉變換

快速傅里葉變換——FFTyL2013-9目前一頁\總數(shù)五十一頁\編于十三點(diǎn)我們關(guān)心的問題快速解決多項(xiàng)式乘法問題衍生問題——高精度乘法目前二頁\總數(shù)五十一頁\編于十三點(diǎn)問題的描述記一個(gè)多項(xiàng)式次數(shù)界為n的多項(xiàng)式A(x)則其中a為每一項(xiàng)的系數(shù)注意最高次系數(shù)為n-1目前三頁\總數(shù)五十一頁\編于十三點(diǎn)問題的描述兩個(gè)多項(xiàng)式相乘我們記一般形式為C的次數(shù)界為A與B次數(shù)界的和普通的時(shí)間復(fù)雜度為O(n^2)目前四頁\總數(shù)五十一頁\編于十三點(diǎn)PART1——中心思想1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前五頁\總數(shù)五十一頁\編于十三點(diǎn)轉(zhuǎn)換思路普通的相乘方法提出概念：點(diǎn)值，插值1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前六頁\總數(shù)五十一頁\編于十三點(diǎn)點(diǎn)值一個(gè)次數(shù)界為n的多項(xiàng)式A(x)的點(diǎn)值表達(dá)就是一個(gè)由n個(gè)點(diǎn)值對所組成的集合：其中每一個(gè)x都不相同，且E.G.多項(xiàng)式的一個(gè)合法點(diǎn)值表達(dá)是1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前七頁\總數(shù)五十一頁\編于十三點(diǎn)插值插值運(yùn)算是點(diǎn)值運(yùn)算的逆運(yùn)算假設(shè)我們得到了一個(gè)有n個(gè)點(diǎn)值對的點(diǎn)值表達(dá)那我們可以確定唯一的一個(gè)次數(shù)界為n的多項(xiàng)式1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前八頁\總數(shù)五十一頁\編于十三點(diǎn)多項(xiàng)式乘法我們來探究一下如何用點(diǎn)值與插值來完成多項(xiàng)式乘法我們確定一組x，求得A與B的點(diǎn)值表達(dá)那我們可以得知C的點(diǎn)值表達(dá)通過插值運(yùn)算，我們可以知道多項(xiàng)式C的系數(shù)1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前九頁\總數(shù)五十一頁\編于十三點(diǎn)注意的地方設(shè)A與B的次數(shù)界均為n則C的次數(shù)界為2n則我們要找出2n個(gè)x來求點(diǎn)值表達(dá)否則不可以進(jìn)行插值運(yùn)算1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前十頁\總數(shù)五十一頁\編于十三點(diǎn)算法流程對于次數(shù)界均為n的多項(xiàng)式A與B1點(diǎn)值運(yùn)算構(gòu)造長度為2n的點(diǎn)值表達(dá)2逐點(diǎn)相乘計(jì)算出C的點(diǎn)值表達(dá)3插值運(yùn)算通過C的點(diǎn)值表達(dá)求出多項(xiàng)式C的每項(xiàng)系數(shù)1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前十一頁\總數(shù)五十一頁\編于十三點(diǎn)時(shí)間復(fù)雜度可以證明，若選取n個(gè)x計(jì)算點(diǎn)值與插值的時(shí)間復(fù)雜度均為O(N^2)本質(zhì)上沒有解決時(shí)間的問題但我們可以巧妙的選擇這些數(shù)來優(yōu)化時(shí)間復(fù)雜度。1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前十二頁\總數(shù)五十一頁\編于十三點(diǎn)PART2——N次單位復(fù)數(shù)根及其性質(zhì)1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前十三頁\總數(shù)五十一頁\編于十三點(diǎn)N次單位復(fù)數(shù)根n次單位復(fù)數(shù)根是滿足的復(fù)數(shù)w。n次單位復(fù)數(shù)根根恰好有n個(gè)，對于k=0,1,...,n-1，這些根是。為了解釋這個(gè)表達(dá)式，我們利用復(fù)數(shù)的指數(shù)形式的定義：下一頁圖說明n個(gè)單位復(fù)數(shù)根均勻地分布在以復(fù)平面的原點(diǎn)為圓心的單位半徑的圓周上。1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前十四頁\總數(shù)五十一頁\編于十三點(diǎn)N次單位復(fù)數(shù)根1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前十五頁\總數(shù)五十一頁\編于十三點(diǎn)性質(zhì)我們需要N次單位復(fù)數(shù)根我們首先來探究這些根的性質(zhì)1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前十六頁\總數(shù)五十一頁\編于十三點(diǎn)性質(zhì)1——主n次單位根我們稱為主n次單位根同時(shí)注意到，n次單位復(fù)數(shù)根都是經(jīng)過旋轉(zhuǎn)而得到的每次旋轉(zhuǎn)都是一定角度n次單位復(fù)數(shù)根可視為公比為主n次單位根的等比數(shù)列1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前十七頁\總數(shù)五十一頁\編于十三點(diǎn)性質(zhì)2——群的性質(zhì)因?yàn)樗酝普?點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前十八頁\總數(shù)五十一頁\編于十三點(diǎn)性質(zhì)3——消去引理&折半引理....消去引理：推論：折半引理：如果n>0為偶數(shù)，那么n次單位復(fù)數(shù)根的平方的集合就是n/2次單位復(fù)數(shù)根的集合。證明：可以知道的平方相同。1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前十九頁\總數(shù)五十一頁\編于十三點(diǎn)性質(zhì)4——求和引理求和引理：對于任意整數(shù)n≥1和不能被n整除的非負(fù)整數(shù)k，有等比數(shù)列求和所以注意k不能是n的倍數(shù)，否則分母為01點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前二十頁\總數(shù)五十一頁\編于十三點(diǎn)PART3——FFT及其關(guān)鍵算法1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前二十一頁\總數(shù)五十一頁\編于十三點(diǎn)DFT——離散傅里葉變換我們希望計(jì)算次數(shù)界為n的多項(xiàng)式在n次單位復(fù)數(shù)根處的值（共n個(gè)）接下來定義結(jié)果yy即為a的離散傅里葉變換（DFT）我們也可記為1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前二十二頁\總數(shù)五十一頁\編于十三點(diǎn)FFT——快速傅里葉變換大前提：n為2的整數(shù)冪（方便計(jì)算）利用復(fù)數(shù)單位根復(fù)數(shù)根的特殊性質(zhì)我們可以在時(shí)間內(nèi)計(jì)算出FFT利用了分治策略1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前二十三頁\總數(shù)五十一頁\編于十三點(diǎn)PART3.1——點(diǎn)值運(yùn)算1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前二十四頁\總數(shù)五十一頁\編于十三點(diǎn)分治策略如何求出單個(gè)數(shù)x的函數(shù)值A(chǔ)(x)？我們定義兩個(gè)新多項(xiàng)式觀察兩個(gè)多項(xiàng)式的特點(diǎn)1分別擁有奇數(shù)下標(biāo)的系數(shù)與偶數(shù)下標(biāo)的系數(shù)2次數(shù)界變?yōu)閚/2（縮小了一半）1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前二十五頁\總數(shù)五十一頁\編于十三點(diǎn)分治策略對于一個(gè)數(shù)x，求A(x)則根據(jù)上兩個(gè)多項(xiàng)式1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前二十六頁\總數(shù)五十一頁\編于十三點(diǎn)分治策略至此我們成功的轉(zhuǎn)換了問題原問題：求一個(gè)多項(xiàng)式A(x)在的函數(shù)值。現(xiàn)問題：求兩個(gè)多項(xiàng)式在的函數(shù)值。1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前二十七頁\總數(shù)五十一頁\編于十三點(diǎn)分治策略現(xiàn)問題：求兩個(gè)多項(xiàng)式在的函數(shù)值。根據(jù)折半引理，并不是n個(gè)不同的值，而是由n/2個(gè)n/2次單位復(fù)數(shù)根組成，而每個(gè)根恰好出現(xiàn)2次于是，我們遞歸地對n/2的多項(xiàng)式A[0](x)與A[1](x)在n/2個(gè)n/2次單位復(fù)數(shù)根進(jìn)行求值1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前二十八頁\總數(shù)五十一頁\編于十三點(diǎn)程序?qū)崿F(xiàn)1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前二十九頁\總數(shù)五十一頁\編于十三點(diǎn)我們關(guān)心的程序?qū)崿F(xiàn)問題1/2：規(guī)定程序遞歸出口3/4/12：定義主n次單位根，更新w值5/6/7/8：定義兩個(gè)多項(xiàng)式并遞歸求解13：返回DFT9/10/11：遞歸結(jié)束后的處理工作1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前三十頁\總數(shù)五十一頁\編于十三點(diǎn)遞歸結(jié)束后的處理工作10：11：1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前三十一頁\總數(shù)五十一頁\編于十三點(diǎn)FFT時(shí)間復(fù)雜度對于運(yùn)行時(shí)間有以下的遞歸式所以采用FFT，我們可以在O(nlgn)時(shí)間內(nèi)實(shí)現(xiàn)點(diǎn)值運(yùn)算（求出次數(shù)界為n的多項(xiàng)式在n次單位復(fù)數(shù)根處的值）。1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前三十二頁\總數(shù)五十一頁\編于十三點(diǎn)PART3.2——插值運(yùn)算1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前三十三頁\總數(shù)五十一頁\編于十三點(diǎn)矩陣乘積我們可以把點(diǎn)值運(yùn)算表示成一個(gè)矩陣方程我們把DFT寫成矩陣乘積y=Vna1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前三十四頁\總數(shù)五十一頁\編于十三點(diǎn)逆矩陣到此為止我們把插值運(yùn)算改寫成y與的逆矩陣的乘積1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前三十五頁\總數(shù)五十一頁\編于十三點(diǎn)逆矩陣定理：證明比較冗長。我們證明，其中In為n*n的單位矩陣?？紤]中(j,j')處的元素：如果j=j'，則此和為1否則根據(jù)求和引理，此和為01點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前三十六頁\總數(shù)五十一頁\編于十三點(diǎn)逆DFT給定矩陣，可以推導(dǎo)出：通過比較DFT的運(yùn)算我們可以看到，對FFT作以下修改就可以計(jì)算出逆DFT：把a(bǔ)與y互換，用，再把計(jì)算結(jié)果的每個(gè)元素除以n。于是我們也可以在O(nlgn)時(shí)間內(nèi)算出逆DFT。1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前三十七頁\總數(shù)五十一頁\編于十三點(diǎn)PART4——標(biāo)程演示1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前三十八頁\總數(shù)五十一頁\編于十三點(diǎn)程序?qū)崿F(xiàn)優(yōu)化因?yàn)橄駛未a中，賦值數(shù)組是不切實(shí)際的但我們發(fā)現(xiàn)，a中的應(yīng)用的系數(shù)是有規(guī)律的所以根據(jù)位移與起始點(diǎn)的不同來優(yōu)化FFT的實(shí)現(xiàn)1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前三十九頁\總數(shù)五十一頁\編于十三點(diǎn)程序?qū)崿F(xiàn)優(yōu)化(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點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前四十頁\總數(shù)五十一頁\編于十三點(diǎn)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;//定義復(fù)數(shù)類型，復(fù)數(shù)運(yùn)算procedureDft(vara:arr;s,t:longint);//a答案數(shù)組，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點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前四十一頁\總數(shù)五十一頁\編于十三點(diǎn)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同學(xué)提供的程序1點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前四十二頁\總數(shù)五十一頁\編于十三點(diǎn)#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點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前四十三頁\總數(shù)五十一頁\編于十三點(diǎn)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(bǔ)[]按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點(diǎn)值與插值2N次單位復(fù)數(shù)根3FFT4演示目前四十四頁\總數(shù)五十一頁\編于十三點(diǎn)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[