快速傅立葉變換

多項式

　　對於多項式$ f\left(x\right)=\sum_{i=0}^{|f|}{f_ix^i} $，其中|f|表示多項式的階數，fi表示多項式f中x^i的係數。

　　多項式的加法定義爲$ c\left(x\right)=a\left(x\right)+b\left(x\right)=\sum_{i=0}^{\max\left(|a|,|b|\right)}{\left(a_i+b_i\right)x^i} $，即$ c_k=a_k+b_k $。

　　多項式的乘法定義爲$ c\left(x\right)=a\left(x\right)\cdot b\left(x\right)=\sum_{k=0}^{|a|+|b|}{\left(\sum_{i+j=k}^{}{a_ib_j}\right)x^k} $，即$ c_k=\sum_{i+j=k}^{}{a_ib_j} $。

　　顯然要計算兩個多項式a(x)，b(x)的乘積，程序的時間複雜度爲O(|a||b|)。

naive(a, b)
    c = 0
    for(i = 0; i <= |a|; i++)
        for(j = 0; j <= |b|; j++)
            c[i + j] = c[i + j] + a[i] * b[j]

　　在多項式階數超過10w的時候，這個方法就徹底頂不住了。不過幸虧還有不少加快多項式乘運算速度的算法，而快速傅立葉變換就是其中之一。

　　先了解一下多項式的其它操做的時間複雜度：多項式的乘法雖然很慢，可是求解一個多項式f在x=x0的時候的取值f(x0)是能夠在O(|f|)時間複雜度內作到的。以及多項式的加法a(x)+b(x)也能夠在O(max(|a|,|b|))的時間複雜度內作到。

　　再瞭解一下多項式的表示方法：多項式的表示方法基本有兩種，一種是經過係數序列(f0,f1,...,fk)來表示一個k階多項式f，這種方法稱爲係數表示法，還有一種就是點值表示法，即用k+1個不一樣的點來表示一個k階多項式。係數表示法你們都很瞭解，下面說一下點值表示法。

　　在代數中，說明過n個不一樣的點能夠惟一肯定一個k階多項式，其中k<n。而若是從n個點還原出原來的多項式，則須要使用到插值算法，著名的插值算法有拉格朗日插值法和牛頓插值法，這裏很少加介紹，兩者的時間複雜度均爲O(n^2)。

　　點值表示法很是適合用於計算多項式乘法，對於多項式乘法c(x)=a(x)*b(x)，假設咱們已經確認了|a|+|b|+1個不一樣的x值x0,x1,...，且分別計算出了a(x0),b(x0),a(x1),b(x1),...，那麼咱們就得知點(xi,a(xi)b(xi))是多項式c上的點，而這組點的數目爲|a|+|b|+1>|c|，故c被這組點惟一確認，在前提下多項式乘法能夠以O(|c|)的時間複雜度運行。

　　固然點值表達式的前提並很差知足，咱們每每須要先經過插值取回多項式，以後再計算在額外的點的多項式值，以後再利用點值乘法算出新的多項式的點值表達式。這整個過程的時間複雜度爲O(|c|^2)。

　　咱們設n爲大於等於c長度的最小2的冪次（即n=2^k>=|c|>2^(k-1)），在運算多項式乘法前，咱們先將a與b經過前面補0將長度擴充到n，以後再運行多項式乘法。下面咱們說明如何在O(nlog2n)的時間複雜度內將長度爲n的多項式從係數表達式轉換爲點值表達式，並在O(nlog2n)的時間複雜度內將n個不一樣的點插值會係數表達式的多項式，而這一算法就是快速傅立葉變換，很顯然這一過程的時間複雜度爲O(nlog2n+n+nlog2n)=O(nlog2n)。

單位複數根

　　首先咱們要謹慎地選取n個點的x值。在複數域中，n次單位複數根是知足w^n=1的全部複數w。由歐拉公式$ e^{iu}=\cos\left(u\right)+i\sin\left(u\right) $可知n次複數根分別爲複數$ e^{\left(2\pi k/n\right)i} $，其中k分別取值0,1,...,n-1的，咱們記爲w(n,0),w(n,1),...,w(n,n-1)。很顯然w(n,i)w(n,j)=w(n,i+j)，因爲w(n,0)=w(n,n)=1，故咱們得知w(n,i)=w(n,n-i)=w(n,i-n)。下面說明這n個複數的有趣性質：

　　消去引理：w(dn,dk)=w(n,k)。

　　證實：略

　　折半引理：2n次單位複數根的平方組成的集合與n次單位複數根組成的集合相同。

　　證實：對於0<=k<n，$ w\left(2n,k\right)^2=\left(e^{\left(2\pi k/2n\right)i}\right)^2=e^{\left(2\pi k/n\right)i}=w\left(n,k\right) $。

　　　　　對於n<=k<2n，因爲w(2n,n)=-1(看歐拉公式），故$ w\left(2n,k\right)^2=\left(-w\left(2n,k-n\right)\right)^2=w\left(2n,k-n\right)^2=w\left(n,k-n\right) $。

　　求和引理：對於任意n>=1和不能被n整除的非負整數k，有$ \sum_{j=0}^{n-1}{w\left(n,k\right)^j}=0 $。

　　證實：$ \sum_{j=0}^{n-1}{w\left(n,k\right)^j}=\frac{w\left(n,k\right)^n-1}{w\left(n,k\right)-1}=\frac{1-1}{w\left(n,k\right)-1}=0 $。

離散傅立葉變換DFT

　　對於一個長度爲n的多項式f(x)，其在n個n次單位複數根w(n,0),w(n,1),...,w(n,n-1)上的取值組成的序列y=(f(w(n,0)),f(w(n,1)),...,f(w(n,n-1)))稱爲f的離散傅立葉變換（DFT），也記做y=DFT(f)。很顯然DFT(f)能夠惟一肯定f。

　　要計算DFT(f)，咱們能夠分治策略。

　　咱們將f切分爲長度爲n/2的兩個多項式even和odd，其中even中僅包含多項式的偶數項係數，而odd中僅包含奇數項係數：

　　even=f₀*x^0+f₂*x^1+f₄*x^2+...+f_n-2x^(n/2-1)

　　odd=f₁*x^0+f₃*x^1+...+f_n-1x^(n/2-1)

　　而很顯然f=even(x^2)+x*odd(x^2)。所以要計算f在w(n,0),w(n,1),...,w(n,n-1)上的取值，只須要計算even和odd在w(n,0)^2,w(n,1)^2,...,w(n,n-1)^2上的取值便可。由折半引理可知{w(n,0)^2,w(n,1)^2,...,w(n,n-1)^2}={w(n/2,0),w(n/2,1),...,w(n/2,n/2-1)}，故咱們實際上要計算的僅爲DFT(even)和DFT(odd)。在獲得DFT(even)和DFT(odd)後，僅需使用O(n)的時間複雜度便可算出DFT(f)。

　　咱們記T(f)表示DFT(f)的時間複雜度，則T(f)=O(n)+2T(f/2)=...=O(kn)+(2^k)*T(f/(2^k))=...=O(nlog2(n))。

　　

DFT(f)
    n = f.length
　  if(n == 1)
　　　　return f[0]

    y = empty-array

    even = 0
    odd = 0
    for(i = 0; i < n / 2; i++)
        even[i] = f[i * 2]
        odd[i] = f[i * 2 + 1]
    yEven = DFT(even)
    yOdd = DFT(odd)
    
    wn1 = e^((2*PI/n)i)
    w = 1
    for(i = 0; i < n / 2; i++)
        y[i]= yEven[i] + w * yOdd[i]
        y[i + n / 2] = yEven[i] - w * yOdd[i]
        w = w * wn1

    return y

　　因爲計算機計算正弦和餘弦函數要花費大量的時間，所以能夠將wn1做爲參數傳入，而在計算DFT(even)時，將wn1*wn1做爲參數傳入便可（w(n,1)^2=w(n/2,1)），這樣能夠節省時間。

離散傅立葉逆變換IDFT

　　離散傅立葉逆變換IDFT將多項式從點值表達式轉換爲係數表達式。即IDFT(DFT(c))=c。觀察下面等式：

$$ \left[\begin{array}{c} y_0\\ y_1\\ \vdots\\ y_{n-1} \end{array}\right]=\left[\begin{matrix} w\left(n,0\right)^0 & w\left(n,0\right)^1 &\cdots & w\left(n,0\right)^{n-1}\\ w\left(n,1\right)^0 & w\left(n,1\right)^1 &\cdots & w\left(n,1\right)^{n-1}\\ \vdots &\vdots &\ddots &\vdots\\ w\left(n,n-1\right)^0 & w\left(n,n-1\right)^1 &\cdots & w\left(n,n-1\right)^{n-1} \end{matrix}\right]\left[\begin{array}{c} c_0\\ c_1\\ \vdots\\ c_{n-1} \end{array}\right] $$

等式左邊爲DFT(c)，而作右邊的向量爲c，方陣爲範德蒙德矩陣，因爲w(n,0),...,w(n,n-1)互不相同（歐拉公式），故方陣可逆。咱們將等式簡記y=Mc。記IM=M^(-1)。

　　定理：IM的j'行j列元素爲IMj'j=w(n,-j'j)/n

　　證實：記P=IM·M，顯然$ P_{j'j}=\sum_{k=0}^{n-1}{w\left(n,-j'k\right)/n\cdot w\left(n,kj\right)}=\frac{1}{n}\sum_{k=0}^{n-1}{w\left(n,k\right)^{j-j'}} $。當j等於j'時，Pj'j=n/n=1，而當j不等於j'時，由求和引理可知Pj'j=0。故P是單位矩陣，所以IM=M^(-1)。

 　　如今咱們能夠利用c=IM·y來計算c了。觀察矩陣下隱含的等式關係：

$$ c_j=\sum_{k=0}^{n-1}{y_k\cdot w\left(n,-jk\right)/n}=\frac{1}{n}\sum_{k=0}^{n-1}{y_kw\left(n,n-j\right)^k} $$

這給了咱們一個啓發，c和DFT(y)/n中包含的值是相同的，只是順序不一樣而已。所以咱們能夠利用DFT以及一些常數時間的操做實現IDFT。

IDFT(y)
    c = 0
    n = y.length
    
    dftY = DFT(y)
    
    c[0] = dftY[0] / n
    for(i = 1; i < n; i++)
        c[i] = dftY[n - i] / n

    return c

　　IDFT和DFT共享相同的時間複雜度O(nlog2n)。

 快速傅立葉變換FFT

　　快速傅立葉變換利用DFT和IDFT計算兩個多項式a(x)和b(x)的乘積。

　　FFT的第一步首先是找到一個合適的二次冪n，並將a和b經過添加0係數項的方式擴展到長度n。以後計算DFT(a)和DFT(b)，在利用點乘計算DFT(c)=DFT(a)·DFT(b)。最後利用IDFT從點值表達式DFT(c)復原出多項式c，即c=FFT(a,b)=IDFT(DFT(a)·DFT(b))。

　　顯然FFT的時間複雜度爲DFT和IDFT的總和，也是O(nlog2n)。而因爲n<2*|c|=2*(|a|+|b|)，所以咱們能夠忽略n與|c|之間的偏差，即時間複雜度能夠寫做O(|c|log2|c|)。是至關優秀的時間複雜度。

FFT(a, b)
    n = 1
    while(n < |a| + |b|)
        n = n * 2
    extend(a, n)
    extend(b, n)
    return IDFT(DFT(a)·DFT(b))

 非遞歸版本DFT

　　實際中，若是你直接使用上面的快速傅立葉變換，你將會發現性能至關不理想。咱們仔細觀察一下FFT的流程，咱們能發現咱們總共作了三次DFT計算（其中一次在調用IDFT中），和一次點乘運算，其中點乘運算是線性時間複雜度，而DFT的時間複雜度爲O(nlogn)。若是咱們能優化DFT，那麼FFT的性能將獲得最大幅度的提高。

　　那麼DFT的缺陷在哪裏，咱們很容易發現每次咱們都須要創建一個多項式y，以保存最終結果，咱們徹底能夠在FFT中複製多項式a，b，以後在DFT中將結果直接寫入到f中並做爲返回值。以後咱們還能夠發現，每次用f調用DFT，須要建立兩個子多項式odd和even，其長度總和爲|f|。所以咱們總共分配的多項式長度爲O(nlog2n)，這是至關大的空間複雜度。而若是咱們能不用分配even和odd，那麼天然能達到優化DFT的結果。

　　觀察對於多項式(a0,a1,a2,a3,a4,a5,a6,a7)調用DFT的過程：

(a0,a1,a2,a3,a4,a5,a6,a7)

|　　　　　　　　　　\

(a0,a2,a4,a8)　　　　(a1,a3,a5,a7)

|　　　　\　　　　　　|　　　　\

(a0,a4)　(a2,a8)　　(a1,a5)　　(a3,a7)

|　　\　　|　　\　　  |　　\　　  |　　\

a0　a4　a2　a8　　a1　a5　　a3　a7

　　若是咱們能將多項式(a0,a1,a2,a3,a4,a5,a6,a7)在DFT的一開始轉換爲(a0,a4,a2,a8,a1,a5,a3,a7)，那麼咱們就能夠用非遞歸的方式實現DFT，一開始將長度爲1的兩個相鄰多項式合併，以後將長度爲2的兩個相鄰多項式合併，再將長度爲4兩個相鄰多項式合併。

　　觀察序列，能發現如下規律，下標二進制最低位爲0的排在二進制最低位爲1的以前。若是最低位相同，則用相同邏輯判斷次低位。容易發現這至關於比較兩個整數二進制(總共log2n位)逆序的大小。所以若是咱們能一開始處理獲得n個數值的二進制序列，以後利用逆序對多項式各位進行重排列，那麼就能實現上述的非遞歸版本的DFT。

　　而要計算全部長度爲m的二進制序列的全部逆序，能夠用下面的方法：

reverse(m)
    n = 2^m
    r = int[n]
    r[0] = 0
    for(i = 1; i < n; i++)
        r[i] = (r[i >> 1] >> 1) | ((1 & i) << (m -1))
    return r

其中咱們使用了一個聰明的想法，i的二進制逆序，咱們能夠先將i右移1位獲得j，因爲j<i，咱們保證r[j]已經獲得，而r[j]做爲j的逆序，等價於將i的逆序r[i]左移了一位，咱們經過對r[j]再右移一位就能夠獲得r[i]的後面m-1位，可是首位始終爲0。考慮到首位正好是i的最低位，咱們能夠將i的最低位向左移動m-1位於r[j]>>1左位或處理，從而獲得r[i]。很容易得出reverse的時間複雜度爲O(n)。

　　最後給出非遞歸版本的DFT：

DFT(p, m)
    r = reverse(m)
    n = 2^m

    for(i = 0; i < n; i++)
        if(r[i] > i)
            swap(p, i, r[i])
    
    for(d = 0; d < m; d++)
        s = 2^d
        w1 = complex(cos(PI/s),sin(PI/s))
        for(i = 0; i < n; i += 2*s)
            w = 1
            for(j = 0; j < s; j++)
                a = i + j
                b = a + s
                t = w *  p[b]
                q = p[a]
                p[b] = q - t
                p[a] = q + t
                w = w * w1

其中r=reverse(m)這個步驟咱們能夠提早執行。以及w1 = complex(cos(PI/s),sin(PI/s))中因爲涉及到了cos和sin的計算，因爲s始終是2的冪，所以可能的值很是少，也能夠一開始就計算好。整個DFT中咱們沒有分配額外的空間，所以空間複雜度能夠認做爲O(n)，而時間複雜度不變爲O(nlog2n)。