FFT入門

FFT(快速傅里葉變換)是上個學期學會的東西，由于接下來要玩母函數，所以現在寫篇博客復習一下。FFT可以在O(nlogn)$O(n/log n)$的時間內完成多項式乘法。

問題

給定兩個十進制數（105$10^5$位），求它們的乘積。

不妨把它歸為一個多項式乘積的問題：每一位都是一個系數，那么1234*2333就變成了：

(x3+2x2+3x+4)?(2x3+3x2+3x+3)$/displaystyle (x^3+2x^2+3x+4)*(2x^3+3x^2+3x+3)$

對于多項式乘積問題，顯然跑O(n2)$O(n^2)$的樸素高精度乘法是過不了的。考慮我們計算a×b$a /times b$的方式：令X$X$為答案，則有：Xn=∑i=0nai×bn?i$/displaystyle X_n = /sum_{i=0}^n a_i/times b_{n-i}$

因此，我們做的事情是直接計算答案的每一位。現在我們換一種思路。

點值表達

之前我們使用的(x5+233x3+x)$(x^{5}+233x^3+x)$這種表達方式，被稱為系數表達。因為它給出了系數向量：(1,0,233,0,1,0)$(1,0,233,0,1,0)$。

但是考慮這樣一個事實：給定n$n$個點，可以唯一確定一個n?1$n-1$次多項式函數。至于如何確定，有高斯消元和拉格朗日插值法。

因此我們擁有了一種全新的表達多項式的方式：點值表達。給出n+1$n+1$個點，可以表達一個多項式。例如：(0,0),(1,1),(2,4)$(0,0),(1,1),(2,4)$是多項式(x2)$(x^2)$的一種點值表達。一個多項式有無數組點值表達。

有什么用呢？ 考慮多項式的加法A+B$A+B$，生成的多項式的點值表達，可以由A$A$和B$B$的點值表達得到。在A$A$和B$B$上取相同的一些x$x$，求出對應的點值表達：A:(x1,ya1),(x2,ya2)?$A:(x_1,y_{a1}),(x_2,y_{a2})/cdots$B:(x1,yb1),(x2,yb2)?$B:(x_1,y_{b1}),(x_2,y_{b2})/cdots$

則A+B:(x1,ya1+yb1),(x2,ya2+yb2)?$A+B :(x_1,y_{a1}+y_{b1}),(x_2,y_{a2}+y_{b2})/cdots$

看圖很好理解：

我們把從點值表達變成系數表達的過程乘坐插值。

那么多項式的乘法也類似。大概是這樣的步驟：

單位復數根

對著數學書腦補一番就解決了。

單位復數根：ωkn$ω_n^k$均勻分布在以(0,0)$(0,0)$為中心的復平面圓上。n=8$n=8$時長成這樣：

ωkn=e2πikn$ω_n^k=e^{2πi/frac{k}{n}}$（注意其中的i$i$是虛數單位）。歐拉告訴我們，eiu=cos(u)+isin(u)$e^{iu}=cos(u)+isin(u)$，故有：ωkn=cos(2πk/n)+isin(2πk/n)$ω_n^k=cos(2πk/n)+isin(2πk/n)$.

看圖應該能腦補出來：盡管ωkn$ω_n^k$有n$n$種取值，(ωkn)2$(ω_n^k)^2$只有n2$/frac{n}{2}$種取值。

在圖中的意義是：第一象限點的平方與第三象限點的平方一致；第二象限點的平方與第四象限點的平方一致。因為(a+b)2=(?a?b)2$(a+b)^2=(-a-b)^2$。

由于這貨的性質，我們選擇單位復數根作為x$x$坐標進行求值。

FFT

關鍵步驟：將A=∑i=0naixi$/displaystyle A=/sum_{i=0}^n a_i x^i$拆分為：A0=a0,a2x,a4x2,?,an?2xn2?1$A_0={a_0,a_2x,a_4x^2,/cdots,a_{n-2}x^{/frac{n}{2}-1}}$A1=a1,a3x,a5x2,…,an?1xn2?1$A_1={a_1,a_3x,a_5x^2,…,a_{n-1}x^{/frac{n}{2}-1}}$ 則有：A(x)=A0(x2)+xA1(x2)$A(x)=A_0(x^2)+xA_1(x^2)$

由于我們使用單位復數根進行求值，則x2$x^2$只有n/2$n/2$種取值。我們把問題規模成功降低了一半！

那么如何插值回來呢？令ωkn=cos(2πk/n)?isin(2πk/n)$ω_n^k=cos(2πk/n)-isin(2πk/n)$即可。（這里變成減號）

上面那段話并不清楚，因為我太弱了，根本講不清楚。要完全理解上面的話，推薦看完代碼，然后啃一啃《導論》。

所以我們就寫出了代碼：uoj #34.多項式乘法

#include <cstdio>#include <cstdlib>#include <cmath>#include <complex>#include <iostream>using namespace std;typedef complex<double> cp;                     //complex庫void fft(cp *a,int n,int flag)                  //作用：求出a的點值表達，存進a{    int i;    cp a0[n/2+1],a1[n/2+1];    if(n==1) return;    cp w_n(cos(2*M_PI/n),sin(flag*2*M_PI/n));   //flag=1:求值  flag=2:插值    cp w(1,0);    for(i=0;i<n/2;i++) a0[i]=a[i*2],a1[i]=a[i*2+1];     //分治    fft(a0,n/2,flag);    fft(a1,n/2,flag);    for(i=0;i<n/2;i++)    {        a[i]=a0[i]+w*a1[i];        a[i+n/2]=a0[i]-w*a1[i];        w=w*w_n;                                //遞推單位復數根    }}cp x[300005]={0},y[300005]={0};int n,m;void init(){    int i;    scanf("%d%d",&n,&m);    for(i=0;i<=n;i++) cin>>x[i].real();    for(i=0;i<=m;i++) cin>>y[i].real();    m+=n;    for(n=1;n<=m;n=n*2);}int main(void){    int i;    init();    fft(x,n,1);                                 //求值    fft(y,n,1);                                 //求值    for(i=0;i<n;i++) x[i]=x[i]*y[i];            //點值乘法    fft(x,n,-1);                                //插值    for(i=0;i<=m;i++)        跑了4460ms。提交記錄  遞歸版的比較慢（廢話），迭代版的跑得比香港記者還快。  
然而我并不會蝴蝶操作那套理論，請看riteme的博客：有關多項式的算法
相關資料
riteme 快速數論變換(NTT) 簡要介紹了快速數論變換。xlightgod UOJ34 多項式乘法適合入門FFT。我就是在那里學習了一個。iamzky 快速傅里葉變換講解很詳細。