1、快速傅里葉變換FFTyL2013-9我們關(guān)心的問題u快速解決多項式乘法問題u衍生問題高精度乘法問題的描述u記一個多項式次數(shù)界為n的多項式A(x)u則u其中a為每一項的系數(shù)u注意最高次系數(shù)為n-110)(njjjxaxA問題的描述u兩個多項式相乘u我們記一般形式為uC的次數(shù)界為A與B次數(shù)界的和u普通的時間復(fù)雜度為O(n2)()()(xBxAxCPART1中心思想2 N次單位復(fù)數(shù)根3 FFT4 演示轉(zhuǎn)換思路u普通的相乘方法u提出概念：點值，插值ijjijibac02 N次單位復(fù)數(shù)根3 FFT4 演示點值u一個次數(shù)界為n的多項式A(x)的點值表達就是一個由n個點值對所組成的集合：u其中每一個x都不相

2、同，且uE.G.u多項式 的一個合法點值表達是),( ,),(),(111100nnyxyxyx)(kkxAy 23)(2xxxA)12, 5(),2 , 0(2 N次單位復(fù)數(shù)根3 FFT4 演示插值u插值運算是點值運算的逆運算u假設(shè)我們得到了一個有n個點值對點值對的點點值表達值表達u那我們可以確定唯一的一個次數(shù)界為n的多項式2 N次單位復(fù)數(shù)根3 FFT4 演示多項式乘法u我們來探究一下如何用點值與插值來完成多項式乘法u我們確定一組x，求得A與B的點值表達u那我們可以得知C的點值表達u通過插值運算，我們可以知道多項式C的系數(shù)),( ,),(),(:111100nnyxyxyxA),( ,),(

3、),(:111100nnyxyxyxB),( ,),(),(:111111000nnnyyxyyxyyxC2 N次單位復(fù)數(shù)根3 FFT4 演示注意的地方u設(shè)A與B的次數(shù)界均為nu則C的次數(shù)界為2nu則我們要找出2n個x來求點值表達u否則不可以進行插值運算2 N次單位復(fù)數(shù)根3 FFT4 演示算法流程u對于次數(shù)界均為n的多項式A與Bu1點值運算u構(gòu)造長度為2n的點值表達u2逐點相乘u計算出C的點值表達u3插值運算u通過C的點值表達求出多項式C的每項系數(shù)2 N次單位復(fù)數(shù)根3 FFT4 演示時間復(fù)雜度u可以證明，若選取n個xu計算點值與插值的時間復(fù)雜度均為O(N2)u本質(zhì)上沒有解決時間的問題u但我們可

4、以巧妙的選擇這些數(shù)來優(yōu)化時間復(fù)雜度。2 N次單位復(fù)數(shù)根3 FFT4 演示PART2N次單位復(fù)數(shù)根及其性質(zhì)1 點值與插值3 FFT4 演示N次單位復(fù)數(shù)根un次單位復(fù)數(shù)根是滿足 的復(fù)數(shù)w。un次單位復(fù)數(shù)根根恰好有n個，對于k=0, 1, . , n-1，這些根是 。為了解釋這個表達式，我們利用復(fù)數(shù)的指數(shù)形式的定義：u下一頁圖說明n個單位復(fù)數(shù)根均勻地分布在以復(fù)平面的原點為圓心的單位半徑的圓周上。1nwnike/2)sin()cos(uiueiu1 點值與插值3 FFT4 演示N次單位復(fù)數(shù)根28w18w08w78w68w58w48wii38w,N1210nnnnnwwww次單位復(fù)數(shù)根為我們記1 點值與

5、插值3 FFT4 演示性質(zhì)u我們需要N次單位復(fù)數(shù)根u我們首先來探究這些根的性質(zhì)1 點值與插值3 FFT4 演示性質(zhì)1主n次單位根u我們稱 為主n次單位根u同時注意到，n次單位復(fù)數(shù)根都是經(jīng)過旋轉(zhuǎn)而得到的u每次旋轉(zhuǎn)都是一定角度un次單位復(fù)數(shù)根可視為公比為主n次單位根的等比數(shù)列1nnww nininwww11 點值與插值3 FFT4 演示性質(zhì)2群的性質(zhì)u因為u所以u推論10nnnwwnkjnkjnknjnwwwwmod)( nknknwwknknww22)(1 點值與插值3 FFT4 演示11nnnww性質(zhì)3消去引理&折半引理.u消去引理：u推論：u折半引理：如果n0為偶數(shù)，那么n次單次單位復(fù)數(shù)根位

6、復(fù)數(shù)根的平方的集合就是n/2次單位復(fù)次單位復(fù)數(shù)根數(shù)根的集合。u證明：可以知道 的平方相同。kndkdnww122/ wwnn)2/(nknknww與222222/)()(knknnnknnknnknwwwwwwknknww2/2)(1 點值與插值3 FFT4 演示性質(zhì)4求和引理u求和引理：對于任意整數(shù)n1和不能被n整除的非負整數(shù)k，有u等比數(shù)列求和u所以u注意k不能是n的倍數(shù)，否則分母為0100)(njjknw011) 1 (11)(11)()(10knkknknnknnknnjjknwwwwww1 點值與插值3 FFT4 演示PART3FFT及其關(guān)鍵算法1 點值與插值2 N次單位復(fù)數(shù)根4 演

7、示DFT離散傅里葉變換u我們希望計算次數(shù)界為n的多項式u在n次單位復(fù)數(shù)根處的值（共n個）u接下來定義結(jié)果yuy即為a的離散傅里葉變換（DFT）u我們也可記為10)(njjjxaxA10)(njkjnjknkwawAy)(aDFTyn1 點值與插值2 N次單位復(fù)數(shù)根4 演示FFT快速傅里葉變換u大前提：n為2的整數(shù)冪（方便計算）u利用復(fù)數(shù)單位根復(fù)數(shù)根的特殊性質(zhì)u我們可以在 時間內(nèi)計算出uFFT利用了分治策略)lg(nnO)(aDFTn1 點值與插值2 N次單位復(fù)數(shù)根4 演示PART3.1點值運算1 點值與插值2 N次單位復(fù)數(shù)根4 演示分治策略u如何求出單個數(shù)x的函數(shù)值A(chǔ)(x)？u我們定義兩個新多

8、項式u觀察兩個多項式的特點u1分別擁有奇數(shù)下標(biāo)的系數(shù)與偶數(shù)下標(biāo)的分別擁有奇數(shù)下標(biāo)的系數(shù)與偶數(shù)下標(biāo)的系數(shù)系數(shù)u2次數(shù)界變?yōu)榇螖?shù)界變?yōu)閚/2（縮小了一半）（縮小了一半）12/12531112/224200)()(nnnnxaxaxaaxAxaxaxaaxA1 點值與插值2 N次單位復(fù)數(shù)根4 演示分治策略u對于一個數(shù)x，求A(x)u則根據(jù)上兩個多項式)()()(2120 xxAxAxAy1 點值與插值2 N次單位復(fù)數(shù)根4 演示分治策略u至此我們成功的轉(zhuǎn)換了問題u原問題：求一個多項式A(x)在 的函數(shù)值。u現(xiàn)問題：求兩個多項式 在 的函數(shù)值。110,nnnnwww)(A)(10 xxA與212120)

9、( ,)( ,)(nnnnwww1 點值與插值2 N次單位復(fù)數(shù)根4 演示分治策略u現(xiàn)問題：求兩個多項式 在 的函數(shù)值。u根據(jù)折半引理， 并不是n個不同的值，而是由n/2個n/2次單位復(fù)數(shù)根組成，而每個根恰好出現(xiàn)2次u于是，我們遞歸地對n/2的多項式A0(x)與A1(x)在n/2個n/2次單位復(fù)數(shù)根進行求值)(A)(10 xxA與212120)( ,)( ,)(nnnnwww212120)( ,)( ,)(nnnnwww1 點值與插值2 N次單位復(fù)數(shù)根4 演示程序?qū)崿F(xiàn)1 點值與插值2 N次單位復(fù)數(shù)根4 演示我們關(guān)心的程序?qū)崿F(xiàn)問題u1/2：規(guī)定程序遞歸出口u3/4/12：定義主n次單位根，更新w值

10、u5/6/7/8：定義兩個多項式并遞歸求解u13：返回DFTu9/10/11：遞歸結(jié)束后的處理工作1 點值與插值2 N次單位復(fù)數(shù)根4 演示遞歸結(jié)束后的處理工作u10：u11：)()()(212010knknknknkknkkwAwAwwAywyy)()()()()()2/(21)2/(2021)2/(201)2/(010)2/(nknnknnknnknknnknknknknkkknknkwAwAwwAwAwwAywyywyy1 點值與插值2 N次單位復(fù)數(shù)根4 演示FFT時間復(fù)雜度u對于運行時間有以下的遞歸式u所以采用FFT，我們可以在O(nlgn)時間內(nèi)實現(xiàn)點值運算（求出次數(shù)界為n的多項式在n

11、次單位復(fù)數(shù)根處的值）。)lg()()2/(2)(nnOnOnTnT1 點值與插值2 N次單位復(fù)數(shù)根4 演示PART3.2插值運算1 點值與插值2 N次單位復(fù)數(shù)根4 演示矩陣乘積u我們可以把點值運算表示成一個矩陣方程u我們把DFT寫成矩陣乘積y=Vna13210)1)(1()1(3)1(21)1(3963)1(2642132113210111111111nnnnnnnnnnnnnnnnnnnnnnnnnnaaaaawwwwwwwwwwwwwwwwyyyyy1 點值與插值2 N次單位復(fù)數(shù)根4 演示逆矩陣u到此為止我們把插值運算改寫成y與 的逆矩陣 的乘積1nnVyaaVy1nVnV1 點值與插值2

12、 N次單位復(fù)數(shù)根4 演示逆矩陣u定理：u證明比較冗長。u我們證明 ，其中In為n*n的單位矩陣。考慮 中(j,j)處的元素：u如果j=j，則此和為1u否則根據(jù)求和引理，此和為0./),(1nwkjVkjnn處元素為的nnnIVV11nnVV10)(101/)(/(nkjjknnkj knkjnj jnnnwwnwVV1 點值與插值2 N次單位復(fù)數(shù)根4 演示逆DFTu給定矩陣 ，可以推導(dǎo)出 ：u通過比較DFT的運算u我們可以看到，對FFT作以下修改就可以計算出逆DFT：u把a與y互換，用 ，再把計算結(jié)果的每個元素除以n。u于是我們也可以在O(nlgn)時間內(nèi)算出逆DFT。1nV)(1yDFTn1

13、01nkkjnkjwyna10njkjnjkwaynnww 代替11 點值與插值2 N次單位復(fù)數(shù)根4 演示PART4標(biāo)程演示1 點值與插值2 N次單位復(fù)數(shù)根3 FFT程序?qū)崿F(xiàn)優(yōu)化u因為像偽代碼中，賦值數(shù)組是不切實際的u但我們發(fā)現(xiàn)，a中的應(yīng)用的系數(shù)是有規(guī)律的u所以根據(jù)位移與起始點的不同u來優(yōu)化FFT的實現(xiàn)1 點值與插值2 N次單位復(fù)數(shù)根3 FFT程序?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 點值與插值2

14、 N次單位復(fù)數(shù)根3 FFTu typeuNode = recordux,y:double;uend;u arr = array0.200000 of Node;u operator + (a,b:Node) c:Node;u begin c.x := a.x + b.x;c.y := a.y + b.y;end;u operator - (a,b:Node) c:Node;u begin c.x := a.x - b.x;c.y := a.y - b.y;end;u operator * (a,b:Node) c:Node;u begin c.x := a.x * b.x - a.y * b.

15、y;c.y := a.x * b.y + a.y * b.x;end;u /定義復(fù)數(shù)類型，復(fù)數(shù)運算u procedure Dft(var a:arr;s,t:longint);/a答案數(shù)組，s起始量，t位移量u beginuif (n t) = 1 then exit;uDft(a,s,t + 1);Dft(a, s + 1 (t + 1) - 1 doubeginuj := s + i (t + 1);uwt := wi (t) * aj + 1 (t + 1) := aj - wt;uend;ufor i := 0 to n t - 1 do as + i 0) and(ansi = 0)

16、 do dec(i);ufor j := i downto 0 douwrite(ansj);uwriteln;uend.u/ 感謝感謝wwt同學(xué)提供的程序同學(xué)提供的程序1 點值與插值2 N次單位復(fù)數(shù)根3 FFTu#includeu#includeu#includeu#includeu#define N 50010/uusing namespace std;uconst double pi=acos(-1);uconst complex I(0,1);uconst double eps=1e-6;uchar aaN,bbN;uint ans4*N;/char ans4*N;!ucomplex

17、a4*N,b4*N;/an-1,an-2,a1,a0uint n;uvoid initial()uu int lenA=strlen(aa),lenB=strlen(bb);u n=max(lenA,lenB);u double t=log2(n);u int tt=int(t);u if(t-tt0)tt+;u n=1(tt+1);/double lengthu int i;u for(i=0; ilenA; i+)ai=aalenA-1-i-0;u while(in)ai+=0;u for(i=0; ilenB; i+)bi=bblenB-1-i-0;u while(in)bi+=0;u1

18、 點值與插值2 N次單位復(fù)數(shù)根3 FFTuvoid bitReverse(complex * a)uu for(int i=1; in-1; i+)u u int j=0;u for(int k=1,tmp=i; kn; j=(j1)|(tmp&1),k=1);u if(ji)swap(ai,aj);u uuvoid iterative_FFT(complex *a,int sig)uu bitReverse(a);u for(int m=2; m=n; m1;u for(int i=0; imh; i+)u u complex wi=exp(i*sig*pi/mh*I);u for(int

19、j=i; jn; j+=m)u u int k=j+mh;u complex u=aj,t=wi*ak;u aj=u+t;u ak=u-t;u u u u if(sig=-1)for(int i=0; in; i+)ai/=n;u1 點值與插值2 N次單位復(fù)數(shù)根3 FFTuvoid convolution()uu iterative_FFT(a,1);u iterative_FFT(b,1);u for(int i=0; in; i+)ai*=bi;/a*bu iterative_FFT(a,-1);uuvoid output()uu int k=0;u ans0=0;u for(int i=0; i=0)coutansk