快速傅裡葉變換（FFT）（下）__快速傅裡葉變換

關於學習FFT演算法的資料個人最推薦的還是演算法導論上的第30章(第三版), 多項式與快速傅裡葉變換, 基礎知識都講得很全面。

FFT演算法基本概念：

FFT(Fast Fourier Transformation)即快速傅裡葉變換, 是離散傅裡葉變換的加速演算法, 可以在O(nlogn)的時間裡完成DFT, 利用相似性也可以在同樣複雜度的時間裡完成逆DFT。
DFT(Discrete Fourier Transform) 即離散傅裡葉變換, 這裡主要就是多項式的係數向量轉換成點值表示的過程。
在ACM-ICPC競賽中, FFT演算法常被用來為多項式乘法加速, 即在O(nlogn)複雜度內完成多項式乘法, 當然實際應用不僅僅限於這些, 時常 出現 需要構造多項式相乘來進行計數的問題, 也需要用FFT演算法來解決, 相關的幾個問題在本文中也會提及。

FFT演算法需要的基礎數學知識：
多項式相關：
多項式相關的定義：一個以x為變數的多項式定義在代數域F上, 將函數A(x)表示為形式和：

稱為係數, 所有係數屬於代數域F, 如果一個多項式的最高非零係數是, 那麼成這個多項式的次數是k, 記做degree(A)=k, 任何一個嚴格大於一個多項式次數的整數都是該多項式的次數界.
關於多項式的加法和乘法相信看到這篇部落格的讀者都會最基本的中學的演算法, 在計算的時候, 如果採用傳統的中學的計算方法, 多項式加法的時間複雜度是O(n), 乘法的時間複雜度是O(n^2) (n是兩個多項式A和B的次數)。

多項式的表示:
在平常的學習中, 最常見的是多項式的係數表達方式, 次數界為n的多項式的係數表達為一個由係數組成的向量。 但是多項式還有一個比較常用的表示方法, 即多項式的點值表達方式。

一個次數界為n的多項式的點值表達就是一個由n個點對組成的集合

使得對任意的整數k=0,1,..n−1，各不相同, 且=A()。

關於點值表達的正確性證明：
對於任意n個點值對組成的集合, 如果存在一個次數界為n的多項式A(x)過這n個點, 那麼：

最左邊這個n*n的矩陣稱為範德蒙德矩陣, 可以用數學歸納法證明它的行列式值為：。

當兩兩不相同時明顯這個行列式的值不為0, 該矩陣可逆, 於是存在唯一解, 所以多項式的點值表達是合理的。
相應的通過n個點的座標直接確定多項式各個係數的值的方法是存在的, 感興趣的讀者可以查詢拉格朗日公式的相關資料, 利用拉格朗日插值公式可以在O(n^2)的時間複雜度內得到多項式的係數表達, 這個也是演算法導論中的一個習題

點值表達方式下多項式的乘法, 不難發現, 如果多項式A(x),B(x)的點值表示分別是：

和

那麼如果多項式C(x)=A(x)B(x), 那麼C(x)的點值表達是：

卷積：
對於兩個多項式的係數向量和, 兩個多項式相乘得到的多項式的係數向量滿足，稱係數向量c是輸入向量a和b的卷積, 記作c=a⊗b。

簡單的多項式乘法的計算方法中, 每一個多項式的係數都通過係數表示方式下卷積的方式來進行計算, 時間複雜度是O(n^2), 但是FFT是先將多項式的從係數標記法轉換成點值標記法(可以在O(nlogn)的時間複雜度下完成, 也就是加速的DFT變換, 然後在點值標記法下進行乘積計算, 在O(n)的時間複雜度內得到結果的點值標記法, 然後進行逆DFT變換, 在O(n*logn)的時間複雜度下完成逆DFT變換得到係數標記法。

而要理解DFT, 則需要一定複數上的數學知識：

複數相關的基礎知識：

單位複數根：n次單位複數根指的是滿足的所有複數ω, n次單位複數根剛好有n個, 他們是, 其中i是複數單位, k=0,1,2...n−1, 在複平面上這n個根均勻的分布在半徑為1的圓上, 關於複數指數的定義如下：

其中倍稱為主n次單位根(這個定義好像接下來沒用到)
關於複數根的幾個定理和引理：

消去引理： 對任何整數 有

證明： 

一個推論： 對任意偶數 n > 0 有

證明：設n = 2*k那麼 

折半引理：如果n > 0是偶數, 那麼n個n次單位複數根的平方的集合就是n/2個n/2次單位複數根的集合

證明：實際上這個引理就是證明了

折半引理對於採用分治對多項式係數向點值表達的轉換有很大作用, 保證了遞迴的子問題是原問題規模的一半。

求和引理：對任意整數n≥1和不能被n整除的非負整數k, 有

這個問題通過等比數列求和公式就可以得到：

DFT與FFT, 以及逆DFT：

DFT:
在DFT變換中, 希望計算多項式A(x)在複數根處的值, 也就是求

稱向量y=(y0,y1,...,yn−1)是係數向量a=(a0,a1,...,an−1)的離散傅裡葉變換, 記為

FFT:

直接計算DFT的複雜度是O(n^2)， 而利用複數根的特殊性質的話, 可以在O(n*logn)的時間內完成, 這個方法就是FFT方法, 在FFT方法中採用分治策略來進行操作, 主要利用了消去引理之後的那個推論。

在FFT的策略中, 多項式的次數是2的整數次冪, 不足的話再前面補0, 每一步將當前的多項式A(x), 次數是2的倍數, 分成兩個部分：


於是就有了

那麼我們如果能求出次數界是n/2的多項式和在n個n次單位複數根的平方處的取值就可以了, 即在

處的值, 那麼根據折半引理, 這n個數其實只有n/2個不同的值, 也就是說, 對於每次分出的兩個次數界n/2的多項式, 只需要求出其n/2個不同的值即可, 那麼問題就遞迴到了原來規模的一半, 也就是說如果知道了兩個子問題的結果, 當前問題可以在兩個子問題次數之和的複雜度內解決, 那麼這樣遞迴問題的複雜度將會是O(nlogn)的, 用a=(a0,a1,...,an−1)表示係數向量, y=(y0,y1,...,yn−1)表示離散變換之後的向量, 這裡給出將算導上的代碼翻譯出來的C++代碼實現(以解決演算法第三版大論第三十章的一個習題, 求(0, 1, 2, 3)的DFT為例)：

#include<bits/stdc++.h>using namespace std;const double eps(1e-8);typedef long long lint; const double PI = acos(-1.0);/* * 這是一個遞迴實現FFT的測試, 測試習題中求DFT(0, 1, 2, 3) */ struct Complex{    double real, image;    Complex(double _real, double _image)    {        real = _real;        image = _image;    }    Complex(){}}; Complex operator + (const Complex &c1, const Complex &c2){    return Complex(c1.real + c2.real, c1.image + c2.image);} Complex operator - (const Complex &c1, const Complex &c2){    return Complex(c1.real - c2.real, c1.image - c2.image);} Complex operator * (const Complex &c1, const Complex &c2){    return Complex(c1.real*c2.real - c1.image*c2.image, c1.real*c2.image + c1.image*c2.real);} Complex* RecursiveFFT(Complex a[], int n)//n表示向量a的維數{    if(n == 1)        return a;    Complex wn = Complex(cos(2*PI/n), sin(2*PI/n));    Complex w = Complex(1, 0);    Complex* a0 = new Complex[n >> 1];    Complex* a1 = new Complex[n >> 1];    for(int i = 0; i < n; i++)        if(i & 1) a1[(i - 1) >> 1] = a[i];        else a0[i >> 1] = a[i];    Complex *y0, *y1;    y0 = RecursiveFFT(a0, n >> 1);    y1 = RecursiveFFT(a1, n >> 1);    Complex* y = new Complex[n];    for(int k = 0; k < (n >> 1); k++)    {        y[k] = y0[k] + w*y1[k];        y[k + (n >> 1)] = y0[k] - w*y1[k];        w = w*wn;    }    return y;} int main(){    Complex* a = new Complex[10];    a[0] = Complex(0, 0);    a[1] = Complex(1, 0);    a[2] = Complex(2, 0);    a[3] = Complex(3, 0);    Complex* ans = new Complex[10];    ans = RecursiveFFT(a, 4);    for(int i = 0; i < 4; i++)        cout<<"("<<ans[i].real<<")"<<"+"<<"("<<ans[i].image<<")"<<"i"<<endl;        return 0;}

可以得到係數向量(0, 1, 2, 3)經過DFT之後是(6, -2-2i, -2, -2+2i)

在上面這個視線中需要注意的就是RecursiveFFT中的y[k] = y0[k] + w*y1[k];和y[k + (n >> 1)] = y0[k] - w*y1[k];的變化, 正是利用了對於偶數的n， 有成立, 得到數組y的所有值

另外有一個概念性的東西： 和都成立, 在這其中將稱為旋轉因子

在進行了DFT操作之後, 將多項式的係數表達成功轉為了點值表達 接下來是如何在O(nlogn)的時間內完成逆DFT的問題, 通過逆DFT將點值表達還原為係數表達

逆DFT:

根據DFT得到的向量y和係數向量a之間的關係, 可以用矩陣乘積的形式來表達他們之間的關係, 即, 也就是

那麼要將DFT變化得到的向量y還原成向量a的話, 只需要用的逆矩陣乘上向量y即可, 這裡需要用到一個定理

在矩陣中, 不難發現對於任意的, 那麼可以找到這樣一個矩陣, 對於任意的,的 ( j , k )出的元素為 ，是矩陣的逆矩陣。

證明如下：要證明這兩個矩陣互逆, 證明其積為單位矩陣即可, 考慮兩個矩陣的乘積在( j , j′ )出的元素, 可以發現這個元素是

當j′=j 時, 這個和是1, 否則根據求和引理, 這個和是0, 故這兩個矩陣互逆

那麼根據這個逆矩陣可以發現要計算的話, 有關係式比較這個式子和之前DFT裡面y和a的關係式子, 可以發現只需要用替換掉即可, 最後結果需要除以n, 所以計算逆DFT的方法和計算DFT和相似, 都可以在O(nlogn)的時間複雜度內解決。

卷積定理：
對任意兩個長度為n的向量a和b, 其中n是2的冪, 有

其中向量a和b用0填充, 使其長度達到2n, 並用⋅表示兩個2n個元素組成向量的點乘(也就是每一維上的數相乘)
這個式子實際上就是多項式的係數表達在乘法時進行的卷積運算得到的結果, 等同於通過將其係數進行DFT變換變成點值表達之後相乘再換回來的過程

關於FFT演算法的迭代實現：

在遞迴實現DFT過程的FFT演算法中, 我們每次將係數向量a分成兩個部分利用折半引理來降低計算的規模, 可以發現在每次分組當中他們滿足這樣一個完全二叉樹的分組(n是2的冪)：

FFT遞迴係數分組

通過上圖的流程可以看出, 最後一層的子節點下標的順序實際上就是其下標轉換成二進位串的倒序的字串按照字典序排列的順序。

比如a0~a7得到的序列a0, a4, a2, a6, a1, a5, a3, a7下標的二進位是000, 100, 010, 110, 001, 101, 011, 111對應的串的倒序是000, 001, 010, 011, 100, 101, 110, 111這個倒序的二進位剛好是0, 1, 2, 3, 4, 5, 6, 7的二進位表示, 於是我們可以在O(nlogn)的複雜度內得到做下面一層的下標順序, 然後可以根據子節點的結果向上迭代得到父親結點的值, 這樣計算的話直接避免了遞迴, 如果直接遞迴的話在一些OJ上可能會造成爆棧的錯誤, 所以還是採用迭代的方式進行比較好。

關於迭代形式的FFT演算法, C++代碼實現如下(同樣以求(0, 1, 2, 3)的DFT變換為例)：

這段代碼的話同時也進行了逆DFT, DFT和逆DFT的過程相似, 加上一個標記判斷當前執行的是哪一種就行了。

#include<bits/stdc++.h>using namespace std;const double eps(1e-8);typedef long long lint; const double PI = acos(-1.0); struct Complex{    double real, image;    Complex(double _real, double _image)    {        real = _real;        image = _image;    }    Complex(){}}; Complex operator + (const Complex &c1, const Complex &c2){    return Complex(c1.real + c2.real, c1.image + c2.image);} Complex operator - (const Complex &c1, const Complex &c2){    return Complex(c1.real - c2.real, c1.image - c2.image);} Complex operator * (const Complex &c1, const Complex &c2){    return Complex(c1.real*c2.real - c1.image*c2.image, c1.real*c2.image + c1.image*c2.real);} int rev(int id, int len){    int ret = 0;    for(int i = 0; (1 << i) < len; i++)    {        ret <<= 1;        if(id & (1 << i)) ret |= 1;    }    return ret;} //當DFT= 1時是DFT, DFT = -1則是逆DFTComplex* IterativeFFT(Complex* a, int len, int DFT)//對長度為len(2的冪)的數組進行DFT變換{    Complex* A = new Complex[len];//用A數組儲存數組a分組之後新的順序    for(int i = 0; i < len; i++)        A[rev(i, len)] = a[i];    for(int s = 1; (1 << s) <= len; s++)    {        int m = (1 << s);        Complex wm = Complex(cos(DFT*2*PI/m), sin(DFT*2*PI/m));        for(int k = 0; k < len; k += m)//這一層結點的包含數組元素個數都是(1 << s)        {            Complex w = Complex(1, 0);            for(int j = 0; j < (m >> 1); j++)//折半引理, 根據兩個子節點計算父親節點            {                Complex t = w*A[k + j + (m >> 1)];                Complex u = A[k + j];                A[k + j] = u + t;                A[k + j + (m >> 1)] = u - t;                w = w*wm;            }        }    }    if(DFT == -1) for(int i = 0; i < len; i++) A[i].real /= len, A[i].image /= len;    return A;} int main(){    Complex* a = new Complex[4];    a[0] = Complex(0, 0);    a[1] = Complex(1, 0);    a[2] = Complex(2, 0);    a[3] = Complex(3, 0);    a = IterativeFFT(a, 4, 1);    cout<<"----------After DFT----------"<<endl;    for(int i = 0; i < 4; i++)        printf("%.9f + (%.9f) i\n", a[i].real, a[i].image);    cout<<"----------After DFT-1----------"<<endl;    a = IterativeFFT(a, 4, -1);    for(int i = 0; i < 4; i++)        printf("%.9f + (%.9f) i\n",a[i].real, a[i].image);    return 0;}

通過上面這個代碼的樣本, FFT演算法的實現基本沒有什麼問題了, 另外演算法導論中的習題有一些很不錯, 便於熟悉這一演算法的很多細節, 這裡就不一一提及了。
經過這些學習之後, 進行一些實戰演練是很有必要的, 接下來是相關習題的練習部分。

HDU 1402 A*B Problem Plus 大整數乘法
HDU 4609 3-idiots FFT計數
UVALive 4671 K-neighbor Substrings FFT算字串Hamming距離 
UVA 12298 Super Poker II FFT計數, long double
URAL 1996 Cipher Massage 3 FFT + KMP 
CodeChef COUNTARI Arithmetic Progressions FFT + 分塊
ZOJ 3856 Goldbach FFT計數
UVALive 6886 Golf Bot FFT模板題
HDU 4093 Xavier is Learning to Count 容斥原理 + FFT,(2011年上海現場賽C題)

HDU 5751 BestCoder Round #84 Eades（線段樹+FFT）

HDU 5730 2016多校1 Shell Necklace （CDQ分治+FFT）

2016 acm香港網路賽 A題 A+B Problem （FFT）

以後還有題目可以看我部落格，可能有。最後感謝Ichimei對FFT的講解。Orz....對理解FFT真的有很大協助。