某知名選手：出多項式題的人就像在販毒，做多項式的人就像在嗑葯。
一直就想寫關於嗑葯的內容了，但是由於嗑葯所需要的時間很久，而且我沒有大塊的時間來寫一篇真正入門的東西，所以一直咕咕咕。
直到現在，為了自我復習整理一遍思路，寫了一篇真正入門的FFT教程。
話不多說，直接進入正題。
一.所需前置芝士
1.多項式是啥？
形如$f(x)=ax^4+bx^3+cx^2+dx+e$這樣的式子，即：$\sum_{i=0}^{n}a_ix^i$叫做多項式。
定義：最高次項為n的多項式叫做n次多項式。
推論：由於常數項的存在，n次多項式最多有(n+1)項。
2.復數是啥？
正常高中所有的知識都是在實數的基礎下盡心運算，但在算法中僅靠實數這些數遠遠不夠。無法達到算法的目的。因此引入數域最大的復數。（注意，復數就是所說的虛數）。
復數的例子：實數中無法表示$\sqrt{-1}$，而這個數在復數中是真實存在的。記作$i$；即：$i^2=-1$
設a,b是實數，那么形如$a+ib$的數叫復數。其中$i$被稱為虛數單位，復數域是目前已知最大的域。
在復平面中，x代表實數，y軸（除原點外的點）代表虛數，從原點(0,0)到(a,b)的向量表示復數$a+bi$
模長：從原點(0,0)到點(a,b)的距離，即$\sqrt{a^2+b^2}$。
幅角：假設以逆時針為正方向，從x軸正半軸到已知向量的轉角的有向角叫做幅角。
復數的運算：
1.加法：實數部相加，虛數部相加。即：$(a+ib)+(c+id)=(a+c)+i(b+d)$
2.減法：實數部相減，虛數部相減。即：$(a+ib)-(c+id)=(a+c)-i(b+d)$
3.乘法：$(a+ib)*(c+id)=ac+iad+ibc+i^2bd=ac-bd+i(ad+bc)$
單位根：在復平面上，以原點為圓心，1為半徑作圓，所得的圓叫單位圓。以圓點為起點，圓的n等分點為終點，做n個向量，設幅角為正且最小的向量對應的復數為$\omega_n$，稱為n次單位根。
注意，上文單位根中所說的n等分中的n必須是2的正整數次冪。
根據復數乘法的運算法則，其余n-1個復數為$\omega_n^2,\omega_n^3,\ldots,\omega_n^n$
單位根的性質：
1.$\omega_n^0=\omega_n^n=1$ 意義：在x實數正半軸上，長度為1。
2.根據復數的定義，我們將復數的實部和虛部作為向量進行運算，也就是說：$\omega_n^k=\cos\theta+i\sin\theta $，而$\theta=2*k\frac{2\pi}{n}$，所以$\omega_n^k=\cos(k\frac{2\pi}{n})+i\sin(k\frac{2\pi}{n})$
3.若z的n次方為1，那么就叫z為n次單位根。
4.$\omega_n^k*\omega_n^k=(\cos(k\frac{2\pi}{n})+i\sin(k\frac{2\pi}{n}))*(\cos(k\frac{2\pi}{n})+i\sin(k\frac{2\pi}{n}))$
$=\cos^2(k\frac{2\pi}{n})-\sin^2(k\frac{2\pi}{n})+i(2\cos(k\frac{2\pi}{n})\sin(k\frac{2\pi}{n}))$
$=\frac{1+\cos (2k\frac{2\pi}{n})}{2}-(\frac{1-\cos (2k\frac{2\pi}{n})}{2})+i\sin(2k\frac{2\pi}{n})$ (根據高中所學的三角降冪公式)
$=\cos (2k\frac{2\pi}{n})+i\sin(2k\frac{2\pi}{n})$
$=\omega_n^{2k}$
5.消去引理：$\omega_n^k=\omega_{2n}^{2k}$
證明：$\omega_{2n}^{2k}=\cos(2k\frac{2\pi}{2n})+i\sin(2k\frac{2\pi}{2n})=\cos(k\frac{2\pi}{n})+i\sin(k\frac{2\pi}{n})=\omega_n^k$
推論：$\omega_n^k=\omega_{dn}{dk}$
6.折半引理：$\omega_{n}^{k+\frac{n}{2}}=-\omega_n^k$
證明：$\omega_n^{\frac{n}{2}}=\cos(\frac{n}{2}*\frac{2\pi}{n})+i\sin(\frac{n}{2}*\frac{2\pi}{n})$
$=\cos \pi+i\sin \pi=-1$
所以$\omega_{n}^{k+\frac{n}{2}}=\omega_n^k *\omega_n^{\frac{n}{2}}=-\omega_n^k$
二.快速傅里葉變換走起
1.點值表示法
我們知道，兩個多項式相乘也就是卷積，正常來算肯定是$O(n^2)$的，而且乍眼一看似乎沒有什么優化方法。但是，就是有人研究出了不損失正確性的$O(nlogn)$的算法，這個研究過程我們稍微提及一下。
首先，我們平常接觸的形如$f(x)=ax^4+bx^3+cx^2+dx+e$的式子表示唯一一個多項式叫做系數表示法。其次，n個點確定唯一一個n-1次多項式，那么我們用n個點也可以唯一表示一個n-1次多項式，這種表示方法叫做點值表示法。
對於兩個用點值表示法表示的多項式如:$A(x)=((x_0,y_a0),(x_1,y_a1),......,(x_n,y_an))$、$B(x)=((x_0,y_b0),(x_1,y_b1),......,(x_n,y_bn))$相乘，那么相乘結果的多項式的點值表示法就是$C(x)=((x_0,y_a0*y_b0),(x_1,y_a1*y_b1),......,(x_n,y_an*y_bn))$，而這樣求得乘積的復雜度是$O(n)$的，原理利用一次函數或者二次函數自己體會就能弄明白。(注意，為了能確定乘積唯一一個多項式，我們定義n為乘積所得多項式$C()$的次數，這樣可以保證得到n+1個點確定唯一一個n次多項式。也就是說，多項式$A()$和$B()$所取的點的個數要相等且等於$(兩個多項式的次數和-2)$)。
那么問題轉換為將多項式系數表示法轉化成點值表示法。
朴素系數轉點值的算法叫DFT（離散傅里葉變換）,優化后為FFT（快速傅里葉變換），點值轉系數的算法叫IDFT（離散傅里葉逆變換），優化后為IFFT(快速傅里葉逆變換)。
對於DFT，想必初中生都會，也就是O(n^2)的暴力取值然后帶入多項式計算。至於IDFT？高斯消元也是可以做的。
那么FFT呢？接下來的部分便主要介紹FFT。
2.FFT快速傅里葉變換
還記得復數嗎？不記得了？趕快去上面翻一翻復數的那些性質，我們在下文會經常的用到。
假設存在一個多項式：$A(x)=a_0x^0+a_1x^1+a_2x^2+......+a_{n-2}x^{n-2}+a_{n-1}x^{n-1}$
然后我們按照下標奇偶性分成兩部分：$A(x)=(a_0x^0+a_2x^2+a_4x^4+......)+(a_1x^1+a_3x^3+a_5x^5+......)$
我們定義兩個多項式，$A_1(x)=(a_0x^0+a_2x^1+a_4x^2+......)$，$A_2(x)=(a_1x^0+a_3x^1+a_5x^2+......)$ （注意：x的次數和$A()$中x的次數不一樣）
我們根據初中知識可以知道$A()$和$A_1()$、$A_2()$的關系：$A(x)=A_1(x^2)+xA_2(x^2)$
我們將單位根$\omega_n^k(k<\frac{n}{2})$代入上面的式子，那么$A(\omega_n^k)=A_1((\omega_n^k)^2)+\omega_n^kA_2((\omega_n^k)^2)$
$=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})$
然后將$\omega_n^{k+\frac{n}{2}}(k<\frac{n}{2})$再一次代入上面的式子，我們可以推導：$A(\omega_n^{k+\frac{n}{2}})=A_1((\omega_n^{k+\frac{n}{2}})^2)+\omega_n^{k+\frac{n}{2}}A_2((\omega_n^{k+\frac{n}{2}})^2)$
$=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(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k}*\omega_n^n)+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k}*\omega_n^n)$
$=A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k})$
$A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k})-\omega_n^{k}A_2(\omega_n^{2k})$
發現了什么？以上兩種取值帶入該多項式后只有常數項不同，而這兩種取值范圍合起來正好是全域且兩個域的大小相同。
也就是說，我們在計算n個不同$A()$的點值表達的時候，我們可以把問題縮小一半，而且這個問題可以遞歸去做（因為要求的是另外兩個多項式）。於是得到了接近於$O(nlogn)$的獲取n個兩兩不同的點值的算法。
等到遞歸到多項式僅僅有一個常數項時，我們無論給這個多項式什么參數，返回值都是這個常數。因此多項式項數為1的時候（次數為0）結束遞歸，返回這個常數。
具體寫法參考一會出現的的代碼。
FFT完結撒花~（逗你的，但FFT真的完了）
3.IFFT(快速傅里葉逆變換)
我們在做題的時候，很少會有人使用點值表示法來表示一個多項式(你可以試試在做二次函數拋物線的時候拋給老師一個點值表示法的多項式)。所以我們還需要一個告訴的算法，把點值表示法轉換成系數表示法。這就是IFFT要做的。
我們假設：$(y_0,y_1,y_2,y_3,......,y_n)$是一個多項式$F()$在$(b_0,b_1,b_2,b_3,......,b_n)$處用FFT求出來的點值表示。其中多項式$F()$就是上文提到的多項式$A()$和多項式$B()$的乘積，而$b_i$其實就是多項式$F()$的n+1個系數。
有些本文上面的內容沒有消化好的會說:$(y_0,y_1,y_2,y_3,......,y_n)$就是該多項式$F()$的系數。
這就大錯特錯了，希望這么想的人一定要好好閱讀一下FFT的內容后再來學習IFFT。
但是----我們不妨就真的把$(y_0,y_1,y_2,y_3,......,y_n)$當作一個n次多項式$G()$的系數。 （注意，這里的多項式$G()$並不是多項式$A()$和多項式$B()$的乘積，而是新設的一個多項式）
我們假設，$G(k)=\sum_{i=0}^{n}y_i(\omega_n^{-k})^i$
$=\sum_{i=0}^{n}(\sum_{j=0}^{n}b_j(\omega_n^i)^j)(\omega_n^{-k})^i$
$=\sum_{i=0}^{n}\sum_{j=0}^{n}b_i(\omega_n^j)^i(\omega_n^{-k})^i$
$=\sum_{j=0}^{n}b_j\sum_{i=0}^{n}(\omega_n^j\omega_n^{-k})^i$
$=\sum_{j=0}^{n}b_j\sum_{i=0}^{n}(\omega_n^{j-k})^i$
發現什么沒有？沒發現？那么請接着推式子。發現了？那么就請跳過這一段。
我們設$S(x)=\sum_{i=0}^{n}x^i$-------------------①式
將方程兩側同時乘x,那么方程變為：$xS(x)=\sum_{i=0}^{n}x^{i+1}$-------------------②式
用②式減①式，得到：$(x-1)S(x)=x^{n+1}-1$
也就是說，我們得到：$S(x)=\frac{x^{n+1}-1}{x-1}$
我們將$\omega_n^k$代入上面的式子，式子變成了$S(\omega_n^k)=\frac{(\omega_n^k)^{n+1}-1}{\omega_n^k-1}$
然后我們發現了一個事情，那就是這個式子如果n+1變成n就會變得更加美妙又間接。
根據點值表示法的定義，我們清楚：如果(n+1+k)個點都在某個n次多項式上，那么這(n+1+k)個點也可以確定一個n次多項式。($k>0$)
因為n是2的正整數冪，所以我們在選取n的值的時候要滿足，n要嚴格大於所求乘積的多項式的次數+1，這樣就可以保證選取n-1個在該多項式上的點就可以唯一確定該多項式。
這么做有什么用呢？我們可以利用數學歸納法的思想把上文中所有提到的n都-1，也就是說，n變成n-1,而n+1可以變成n。
再來觀察這個式子：$S(\omega_n^k)=\frac{(\omega_n^k)^{n}-1}{\omega_n^k-1}$
$=\frac{(\omega_n^n)^k-1}{\omega_n^k-1}$
$=\frac{1-1}{\omega_n^k-1}$
因為$\omega_n^k$中的k原來取遍0~n,那么現在k只可以取遍0~n-1。
所以當k不等於0的時候$S(\omega_n^k)=\frac{1-1}{\omega_n^k-1}=0$
那么當k等於0的時候呢？顯然，$S(\omega_n^k)=n$
現在再來看這個式子：$G(k)=\sum_{j=0}^{n}a_j\sum_{i=0}^{n}(\omega_n^{j-k})^i$
當$j=k$的時候，后半部分的值為n，而當$j!=k$的時候，后半部分的值為0.
因此可以得到：$G(k)=nb_k$
所以：$b_k=\frac{G(k)}{n}$
至於$G(k)$怎么求呢？別忘了，我們把點值$(y_0,y_1,y_2,y_3,......,y_n)$當作了一個n次多項式$G()$的系數。而$G_k$其實就是多項式$G(k)$的點值。所以求一個已知系數的n-1次多項式的n點值我們用什么？FFT！FFT！FFT！
因此我們再用一邊快速傅里葉變換把多項式$A()$、$B()$快速傅里葉變換后相乘所得多項式$F()$的點值當作另一個多項式$G()$的系數求出$G()$的點值表達式。然后$G()$的n個點值各自除n就是$F()$的n個系數系數。
至此，IFFT完結撒花~
這些推導一定要理解不要死背，要不然只會做FFT的板子啊~。
四.實踐中創新
1.由於c++自帶的復數庫complex太慢，因此我們自己定義復數類：
struct complex { double x,y; complex (double xx=0,double yy=0){x=xx,y=yy;} }a[200010],b[200010]; complex operator + (complex a,complex b){ return complex(a.x+b.x , a.y+b.y);} complex operator - (complex a,complex b){ return complex(a.x-b.x , a.y-b.y);} complex operator * (complex a,complex b){ return complex(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x);}
2.利用FFT運算兩個多項式的卷積：
#include <iostream> #include <cstdio> #include <cmath> #define inc(i,a,b) for(register int i=a;i<=b;i++) using namespace std; const double Pi=acos(-1.0); struct complex { double x,y; complex (double xx=0,double yy=0){x=xx,y=yy;} }a[200010],b[200010]; complex operator + (complex a,complex b){ return complex(a.x+b.x , a.y+b.y);} complex operator - (complex a,complex b){ return complex(a.x-b.x , a.y-b.y);} complex operator * (complex a,complex b){ return complex(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x);} void fft(int nowlimit,complex *now,int type){ if(nowlimit==1) return; complex a1[nowlimit>>1],a2[nowlimit>>1]; for(int i=0;i<=nowlimit;i+=2){ a1[i>>1]=now[i]; a2[i>>1]=now[i+1]; } fft(nowlimit>>1,a1,type); fft(nowlimit>>1,a2,type); complex wn=complex(cos(2.0*Pi/nowlimit),type*sin(2.0*Pi/nowlimit)),w=complex(1,0); for(int i=0;i<(nowlimit>>1);i++,w=w*wn){ now[i]=a1[i]+w*a2[i]; now[i+(nowlimit>>1)]=a1[i]-w*a2[i]; } } int main(){ int n,m; cin>>n>>m; inc(i,0,n) cin>>a[i].x; inc(i,0,m) cin>>b[i].x; int limit=1; while(n+m>=limit) limit<<=1; fft(limit,a,1); fft(limit,b,1); inc(i,0,limit) a[i]=a[i]*b[i]; fft(limit,a,-1); inc(i,0,n+m) printf("%d ",(int)(a[i].x/limit+0.5)); }
我們來觀察一下上面的代碼，不難看出，這是一個遞歸版FFT，因此這個FFT慢到家了，連模板都無法AC，甚至比n^2跑的還慢，但是其中有許多技巧我需要簡單說一說。
2.1.我們來看這段代碼：
complex wn=complex(cos(2.0*Pi/nowlimit),type*sin(2.0*Pi/nowlimit)),w=complex(1,0);
其中注意，wn表示的是在當前區間長度nowlimit下的nowlimit次單位根。也就是說，$wn^{nowlimit}=1$，這意味着把復平面均分成nowlimit份。而nowlimit一定是2的正整數冪。而w表示的就是$\omega_{nowlimit}^0$。
2.2.我們再來看這段代碼：
for(int i=0;i<(nowlimit>>1);i++,w=w*wn){ now[i]=a1[i]+w*a2[i]; now[i+(nowlimit>>1)]=a1[i]-w*a2[i]; }
請注意，由於我們的數組now是一個指針，所以遞歸時改變該層的now數組其實就是改變遞歸上一層時的a1數組或者a2數組。（因為遞歸時我們調用了fft(nowlimit,a1,type)和fft(nowlimit,a2,type)）；
而now數組在第一次遞歸時指向的是原數組本身(系數數組)，所以在fft后，原系數數組就變成了點值表達式數組。
有些人可能會問：我們推出來的式子不是$A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})$嗎？怎么到代碼里就變成$now[i]=a_1[i]+w*a_2[i]$了呢？說好的平方呢？
其實代碼沒有錯，因為我們wn表示的是在當前遞歸層區間長度為nowlimit下的nowlimit次單位根，而在上一遞歸層中的nowlimit是這一層的2倍。因此上一層的wn正好是這一層wn的平方根。也就是說：上一層wn的平方正好是這一層的wn。再換句話說，就是上一層選取的單位根的個數正好是這一層選取單位根個數的兩倍，而這等價於上一層所得到的點值的個數是這一層所得點值個數的兩倍。
而至於$now[i]=a_1[i]+w*a_2[i]$中的i是什么呢？其實i只是一個尋址符。我們把復平面均分成(1<<n)份，從x軸正半軸開始，以逆時針為正方向。我們在每層遞歸時按照正方向的順序依次選取單位根，其中第i個選取的單位根代入多項式所得的答案(重點！)就是$a_1[i]$。$a_2[i]$同理，單位根所選取的數都是一樣的，但答案與$a_1[i]$不同。因為雖然單位根的值相同但這次是代入多項式$a2()$時所得的答案。
2.3.在主函數中，我們進行了一次fft(limit,a,-1)。這是在干什么呢？
之前說過，IFFT中把多項式$A()$、$B()$快速傅里葉變換后相乘所得多項式$F()$的點值當作另一個多項式$G()$的系數求出$G()$的點值表達式。然后$G()$的n個點值各自除n就是$F()$的n個系數系數。而$G(k)=\sum_{i=0}^{n}y_i(\omega_n^{-k})^i$，次數是負數（感性理解一下），因此我們這里的type代入-1。
updata:對於上面一行中次數是負數這一點，我們可以想一想三角函數。因為在x軸上方的角的sin值都為正，x軸下方的角的sin值都為負，y軸右側的cos值都為正，y軸左側的cos值都為負。而且在單位根次數取負數時就相當於取以x軸為對稱軸的那個單位根。所以sin值變為原來的相反數，cos值不變。
3.關於常數優化
常數什么的才不是關鍵呢？（啊呸！這卡的比$n^2$還慢而且還爆棧）
常數優化什么的，如果是寫的如此糟糕的fft的話就一定能做到的吧。
3.1沒錯，我們來看一個聽起來很nb的操作：蝴蝶變換：
for(int i=0;i<(limit>>1);i++,w=w*Wn) { complex t=w*a2[i]; a[i]=a1[i]+t, a[i+(limit>>1)]=a1[i]-t; }
我們發現了什么？沒錯，你沒看錯。僅僅是把$w*a_2[i]$從算兩次變成了算一次，這樣就優化掉了一個大大的復數乘法啦~
3.2還有什么優化呢？比如說遞歸變成循環模擬？
你沒想錯，我們手動用循環模擬遞歸過程，這樣在優化常數的時候同時解決了爆棧這個問題。
但如果我們想要不遞歸，就要提前知道遞歸最底層時每一層的系數是多少，而要知道這個，似乎只有遞歸一條路。是否有其他辦法呢？
我們打表觀察：
發現了啥？沒錯，原序列與后序列在二進制表示下的差別只是翻轉原序列過得到的。
這樣就可以完全避免遞歸。
我們先上代碼，解釋在代碼之后。
#include <iostream> #include <cstdio> #include <cmath> #define inc(i,a,b) for(register int i=a;i<=b;i++) using namespace std; const double Pi=acos(-1.0); struct complex { double x,y; complex (double xx=0,double yy=0){x=xx,y=yy;} }a[400010],b[400010]; complex operator + (complex a,complex b){ return complex(a.x+b.x , a.y+b.y);} complex operator - (complex a,complex b){ return complex(a.x-b.x , a.y-b.y);} complex operator * (complex a,complex b){ return complex(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x);} int limit=1,num=0; int rev[400010]; void fft(complex *now,int type){ inc(i,0,limit-1) if(i<rev[i]) swap(now[i],now[rev[i]]); for(int mid=1;mid<limit;mid<<=1){ //枚舉待合並區間的中點 complex wn=complex(cos(Pi/mid),type*sin(Pi/mid)); for(int r=mid<<1,j=0;j<limit;j+=r){ //r是當前區間的大小，j是當前區間的左端點 complex w=complex(1,0); for(int k=0;k<mid;k++,w=w*wn){ //k是當前在區間的什么位置，w是當前的單位根。 complex x=now[j+k],y=now[j+k+mid]*w; now[j+k]=x+y; now[j+k+mid]=x-y; } } } } int main(){ int n,m; cin>>n>>m; inc(i,0,n) cin>>a[i].x; inc(i,0,m) cin>>b[i].x; while(n+m>=limit) limit<<=1,num++; inc(i,0,limit-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(num-1)); //對於這里的解釋，最好的理解方法就是對着剛才的那張圖自己模擬遞推過程。模擬一遍后就知道遠原理了。 fft(a,1); fft(b,1); inc(i,0,limit) a[i]=a[i]*b[i]; fft(a,-1); inc(i,0,n+m) printf("%d ",(int)(a[i].x/limit+0.5)); }
3.2.1 關於rev[i]
rev[i]表示原序列第i個元素在后序列的位置是rev[i]。而rev[i]的計算自己模擬一遍就能明白
3.2.2 關於得到后序列
我們在得到后序列的時候加了if(i<rev[i])，這是為了交換就交換一遍。如果不寫這段代碼，那么交換后的序列和原序列一模一樣。
3.2.3. 關於$complex wn=complex(cos(Pi/mid),type*sin(Pi/mid));$
我們會發現，為什么Pi不用*2了呢？這是因為mid指的是區間的一半，而我們要除的是整個區間的長度，所以分母的2和分子的2抵消了。
3.2.4. 關於mid
mid指的是當前區間的中點，但要注意，mid要向上取整。也就是說，mid同樣是該區間右半部分的第一個元素。