本文只討論FFT在信息學奧賽中的應用
文中內容均為個人理解，如有錯誤請指出，不勝感激
前言
先解釋幾個比較容易混淆的縮寫吧
DFT：離散傅里葉變換—>$O(n^2)$計算多項式乘法
FFT：快速傅里葉變換—>$O(n*\log(n)$計算多項式乘法
FNTT/NTT：快速傅里葉變換的優化版—>優化常數及誤差
FWT：快速沃爾什變換—>利用類似FFT的東西解決一類卷積問題
MTT：毛爺爺的FFT—>非常nb/任意模數
FMT 快速莫比烏斯變化—>感謝stump提供
多項式
系數表示法
設$A(x)$表示一個$n-1$次多項式
則$A(x)=\sum_{i=0}^{n} a_i * x^i$
例如：$A(3)=2+3*x+x^2$
利用這種方法計算多項式乘法復雜度為$O(n^2)$
（第一個多項式中每個系數都需要與第二個多項式的每個系數相乘）
點值表示法
將$n$互不相同的$x$帶入多項式，會得到$n$個不同的取值$y$
則該多項式被這$n$個點$(x_1,y_1),(x_2,y_2),\dots,(x_n,y_n)$唯一確定
其中$y_i=\sum_{j=0}^{n-1} a_j*x_i^j$
例如：上面的例子用點值表示法可以為$(0,2),(1,6),(2,12)$
利用這種方法計算多項式乘法的時間復雜度仍然為$O(n^2)$
（選點$O(n)$，每次計算$O(n)$）
我們可以看到，兩種方法的時間復雜度都為$O(n^2)$，我們考慮對其進行優化
對於第一種方法，由於每個點的系數都是固定的，想要優化比較困難
對於第二種方法，貌似也沒有什么好的優化方法，不過當你看完下面的知識，或許就不這么想了
復數
在介紹復數之前，首先介紹一些可能會用到的東西
向量
同時具有大小和方向的量
在幾何中通常用帶有箭頭的線段表示
圓的弧度制
等於半徑長的圓弧所對的圓心角叫做1弧度的角，用符號rad表示，讀作弧度。用弧度作單位來度量角的制度叫做弧度制
公式:
$1^{\circ }=\dfrac{\pi}{180}rad$
$180^{\circ }=\pi rad$
平行四邊形定則
(好像畫的不是很標准。。)
平行四邊形定則：AB+AD=AC
復數
定義
設$a,b$為實數，$i^2=-1$，形如$a+bi$的數叫復數，其中$i$被稱為虛數單位，復數域是目前已知最大的域
在復平面中，$x$代表實數，$y$軸（除原點外的點）代表虛數，從原點$(0,0)$到$(a,b)$的向量表示復數$a+bi$
模長：從原點$(0,0)$到點$(a,b)$的距離，即$\sqrt{a^2+b^2}$
幅角：假設以逆時針為正方向，從$x$軸正半軸到已知向量的轉角的有向角叫做幅角
運算法則
加法：
因為在復平面中，復數可以被表示為向量，因此復數的加法與向量的加法相同，都滿足平行四邊形定則（就是上面那個）
乘法：
幾何定義：復數相乘，模長相乘，幅角相加
代數定義：$$(a+bi)*(c+di)$$
$$=ac+adi+bci+bdi^2$$
$$=ac+adi+bci-bd$$
$$=(ac-bd)+(bc+ad)i$$
單位根
下文中，默認$n$為$2$的正整數次冪
在復平面上，以原點為圓心，$1$為半徑作圓，所得的圓叫單位圓。以圓點為起點，圓的$n$等分點為終點，做$n$個向量，設幅角為正且最小的向量對應的復數為$\omega_n$，稱為$n$次單位根。
根據復數乘法的運算法則，其余$n-1$個復數為$\omega_n^2,\omega_n^3,\ldots,\omega_n^n$
注意$\omega_n^0=\omega_n^n=1$（對應復平面上以$x$軸為正方向的向量）
那么如何計算它們的值呢？這個問題可以由歐拉公式解決$$\omega_{n}^{k}=\cos\ k *\frac{2\pi}{n}+i\sin k*\frac{2\pi}{n}$$
例如
圖中向量$AB$表示的復數為$8$次單位根
單位根的幅角為周角的$\frac{1}{n}$
在代數中，若$z^n=1$，我們把$z$稱為$n$次單位根
單位根的性質
- $\omega _{n}^{k}=\cos k\dfrac{2\pi}{n}+i\sin k\dfrac {2\pi }{n}$（即上面的公式）
- $\omega _{2n}^{2k}=\omega _{n}^{k}$
證明：
$$\omega _{2n}^{2k}=\cos 2k*\frac{2\pi}{2n}+i\sin2k*\frac{2\pi}{2n}$$
$$=\omega _{n}^{k}$$
- $\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}^{0}=\omega _{n}^{n}=1$
講了這么多，貌似跟我們的正題沒啥關系啊。。
OK！各位坐穩了，前方高能！
快速傅里葉變換
我們前面提到過，一個$n$次多項式可以被$n$個點唯一確定。
那么我們可以把單位根的$0$到$n-1$次冪帶入，這樣也可以把這個多項式確定出來。但是這樣仍然是$O(n^2)$的呀！
我們設多項式$A(x)$的系數為$(a_o,a_1,a_2,\ldots,a_{n-1})$
那么$$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)$$
我們將$\omega_n^k (k<\frac{n}{2}) $代入得
$$A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})$$
$$=A_1(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_2(\omega_{\frac{n}{2}}^{k})$$
同理，將$\omega_n^{k+\frac{n}{2}}$代入得
$$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})$$
大家有沒有發現什么規律？
沒錯！這兩個式子只有一個常數項不同！
那么當我們在枚舉第一個式子的時候，我們可以$O(1)$的得到第二個式子的值
又因為第一個式子的$k$在取遍$[0,\frac{n}{2}-1]$時，$k+\frac{n}{2}$取遍了$[\frac{n}{2},n-1]$
所以我們將原來的問題縮小了一半！
而縮小后的問題仍然滿足原問題的性質，所以我們可以遞歸的去搞這件事情！
直到多項式僅剩一個常數項，這時候我們直接返回就好啦
時間復雜度：
不難看出FFT是類似於線段樹一樣的分治算法。
因此它的時間復雜度為$O(nlogn)$
快速傅里葉逆變換
不要以為FFT到這里就結束了。
我們上面的討論是基於點值表示法的。
但是在平常的學習和研究中很少用點值表示法來表示一個多項式。
所以我們要考慮如何把點值表示法轉換為系數表示法，這個過程叫做傅里葉逆變換
$(y_0,y_1,y_2,\dots,y_{n-1})$為$(a_0,a_1,a_2,\dots,a_{n-1})$的傅里葉變換（即點值表示）
設有另一個向量$(c_0,c_1,c_2,\dots,c_{n-1})$滿足
$$c_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i$$
即多項式$B(x)=y_0,y_1x,y_2x^2,\dots,y_{n-1}x^{n-1}$在$\omega_n^{0},\omega_n^{-1},\omega_n^{-2},\dots,\omega_{n-1}^{-(n-1)}$處的點值表示
emmmm又到推公式時間啦
$(c_0,c_1,c_2,\dots,c_{n-1})$滿足
$$c_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i$$
$$=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i$$
$$=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^j)^i)(\omega_n^{-k})^i$$
$$=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^j)^i(\omega_n^{-k})^i)$$
$$=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^j)^i(\omega_n^{-k})^i$$
$$=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^{j-k})^i$$
$$=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i)$$
設$S(x)=\sum_{i=0}^{n-1}x^i$
將$\omega_n^k$代入得
$$S(\omega_n^k)=1+(\omega_n^k)+(\omega_n^k)^2+\dots(\omega_n^k)^{n-1}$$
當$k!=0$時
等式兩邊同乘$\omega_n^k$得
$$\omega_n^kS(\omega_n^k)=\omega_n^k+(\omega_n^k)^2+(\omega_n^k)^3+\dots(\omega_n^k)^{n}$$
兩式相減得
$$\omega_n^kS(\omega_n^k)-S(\omega_n^k)=(\omega_n^k)^{n}-1$$
$$S(\omega_n^k)=\frac{(\omega_n^k)^{n}-1}{\omega_n^k-1}$$
$$S(\omega_n^k)=\frac{(\omega_n^n)^{k}-1}{\omega_n^k-1}$$
$$S(\omega_n^k)=\frac{1-1}{\omega_n^k-1}$$
觀察這個式子，不難看出它分母不為0，但是分子為0
因此，當$K!=0$時，$S(\omega^{k}_{n})=0$
那當$k=0$時呢？
很顯然，$S(\omega^{0}_{n})=n$
繼續考慮剛剛的式子
$$c_k=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i)$$
當$j \neq k$時，值為$0$
當$j=k$時，值為$n$
因此，
$$c_k=na_k$$
$$a_k=\frac{c_k}{n}$$
這樣我們就得到點值與系數之間的表示啦
理論總結
至此，FFT的基礎理論部分就結束了。
我們來小結一下FFT是怎么成功實現的
首先，人們在用系數表示法研究多項式的時候遇阻
於是開始考慮能否用點值表示法優化這個東西。
然后根據復數的兩條性質（這個思維跨度比較大）得到了一種分治算法。
最后又推了一波公式，找到了點值表示法與系數表示法之間轉換關系。
emmmm
其實FFT的實現思路大概就是
系數表示法—>點值表示法—>系數表示法
引用一下遠航之曲大佬的圖
當然，在實現的過程中還有很多技巧
我們根據代碼來理解一下
遞歸實現
遞歸實現的方法比較簡單。
就是按找我們上面說的過程，不斷把要求的序列分成兩部分，再進行合並
在c++的STL中提供了現成的complex類，但是我不建議大家用，畢竟手寫也就那么幾行，而且萬一某個毒瘤卡STL那豈不是很ＧＧ？
// luogu-judger-enable-o2 // luogu-judger-enable-o2 #include<iostream> #include<cstdio> #include<cmath> using namespace std; const int MAXN = 4 * 1e6 + 10; inline int read() { char c = getchar(); int x = 0, f = 1; while (c < '0' || c > '9') {if (c == '-')f = -1; c = getchar();} while (c >= '0' && c <= '9') {x = x * 10 + c - '0'; c = getchar();} return x * f; } const double Pi = acos(-1.0); struct complex { double x, y; complex (double xx = 0, double yy = 0) {x = xx, y = yy;} } a[MAXN], b[MAXN]; 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 fast_fast_tle(int limit, complex *a, int type) { if (limit == 1) return ; //只有一個常數項 complex a1[limit >> 1], a2[limit >> 1]; for (int i = 0; i <= limit; i += 2) //根據下標的奇偶性分類 a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1]; fast_fast_tle(limit >> 1, a1, type); fast_fast_tle(limit >> 1, a2, type); complex Wn = complex(cos(2.0 * Pi / limit) , type * sin(2.0 * Pi / limit)), w = complex(1, 0); //Wn為單位根，w表示冪 for (int i = 0; i < (limit >> 1); i++, w = w * Wn) //這里的w相當於公式中的k a[i] = a1[i] + w * a2[i], a[i + (limit >> 1)] = a1[i] - w * a2[i]; //利用單位根的性質，O(1)得到另一部分 } int main() { int N = read(), M = read(); for (int i = 0; i <= N; i++) a[i].x = read(); for (int i = 0; i <= M; i++) b[i].x = read(); int limit = 1; while (limit <= N + M) limit <<= 1; fast_fast_tle(limit, a, 1); fast_fast_tle(limit, b, 1); //后面的1表示要進行的變換是什么類型 //1表示從系數變為點值 //-1表示從點值變為系數 //至於為什么這樣是對的，可以參考一下c向量的推導過程， for (int i = 0; i <= limit; i++) a[i] = a[i] * b[i]; fast_fast_tle(limit, a, -1); for (int i = 0; i <= N + M; i++) printf("%d ", (int)(a[i].x / limit + 0.5)); //按照我們推倒的公式，這里還要除以n return 0; }
update:遞歸版我本地是可以AC的，只要開大數組就可以了，但是交到洛谷上會WA0
woc? 臉好疼。。。。。。
咳咳。
速度什么的才不是關鍵呢？
關鍵是我們AC不了啊啊啊
表着急，AC不了不代表咱們的算法不對，只能說這種實現方法太low了
下面介紹一種更高效的方法
迭代實現
再盜一下那位大佬的圖
觀察一下原序列和反轉后的序列？
聰明的你有沒有看出什么顯而易見的性質？
沒錯！
我們需要求的序列實際是原序列下標的二進制反轉！
因此我們對序列按照下標的奇偶性分類的過程其實是沒有必要的
這樣我們可以$O(n)$的利用某種操作得到我們要求的序列，然后不斷向上合並就好了
// luogu-judger-enable-o2 #include<iostream> #include<cstdio> #include<cmath> using namespace std; const int MAXN = 1e7 + 10; inline int read() { char c = getchar(); int x = 0, f = 1; while (c < '0' || c > '9') {if (c == '-')f = -1; c = getchar();} while (c >= '0' && c <= '9') {x = x * 10 + c - '0'; c = getchar();} return x * f; } const double Pi = acos(-1.0); struct complex { double x, y; complex (double xx = 0, double yy = 0) {x = xx, y = yy;} } a[MAXN], b[MAXN]; 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 N, M; int l, r[MAXN]; int limit = 1; void fast_fast_tle(complex *A, int type) { for (int i = 0; i < limit; i++) if (i < r[i]) swap(A[i], A[r[i]]); //求出要迭代的序列 for (int mid = 1; mid < limit; mid <<= 1) { //待合並區間的長度的一半 complex Wn( cos(Pi / mid) , type * sin(Pi / mid) ); //單位根 for (int R = mid << 1, j = 0; j < limit; j += R) { //R是區間的長度，j表示前已經到哪個位置了 complex w(1, 0); //冪 for (int k = 0; k < mid; k++, w = w * Wn) { //枚舉左半部分 complex x = A[j + k], y = w * A[j + mid + k]; //蝴蝶效應 A[j + k] = x + y; A[j + mid + k] = x - y; } } } } int main() { int N = read(), M = read(); for (int i = 0; i <= N; i++) a[i].x = read(); for (int i = 0; i <= M; i++) b[i].x = read(); while (limit <= N + M) limit <<= 1, l++; for (int i = 0; i < limit; i++) r[i] = ( r[i >> 1] >> 1 ) | ( (i & 1) << (l - 1) ) ; // 在原序列中 i 與 i/2 的關系是 ： i可以看做是i/2的二進制上的每一位左移一位得來 // 那么在反轉后的數組中就需要右移一位，同時特殊處理一下奇數 fast_fast_tle(a, 1); fast_fast_tle(b, 1); for (int i = 0; i <= limit; i++) a[i] = a[i] * b[i]; fast_fast_tle(a, -1); for (int i = 0; i <= N + M; i++) printf("%d ", (int)(a[i].x / limit + 0.5)); return 0; }