天涯明月刀ol|天涯明月刀迅雷下载
  1. 當前所在位置:
  2. 首頁
  3. 深海捕魚

FFT詳解 什么是FFT

2018-12-24 admin
最新更新
2018.6.29 新增了“關于DFT的再次思考” 和 “有關FFT算法實現機理的再討論”兩個小節,希望能對沒學明白的同學有所幫助。
 
2018.6.23 回復了評論區的兩個評論,如果同學們覺得我的講解有錯誤,或者是有講得不夠清晰的地方,請在評論區評論評論,我會虛心采納的,謝謝。
 
2017.12.8 利用LATEXLATEX對公式進行重新整理,用以輔助原來的圖片公式。但是為了體現出這篇博客的”歷史滄桑感”,我保留了原有的所有圖片(暫時只完成了一部分的公式轉換,公式實在是太多了!)。—— GGN
 
序言
(來學習的同學們請直接跳過這里就好了…)
 
非常感謝Leo學長為我耐心的講解,在此對學長表示誠摯的敬意。停課學習競賽的那一個月里,學長們不惜耽誤自己的時間為我耐心講解算法,對我的幫助是我永難忘記的。最后衷心祝愿高三學長們在高考中取得優異的成績!
 
同時我也要向學長們致歉,感謝學長們曾經對我的包容。
 
2017.11.19 GGN 補加序言
 
前言
一直想學FFT,當時由于數學基礎太差,導致啥都學不懂。請教了機房里的幾位學長大神,結果還是沒太明白。因此下定決心寫一篇關于“FFT”的文章,一篇起碼我能看得懂的“FFT”。不過這好像并不是一件簡單的事,因為首先我要學會“FFT”的原理。
 
建議同學們先自學一下“復數(虛數)”的性質、運算等知識,不然看這篇文章有很大概率看不懂。最后你就會發現,這不是一個“算法”的博客,而是一個數學的博客。
 
這是當時教我FFT的機房學長的博客:
 
Leo 的 《FFT與多項式乘法》
 
(溫馨提示:整篇文章都不是很好理解,提醒同學們保持清醒的頭腦。)
 
1.什么是FFT?
FFT,即為快速傅氏變換,是離散傅氏變換的快速算法,它是根據離散傅氏變換的奇、偶、虛、實等特性,對離散傅立葉變換的算法進行改進獲得的。它對傅氏變換的理論并沒有新的發現,但是對于在計算機系統或者說數字系統中應用離散傅立葉變換,可以說是進了一大步。——360百科
 
FFT(Fast Fourier Transformation)就是“快速傅里葉變換”的意思,它是一種用來計算DFT(離散傅里葉變換)和IDFT(離散傅里葉反變換)的一種快速算法。這種算法運用了一種高深的數學方式、把原來復雜度為O(n2)O(n2)的樸素多項式乘法轉化為了O(nlogn)O(nlogn)的算法。
 
2.多項式乘法的樸素算法
假設我們現在不會FFT,我們是又是如何計算多項式乘法的呢?現在我們有兩個關于x的二次多項式:
 
f(x)=a1x2+b1x+c1
f(x)=a1x2+b1x+c1
 
g(x)=a2x2+b2x+c2
g(x)=a2x2+b2x+c2
 
 
我命令K(x)為他們的乘積,則有:
 
K(x)=f(x)×g(x)=a1a2x4+(a1b2+a2b1)x3+(a1c2+a2c1+b1b2)x2+(b1c2+b2c1)x+c1c2
K(x)=f(x)×g(x)=a1a2x4+(a1b2+a2b1)x3+(a1c2+a2c1+b1b2)x2+(b1c2+b2c1)x+c1c2
 
 
如果在程序中,我們用一個數組來儲存一個多項式的各個項的系數,我們如何去做這樣一個復雜的乘法呢?
 
#include<iostream>
#include<vector>
#include<cstdlib>
using namespace std;
 
vector<double>ForceMul(vector<double>A,vector<double>B)//表示A,B兩個多項式相乘的結果
{
    vector<double>ans;
    int aLen=A.size();//A的元素個數
    int bLen=B.size();//B的元素個數
    int ansLen=aLen+bLen-1;//ans的元素個數=A的元素個數+B的元素個數-1
    for(int i=1;i<=ansLen;i++)//初始化ans
        ans.push_back(0);
    for(int i=0;i<aLen;i++)
        for(int j=0;j<bLen;j++)
            ans[i+j]+=A[i]*B[j];//A的i次項 與 B的j次項 相乘的結果 累加到ans的[i+j]次位
    return ans;//返回ans
}
 
int main()
{
    vector<double>A,B;
    cout<<"input A:";
    for(int i=0;i<3;i++)//從0次項開始輸入A的各項系數
    {
        int x;
        cin>>x;
        A.push_back(x);
    }
    cout<<"input B:";
    for(int i=0;i<3;i++)//從0次項開始輸入B的各項系數
    {
        int x;
        cin>>x;
        B.push_back(x);
    }
    vector<double>C=ForceMul(A,B);//C=A與B暴力相乘
    cout<<"output C:";
    for(int i=0;i<5;i++)//從0次項開始輸出C的各項系數
        cout<<C[i]<<" ";
    cout<<endl;
    system("pause");
    return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
這就是樸素算法,它的復雜度為O(lenA*lenB)。如果lenA=lenB=10^5,程序時間就會爆掉,那么我們如何進行優化呢?
 
3.系數表示法與點值表示法
系數表表示法,就是用一個多項式的各個項的系數表示這個多項式,也就是我們平時所用的表示法。例如,我們可以這樣表示:
 
f(x)=a0+a1x+a2x2+..+anxn⇔f(x)={a0,a1,a2,..,an}
f(x)=a0+a1x+a2x2+..+anxn⇔f(x)={a0,a1,a2,..,an}
 
 
點值表示法,就是把這個多項式理解成一個函數,用這個函數上的若干個點的坐標來描述這個多項式。
 
(兩點確定一條直線,三點確定一條拋物線…同理n+1個點確定一個n次函數,其原理來自于“高斯消元”,下文會有介紹。)
 
因此表示成這樣:(注意:x[0]->x[n]是n+1個點)
 
f(x)=a0+a1x+a2x2+..+anxn⇔f(x)={(x0,y0),(x1,y1),(x2,y2),..,(xn,yn)}
f(x)=a0+a1x+a2x2+..+anxn⇔f(x)={(x0,y0),(x1,y1),(x2,y2),..,(xn,yn)}
 
 
為什么n+1個確定的點能確定一個唯一的多項式呢?你可以嘗試著把這n+1個點的值分別代入多項式中:
 
f(x0)=y0=a0+a1x0+a2x20+..+anxn0
f(x0)=y0=a0+a1x0+a2x02+..+anx0n
f(x1)=y1=a0+a1x1+a2x21+..+anxn1
f(x1)=y1=a0+a1x1+a2x12+..+anx1n
f(x2)=y2=a0+a1x2+a2x22+..+anxn2
f(x2)=y2=a0+a1x2+a2x22+..+anx2n
...
...
f(xn)=yn=a0+a1xn+a2x2n+..+anxnn
f(xn)=yn=a0+a1xn+a2xn2+..+anxnn
 
 
你會得到n+1個方程,其中x[0~n]和y[0~n]是已知的,a[0~n]是未知的。n+1的未知數,n+1個方程所組成的方程組為“n+1元一次方程”,因為它是“一次方程”,所以(一般情況下,不考慮無解和無數解)可以通過“高斯消元”解得所有未知數唯一確定的值。也就是說,用點知表示法可以確定出 唯一確定的 系數表示法中的 每一位的 系數。
 
這種把一個多項式轉化成“離散的點”表示的方法就叫做“DFT”(離散傅里葉變換)。
 
把“離散的點”還原回多項式的方法就叫做“IDFT”(離散傅里葉反變換)。
 
(請同學們重視上面的這兩句話,因為這是我能想到的最好理解的解釋方法了。)
 
有興趣的可以看一下360百科給出的DFT定義:
 
離散傅里葉變換(Discrete Fourier Transform,縮寫為DFT),是傅里葉變換在時域和頻域上都呈離散的形式,將信號的時域采樣變換為其DTFT的頻域采樣。在形式上,變換兩端(時域和頻域上)的序列是有限長的,而實際上這兩組序列都應當被認為是離散周期信號的主值序列。即使對有限長的離散信號作DFT,也應當將其看作其周期延拓的變換。在實際應用中通常采用快速傅里葉變換計算DFT。——360百科
 
4.“復數”的引入
初中我們可能學過一些定理:比如−1−−−√−1不存在。但是,當時我們的“數域”僅僅止步于“實數”,而現在我們要介紹一個新的數域——“虛數”。我們令:i=−1−−−√i=−1表示這個“虛數單位”,“ii”對于虛數的意義就相當于是“數字1”對于實數的意義。
 
這是百科對“虛數”的定義:
 
虛數可以指不實的數字或并非表明具體數量的數字。——360百科
 
在數學中,虛數就是形如a+b*i的數,其中a,b是實數,且b≠0,i² = - 1。虛數這個名詞是17世紀著名數學家笛卡爾創立,因為當時的觀念認為這是真實不存在的數字。后來發現虛數a+b*i的實部a可對應平面上的橫軸虛部b與對應平面上的縱軸,這樣虛數a+b*i可與平面內的點(a,b)對應。——360百科
 
然后是百科對“復數”的定義:
 
復數x被定義為二元有序實數對(a,b) ,記為z=a+bi,這里a和b是實數,i是虛數單位。在復數a+bi中,a=Re(z)稱為實部,b=Im(z)稱為虛部。當虛部等于零時,這個復數可以視為實數;當z的虛部不等于零時,實部等于零時,常稱z為純虛數。復數域是實數域的代數閉包,也即任何復系數多項式在復數域中總有根。 復數是由意大利米蘭學者卡當在十六世紀首次引入,經過達朗貝爾、棣莫弗、歐拉、高斯等人的工作,此概念逐漸為數學家所接受。
 
正如上文所說,一個復數可以看成是“復平面”上的一個點。復平面就是以實數部為x軸,以虛部為y軸所組成的類似“直角坐標系”的一個平面坐標體系。同樣,我們也可以用“極坐標”來表示一個平面中的點。
 
 
 
上圖就是在這個復平面上的一個點的三種表示方法。
 
然后,思考一個簡單的問題:兩個復數的乘法有沒有某種特定的幾何意義?(只是一個數學性質,在此不進一步深究,可用三角函數證明。)
 
 
 
兩復數相乘,“長度相乘,極角相加”。不難想象如果兩個復數它們到“坐標原點”的距離都是1,那么它們的乘積到坐標原點的距離還是一,只不過是繞著原點進行了旋轉。
 
5.單位復根
現在,回到我們剛才講到的“點值表示法”的問題,考慮這樣一個問題,如果我有兩個用點值表示的多項式,如何表示它們兩個多項式的乘積呢?(假設這兩個多項式選取的所有點的x值恰好相同。)
 
f(x)={(x0,f(x0)),(x1,f(x1)),(x2,f(x2)),..,(xn,f(xn))}
f(x)={(x0,f(x0)),(x1,f(x1)),(x2,f(x2)),..,(xn,f(xn))}
g(x)={(x0,g(x0)),(x1,g(x1)),(x2,g(x2)),..,(xn,g(xn))}
g(x)={(x0,g(x0)),(x1,g(x1)),(x2,g(x2)),..,(xn,g(xn))}
 
 
如果F(x)=f(x)×g(x)F(x)=f(x)×g(x),那么就有F(x0)=f(x0)∗g(x0)F(x0)=f(x0)∗g(x0)(x0x0為任意數)。也就是說,如果我把兩個函數的點值表示法中的xx值相同的點的y值乘在一起就是它們的乘積(新函數)的點值表示。(這是一個O(n)O(n)的操作。)
 
f(x)×g(x)={(x0,f(x0)g(x0)),(x1,f(x1)g(x1)),(x0,f(x2)g(x2)),..,(xn,f(xn)g(xn))}
f(x)×g(x)={(x0,f(x0)g(x0)),(x1,f(x1)g(x1)),(x0,f(x2)g(x2)),..,(xn,f(xn)g(xn))}
 
 
但是我們要的是系數表達式,而不是點制表達式,如果用高斯消元暴力地去解一個“n+1元方程組”恐怕就要把時間復雜度拉回O(n2)O(n2)(甚至更高)。為什么呢?因為當我們計算x0,x20,...,xn0x0,x02,...,x0n時會浪費大量的時間。這個數學運算看似是沒有辦法加速的,而實際上我們可以找到一種神奇的“x值”,帶進去之后不用反復地去做無用n次方操作。我能想到的第一個數就是“1”,因為1的幾次方都是1。然后我能想到的數就是“-1”,因為“-1”的平方是1。把這樣的數帶進去就可以減少我們的運算次數。
 
但是問題又來了,我們只有“-1”和“1”兩個數,但是我們要至少帶進去n+1個不同的數才能進行系數表示。考慮一個問題:也許“虛數”可能會幫上我們大忙。我們需要的是滿足“ωk=1ωk=1”的數(k為整數),然后我們就會發現“i”好像也滿足這個條件:i*i=-1,(i*i)*(i*i)=1=i^4,當然“-i”也有這個性質。然而僅僅4個數還是不能滿足我們的需求。
 
 
 
看上圖中的紅圈,紅圈上的每一個點距原點的距離都是1個單位長度,所以說如果說對這些點做k次方運算,它們始終不會脫離這個紅圈。因為它們在相乘的時候r始終=1,只是θ的大小在發生改變。而這些點中有無數個點經過若干次方之后可以回到“1”。因此,我們可以把這樣的一組神奇的x帶入函數求值。
 
像這種能夠通過k次方運算回到“1”的數,我們叫它“復根”用“ω”表示。如果ωk=1ωk=1,那么我們稱“ωω為1的k次復根”計做“ωnkωkn”:
 
 
 
其中“n”就是一個序號數,我們把所有的“負根”按照極角大小逆時針排序從零開始編號。以“四次負根”“ωn4ω4n”為例:
 
 
 
你會發現:其實k次負根就相當于是給圖中的圓周平均分成k個弧,弧與弧之間的端點就是“復根”,另外ω24=−1=i2=(ω14)2,ω34=−i=i3=(ω14)3ω42=−1=i2=(ω41)2,ω43=−i=i3=(ω41)3,ω04ω40是這個圓與“Real”軸正半軸的交點,所以無論k取多少,ω04ω40始終=1。我們只需要知道ω14ω41,就能求出ωnkωkn,所以我們稱“ω1kωk1”為“單位復根”。
 
其實,我們用“ωkωk”表示單位復根,ω1kωk1表示的是“單位復根”的“1次方”也就是它本身,其他的就叫做k次單位復根的n次方。
 
 
 
6.FFT的主要流程之DFT
說了這么多了,終于說到FFT了。FFT運用到了一種分治的思想,分治地去求當x=ω(k)[n]時整個多項式的值(永遠都不要忘了每一個步驟的目的是什么)。你可以把一個多項式分成奇數次數項,和偶數次數項兩部分,然后再用分治的思想去處理它的“奇數次項”和“偶數次項”。
 
 
 
我們用DFT(F(x))[k]表示當x=ω^k時F(x)的值,所以有:DFT(F(x))[k]=DFT(G(x^2))[k]+ω^k*DFT(H(x^2)),也就是:
 
 
 
(把當前單位復根的平方分別以DFT的方式帶入G函數和H函數求值。)
 
但是這個二分最大的局限就是只能處理長度為2的整數次冪的多項式,因為如果長度不為的整數次冪,二分到最后就會出現左半部分右半部分長度不一致的情況(導致程序取不到系數而爆炸),所以我們在做第一次DFT之前一定要把這個多項式補成長度為2^n(n為整數)的多項式(補上去的高次項系數為“0”),長度為2^n的多項式的最高次項為2^n-1次項。當我們向這個式子中“帶入數值”的時候,一定要保證我帶入的每個數都是不一樣的,所以要帶入1的2^n的單位復根的各個冪(因為1的2^n復根恰好有2^n個)。
 
【2018.6.30】 我覺得上面這一段內容講得還不夠清晰,如果同學們沒明白這段話的意思,可以先跳到后面“關于DFT的再次思考”,看看那里的講解。看完以后再回到這里,繼續閱讀。
 
但是這個算法還需要從“分治”的角度繼續優化。我們每一次都會把整個多項式的奇數次項和偶數次項系數分開,一只分到只剩下一個系數。但是,這個遞歸的過程需要更多的內存。因此,我們可以先“模仿遞歸”把這些系數在原數組中“拆分”,然后再“倍增”地去合并這些算出來的值。然而我們又要如何去拆分這些數呢?
 
 
 
貌似并沒有什么規律可循,但實際上你可要看仔細了,把這些下標都轉化成二進制看看吧:
 
 
 
是不是發現了一些神奇的特點:“拆分”之后的序列的下標恰好為長度為3位的二進制數的翻轉。也就是說我們對原來的每個數的下標進行長度為三的二進制翻轉就是新的下標。而為什么是長度為3呢?因為8=2^3。為了證明這一點,我們可以再舉一個簡單例子(它還是成立的):
 
 
 
(一個數學性質,在此不再證明。我感覺這個原理有點像是“基數排序”,感興趣的同學可以去看看。)
 
關于二進制翻轉是如何實現的,本文中并沒有介紹,強烈建議看一下這篇文章:
 
學習鏈接:補充——FFT中的二進制翻轉問題
 
7.FFT的主要流程之IDFT
我們先整理一下思路,IDFT是做什么的?IDFT(傅里葉反變換)就是把一個用“點值表示法”表示的多項式,轉化成一個用“系數表示法”表示的多項式,但是這似乎并不是很容易。然而其實我們剛剛恰好做了一些非常機智的事情——把“單位復根”的若干次方帶入了原多項式。我們可以表示一下這些多項式(這里使用一個矩陣表示,不會的建議自學)。
 
 
 
如果我們想把這個表達式還原成只含有“a系數”的矩陣,那么就要在中間那個“巨大的矩陣”身上乘上一個它的“反矩陣”(反對稱矩陣)就可以了。這個矩陣的中有一種非常特殊的性質,對該矩陣的每一項取倒數,再除以n就可以得到該矩陣的反矩陣。而如何改變我們的操作才能使計算的結果文原來的倒數呢?這就要看我們求“單位復根的過程了”:根據“歐拉函數”e^iπ=-1,我么可以得到e^2πi=1。如果我要找到一個數,它的k次方=1,那么這個數ω[k]=e^(2πi/k)(因為(e^(2πi/k))^k=e^(2πi)=1)。而如果我要使這個數值變成“1/ω[k]”也就是“(ω[k])^-1”,我們可以嘗試著把π取成-3.14159…,這樣我們的計算結果就會變成原來的倒數,而其它的操作過程與DFT是完全相同的(這真是極好的)。我們可以定義一個函數,向里面摻一個參數“1”或者是“-1”,然后把它乘到“π”的身上。傳入“1”就是DFT,傳入“-1”就是IDFT,十分的智能。
 
機房學長的代碼寫得實在是太好了Orz,忍不住引用了一下:
 
下面的代碼來自 Leo的 《FFT與多項式乘法》
 
歡迎同學們一同Orz ↑
 
這是整個代碼的一個局部(出于禮貌,保留了原有的代碼格式,但是加了一些注釋):
 
int rev[maxl];
void get_rev(int bit)//bit表示二進制位數,計算一個數在二進制翻轉之后形成的新數
{
    for(int i=0;i<(1<<bit);i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}
void fft(cd *a,int n,int dft)//n表示我的多項式位數
{
    for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
    //中間的那個if保證了每個數做多只被交換了1次
    //如果不寫那么會有一些數被交換兩次,導致最終的位置沒有變
    for(int step=1;step<n;step<<=1)//模擬一個合并的過程
    {
        cd wn=exp(cd(0,dft*PI/step));//計算當前單位復根
        for(int j=0;j<n;j+=step<<1)
        {
            cd wnk(1,0);//計算當前單位復根
            for(int k=j;k<j+step;k++)
            {//蝴蝶操作
                cd x=a[k];
                cd y=wnk*a[k+step];
                a[k]=x+y;//這就是上文中F(x)=G(x)+ωH(x)的體現
                a[k+step]=x-y;
                    //后半個“step”中的ω一定和“前半個”中的成相反數
                    //“紅圈”上的點轉一整圈“轉回來”,轉半圈正好轉成相反數
                wnk*=wn;
            }
        }
    }
    if(dft==-1) for(int i=0;i<n;i++) a[i]/=n;
    //考慮到如果是IDFT操作,整個矩陣中的內容還要乘上1/n  
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
關于DFT的再次思考
2018.6.29 我最近有思考了一些關于DFT原理的問題,覺得之前講得有些草率,決定再寫一個表述更清晰一些的版本。每個人都有適合的自己講課風格,如果沒看懂上面的講解,不妨試試看看這里。
 
說到DFT(離散傅里葉變換)這個環節,就必須先說這個環節的目的。做音頻處理的人用這個過程來實現值域與頻域的轉換,而我們OIer主要是為了把一個用系數表示法表示的函數轉化為一個用點值表示法表示的函數。
 
我在前文中提到了一類神奇的數——單位復根,并且了解了關于它的一些性質: 
1.我們通常用 ωn=e2πinωn=e2πin表示n次單位復根(的一次方)
 
2.n次單位復根的值隨指數變化而循環,有 ωn+kn=ωknωnn+k=ωnk,(也就是說k為整數時,ωknωnk恰有n種不同的取值,分別為ω0n,ω1n,ω2n,...,ωn−1nωn0,ωn1,ωn2,...,ωnn−1,k等于其余整數值時,都可以通過幾次+n或者-n 平移 到這個區間里)
 
3.折半引理: ωkn=e2πin⋅k=e2πin2⋅k2=ωk2n2ωnk=e2πin⋅k=e2πin2⋅k2=ωn2k2,這個性質一會兒會用到,一定要理解透徹。
 
一個n-1次多項式可以用n個系數去描述,也可以用函數上的n個點去描述,n+1個點當然也可以,但是多余。根據上述2號結論,我們知道n次復根恰好有n種不同的取值,那我們不妨就把這n個不同的值帶入多項式。
 
這樣實際上,問題就轉化成——把一個長度為n的向量A,變換成另一個長度為n的向量B。(我們記A與B的下標從0開始,到n-1結束,共n個位置。)
 
其中: 
A 為系數向量,AkAk 中儲存原多項式中xkxk項前的系數。 
B 為點值向量,BkBk 中儲存 x=ωknx=ωnk 時原多項式的值。
 
為了描述方便,我們稱 B向量 為多項式 A0+A1x+A2x2+...+An−1xn−1A0+A1x+A2x2+...+An−1xn−1的變換結果。
 
根據這個定義,我們可以表示出B向量某一個位置k上的值BkBk:
 
Bk=A0+A1⋅(ωkn)+A2⋅(ωkn)2+A3⋅(ωkn)3+...+An−1⋅(ωkn)n−1Bk=A0+A1⋅(ωnk)+A2⋅(ωnk)2+A3⋅(ωnk)3+...+An−1⋅(ωnk)n−1
我們不妨假設n為偶數(因為多向式可以在最后補上一些系數為0的高次項,所以這一假設對運算結果沒有影響),我們把這個多項式的奇數次項和偶數次項分開,看看會不會得到一些驚人的結論。
 
推導(因為n為偶數,所以n-2為偶數、n-1為奇數):
 
Bk=(A0+A2⋅(ωkn)2+A4⋅(ωkn)4...+An−2⋅(ωkn)n−2)+(A1⋅(ωkn)+A3⋅(ωkn)3+A5⋅(ωkn)5+...+An−1⋅(ωkn)n−1)Bk=(A0+A2⋅(ωnk)2+A4⋅(ωnk)4...+An−2⋅(ωnk)n−2)+(A1⋅(ωnk)+A3⋅(ωnk)3+A5⋅(ωnk)5+...+An−1⋅(ωnk)n−1)
看著不“對稱”,把右半部分(也就是所有奇數項)都提出一個ωknωnk來:
 
Bk=(A0+A2⋅(ωkn)2+A4⋅(ωkn)4...+An−2⋅(ωkn)n−2)+ωkn⋅(A1+A3⋅(ωkn)2+A5⋅(ωkn)4+...+An−1⋅(ωkn)n−2)Bk=(A0+A2⋅(ωnk)2+A4⋅(ωnk)4...+An−2⋅(ωnk)n−2)+ωnk⋅(A1+A3⋅(ωnk)2+A5⋅(ωnk)4+...+An−1⋅(ωnk)n−2)
這次,左半部分和右半部分的指數得到了統一,都是0次到n-2次,也就是都是偶數了(后面還會用到這個偶數指數的性質)。
 
好像可以用一次折半引理,不用白不用,試試看:
 
(ωkn)2=ω2kn=ωkn2(ωnk)2=ωn2k=ωn2k
所以有:
 
Bk=(A0+A2⋅(ωkn2)+A4⋅(ωkn2)2...+An−2⋅(ωkn2)n2−1)+ωkn⋅(A1+A3⋅(ωkn2)1+A5⋅(ωkn2)2+...+An−1⋅(ωkn2)n2−1)Bk=(A0+A2⋅(ωn2k)+A4⋅(ωn2k)2...+An−2⋅(ωn2k)n2−1)+ωnk⋅(A1+A3⋅(ωn2k)1+A5⋅(ωn2k)2+...+An−1⋅(ωn2k)n2−1)
誒?好像…A0+A2⋅(ωkn2)+A4⋅(ωkn2)2...+An−2⋅(ωkn2)n2−1A0+A2⋅(ωn2k)+A4⋅(ωn2k)2...+An−2⋅(ωn2k)n2−1這個式子不就是多項式A0+A2⋅x+A4⋅x2...+An−2⋅xn2−1A0+A2⋅x+A4⋅x2...+An−2⋅xn2−1 在 x=ωkn2x=ωn2k處的點值嗎?右半部分?右半部分好像也很像,像是多項式A1+A3⋅x+A5⋅x2...+An−1⋅xn2−1A1+A3⋅x+A5⋅x2...+An−1⋅xn2−1 在 x=ωkn2x=ωn2k處的點值。而且這兩個多項式都恰好有 n2n2 項。
 
把n個n次復根帶入一個長度為n的多項式,是我們要做的DFT過程,那把n/2個 n/2次復根帶入一個長度為n/2的多項式,是不是也可以看成是DFT?
 
假如我們已經知道了多項式A0+A2⋅x+A4⋅x2...+An−2⋅xn2−1A0+A2⋅x+A4⋅x2...+An−2⋅xn2−1的變換結果為長度為n/2的向量L,多項式A1+A3⋅x+A5⋅x2...+An−1⋅xn2−1A1+A3⋅x+A5⋅x2...+An−1⋅xn2−1的9 變換結果 (變換結果一詞在前文中有定義)為長度為n/2的向量R, 實際上B向量中的任意一個位置就可以O(1)算出了(接下來我們說具體怎么 O(1) 求)。
 
因為根據上文的定義有:Bk=Lk+ωkn⋅RkBk=Lk+ωnk⋅Rk,其中0≤k<n20≤k<n2,因為L和R的長度只有n2n2。這種方法可以求出B向量的左半部分。那右半部分該怎么求呢?
 
其實還有這樣一個式子,Bn2+k=Lk−ωkn⋅RkBn2+k=Lk−ωnk⋅Rk,其中0≤k<n20≤k<n2。
 
我們知道ωn2+kn=ωkn⋅ωn2n=ωkn⋅e2πin⋅n2=ωkn⋅eπi=−ωknωnn2+k=ωnk⋅ωnn2=ωnk⋅e2πin⋅n2=ωnk⋅eπi=−ωnk。
 
前文中我們提過一嘴,Bk=(A0+A2⋅(ωkn)2+A4⋅(ωkn)4...+An−2⋅(ωkn)n−2)+ωkn⋅(A1+A3⋅(ωkn)2+A5⋅(ωkn)4+...+An−1⋅(ωkn)n−2)Bk=(A0+A2⋅(ωnk)2+A4⋅(ωnk)4...+An−2⋅(ωnk)n−2)+ωnk⋅(A1+A3⋅(ωnk)2+A5⋅(ωnk)4+...+An−1⋅(ωnk)n−2)
這個式子中,除了奇數次項提出的一個ωknωnk外,其余所有關于ωknωnk的項都是偶數次項。一個數的平方與其相反數的平方一定相等,一個數的偶次方與其相反數的偶次方也一定相等。ωn2+kn=−ωknωnn2+k=−ωnk而,這就說明,在對B向量右半部分的求解中,我們仍然可以在L向量和R向量中得到我們想要的結果。
 
綜上,有:
 
Bk=Lk+ωkn⋅RkBk=Lk+ωnk⋅Rk,其中0≤k<n20≤k<n2 
Bn2+k=Lk−ωkn⋅RkBn2+k=Lk−ωnk⋅Rk,其中0≤k<n20≤k<n2
我們把這個 通關過 L和R 求解 B 的方法 稱為 信息合并 。
 
同理,我們想要求解L向量和R向量,繼續遞歸即可,把一個長度為n2n2的多項式分成兩個長度為n4n4的多項式,然后進行DFT,再進行信息合并。遞歸直到到多項式長度為1為止(因為長度為1的多項式只有常數項,變換結果就是其本身。)
 
這就是這個算法的主要思想,但是思想與代碼之間還是有很大的差距的,我知道很多同學看不懂上面的代碼,我在這里解釋一下這個代碼的實現機理。
 
有關FFT算法實現機理的再討論
每一次遞歸,我們都要求被處理的多項式的項數是偶數(或者是1),那么我們不妨先在多項式系數表示法尾部補0,把它的長度補成一個2的整數次冪,這樣每次都可以繼續分治(因為整個多項式的長度不會超過原多項式的二倍,所以時間復雜度不變)。
 
為了節省空間,我們從頭到尾只使用一個長度為n的數組arr(下標也是從0到n-1)來實現整個DFT的過程,最開始的時候我們用arr儲存A。
 
接下來的內容表示arr數組中內容的“變遷”。
 
arr={A0,A1,...,An2−1,An2,An2+1,...,An−1}
arr={A0,A1,...,An2−1,An2,An2+1,...,An−1}
奇數偶數分項:
 
arr={A0,A2,...,An−2,A1,A3,...,An−1}
arr={A0,A2,...,An−2,A1,A3,...,An−1}
遞歸計算,直接把計算結果L和R覆蓋到原數組對應的位置上
 
arr={L0,L1,...,Ln2−1,R0,R1,...,Rn2−1}
arr={L0,L1,...,Ln2−1,R0,R1,...,Rn2−1}
最后信息合并,得到:
 
arr={B0,B1,...,Bn2−1,Bn2,Bn2+1,...,Bn−1}
arr={B0,B1,...,Bn2−1,Bn2,Bn2+1,...,Bn−1}
信息合并,又叫做蝴蝶操作,它為什么會有這么美麗的一個名字呢,我想用一個圖來談談我的見解。
 
 
 
我們根據上文關于信息合并的推導可知,L0L0和R0R0只決定B0B0和Bn2Bn2, L1L1和R1R1只決定B1B1和Bn2+1Bn2+1,…,LkLk和RkRk只決定BkBk和Bn2+kBn2+k。所以,我們可以直接把BkBk和Bn2+kBn2+k的值覆蓋到LkLk和RkRk的位置上(而且發現它們恰好在數組上的位置也是對應的)。
 
我們用一個變量k從0循環到n2−1n2−1,每次計算如下兩個值:
 
x=Lk,y=ωkn⋅Rk
x=Lk,y=ωnk⋅Rk
這樣就有:
 
Bk=x+y,Bn2+k=x−y
Bk=x+y,Bn2+k=x−y
而ωkn=ωk−1n⋅ωnωnk=ωnk−1⋅ωn,是可以在循環的時候遞推出來的。
 
而事實上,我們前面也討論過,分治的過程是可以利用循環來實現的。我們先只奇數偶數分項,然后得到了下標二進制反轉后的序列。(這兩個詞我沒給定義,關于下表二進制反轉的問題,我覺得第6節“FFT的主要流程之DFT”講得挺明白的。)
 
然后,逐層實現信息的合并,第一層,把相鄰的兩個長度為1的變換結果,合并成一個長度為2的變換結果,信息塊個數除以2。第二層,把相鄰的兩個長度為2的的變換結果,合并成一個長度為4的變換結果…這樣我們只需要O(nlogn)的時間復雜度就可以求出整個多項式的變換結果(每層合并都需要花費O(n)的總時間代價,一共有O(logn)層合并)。
 
IDFT的實現原理不再贅述,見上文第7節“FFT的主要流程之IDFT”。
 
現在我把代碼再貼一次:
 
int rev[maxl];
void get_rev(int bit)//bit表示二進制位數,計算一個數在二進制翻轉之后形成的新數
{
    for(int i=0;i<(1<<bit);i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}
void fft(cd *a,int n,int dft)//n表示我的多項式位數
{
    for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
    //中間的那個if保證了每個數做多只被交換了1次
    //如果不寫那么會有一些數被交換兩次,導致最終的位置沒有變
    for(int step=1;step<n;step<<=1)//模擬一個合并的過程
    {
        cd wn=exp(cd(0,dft*PI/step));//計算當前單位復根
        for(int j=0;j<n;j+=step<<1)
        {
            cd wnk(1,0);//計算當前單位復根
            for(int k=j;k<j+step;k++)
            {//蝴蝶操作
                cd x=a[k];
                cd y=wnk*a[k+step];
                a[k]=x+y;//這就是上文中F(x)=G(x)+ωH(x)的體現
                a[k+step]=x-y;
                    //后半個“step”中的ω一定和“前半個”中的成相反數
                    //“紅圈”上的點轉一整圈“轉回來”,轉半圈正好轉成相反數
                    //一個數相反數的平方與這個數自身的平方相等..
                wnk*=wn;
            }
        }
    }
    if(dft==-1) for(int i=0;i<n;i++) a[i]/=n;
    //考慮到如果是IDFT操作,整個矩陣中的內容還要乘上1/n  
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
感覺我一年前寫的注釋雖然不是很清晰,但是還是蠻正確的,就沒有改動。希望同學們能有一個更深入的理解~
 
8.后記
最后總結一下FFT的優化思想:
 
 
 
然后是FFT的優化理念:
 
 
 
(完全自創,如有雷同,純屬巧合,盡管并沒什么用。)
 
FFT其實還可以用來計算BIGNUM乘法,因為我們可以把一個長整數理解成a[0]+a[1]*10+a[2]*10^2+…+a[n]*10^n。把“10”當成未知數,這個多項式每一個次方項的系數就是BIGNUM每一數位上的數。而這時,數組長度“n”就不能單純的取這個十進制數的長度,而要取大于等于兩個十進制數長度加和的最小的2的正整數次冪。因為我們要保證DFT得到的離散點的個數足夠表示我最終生成的新多項式(也就是取的點的個數要大于等于這個結果多項式的長度)。
 
總之,這真的是一個很好的算法,但是要注意,除題中的數據范圍(多項式長度)超過了10^4,否則暴力是可以解決的。而且FFT常數巨大,請同學們一定要慎用!慎用!慎用!
 
寫了一整天,可算是把這篇長篇大論的博客寫完了。FFT就是一個典型的用數學方法對問題實現優化的方法。DFT對數學基礎要求很高,但是信息學與數學的聯系是相當緊密的,所以這種情況時常會出現。數學不好的同志們一定要加油哦~
 
趕稿匆忙,如有謬誤,望同學們諒解。
 
2017.11.19【補】第一次自己手寫FFT,感覺很開心
 
#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<complex>
using namespace std;
 
typedef complex<double> cd;//復數類的定義
const int maxl=2094153;//nlogn的最大長度(來自leo學長的博客)
const double PI=3.14159265358979;//圓周率,不解釋
 
cd a[maxl],b[maxl];//用于儲存變換的中間結果
int rev[maxl];//用于儲存二進制反轉的結果
void getrev(int bit){
    for(int i=0;i<(1<<bit);i++){//高位決定二進制數的大小
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }//能保證(x>>1)<x,滿足遞推性質
}
 
void fft(cd* a,int n,int dft){//變換主要過程
    for(int i=0;i<n;i++){//按照二進制反轉
        if(i<rev[i])//保證只把前面的數和后面的數交換,(否則數組會被翻回來)
            swap(a[i],a[rev[i]]);
    }
    for(int step=1;step<n;step<<=1){//枚舉步長的一半
        cd wn=exp(cd(0,dft*PI/step));//計算單位復根
        for(int j=0;j<n;j+=step<<1){//對于每一塊
            cd wnk(1,0);//!!每一塊都是一個獨立序列,都是以零次方位為起始的
            for(int k=j;k<j+step;k++){//蝴蝶操作處理這一塊
                cd x=a[k];
                cd y=wnk*a[k+step];
                a[k]=x+y;
                a[k+step]=x-y;
                wnk*=wn;//計算下一次的復根
            }
        }
    }
    if(dft==-1){//如果是反變換,則要將序列除以n
        for(int i=0;i<n;i++)
            a[i]/=n;
    }
}
 
int output[maxl];
char s1[maxl],s2[maxl];
int main(){
    scanf("%s%s",s1,s2);//讀入兩個數
    int l1=strlen(s1),l2=strlen(s2);//就算"次數界"
    int bit=1,s=2;//s表示分割之前整塊的長度
    for(bit=1;(1<<bit)<l1+l2-1;bit++){
        s<<=1;//找到第一個二的整數次冪使得其可以容納這兩個數的乘積
    }
    for(int i=0;i<l1;i++){//第一個數裝入a
        a[i]=(double)(s1[l1-i-1]-'0');
    }
    for(int i=0;i<l2;i++){//第二個數裝入b
        b[i]=(double)(s2[l2-i-1]-'0');
    }
    getrev(bit);fft(a,s,1);fft(b,s,1);//dft
    for(int i=0;i<s;i++)a[i]*=b[i];//對應相乘
    fft(a,s,-1);//idft
    for(int i=0;i<s;i++){//還原成十進制數
        output[i]+=(int)(a[i].real()+0.5);//注意精度誤差
        output[i+1]+=output[i]/10;
        output[i]%=10;
    }
    int i;  
    for(i=l1+l2;!output[i]&&i>=0;i--);//去掉前導零
    if(i==-1)printf("0");//特判長度為0的情況
    for(;i>=0;i--){//輸出這個十進制數
        printf("%d",output[i]);
    }
    putchar('\n');
    return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
[2017.12.1]在此補充一下NTT(快速數論變換)的多項式乘法的代碼。
 
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
//#include<complex>
using namespace std;
//typedef complex<double> cd;
 
typedef long long LL;
 
void exgcd(int a,int b,int& x,int& y){
    if(b==0){
        x=1;
        y=0;
        return;
    }
    int x0,y0;
    exgcd(b,a%b,x0,y0);
    x=y0;y=x0-int(a/b)*y0;
}
 
int Inv(int a,int p){
    int x,y;
    exgcd(a,p,x,y);
    x%=p;
    while(x<0)x+=p;
    return x;
}
 
int qpow(int a,int b,int p){
    if(b<0){
        b=-b;
        a=Inv(a,p);
    }
    LL ans=1,mul=a%p;
    while(b){
        if(b&1)ans=ans*mul%p;
        mul=mul*mul%p;
        b>>=1;
    }
    return ans;
}
 
#define maxn (65537*2)
const int MOD=479*(1<<21)+1,G=3;
 
int rev[maxn];
void get_rev(int bit){
    for(int i=0;i<(1<<bit);i++){
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    }
}
 
//from internet
//for(int i=0; i<NUM; i++)  
//{  
//    int t = 1 << i;  
//    wn[i] = quick_mod(G, (P - 1) / t, P);  
//}  
 
LL a[maxn],b[maxn];
void ntt(LL* a,int n,int dft){
    for(int i=0;i<n;i++){
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    }
    for(int step=1;step<n;step<<=1){
        LL wn;
        wn=qpow(G,dft*(MOD-1)/(step*2),MOD);
        for(int j=0;j<n;j+=step<<1){
            LL wnk=1;//這里一定要用long long不然會迷之溢出
            for(int k=j;k<j+step;k++){
                LL x=a[k]%MOD,y=(wnk*a[k+step])%MOD;//這里也要用long long
                a[k]=(x+y)%MOD;
                a[k+step]=((x-y)%MOD+MOD)%MOD;
                wnk=(wnk*wn)%MOD;
            }
        }
    }
    if(dft==-1){
        int nI=Inv(n,MOD);
        for(int i=0;i<n;i++)
            a[i]=a[i]*nI%MOD;
    }
}
 
#include<cstring>
char s1[maxn],s2[maxn];
 
int main(){
    //scanf("%*d");
    scanf("%s%s",s1,s2);
    int l1=strlen(s1),l2=strlen(s2);
    for(int i=0;i<l1;i++)a[i]=s1[l1-i-1]-'0';
    for(int i=0;i<l2;i++)b[i]=s2[l2-i-1]-'0';
    int bit,s=2;
    for(bit=1;(1<<bit)<(l1+l2-1);bit++)s<<=1;
    get_rev(bit);ntt(a,s,1);ntt(b,s,1);
    for(int i=0;i<s;i++)
        a[i]=a[i]*b[i]%MOD;
    ntt(a,s,-1);
    for(int i=0;i<s;i++){
        a[i+1]+=a[i]/10;
        a[i]%=10;
    }
    int cnt=s;
    while(cnt>=0 && a[cnt]==0)cnt--;
    if(cnt==-1)printf("0");
    for(int i=cnt;i>=0;i--){
        printf("%d",a[i]);
    }
    putchar('\n');
    return 0;
--------------------- 
 
捕魚駕到 天涯明月刀ol