在學習一項算法以前,咱們總該關心這個算法到底是爲了幹什麼。c++
(如下應用只針對OI)git
一句話:求多項式乘法(固然它的實際用處不少)算法
設多項式函數
\(A(x)=a_0+a_1x+a_2x^2+\ldots+a_nx^n\)學習
\(B(x)=b_0+b_1x+b_2x^2+\ldots+b_mx^m\)優化
咱們想求 \(F(x)=A(x)B(x)=\sum\limits_{i=0}^n\sum\limits_{j=0}^ma_ib_jx^{i+j}\)spa
肉眼可見的複雜度\(O(nm)\).net
而使用通常的FFT則複雜度可以在\(O(nlogn)\) (實則帶比較大的常數,但數據範圍大時確定比\(O(nm)\)好)code
寫這篇博客是爲了讓本身和看到的人可以梳理一遍FFT這個算法。orm
如下FFT的實現使用比較常見的複數單位根作法。
首先要知道複數是什麼,高中選修即有,網上資料也不少,不贅述,這裏只講FFT直接相關。
定義: \(z^n=1\) 那麼 \(z\) 就是n次的單位根
對於一個\(n\)次單位根,就是將複數單位圓等分紅\(n\)份,每一個數 \(w_n^k (k\in[0,n-1) )\) 就被稱做n次的第k個單位根
如圖
圖中表示了全部8次單位根,其中\(w_8^0=1,w_8^1=1+i,w_8^2=i,w_8^3=-1+i\ldots\) 所謂單位根就是這些複數,它們的平方均爲1。
雖然單位根的次數\(n\)在定義上能夠取全體正整數,可是咱們在FFT中只須要考慮n爲2的整數次冪,以方便咱們利用性質。
單位根有一些性質,咱們先列出來看看 (下邊的n均爲2的整次冪!)
\(w_n^i\times w_n^j=w_n^{i+j}\),這個性質實際上不能叫性質,從複數乘法,輻角的概念出發即可以得證。
\(w_{2n}^{2k}=w_n^k\)
理解,能夠看到 \(2k\) 和 \(k\) 在 \(2n\) 和 \(n\) 中佔的比例是相同的,聯想單位圓便可(\(w_8^2\)顯然等於\(w_4^1\),均爲1)
我看了好幾篇博客,包括知乎,對摺半引理的說法是: 對於 \(n\) 等於 \(2\) 的整次冪, \(w_n^k\)通過平方能夠獲得全部\(w_{\frac{n}{2}}^k\)(單純來講只須要\(n\)是二的倍數就好)
兩個小性質均可以證實
\(w_n^{2k}=w_{\frac{n}{2}}^k\)
由消去引理,這個小性質是顯然的。從這個性質出發,若是 \(k\) 取遍\([0,n-1]\),那麼首先右邊是全了,對於左邊,\(2k\in[0,2n-2]\),經歷了一次循環, (按道理這裏須要指出一個循環性質,但我以爲不必) \(w_n^{k+n}=w_n^{k}\times w_n^{n}\), 而\(w_n^n\)即\(w_n^0=1\),因此 \(2k\) 過了\(n-1\)以後會通過\(\frac{n}{2}\)個相同的點
\(w_n^k=-w_n^{k+\frac{n}{2}}\)
旋轉\(180^\circ\) 事情變得顯然;從複數乘法的角度來講也能夠:\(w_n^{k+\frac{n}{2}}=w_n^k\times w_n^{\frac{n}{2}}\),而\(w_n^{\frac{n}{2}}\)剛好爲\(-1\)。
那麼這個性質告訴咱們前一半\(k\in[0,\frac{n}{2}-1]\)和後一半 \(k\in[\frac{n}{2},n-1]\)的平方是同樣的,因此同理,兩半的平方值是同樣的,也就是兩部分平方只肯定了\(\frac{n}{2}\)個值,那麼利用前邊第一個小性質,咱們不須要循環性質就可以證實折半引理。
折半引理告訴咱們,若是知道了n的全部單位根,那麼\(\frac{n}{2}\)的單位根咱們能夠知道,反之,咱們利用第二個性質也能夠從\(\frac{n}{2}\)的全部單位根,加一個負號,就能推出全部的 \(n\) 次單位根。
\(\sum\limits_{i=0}^{n-1}w_n^k=0\)
從折半引理看,咱們把 \(w_n^k\) 分紅兩部分\(w_n^a (a\in[0,\frac{n}{2}-1] ),w_n^b(b\in[\frac{n}{2},n-1])\),也就是圓的上半部分和下半部分,對應有每一個\(w_n^a+w_n^b=0\),其實是\(w_n^a+w_n^{a+\frac{n}{2}}=0\),因此兩部分的累加和就是0(咱們上邊已經假設n是2的整次冪了,因此必定是相同大小的兩部分)
首先簡單說一下輻角,在複數平面上,從實數軸正半軸逆時針旋轉到複數\(z\)到原點的連線,中間通過的角度就叫作輻角(大小在\(2\pi\)以內)
從圓的知識類比一下,咱們能夠知道
在敘述 \(n\) 次單位根這件事情的時候,咱們是將複數單位圓 \(n\) 等分以後獲得單位根,這樣來看的話,顯然\(w_n^1\)的輻角是 \(\theta=\frac{2\pi}{n}\), 那麼\(w_n^k=(w_n^1)^k\) ,那麼\(w_n^k\)的輻角就是 \(k\theta\),咱們知道在普通二維平面上單位圓上的點表示爲\((cos\theta,sin\theta)\),類比到複平面,就會獲得 \(w_n^k=cos\theta+sin\theta\ i \ \ \ \ (\theta=\frac{2\pi}n)\)
至此,咱們已經瞭解單位根的性質及其表示,然而咱們還不知道爲何要介紹這麼一個東西。下面咱們會了解的。
一開始,咱們作介紹時設了兩個多項式
\(A(x)=a_0+a_1x+a_2x^2+\ldots+a_nx^n\)
\(B(x)=b_0+b_1x+b_2x^2+\ldots+b_mx^m\)
仔細看的話,咱們發現這兩個多項式像什麼?沒錯,像函數,實際上它徹底能夠當成一個函數表達式。
即\(y=A(x)\)
咱們知道,求解一次函數只須要知道兩個點的座標,求解二次函數只須要知道三個點的座標,那麼咱們能夠換個說法,兩個點肯定了一個一次函數,三個點肯定了一個二次函數,一樣的,咱們能夠說 \(n+1\) 個點肯定了一個 \(n\) 次函數,肯定的方法很簡單,咱們帶入聯當即可(想一想咱們是怎麼經過兩個點肯定一條直線的)。
如今咱們把名詞換一下,」函數「變成「多項式」,就能夠改爲說 \(n+1\) 個不一樣點肯定一個 \(n\) 次多項式。
即點集\(((x_0,y_0),(x_1,y_1),\ldots ,(x_n,y_n))\) 可以肯定一個 \(n\) 次多項式。
代入\(y=A(x)\)看的形象一點 ,咱們解下邊的方程組(點是咱們已知且在函數上的),就可以求出來全部的係數 \(a\) (就是和求函數解析式同樣的概念),若是用高斯消元來解這樣一個方程組,複雜度是 \(O(n^3)\)的,這個過程被稱做離散傅里葉逆/反變換(IDFT) 。
同時,若是咱們知道一個多項式的係數表達式,咱們只須要隨便找\(n+1\)個點 \(x\) 代入到\(A(x)\)中,就能夠獲得全部的\(y\)。
這樣咱們就能夠將一個多項式從係數表達式 \(O(n^2)\) 的轉換成點值表達式。
把多項式的係數表示換成點值表示的變換被稱做離散傅里葉變換(DFT)。(不當心先介紹了逆變換)
點值表達式之間想要肯定新的點值表達式,只用將y相乘便可(\(x\)對應相等),好比說
\(A=((x_0,y_0),(x_1,y_1)\ldots,(x_n,y_n))\)
\(B=((x_0,y_1'),(x_1,y_2')\ldots,(x_n,y_n') )\)
\(F=AB=((x_0,y_1y_1'),(x_1,y_2y_2'),(x_2,y_2y_2')\ldots,(x_n,y_ny_n') )\)
由於咱們的每一個\(y_iy_i'\)都是一次多項式乘法的結果,按照定義,\(F(x_k)=A(x_k)B(x_k)=y_ky_k'=\sum\limits_{i=0}^n\sum\limits_{j=0}^ma_ib_jx_k^{i+j}\)咱們的多項式乘法自己是隻關心繫數的,和 \(x\) 無關(就像函數表達式\(y=f(x)\)只和各項係數有關而與具體的點無關同樣),而根據咱們點值表示法的定義 \(y_k=F(x_k)\),因此\((x_0,y_1y_1')\)就是\(F(x)\)上的一個點 \((x_k,F(x_k))\)
一個小問題 一開始咱們知道,多項式\(A(x)B(x)\)的結果理應獲得最高次項爲\(n+m\)的\(F(x)\),而A,B的直接轉化的點值表示法分別只有\(n\)個點和\(m\)個點,最後要真正肯定\(F(x)\),須要 \(n+m+1\)個點怎麼辦?好辦,咱們一開始把A和B多算一些y,補若干個點補到\(n+m+1\)個點,多一些點並不會影響什麼,而這樣就可以創造出\(n+m+1\)個可以肯定\(F(x)\)的點集了。
AB的運算獲得的新點集,就可以肯定新的 多項式 \(F(x)\) 。而在已知 \(y_i\)和 \(y_i'\) 的狀況下,運算AB須要的複雜度是 \(O(n)\)的!
從前置知識點值表示法的內容中,咱們能夠獲得一種新的求多項式 \(F(x)=A(x)B(x)\) 的方法
DFT,將\(A(x)\) \(B(x)\) 經過離散傅里葉變換 (DFT)轉化爲點值表示法,爲了可以真正肯定\(F(x)\),A和B都要轉化爲\(n+m+1\)個點以肯定最高次項爲\(n+m\)的\(F(x)\),複雜度 \(O(n^2)\)
將\(A(x)\)的全部點值和\(B(x)\)的全部點值對應相乘,獲得\(F(x)\)的點值表示法。複雜度 \(O(n)\)
IDFT,以全部 \(x\) 的各個次項作係數矩陣(以\(x_i^1,x_i^2\ldots\)作係數矩陣),進行高斯消元,獲得\(F(x)\)的全部係數。複雜度 \(O(n^3)\)
很是漂亮啊!!咱們一通操做,直接獲得一個O(n³)的作法可以進行多項式乘法!!秒!
咱們進行這麼多操做,研究了點值表示法可以對應多項式的係數表示法的事實,固然不是爲了花裏胡哨獲得一個\(O(n^3)\)的作法
這個時候,咱們就要考慮優化了,優化什麼呢?咱們不忍心點值表示運算明明已經達到了極優秀的\(O(n)\),而DFT和IDFT卻嚴重拖慢運行。
這個時候,咱們就要用到複數單位根了。
平時咱們作函數之類的事情,一般規定 \(x\in R\) 實際上咱們能夠擴展數域, 使得 \(x\in C\) 即讓\(x\)從實數擴展到整個複數域,做爲一個橫座標來使用。而後咱們考慮將DFT的過程當中任選的\(x\),變成有規律的複數單位根,爲了使得這些 \(x\) 有規律咱們還規定,對於一個\(n-1\)次多項式,咱們向多項式 \(A(x)\) 和\(B(x)\) 代入全部的單位根,即代入全部的 \(w_n^k\) ,而全部的單位根是不一樣的複數,因此點集\((w_n^k,A(w_n^k))\) 就可以做爲\(A(x)\)的點值表示法,\(B(x)\)同理。
接下來咱們將 \(n\) 次單位根代入到咱們的多項式中,爲了使用咱們上邊指出過的複數單位根的性質,以及爲了可以惟一肯定咱們的目標多項式,咱們選擇的 \(n\) 應當是2的整次冪,而且它要大於等於n+m+1。
代入以後 \(A(x)=((w_n^0,A(w_n^0),),(w_n^1,A(w_n^1))\ldots,(w_n^{n-1},A(w_n^{n-1}))\),考慮如何快速求出每一個\(A(w_n^k)\),咱們從消去引理知道,只要知道偶數項的\(w_n^k=w_\frac{n}2^{\frac{k}2}\),這提示咱們須要將原來的表達式進行一波奇偶分組,從係數表達式入手,如今\(A(x)\)是一個\(n-1\)次多項式。
設\(A_1(x)=(a_0+a_2x+a_4x^2+\ldots+a_{n-2}x^{\frac{n}2-1}),A_2(x)=(a_1+a_3x+\ldots+a_{n-1}x^{\frac{n}2-1})\)
咱們會獲得\(A(x)=A_1(x^2)+xA_2(x^2)\)
將全部\(x\) 寫成單位根的樣子
對於\(k\in[0,\frac{n}2-1]\),\(A(w_n^k)=A_1(w_\frac{n}2^k)+w_n^kA_2(w_\frac{n}2^k)\)
對於\(k\in[\frac{n}2,n-1]\), 從\(w_n^k=-w_n^{k+\frac{n}{2}},w_n^k=w_\frac{n}2^\frac{k}2\)因爲\(k\)已經大於\(\frac{n}2\),因此這裏全部的\(A_1(w_\frac{n}2^k)=A_1(w_\frac{n}2^{k-\frac{n}2})\),\(A_2\)同理
因此\(A(w_n^k)=A_1(w_{\frac{n}2}^{k-\frac{n}2})-w_n^{k-\frac{n}2}A_2(w_{\frac{n}2}^{k-\frac{n}2})\) 也就是咱們能夠直接引用在前一半已經算出來過的\(A_1A_2\),來計算出後一半的全部\(A(w_n^k)\) ,這樣咱們把計算的規模縮小了一半,同理,咱們要知曉前一半\(A_1\)和\(A_2\)的值,也能夠從一半的前一半推之,即遞歸即可以求解,最後須要進行log層計算,而後只須要在遞歸的底層對一個點進行 \(O(n)\) 的多項式運算而後每層\(O(n)\)的向上合併,最終的複雜度是 \(T(n)=2T(\frac{n}2)+O(n)\) 其最終複雜度是\(O(nlogn)\)的。這樣咱們經過將多項式拆分紅兩部分而後進行類似性的遞歸,使得咱們的DFT可以在\(O(nlogn)\) 的複雜度下獲得多項式的點值表達式,這個過程就被稱做FFT。
如何快速的把獲得的點值表示法的多項式變成係數表示法,先放一個結論吧
也就是說,若是咱們知道了一個多項式,能夠用FFT化做點值表示,而若是咱們知道了全部點值反過來求全部係數,咱們發現求每一個係數的形式和求每一個點值的形式是相同的不一樣的是多了一個係數\(\frac{1}n\)以及在\(w\)上標多了一個負號,而這是不影響單位根的三項引理的,也就說,這個反向的過程依然能夠用FFT解決。如此一來,咱們就能夠一樣的在\(O(nlogn)\) 的複雜將點值表示轉化爲係數表示。
關於IFFT的證實:
咱們把每一個點展開成多項式,即把\((w_n^k,A(w_n^k) )\)還原成多項式的模樣而後以多項式的係數做爲列向量相乘,就可以獲得全部的\(y\)
咱們的點值表示法所獲得的矩陣是一個範德蒙矩陣,而咱們進行IFFT的過程就是已知向量 \(y\) 以及點矩陣\(W\) 求向量\(a\)
那麼咱們只須要計算 \(yW^{-1}\) 便可,對於這個範德蒙矩陣求逆。
直接構造矩陣S
根據求和引理,\(c_{i,j}=w_{i,k}s_{k,j}\)兩個矩陣相乘就會獲得矩陣C
這個矩陣乘上一個\(\frac{1}n\)就是單位矩陣,因此以\(w_n^{-k}\)做爲點的橫座標依次代入矩陣再乘一個\(\frac{1}n\)就是咱們須要的\(W^{-1}\),因此只須要咱們再用\(w_n^{-k}\)作一遍FFT,將結果乘\(\frac{1}n\)便可。
這樣一來,咱們用點值表示法就可以\(O(nlogn)\)快速求出係數表示了。
蝴蝶變換是FFT的一個優化,之因此要優化,是由於純遞歸求解FFT效率並不高,因此考慮如何迭代求解。
每次奇偶分組的時候,咱們都會挑一個當前所在位置編號\(i\) 爲偶數的分到左邊,爲奇數的分到右邊。
偶數的偶數的偶數次項到最後,就會是0,多一個奇數就會在頭多1...(感性理解)
總之從變換的規律來看,咱們遞歸到最後,係數的分組是原來係數的二進制反轉
好比咱們有一個7次多項式,有8項係數\(a_0-a_7\)
分組獲得的就是
你看最後一組的分配,就是每一個數二進制反轉的結果。
換言之,咱們一開始若是將係數二進制反轉,而後從下往上推,每次都能推出一個長度更長的多項式,這樣子問題就能夠一步步向上推出全局的答案。
終於,咱們獲得了整個FFT的算法(typroa顯示共4862詞(不包括下邊的代碼),碼了一天,淦)!下面會在例題裏貼上一個模板,而後隨着作題會更新一些題目總結吧。
【學習筆記】超簡單的快速傅里葉變換(FFT)(含全套證實)
單位根-百度百科
OIwiki(裏邊還有較詳細的二進制反轉證實,我懶了不想寫了)
百度上的傅里葉變換(主要看嚴謹概念來的,然鵝看不大懂,而後題目就多了個【OI向】)
[學校的非公開資料]這個就莫得鏈接
代碼以下:
#include<bits/stdc++.h> #define reg register typedef long long ll; using namespace std; inline int qr(){ int x=0,f=0;char ch=0; while(!isd igit(ch)){f|=ch=='-';ch=getchar();} while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();} return f?-x:x; } const int maxn=3e6+199; const double PI=3.14159265358979; int n,m; struct Complex{ double x,y; Complex operator+(const Complex &t)const { return (Complex){x+t.x,y+t.y}; } Complex operator-(const Complex &t)const{ return (Complex){x-t.x,y-t.y}; } Complex operator* (const Complex &t)const { return (Complex){x*t.x-y*t.y,x*t.y+y*t.x}; } }a[maxn],b[maxn]; int rev[maxn],bit,tot;//反轉的二進制數 void fft(Complex c[],int inv){//inv表示的是單位根轉動的方向 for(reg int i=0;i<tot;i++){//蝴蝶變換! if(i<rev[i]){//小於rev[i]時交換一次就好 swap(c[i],c[rev[i]]); } } for(reg int mid=1;mid<tot;mid<<=1){//從底層開始枚舉,即從最小的子問題開始 Complex w1=(Complex){cos(PI/mid),inv*sin(PI/mid)};//求得單位根 for(reg int i=0;i<tot;i+=mid<<1){//枚舉每個子問題 Complex wk=(Complex){1,0}; for(reg int j=0;j<mid;j++,wk=wk*w1){//掃左半邊,此時大項目A的第i+j項就是原來處理出的左半邊和右半邊的結合,即以前的第i+j項和i+j+mid的組合,模擬遞歸 Complex x=c[i+j],y=wk*c[i+j+mid];//A1和A2 c[i+j]=x+y,c[i+j+mid]=x-y;//左半邊和右半邊的處理 }//c存儲點值A(wk)=A1*(wk)+wk*A2(wk) } } } int main(){ // freopen(".in","r",stdin); // freopen(".out","w",stdout); n=qr();m=qr(); for(reg int i=0;i<=n;i++) a[i].x=qr(); for(reg int j=0;j<=m;j++) b[j].x=qr(); while((1<<bit)<n+m+1) bit++;//找到第一個大於n+m+1的2的整次冪 tot=1<<bit;//擴展多項式長度 for(reg int i=0;i<tot;i++){ rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));//記錄每一個數的二進制反轉 } //rev記錄的是每一個二進制數的二進制逆轉,咱們這個遞推是求出來除去第0位的逆再把第0位放到最高位,就可以獲得這個數的逆 fft(a,1),fft(b,1);//轉化爲點 for(reg int i=0;i<tot;i++) a[i]=a[i]*b[i];//求出目標多項式的點值表示 fft(a,-1);//求逆轉化 for(reg int i=0;i<=n+m;i++){ printf("%d ",(int)(a[i].x/tot+0.5));//除去tot,就像咱們以前分析的 ,+0.5來以防精度出問題 } return 0; }