【Learning】多項式乘法與快速傅裏葉變換(FFT)
簡介:
FFT主要運用於快速卷積,其中一個例子就是如何將兩個多項式相乘,或者高精度乘高精度的操作。
顯然暴搞是$O(n^2)$的復雜度,然而FFT可以將其將為$O(n lg n)$。
這看起來十分玄學,因為怎麽看它們的相乘操作都逃不過$O(n^2)$,FFT是如何再減少復雜度的呢?
講到FFT就不可避免地出現公式,但實際上它們都是比較容易理解的。
全局思路
設兩個次數界均為$n$的多項式$\begin{aligned}A(x)&=a_0x^0+a_1x^1+a_2x^2+...+a_{n-1}x^{n-1}\\B(x)&=b_0x^0+b_1x^1+b_2x^2+...+b_{n-1}x^{n-1}\end{aligned}$
那麽我們要求$C=A*B$。
我們把$A$、$B$和$C$的函數圖像畫出來,考慮每一個取值點:$C(x_1)=A(x_1)*B(x_1)$,$C(x_2)=A(x_2)*B(x_2)$,$C(x_3)=A(x_3)*B(x_3)$......
如果我們能對$A$和$B$求出它們在$x_1,x_2,...x_{2n-1}$處的值$A(x_1),A(x_2),...,A(x_{2n-1})$與$B(x_1),B(x_2),...,B(x_{2n-1})$。
那麽將它們一一對應地相乘,就可以得到$C(x_1),C(x_2),...,C(x_{2n-1})$。其中,$C(x)=A(x)*B(x)$。
接著,我們可以利用$C$的這$2n-1$個點,對$C$進行插值,求出$C$的解析式。由於兩個次數界為$n$的多項式相乘後是一個次數界為$2n-1$的多項式,因此我們需要$2n-1$個點,才能對$C$進行準確插值。
問題是,上面三步的時間復雜度分別為$O(n^2)$、$O(n)$和$O(n^2)$,還是沒有什麽改進。它們的名字分別是:DFT,點值乘法,逆DFT。
改進以後,它們分別是FFT,點值乘法,逆FFT。時間復雜度分別為$O(nlgn)$,$O(n)$,$O(nlgn)$。
那就來改進吧(DFT $O(n^2)$ ---> FFT $O(nlgn)$)
本質老是沒有飛躍,多半是廢了
但如果用復數呢?
定義$n$次單位復數根為滿足$\omega^n=1$的復數$\omega$。$n$次單位復數根恰好有$n$個:$\omega_n^0,\omega_n^1,...,\omega_n^{n-1}$,它們的$n$次方都為1。
其中,$\omega_n^0=e^{2\pi i/n}$, $\omega_n^x=(\omega_n^0)^x$, $\omega_n^{n/2}=(e^{2\pi i/n})^{n/2}=e^{\pi i}=-1$。
(以下用$n$來代替之前提到的$2n-1$;用$A$表示一個次數界為$n$的多項式,上文提到的"$A$"和"$B$"的操作都是同理的)
我們把這$n$個$n$次單位復數根作為取值點,求出$A(\omega_n^0),A(\omega_n^1),...,A(\omega_n^{n-1})$。
這一步叫做離散傅裏葉變換(DFT)。對於$k=0,1,...,n-1$,它要求$y_k=A(\omega_n^k)=\sum\limits_{j=0}^{n-1}a_j\omega_n^{kj}$
然而若不利用單位復數根的性質,復雜度仍然是$O(n^2)$的。
單位復數根的性質
這$n$個復數有神秘性質,主要用到三個:
1. $\omega_n^{k+n/2}=\omega_n^k*\omega_n^{n/2}=\omega_n^k*-1=-\omega_n^k$,
什麽意思呢:比如$n=8$時,$n$個單位負數根,$\omega_n^0$和$\omega_n^4$互為相反數,$\omega_n^1$和$\omega_n^5$互為相反數.....也就是$[0,n/2-1)$與$[n/2-1,n)$對應的單位復數根互為相反數。
2.消去引理:$\omega_{an}^{ak}=(e^{2\pi i/an})^{ak}=(e^{2\pi i/n})^k=\omega_n^k$,類似於分數約分。
3.折半引理:$$\begin{aligned}(\omega_n^k)^2&=\omega_n^{2k}=\omega_{n/2}^k\\(\omega_n^{k+n/2})^2=\omega_n^{2k+n}=\omega_n^{2k}*\omega_n^n&=\omega_n^{2k}=\omega_{n/2}^k\end{aligned}$$
什麽意思呢?如果把$n$個單位根分成兩組$\omega_n^0...\omega_n^{n/2-1}$ 和 $\omega_n^{n/2}...\omega_n^{n-1}$,兩兩對應位置的單位根的平方是相同的。
如$n==8$時:
$(\omega_8^0)^2=(\omega_8^4)^2=\omega_{4}^0\\(\omega_8^1)^2=(\omega_8^5)^2=\omega_{4}^1\\(\omega_8^2)^2=(\omega_8^6)^2=\omega_{4}^2\\(\omega_8^3)^2=(\omega_8^7)^2=\omega_{4}^3$
也就是$n$個$n$次單位復數根的平方的集合,等於$n/2$個$n/2$次單位復數根的集合。
多項式的拆分
我們回來看一下$A$可以如何拆分:記$A$的系數為$a_0,a_1,...,a_{n-1}$。
如果我們設$\begin{aligned}A_0(x)&=a_0x^0+a_2x^1+a_4x^2+...+a_{n-2}x^{n/2}\\A_1(x)&=a_1x^0+a_3x^1+a_5x^2+...+a_{n-1}x^{n/2}\end{aligned}$,也就是將$A$的系數奇偶分組,成為兩個次數界為$n/2$的多項式。
那麽有$$A(x)=A_0(x^2)+x*A_1(x^2)$$
我們求的是$A(\omega_n^0),A(\omega_n^1),...,A(\omega_n^{n-1})$,那麽轉換一下就變成求
$$\begin{aligned}
A(\omega_n^0)&=A_0((\omega_n^0)^2)+\omega_n^0*A_1((\omega_n^0)^2)\\
A(\omega_n^1)&=A_0((\omega_n^1)^2)+\omega_n^1*A_1((\omega_n^1)^2)\\
&...\\
A(\omega_n^{n-1})&=A_0((\omega_n^{n-1})^2)+\omega_n^{n-1}*A_1((\omega_n^{n-1})^2)\\
\end{aligned}$$
求解$A_0$和$A_1$在$n$個單位復數根,我們用遞歸實現。
我們發現代入$A_0$和$A_1$的參數是一個單位復數根的平方,這意味著代入$A_0$和$A_1$的單位復數根並沒有$n$個。根據折半引理,代入$A_0$和$A_1$的總共只有$n/2$個不同的數:$\omega_{n/2}^0,\omega_{n/2}^1,...,\omega_{n/2}^{n/2-1}$,因為$(\omega_n^k)^2=(\omega_n^{k+n/2})^2$。
我們像上面把單位復數根分為$[0,n/2)$和$[n/2,n)$兩組,觀察$A(\omega_n^k)$和$A(\omega_n^{k+n/2})$,也就是相對的兩個單位復數根的代入:
\begin{aligned}
A(\omega_n^k)&=A_0(\omega_{n/2}^k)+\omega_n^k*A_1(\omega_{n/2}^k)\\
A(\omega_n^{k+n/2})&=A_0(\omega_{n/2}^k)+\omega_n^{k+n/2}*A_1(\omega_{n/2}^k)\\
&=A_0(\omega_{n/2}^k)-\omega_n^k*A_1(\omega_{n/2}^k)
\end{aligned}
它們長得好像!
這下可好,我們只需要遞歸求解$A_0(\omega_n^0...\omega_n^{n/2-1})$和$A_1(\omega_n^0...\omega_n^{n/2-1})$,就可以求出$A(\omega_n^0...\omega_n^{n-1})$了。
時間復雜度下降的原因就在於,用$n/2$次的遞歸得到的數據,可以求出右半邊的數值。
點值乘法 (呵呵 $O(n)$ 不優化了吧這個)
對於$A$和$B$都進行DFT後,我們對$n$個點值直接相乘,得到$C$的$n$個點值。
IFFT (IDFT $O(n^2)$ ---> IFFT $O(nlgn)$)
如果我們知道$C$的$n$個點值,如何知道$C$的解析式呢?
我們看一下DFT的矩陣形式:$y=V_na$,分別與下式對應:
$$\begin{bmatrix}
y_0\\y_1\\y_2\\y_3\\.\\.\\y_{n-1}
\end{bmatrix}
=
\begin{bmatrix}
1&1&1&1&...&1\\
1&\omega_n&\omega_n^2&\omega_n^3&...&\omega_n^{n-1}\\
1&\omega_n^2&\omega_n^4&\omega_n^6&...&\omega_n^{2(n-1)}\\
1&\omega_n^3&\omega_n^6&\omega_n^9&...&\omega_n^{3(n-1)}\\
...&...&...&...&...&...\\
1&\omega_n^{n-1}&\omega_n^{2(n-1)}&\omega_n^{3(n-1)}&...&\omega_n^{(n-1)(n-1)}
\end{bmatrix}
*
\begin{bmatrix}
a_0\\a_1\\a_2\\a_3\\.\\.\\a_{n-1}
\end{bmatrix}$$
我們所求的是$a$,而$a=yV_n^{-1}$,求出$V_n$的逆矩陣即萬事大吉了。
定理:對於$V_n$,$(k,j)$處的元素為$\omega_n^{kj}$。
而對於$V_n^{-1}$,$(k,j)$處的元素為$\omega_n^{-kj}/n$。
如果想簡單證明的話,將$V_n^{-1}$寫出來算一算就好。(可以參見算導)
那麽$a_j=\frac{1}{n}\sum\limits_{k=0}^{n-1}y_k\omega_n^{-kj}$。
看回上面DFT的算式表達,我們發現它們長得幾乎一樣:IFFT的表達,僅僅是多了一個$\frac{1}{n}$,以及單位復數根的指數取負數。
這就非常棒了:IFFT的程序其實和FFT一樣,只不過單位復數根替換一下,算完以後,每一個數值都除去$n$即可,具體參見代碼解釋。
END
FFT的應用,主要是將問題轉化成如DFT式子的形式,用FFT來進行加速或計算的操作。
附上遞歸版代碼和非遞歸版代碼:
#include <cstdio> #include <vector> #include <cmath> #define max(a,b) ((a)>(b)?(a):(b)) using namespace std; const int N=50010; const double Pi=3.14159265358979323846; struct Comp{//手寫了一個復數類 double a,b; Comp(){a=b=0.0;} Comp(double x,double y){a=x;b=y;} friend Comp operator + (Comp x,Comp y){ return Comp(x.a+y.a,x.b+y.b); } friend Comp operator - (Comp x,Comp y){ return Comp(x.a-y.a,x.b-y.b); } friend Comp operator * (Comp x,Comp y){ return Comp(x.a*y.a-x.b*y.b,x.b*y.a+x.a*y.b); } }; typedef vector<Comp> vc; int A,B,type,len; vc a,b,c; vc fft(vc u,int flag){//flag標識是否為逆FFT int n=u.size(); if(n==1) return u;//規模為1時,只有一個常數項的多項式的FFT就為這個常數,可以直接返回了 Comp w_n=Comp(cos(2*Pi/n),sin(2*Pi/n)),w=Comp(1,0);//算出單位復數根的底w_n;w是用來叠代的,減少計算次數 if(flag) w_n.b*=-1.0;//逆FFT與FFT的不同 vc a0,a1,v; a0.clear(); a1.clear(); v.clear(); for(int i=0;i<n;i++){//系數按奇偶分組 if(i&1) a1.push_back(u[i]); else a0.push_back(u[i]); v.push_back(Comp(0,0)); } //遞歸求解A0和A1 a0=fft(a0,flag); a1=fft(a1,flag); //用一半的數據,綜合算出全部的結果,w在此處不斷乘上w_n,就保證它是w_n的k次方 for(int k=0;k<=n/2-1;k++){ v[k]=a0[k]+w*a1[k]; v[k+n/2]=a0[k]-w*a1[k]; w=w*w_n; } return v; } int main(){ //原題:求兩個多項式相乘後的系數(系數都為整數) scanf("%d%d%d",&A,&B,&type); A++; B++; for(int i=0,x;i<A;i++) scanf("%d",&x),a.push_back(Comp(x,0)); for(int i=0,x;i<B;i++) scanf("%d",&x),b.push_back(Comp(x,0)); len=1;//算出高位補齊len(上文提到的至少需要2n-1個點),並把兩個多項式的次數都擴展到len //代碼裏的len指的是上文提到的n while(len<(max(A,B)*2)) len<<=1; for(int i=A;i<len;i++) a.push_back(Comp(0,0)); for(int i=B;i<len;i++) b.push_back(Comp(0,0)); //求兩個多項式在n個單位復數根的值O(nlgn) a=fft(a,0); b=fft(b,0); //點值乘法 O(n) for(int i=0;i<len;i++) c.push_back(a[i]*b[i]); //對點值乘法的結果進行逆FFT O(nlgn) c=fft(c,1); for(int i=0;i<A+B-1;i++) printf("%d ",(int)(c[i].a/len+0.5));//除去len,四舍五入(這題是整數) return 0; }遞歸版FFT
#include <cstdio> #include <iostream> #include <cmath> #define max(a,b) ((a)>(b)?(a):(b)) using namespace std; const int N=50010; const double Pi=3.14159265358979323846; struct Comp{ double a,b; Comp(){a=b=0.0;} Comp(double x,double y){a=x;b=y;} friend Comp operator + (Comp x,Comp y){return Comp(x.a+y.a,x.b+y.b);} friend Comp operator - (Comp x,Comp y){return Comp(x.a-y.a,x.b-y.b);} friend Comp operator * (Comp x,Comp y){return Comp(x.a*y.a-x.b*y.b,x.b*y.a+x.a*y.b);} }a[N*4],b[N*4]; int A,B,type,n; inline int rev(int x){ int ret=0; for(int i=1;i<n;i<<=1,x>>=1) ret=(ret<<1|(x&1)); return ret; } void fft(Comp *a,int f){ int lg=log2(n),len; Comp w,w_n,u,v; for(int i=0,t;i<n;i++) if(i<(t=rev(i))) swap(a[i],a[t]); for(int i=1;i<=lg;i++){ len=1<<i; w_n=Comp(cos(2*Pi/len),sin(2*Pi/len)*f); for(int j=0;j<n;j+=len){ w=Comp(1,0); for(int k=0;k<=len/2-1;k++){ u=a[j+k]; v=w*a[j+k+len/2]; a[j+k]=u+v; a[j+k+len/2]=u-v; w=w*w_n; } } } } int main(){ scanf("%d%d%d",&A,&B,&type); for(int i=0,x;i<A;i++) scanf("%lf",&a[i].a); for(int i=0,x;i<B;i++) scanf("%lf",&b[i].a); for(n=1;n<A+B;n<<=1); fft(a,1); fft(b,1); for(int i=0;i<n;i++) a[i]=a[i]*b[i]; fft(a,-1); for(int i=0;i<A+B-1;i++) printf("%d\n",(int)(a[i].a/n+0.5)); return 0; }非遞歸版FFT(常數小)
【Learning】多項式乘法與快速傅裏葉變換(FFT)