FFT

首先要說明一個誤區，很多人認為FFT只是用來處理多項式乘的，其實FFT是用來實現多項式的係數表示法和點值表示法的快速轉換的，所以FFT的用處遠不止多項式乘。

FFT的前置知識：點值表示法，複數運算，三角函式。

多項式的係數表示法和點值表示法

係數表示法

\[A(x)=\sum_{i=0}^n a_i*x^i \]

點值表示法

不妨將A視為關於x的函式，點值表示法就是在A的影象上取n個點，則該多項式可以被這n個點唯一確定。

\[a_i=A(x_i)\\ A(x)=\{(x_1,a_1),(x_2,a_2),\dots,(x_n,a_n)\} \]

點值表示法有什麼好處呢？

我們知道係數表示法下兩多項式相乘是\(O(n^2)\)

，但在點值表示法下奇蹟出現了：

\[A(x)=\{(x_1,a_1),(x_2,a_2),\dots,(x_n,a_n)\} \\ B(x)=\{(x_1,b_1),(x_2,b_2),\dots,(x_n,b_n)\} \\ A(x)*B(x)=\{(x_1,a_1*b_1),(x_2,a_2*b_2),\dots,(x_n,a_n*b_n)\} \]

顯然這個是可以O(n)實現的。雖然但是，我們幾乎不會在計算中用到點值表示法，但這也給了我們一個解決多項式乘的思路。係數轉點值，相乘，點值轉系數。

又很顯然，我們可以隨便取n個數往函式裡帶，可惜這樣又使複雜度回到了\(O(n^2)\)。

於是FFT出現了，FFT使我們可以用\(O(n\log n)\)

的複雜度將係數轉換成一組特殊的點值，並再把點值轉回係數。

複數的計算

簡單的理解，複數就是實數加虛數，多少都知道點虛數吧，沒錯，知道點就夠了。

\[i=\sqrt{-1}\\ z_1=a+bi\\ z_2=c+di\\ z_1+z_2=a+c+(b+d)i\\ z_1-z_2=a-c+(b-d)i\\ z_1*z_2=ac-bd+(da+bc)i \]

單位根

百度一下

記住以下性質（感興趣可以自己推推，就是基礎的三角函式）

\[\omega_n^k=cosk\frac{2\pi}{n}+sink\frac{2\pi}{n} \\ \omega_n^0=\omega_n^n=1\\ \omega_{2n}^{2k}=\omega_n^k \]

好了，現在你已經掌握了所有FFT的前置知識了，自己來推推FFT吧

。

正式開始FFT

將\(\omega_n^0,\dots,\omega_n^n-1\)這n個數帶入得到點值表示，於是：

\[A(x)=a_0+a_1*x+a_2*{x^2}+a_3*{x^3}+a_4*{x^4}+a_5*{x^5}+ \dots+a_{n-2}*x^{n-2}+a_{n-1}*x^{n-1}\\ A(x)=(a_0+a_2*{x^2}+a_4*{x^4}+\dots+a_{n-2}*x^{n-2})+(a_1*x+a_3*{x^3}+a_5*{x^5}+ \dots+a_{n-1}*x^{n-1})\\ A_1(x)=a_0+a_2*{x}+a_4*{x^2}+\dots+a_{n-2}*x^{\frac{n}{2}-1}\\ A_2(x)=a_1+a_3*{x}+a_5*{x^2}+ \dots+a_{n-1}*x^{\frac{n}{2}-1}\\ A(x)=A_1(x^2)+xA_2(x^2)\\ \]

我們將ωkn(k<n2)ωnk(k<n2)代入得

\[A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})\\ A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}(\omega_n^{2k+n})=A_1(\omega_n^{2k}*\omega_n^n)-\omega_n^kA_2(\omega_n^{2k}*\omega_n^n)=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k}) \]

摘自attak的blog

觀察這兩個式子

\[A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})\\\ A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k})\\ \]

顯然（怎麼又是顯然）只要求出\(A_1(\omega_n^{2k})\)和\(\omega_n^kA_2(\omega_n^{2k})\)就可以得出\(A(\omega_n^k)\)和\(A(\omega_n^{k+\frac{n}{2}})\)，而\(A_1(\omega_n^{2k})\)和\(\omega_n^kA_2(\omega_n^{2k})\)又可以進一步遞迴，正是這一點性質使FFT的複雜度達到了優秀的\(O(n\log n)\)

實現

個人認為這是所有演算法最難也最重要的部分，然而很多blog都是將這部分一筆帶過，所以我決定來詳細的講講。

遞迴寫法

竟然是遞迴，那必然有遞迴極限。每次遞迴多項式的項數剩下一半，只剩一項時\(A(x)=a_1\)與\(x\)無關，所以直接返回就好了。

在求出了\(A_1(\omega_n^{2k})\)和\(\omega_n^kA_2(\omega_n^{2k})\)後做兩次多項式加減可以得出\(A(\omega_n^k)\)和\(A(\omega_n^{k+\frac{n}{2}})\)

由於多次的遞迴和陣列新建和賦值使遞迴寫法的常數大的出奇，所以我們需要更好的寫法。

重難點來了！

遞推做法

首先觀察這張圖。

我們的遞迴做法就是從上向下將原數列對半拆，再合併。

我們的合併順序就是從下到上合併，詳細的說先合併\(a_0和a_4,a_2和a_6,a1和a_5,a_3和a_7,然後a_0,a_4合併好的整體與a_2和a_6合併好的整體再合併\dots\)

觀察最上和最下層的數列的二進位制表示，發現就是將二進位制翻轉了。我們有這個神仙操作可以快速的翻轉二進位制。

rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); ///rev[i]儲存i二進位制翻轉後的數

於是我們可以\(O(n)\)得到最下層的數列，再依次往上推，即不用遞迴也不用開大量陣列，讓程式碼變得飛快！

void fft(cp *a,int n,int inv) 
    int bit=0;
    while((1<<bit)<n)bit++;
    for(int i=0;i<=n-1;i++)
    {
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));	///翻轉二進位制 3(2)=011 to 6(2)=011
        if(i<rev[i]) swap(a[i],a[rev[i]]); 
	}
    for(int mid=1;mid<n;mid*=2)	///mid表示當前合併的行中每一塊的長度的一半 
    {
        cp ur(cos(pi/mid),inv*sin(pi/mid));	///ur代表單位根 
        for(int i=0;i<n;i+=mid*2) ///i表示當前列舉到這一橫行的哪一塊的開頭下標 
        {
            cp tmp(1,0);
            for(int j=0;j<mid;j++,tmp*=ur) ///與遞迴的寫法一樣，合併這一塊和後面一塊。 
            {
                cp x=a[i+j],y=tmp*a[i+j+mid];
                a[i+j]=x+y,a[i+j+mid] 
            }
		}
    }
}