快速傅里葉變換(Fast Fourier Transform,FFT)

前言 關(guān)于為什么又開始在這個(gè)博客上寫算法

FFT是我高中學(xué)競賽階段接觸到的最后一個(gè)算法，也是我高中一直沒有學(xué)會的算法之一，在高中畢業(yè)的時(shí)候認(rèn)為自己大學(xué)本科一定學(xué)的是CS，所以想要把這些自己之前沒有學(xué)會的算法都再琢磨一下，可是基于當(dāng)時(shí)的知識儲備我還完全不知道復(fù)數(shù)域和復(fù)平面以及多項(xiàng)式的各種定理，于是半途而廢了，F(xiàn)FT也成了在算法學(xué)習(xí)上我的一個(gè)心結(jié)之一。
但是巧合的是，我大學(xué)本科并沒有學(xué)CS，而是學(xué)了數(shù)學(xué)，這些我一直未曾學(xué)會的算法也就被一直擱置了兩年，直到昨天我在上數(shù)據(jù)結(jié)構(gòu)課時(shí)，數(shù)學(xué)院數(shù)據(jù)結(jié)構(gòu)課所學(xué)內(nèi)容跟OI比起來簡直不值一提，于是無聊的我突發(fā)奇想翻開了自己之前洛谷的提交記錄，發(fā)現(xiàn)FFT的版子還一直放在我的題單，于是就想利用這些時(shí)間把這些之前沒學(xué)懂的算法解決掉好了，也為明年考研跨考CS打下基礎(chǔ)。
這次在數(shù)學(xué)院歷練了兩年多的我拿著高等代數(shù)數(shù)學(xué)分析復(fù)變函數(shù)的知識再來看FFT，頓時(shí)覺得這些知識難度與我兩年前看的時(shí)候相比要簡單了許多，在昨天晚上解決了自己所有的疑難點(diǎn)實(shí)現(xiàn)了代碼拿到了AC。
看到我的這個(gè)博客之前發(fā)布的線段樹網(wǎng)絡(luò)流等等算法講解總瀏覽量都接近一萬了，那就繼續(xù)吧，讓這個(gè)斷更了4年的博客重新活起來！


奈芙蓮封面是這個(gè)博客的靈魂！！

問題背景

多項(xiàng)式乘法
給定一個(gè)n次多項(xiàng)式A(x)和m次多項(xiàng)式B(n)，求出F(x)與G(x)的卷積。

輸入格式

第一行兩個(gè)整數(shù)n,mn,m。
接下來一行n+1個(gè)數(shù)字，從低到高表示F(x)的系數(shù)。
接下來一行m+1個(gè)數(shù)字，從低到高表示G(x)的系數(shù)。

輸出格式

一行n+m+1個(gè)數(shù)字，從低到高表示F(x) \cdot G(x)的系數(shù)。

輸入輸出樣例

輸入樣例

1 2
1 2
1 2 1

輸出樣例

1 4 5 2

問題分析

假設(shè)兩個(gè)多項(xiàng)式A(x)=a_0+a_1x+a_2x^2+\cdots +a_{n-1}x_{n-1}和B(x)=b_0+b_1x+b_2x^2+\cdots +b_{n-1}x^{n-1}，兩個(gè)多項(xiàng)式可以寫作：
\begin{equation} A(x)=\sum_{i=0}^{n-1}a_ix^i\\ B(x)=\sum_{i=0}^{n-1}a_ix^i \end{equation}
傳統(tǒng)方法是利用兩個(gè)多項(xiàng)式的系數(shù)進(jìn)行卷積運(yùn)算，得到一個(gè)2n-2次多項(xiàng)式C(x)：
\begin{equation} C(x)=\sum_{i=0}^{2n-2}c_ix_i \end{equation}
這種卷積運(yùn)算的時(shí)間復(fù)雜度為O(n^2)，顯然在數(shù)據(jù)范圍較大的情況下難以承受，而利用快速傅里葉變換可將時(shí)間復(fù)雜度降為O(nlogn)。

FFT介紹

快速傅里葉變換 (fast Fourier transform),即利用計(jì)算機(jī)計(jì)算離散傅里葉變換（DFT)的高效、快速計(jì)算方法的統(tǒng)稱，簡稱FFT?？焖俑道锶~變換是1965年由J.W.庫利和T.W.圖基提出的。采用這種算法能使計(jì)算機(jī)計(jì)算離散傅里葉變換所需要的乘法次數(shù)大為減少，特別是被變換的抽樣點(diǎn)數(shù)N越多，F(xiàn)FT算法計(jì)算量的節(jié)省就越顯著。
FFT（Fast Fourier Transformation） 是離散傅氏變換（DFT）的快速算法。即為快速傅氏變換。它是根據(jù)離散傅氏變換的奇、偶、虛、實(shí)等特性，對離散傅立葉變換的算法進(jìn)行改進(jìn)獲得的。
——百度百科

要解決的問題是多項(xiàng)式的乘法，而一個(gè)多項(xiàng)式的表示方法并不唯一，傳統(tǒng)意義上多項(xiàng)式的表示利用的是系數(shù)表示法，即對于一個(gè)多項(xiàng)式A(x)=a_0+a_1x+a_2x^2+\cdots +a_{n-1}x_{n-1}可由一個(gè)系數(shù)向量(a_0,a_1,a_2,\cdots,a_{n-1})唯一表示。

而除了系數(shù)表示法之外，多項(xiàng)式也可以利用點(diǎn)值表示，對于多項(xiàng)式A(x)，選定n個(gè)x值(x_0,x_1,x_2,\cdots,x_{n-1})帶入多項(xiàng)式進(jìn)行計(jì)算，得到n個(gè)點(diǎn)值\{ (x_0,A(x_0)),(x_1,A(x_1)),(x_2,A(x_2)),\cdots,(x_{n-1},A(x_{n-1})) \}這n個(gè)點(diǎn)值也可唯一表示多項(xiàng)式A(x)。

因此當(dāng)我們利用同一個(gè)向量x得到了兩個(gè)兩個(gè)多項(xiàng)式的點(diǎn)值表示法后，用對應(yīng)點(diǎn)值相乘，得到(A(x_0)B(x_0),A(x_1)B(x_1),A(x_2)B(x_2),\cdots,A(x_{n-1})B(x_{n-1}))即為兩個(gè)多項(xiàng)式的乘積多項(xiàng)式的點(diǎn)值表示，這個(gè)過程的時(shí)間復(fù)雜度為O(n)。

需要注意的一點(diǎn)是，當(dāng)兩個(gè)多項(xiàng)式次數(shù)為n,m時(shí)，他們的乘積多項(xiàng)式次數(shù)為n+m，因此利用點(diǎn)值表示計(jì)算時(shí)，計(jì)算兩個(gè)多項(xiàng)式的點(diǎn)值表示時(shí)應(yīng)選用n+m+1個(gè)變量，才能使得到的結(jié)果唯一表示乘積多項(xiàng)式C(x)。

然而對于一個(gè)多項(xiàng)式A(x)來說，代入任意選定的變量(x_0,x_1,x_2,\cdots,x_{n-1})計(jì)算他的點(diǎn)值表示法時(shí)間復(fù)雜度依然是O(n^2)，并沒有起到優(yōu)化的效果，而快速傅里葉變換解決了這個(gè)問題，使得系數(shù)表示法轉(zhuǎn)化為點(diǎn)值表示法的時(shí)間復(fù)雜度降低為O(nlogn)。

快速傅里葉變換FFT

FFT在計(jì)算多項(xiàng)式的系數(shù)表示法變換為點(diǎn)值表示法時(shí)，選定復(fù)平面上單位圓上的單位復(fù)根\omega_n^k作為變量計(jì)算多項(xiàng)式的點(diǎn)值，在這里單位根\omega_n^k滿足以下的一些性質(zhì)（如果有不理解的可以自行查閱復(fù)數(shù)的一些相關(guān)知識）：
\begin{aligned} \omega_n^k&=e^{\frac{2k\pi}{n}i}=\cos\frac{2k\pi}{n}+i\sin\frac{2k\pi}{n}\\ n&=2^i,i\in N^+\quad k\in [0,n)\\ \omega_n^0&=\omega_n^n=1\\ \omega_n^{\frac{n}{4}}&=i\\ \omega_n^{\frac{3n}{4}}&=-i\\ \omega_{2n}^{2k}&=\omega_{n}^k\\ \omega_{n}^{k+\frac{n}{2}}&=-\omega_n^k\\ \end{aligned}
以上的這些性質(zhì)都可以由Euler公式得到，推導(dǎo)過程并非FFT重點(diǎn)這里就省略了。

對于一個(gè)多項(xiàng)式A(x)=a_0+a_1x+a_2x^2+\cdots +a_{n-1}x_{n-1}，我們可以對其進(jìn)行劃分，將偶數(shù)次項(xiàng)與奇數(shù)次項(xiàng)分開，在這里假設(shè)n-1為奇數(shù)，得到：
A(x)=(a_0+a_2x^2+\cdots +a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots +a_{n-1}x^{n-2})
我們分別定義兩個(gè)多項(xiàng)式A_1(x),A_2(x)：
\begin{aligned} &A_1(x)=a_0+a_2x^1+\cdots +a_{n-2}x^{\frac{n}{2}-1}\\ &A_2(x)=a_1+a_3x^2+\cdots +a_{n-1}x^{\frac{n}{2}-1}\\ \end{aligned}
那么原多項(xiàng)式A(x)就可以表示為：
A(x)=A_1(x^2)+xA_2(x^2)
將\omega_n^k代入上式：
\begin{aligned} A(\omega_n^k)&=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})\\ &=A_1(\omega_{\frac{n}{2}}^k)+\omega_n^kA_2(\omega_{\frac{n}{2}}^k) \end{aligned}
將\omega_n^{k+\frac{n}{2}}代入上式：
\begin{aligned} A(\omega_n^{k+\frac{n}{2}})&=A_1(\omega_{n}^{2k+n})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n})\\ &=A_1(\omega_n^n\times \omega_n^{2k})-\omega_n^kA_2(\omega_n^n\times \omega_n^{2k})\\ &=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k})\\ &=A_1(\omega_{\frac{n}{2}}^k)-\omega_n^kA_2(\omega_{\frac{n}{2}}^k) \end{aligned}
由此可以發(fā)現(xiàn)，A(\omega_n^k)和A(\omega_n^{k+\frac{n}{2}} )在計(jì)算的過程中只有一個(gè)符號不同，因此在進(jìn)行枚舉計(jì)算A(\omega_n^k)時(shí)即可直接得到A(\omega_n^{k+\frac{n}{2}} )的值，利用這種方法進(jìn)行分治，便可以將復(fù)雜度降至O(nlogn)。

逆傅里葉變換IFFT

利用上述的方法得到了乘積多項(xiàng)式的點(diǎn)值表示法，那么現(xiàn)在需要解決的問題時(shí)如何將點(diǎn)值表示法再轉(zhuǎn)換回系數(shù)表示法。

假設(shè)得到多項(xiàng)式的FFT點(diǎn)值表示為(y_0,y_1,y_2,\cdots ,y_{n-1})，其系數(shù)表示為(a_0,a_1,a_2,\cdots ,a_{n-1})，根據(jù)FFT原理，y_k可如下表示：
y_k=\sum_{i=0}^{n-1}a_i(\omega_n^k)^i
取\omega_n^k的共軛復(fù)數(shù)\omega_n^{-k}，如下定義向量(c_0,c_1,c_2,\cdots ,c_{n-1})：
c_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i
那么由定義可以推導(dǎo)出如下的公式：
\begin{aligned} c_k&=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i\\ &=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^{j-k})^i\\ &=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i \end{aligned}
對于復(fù)平面上的單位根\omega_n^k，有如下的性質(zhì)：
\begin{aligned} \sum_{i=0}^{n-1}(\omega_n^{j-k})^i&=0\quad(j\neq k)\\ \sum_{i=0}^{n-1}(\omega_n^{j-k})^i& =\omega_n^ 0=1 \quad (j= k ) \end{aligned}
因此可以得到：
\begin{aligned} &c_k=na_k\\ &a_k=\frac{c_k}{n} \end{aligned}
因此，利用FFT得到了多項(xiàng)式的點(diǎn)值表示后只需要將變量換為原本選定的單位根的共軛復(fù)數(shù)再進(jìn)行一次FFT就能得到多項(xiàng)式的系數(shù)表示。

這里需要說明一個(gè)問題，我們以上的討論都是建立在n為2的冪次的條件下的，那么當(dāng)n不是2的冪次時(shí)需要將n擴(kuò)大為大于n的最小的2的冪次，在進(jìn)行逆傅里葉變換時(shí)，通過上述的推導(dǎo)可以發(fā)現(xiàn)，當(dāng)我們計(jì)算的c_k中k的值大于原本的n時(shí)，就不存在j=k的情況了，因此求得的a_k為0，表示該多項(xiàng)式的k次項(xiàng)系數(shù)為0。

代碼實(shí)現(xiàn)

下面是利用遞歸實(shí)現(xiàn)FFT的代碼，具體的解釋見代碼注釋好了。

#include<bits/stdc++.h>
#define rep(i,a,b) for(register int i=a;i<=b;++i)
using namespace std;

const int N=1e6+10;
const double Pi=acos(-1.0);
int n,m,l=1;

inline int read(){
    int x=0,f=1;
    char c=getchar();
    while(c<'0' || c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c<='9' && c>='0'){x=(x<<3)+(x<<1)+(int)(c-'0');c=getchar();}
    return x*f;
}

struct Complex{
    double re,im;
    Complex (double xx=0,double yy=0){re=xx,im=yy;}
};

Complex operator + (Complex a,Complex b){
    return Complex(a.re+b.re,a.im+b.im);
}

Complex operator - (Complex a,Complex b){
    return Complex(a.re-b.re,a.im-b.im);
}

Complex operator * (Complex a,Complex b){
    return Complex(a.re*b.re-a.im*b.im,a.im*b.re+a.re*b.im);
}

Complex a[N<<1],b[N<<1];

inline void FFT(Complex *a,int t,int inv){//inv表示當(dāng)前進(jìn)行FFT還是IFFT
    if(t==1)return;//計(jì)算長度為1時(shí)回溯
    int mid=t>>1;
    static Complex tmp[N];
    rep(i,0,mid-1)
        tmp[i]=a[2*i],tmp[i+mid]=a[2*i+1];//將多項(xiàng)式按照次數(shù)奇偶進(jìn)行二分
    rep(i,0,t-1)a[i]=tmp[i];
    FFT(a,mid,inv);
    FFT(a+mid,mid,inv);
    Complex wn(cos(2.0*Pi/t),inv*sin(2.0*Pi/t)),w(1,0);//單位根
    rep(i,0,mid-1){
        tmp[i]=a[i]+w*a[i+mid];
        tmp[i+mid]=a[i]-w*a[i+mid];
        w=w*wn;
    }
    rep(i,0,t-1)a[i]=tmp[i];
}

int main(){
    n=read(),m=read();
    rep(i,0,n) a[i].re=read();
    rep(i,0,m) b[i].re=read();
    while(l<=(n+m))l<<=1;//找到大于等于n+m的最小2的冪次
    FFT(a,l,1);
    FFT(b,l,1);
    rep(i,0,l)
        a[i]=a[i]*b[i];
    FFT(a,l,-1);
    rep(i,0,n+m)
        printf("%d ",(int)(a[i].re/l+0.5));//輸出整數(shù)，需要注意精度
    return 0;
}

但是問題到這里并沒有結(jié)束

當(dāng)有的毒瘤數(shù)據(jù)范圍非常大的時(shí)候，用遞歸進(jìn)行計(jì)算時(shí)，大量的遞歸會造成棧溢出，那么是否有不用遞歸的做法？

FFT的迭代實(shí)現(xiàn)

對于這樣一個(gè)序列(a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7)，我們觀察對其進(jìn)行二分的過程：


我們發(fā)現(xiàn)了一個(gè)神奇的性質(zhì)，在對這個(gè)序列進(jìn)行二分以后的序列的二進(jìn)制可以由原序列的二進(jìn)制進(jìn)行翻轉(zhuǎn)得到，那么我們可以利用這個(gè)性質(zhì)，用一個(gè)O(n)的方法可以直接得到最終的序列，從而省去了遞歸的過程，用最終的序列反向遞推實(shí)現(xiàn)即可。

Rader算法

Rader算法即為實(shí)現(xiàn)上述操作的一種算法，對于N個(gè)數(shù)，我們把遞增自然數(shù)(0,1,2,3,\cdots)稱為順序數(shù)列；對順序數(shù)列中的每一個(gè)數(shù)，將其二進(jìn)制倒序后轉(zhuǎn)化為十進(jìn)制，稱為倒序數(shù)列。
對于一個(gè)順序數(shù)列，第i個(gè)數(shù)的二進(jìn)制可以視為將第i/2（這里是整除）個(gè)數(shù)的二進(jìn)制左移一位，再根據(jù)i的奇偶性對其末尾加1或者不加1。
那么要得到它的倒序數(shù)列，只需要將這個(gè)操作反向進(jìn)行即可，即第i個(gè)數(shù)的二進(jìn)制可以視為將第i/2個(gè)數(shù)的二進(jìn)制右移一位，再根據(jù)i的奇偶性對其最高加1或者不加1，這里最高位即為第\log_{2}n位。

迭代進(jìn)行FFT（蝴蝶變換）

利用Rader算法求得了遞推序列，那么如何通過迭代得到最終的答案？
這其實(shí)跟迭代實(shí)現(xiàn)01背包的做法思路差不多，對于求A(x)在n次單位根的各冪次的點(diǎn)值時(shí)，m=n/2次單位根的各冪次在A_1或A_2處的點(diǎn)值已經(jīng)被計(jì)算并且儲存在了A數(shù)組中，那么在下一層的迭代過程中直接使用A數(shù)組存儲的答案繼續(xù)進(jìn)行迭代計(jì)算即可。

迭代優(yōu)化FFT代碼實(shí)現(xiàn)

詳細(xì)解釋依舊見代碼注釋。

#include<bits/stdc++.h>
#define rep(i,a,b) for(register int i=a;i<=b;++i)
using namespace std;

const int N=1e7+10;
const double Pi=acos(-1.0);
int n,m,l=1,t=0;
int bin[N<<1];

inline int read(){
    int x=0,f=1;
    char c=getchar();
    while(c<'0' || c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c<='9' && c>='0'){x=(x<<3)+(x<<1)+(int)(c-'0');c=getchar();}
    return x*f;
}

struct Complex{
    double re,im;
    Complex (double xx=0,double yy=0){re=xx,im=yy;}
};

Complex operator + (Complex a,Complex b){
    return Complex(a.re+b.re,a.im+b.im);
}

Complex operator - (Complex a,Complex b){
    return Complex(a.re-b.re,a.im-b.im);
}

Complex operator * (Complex a,Complex b){
    return Complex(a.re*b.re-a.im*b.im,a.im*b.re+a.re*b.im);
}

Complex a[N<<1],b[N<<1];

inline void FFT(Complex *A,int inv){
    rep(i,0,l-1)
        if(i<bin[i]) swap(A[i],A[bin[i]]);//得到倒序數(shù)列
    for(int mid=1;mid<l;mid<<=1){//枚舉當(dāng)前迭代區(qū)間的中點(diǎn)
        Complex wn(cos(Pi/mid),inv*sin(Pi/mid));//單位根
        for(int r=mid<<1,j=0;j<l;j+=r){//j枚舉迭代區(qū)間位置，用r來得到區(qū)間右端點(diǎn)
            Complex w(1,0);
            rep(k,0,mid-1){//k枚舉當(dāng)前迭代區(qū)間
                Complex x=A[k+j],y=w*A[k+j+mid];
                A[k+j]=x+y;
                A[k+j+mid]=x-y;
                w=w*wn;
            }
        }
    }
}

int main(){
    n=read(),m=read();
    rep(i,0,n) a[i].re=read();
    rep(i,0,m) b[i].re=read();
    while(l<=(n+m)){
        l<<=1,t++;
    }
    rep(i,0,l)//Rader算法
        bin[i]=(bin[i>>1]>>1)|((i&1)<<(t-1));
    FFT(a,1);FFT(b,1);
    rep(i,0,l)
        a[i]=a[i]*b[i];
    FFT(a,-1);
    rep(i,0,n+m)
        printf("%d ",(int)(a[i].re/l+0.5));
    return 0;
}
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

  • 快速傅里葉變換(Fast Fourier Transform,FFT)用來計(jì)算離散傅里葉變換（Discrete F...
    簡為2016閱讀 3,910評論 0贊 3
  • 本章涉及知識點(diǎn)：1、多項(xiàng)式乘法的時(shí)間復(fù)雜度2、多項(xiàng)式的表示：系數(shù)3、多項(xiàng)式的表示：點(diǎn)值4、復(fù)數(shù)的表示5、單位復(fù)數(shù)根...
    PrivateEye_zzy閱讀 10,166評論 0贊 11
  • 【注】參考自算法導(dǎo)論。 1. 簡介 快速傅里葉變換（FFT）是實(shí)現(xiàn)離散傅里葉變換（DFT）和離散傅里葉逆變換（ID...
    BlueHeart0621閱讀 507評論 0贊 1
  • 前言 ??本文是個(gè)人學(xué)習(xí)心得的分享，希望大家在閱讀文章后能在評論中一起學(xué)習(xí)交流！另外還可以訪問我的HelloCUD...
    liadrinz閱讀 7,878評論 0贊 4
  • 概述 ??希爾伯特空間是一個(gè)完備的內(nèi)積空間，其標(biāo)準(zhǔn)正交函數(shù)系，直觀來看就是向量空間中基的延伸。其為基于任意正交系上...
    殉道者之花火閱讀 2,096評論 1贊 5

友情鏈接更多精彩內(nèi)容