我們知道,有兩種表示多項式的方法,即係數表示法和點值表示法。什麼是係數表示法?所謂的係數表示法,舉個例子如下圖所示,A(x)=6x^3 + 7x^2 - 10x + 9,B(x)=-2x^3+4x-5,則C(x)=A(x)*B(x)就是普通的多項式相乘的算法,係數與係數相乘,這就是所謂的係數表示法。
那麼,何謂點值表示法?。一個次數為n次的多項式A(x)的點值表示就是n個點值所形成的集合:{(x0,y0), (x1,y1),..., (xn-1,yn-1)}。其中所有xk各不相同,且當k=0,1,……,n-1時,有y(k)=A(xk)。
一個多項式可以由很多不同的點值表示,這是由於任意n個相異點x0,x1,...,xn-1組成的集合,都可以看做是這種點值法的表示方法。對於一個用係數形式表示的多項式來說,在原則上計算其點值表示是簡單易行的,我們只需要選取n個相異點x0,x1,...,xn-1,然後對k=0,1,...,n-1,求出A(xk),然後根據霍納法則,求出這n個點的值所需要的時間為O(n^2)。
稍後,你將看到,如果巧妙的選取xk的話,適當的利用點值表示可以使多項式的乘法可以在線性時間內完成。所以,如果我們能恰到好處的利用係數表示法與點值表示法的相互轉化,那麼我們可以加速多項式乘法(下面,將詳細闡述這個過程),達到O(n*logn)的最佳時間複雜度。
前面說了,當用係數表示法,即用一般的算法表示多項式時,兩個 n次多項式的乘法需要用 O(n^2)的時間才能完成。但採用點值表示法時,多項式乘法的運行時間僅為O(n)。接下來,咱們要做的是,如何利用係數表示法與點值表示法的相互轉化來實現多項式係數表示法的O(n*logn)的快速乘法。
第二節、多項式的快速乘法在此之前,我們得明白兩個概念,求值與插值。通俗的講,所謂求值(係數->點值),是指係數形式的多項式轉換為點值形式的表示方式。而插值(點值->係數)正好是求值的逆過程,即反過來,它是已知點值表示法,而要你轉換其為多項式的係數表示法(用n個點值對也可以唯一確定一個不超過n-1次多項式,這個過程稱之為插值)。而這個係數表示法與點值表示法的相互轉化,即無論是從係數表示法轉化到點值表示法所謂的求值,還是從點值表示法轉化到係數表示法所謂的插值,求值和插值兩種過程的時間複雜度都是O(n*logn)。
前面,我們已經說了,在多項式乘法中,如果用係數表示法表示多項式,那麼多項式乘法的複雜度將為O(n^2),而用點值表示法表示的多項式,那麼多項式的乘法的複雜度將是線性複雜度為O(n),即: 適當的利用點值表示可以使多項式的乘法可以在線性時間內完成。此時讀者可以發揮你的想像,假設我們以下面這樣一種過程來計算多項式的乘法(不過在此之前,你得先把兩個要相乘的多項式A(x)和B(x)擴充為次數或者說係數為2n次的多項式),我們將會得到我們想要的結果:
係數表示法轉化為點值表示法。先用係數表示法表示一個多項式,然後對這個多項式進行求值操作,即多項式從係數表示法變成了點值表示法,時間複雜度為O(n*logn);
點值表示法計算多項式乘法。用點值表示法表示多項式後,計算多項式的乘法,線性複雜度為O(n);
點值表示法轉化為係數表示法。最終再次將點值表示法表示的多項式轉變為係數表示法表示的多項式,這一過程,為O(n*logn)。
那麼,綜上上述三項操作,係數表示法表示的多項式乘法總的時間複雜度已被我們降到了O(n*logn+n+n*logn)=O(N*logN)。
如下圖所示,從左至右,看過去,這個過程即是完成多項式乘法的計算過程。不過,完成這個過程有兩種方法,一種就是前面第一節中所說的普通方法,即直接對係數表示法表示的多項式進行乘法運算,複雜度為O(n^2),它體現在下圖中得Ordinary multiplication過程。還有一種就是本節上文處所述的三個步驟:1、將多項式由係數表示法轉化為點值表示法(點值過程);2、 利用點值表示法完成多項式乘法;3、將點值表示法再轉化為係數表示法(插值過程)。三個步驟下來,總的時間複雜度為O(N*logN)。它體現在下圖中的Evaluation + Pointwise multiplication + Interpolation 三個合過程。
由上圖中,你也已經看到了。我們巧妙的利用了關於點值形式的多項式的線性時間乘法算法,來加快了係數形式表示的多項式乘法的運算速度。在此過程中,我們的工作最關鍵的就在於以O(n*logn)的時間快速把一個多項式從係數形式轉換為點值形式(求值,我們即稱之為FFT),然後再以O(n*logn)的時間從點值形式轉換為係數形式(插值,我們即稱之為FFT逆)。最終達到了我們以最終的係數形式表示的多項式的快速乘法為O(N*logN)的時間複雜度。好不令人心生快意。
對上圖,有一點必須說明,項w2n是2n次單位復根。且不容忽視的是,在上述兩種表示法即係數表示法和點值表示法相互轉換的過程中, 都使用了單位復根,才得以在O(n*logn)的時間內完成求值與插值運算(至於何謂單位復根,請參考相關資料。因為為了保證本文的通俗易懂性,無意引出一大堆公式或定理)。
第三節、FFT註:本節有相當部分文字來自算法導論中文版第二版以及維基百科。
在具體介紹FFT之前,我們得熟悉知道折半定理是怎麼一回事兒:如果n>0且n為偶數,那麼,n個n次單位復根的平方等於n/2個n/2個單位復根。我們已經知道,通過使用一種稱為快速傅立葉變換(FFT)的方法,可以在O(n*logn)的時間內計算出DFTn(a),而若採用直接的方法複雜度為O(n^2)。FFT就是利用了單位復根的特殊性質。
你將看到,FFT方法運用的策略為分治策略,所謂分治,即分而治之。兩個多項式A(x)與B(x)相乘的過程中,FFT用A(x)中偶數下標的係數與奇數下標的係數,分別定義了兩個新的次數界為n/2的多項式A[0](x)和A[1](x):
A[0](x) =a0 +a2x +a4x2 +··· +an-2xn/2-1,
A[1](x) =a1 +a3x +a5x2 +··· +an-1xn/2-1.
注意,A[0]包含了A中所有偶數下標的係數(下標的相應二進位數的最後一位為0),而A[1]包含A中所有奇數下標的係數(下標的相應二進位數的最後一位為1)。有下式:
這樣,求A(x)在處的問題就轉換為:
1)求次數界為n/2的多項式A[0]與A[1]在點
的值
2)將上述結果進行組合。
下面,我們用N次單位根WN來表示。
為了簡單起見,我們下面設待變換序列長度n = 2r。根據上面單位根的對稱性,求級數時,可以將求和區間分為兩部分:
Fodd(k) 和 Feven(k)是兩個分別關於序列
奇數號和偶數號序列N/2點變換。由此式只能計算出yk的前N/2個點,對於後N/2個點,注意Fodd(k) 和 Feven(k) 都是周期為N/2的函數,由單位根的對稱性,於是有以下變換公式:
。
這樣,一個N點變換就分解成了兩個N/2點變換,照這樣可繼續分解下去。這就是庫利-圖基快速傅立葉變換算法的基本原理。此時,我們已經不難分析出此時算法的時間複雜度將為O(NlogN)。
ok,沒想到,本文之中還是出現了這麼多的公式(沒辦法,搞這個FFT就是個純數學活兒。之前買的一本小波與傅立葉分析基礎也是如此,就是一本數學公式和定理的聚集書。不過,能看懂更好,實在無法弄懂也只權當做個稍稍了解)。
好了,下面,咱們來簡單理解下FFT中的蝶形算法,本文將告結束。如下圖所示:
有人這樣解釋蝶形算法:對於N(2的x次方)點的離散信號,把它按索引位置分成兩個序列,分別為0,2,4,....,2K(記為A)和1,3,5,....,2K-1(記為B),由傅立葉變換可以推出N點的傅立葉變換前半部分的結果為A+B*旋轉因子,後半部分的結果為A-B*旋轉因子。於是求N點的傅立葉變換就變成分別求兩個N/2點序列的傅立葉變換,對每一個N/2點的序列,遞歸前面的步驟,直到只有兩點的序列,就只變成簡單的加減關係了。把這些點的加減關係用線連接,看上去就是個蝶形。ok,更多可參考算法導論第30章。
舉一個例子,我們知道,4X4的矩陣運算如果按常規算法的話仍要進行64次乘法運算和48次加法,這將耗費較多的時間,於是在H.264中,有一種改進的算法(蝶形算法)可以減少運算的次數。這種矩陣運算算法構造非常巧妙,利用構造的矩陣的整數性質和對稱性,可完全將乘法運算轉化為加法運算。如下圖所示: