FFT求卷積（多項式乘法）
卷積
如果有兩個無限序列a和b，那么它們卷積的結果是：\(y_n=\sum_{i=-\infty}^\infty a_ib_{n-i}\)。如果a和b是有限序列，a最低的項為a0，最高的項為an，b同理，我們可以把a和b超出范圍的項都設置成0。那么可以得出：y0=a0b0，y1=a1b0+a0b1，y2=a0b2+a1b1+a2b0……，y(n+m)=a(n)b(m)。
構造兩個多項式A(x)和B(x)：
\(A=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}+a_nx^n\)，
\(B=b_0+b_1x+b_2x^2+...+b_{m-1}x^{m-1}+b_mx^m\)。
那么\(A(x)*B(x)=C(x)=a_0b_0+(a_0b_1+a_1b_0)x+...+a_nb_mx^{n+m}\)，把系數提取出來，可以發現兩序列卷積可以轉換為用序列作系數進行多項式乘法。
多項式
一個多項式既可以用系數表示，也可以用點值表示。n個點可以表示一個n-1次多項式。
如果用系數表示法來多項式乘法，時間復雜度是\(O(n^2)\)的，而用點值表示法只需要\(O(n)\)的時間。然而我們需要的是系數表示法。所以我們需要找到一個優秀的算法將它們兩者轉換，這就是（我們眼中的）FFT。
復數
設\(i^2=-1\)，a，b為實數，形如\(a+bi\)的數叫做復數。
用x軸表示a的大小，y軸表示b的大小，構造出的平面直角坐標系叫做復平面。復數的模長是原點到\((a, b)\)的距離，即\(\sqrt{a^2+b^2}\)。復數的輻角即為以逆時針為正方向，從x軸正半軸到已知向量的轉角。
復數的加減法則是顯然的，可以看作向量的加減。
復數可以寫成\(N(cos\alpha+isin\alpha )\)，\(\alpha\)表示復數的輻角。設\(z_1=A(cos\alpha + isin\alpha)\)，\(z_2=B(cos\beta + isin\beta)\)，那么\(z_1z_2=AB[(cos\alpha cos\beta-sin\alpha sin\beta)+i(sin\alpha cos\beta+cos\alpha sin\beta)]=AB[cos(\alpha+\beta)+isin(\alpha+\beta)]\)。也就是說，兩復數相乘，模長相乘，輻角相加。如果寫成普通形式的話，就是\((a+bi)(c+di)=(ac-bd)+(bc+ad)i\)。
單位根
在復平面上，以原點為圓心，1為半徑作圓，所得得圓為單位圓。從x軸正半軸開始將圓n等分，聯向第一個等分點所代表的復數\(\omega_n\)叫做n次單位根，意思是說\(w_n\)的n次方為1（根據復數的乘法運算法則）。可以推得，其他等分點代表的向量為\(\omega_n^1\),\(\omega_n^2\)……，一直到\(\omega_n^n = \omega_n^0=1\)。顯然\(\omega_n^k=cosk*\frac{2\pi}{n}+isink*\frac{2\pi}{n}\)。
單位根有幾個性質：
- 消去引理：\(\omega_{2n}^{k+n}=-\omega_{2n}^k\)。這是最重要的性質，使得分叉個數為2。
- 折半引理：\({(\omega_n^k)^2}={\omega_{n/2}^k}\)。這保證了FFT中子問題和原問題的規模都是n。
- 求和引理：\(\sum_{i=0}^{n-1}(\omega_n^k)^i=\left\{ \begin{aligned} 0, n\nmid k \\ n, n\mid k \end{aligned} \right.\)這是用來證明逆變換的。
DFT
前面說過，DFT是要把多項式的系數表達轉成點值表達。設多項式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)\)
根據單位根的性質，將前面一半的值帶入可得：
\(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})\)(折半引理的作用：將問題分解成條件完全相同的子問題)
同理帶入后面的值：
\(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^kA_2(\omega_n^{2k})\)(消去引理：使得分叉個數為2。如果沒有這個引理的話，就必須再去算一遍\(A(\omega_n^{k+\frac{n}{2}})\)的值，分叉個數變成4了。)
由於這兩個式子只有加號減號不同，我們只需計算前面一半的點值即可。這樣就將問題規模縮小了一半。當n=1時，點值是一個常數，直接返回即可。不難看出這是一個分治算法，時間復雜度為\(O(nlogn)\)。
為IDFT作准備
我們發現，FFT其實是在求下圖的\(y_i\)：（實在打不出來qwq）
那么現在的問題是，已知\(y_i\)，如何推回\(a_i\)？
由於我太弱了，繼續上圖吧。。
（補充一下，只要求出那個范德蒙德行列式的逆矩陣，乘在等式兩邊，那么就可以通過\(y_i\)推出\(a_i\)）
怎么構造\(V^{-1}\)，使得\(v_i^Tv_j^{-1}=\left\{ \begin{aligned} 1, i=j \\ 0, i\ne j\end{aligned} \right.\)呢？
是不是很神？這樣我們就構造出了\(V^{-1}\)：
現在，問題就變成了用\(\omega_n^{-1}\)為本原單位根，對y向量作FFT以后除以n。PPT里說的吼啊，稍微修改一下代碼就行了。
遞歸實現FFT
#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
const int maxn=2e6+5;
const double Pi=3.1415926535898;
int t, n, m, len=1;
struct Cpx{ //復數
double x, y;
Cpx (double t1=0, double t2=0){ x=t1, y=t2; }
}A[maxn*2], B[maxn*2], C[maxn*2];
Cpx operator +(Cpx a, Cpx b){ return Cpx(a.x+b.x, a.y+b.y); }
Cpx operator -(Cpx a, Cpx b){ return Cpx(a.x-b.x, a.y-b.y); }
Cpx operator *(Cpx a, Cpx b){ return Cpx(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x); }
void fdft(Cpx *a, int n, int flag){ //快速將當前多項式從系數表達轉換為點值表達
if (n==1) return; //如果只有1項系數為k，唯一的點值就是(w[1,1],k*w[1,1])=(1, k)
Cpx a1[(n>>1)+1], a2[(n>>1)+1];
for (int i=0; i<(n>>1); ++i) a1[i]=a[i<<1], a2[i]=a[i<<1|1];
fdft(a1, n>>1, flag); fdft(a2, n>>1, flag);
Cpx w1(cos(2*Pi/n), flag*sin(2*Pi/n)), w(1, 0); //idft用的負根
for (int i=0; i<(n>>1); ++i, w=w*w1){
a[i]=a1[i]+w*a2[i];
a[i+(n>>1)]=a1[i]-w*a2[i];
}
}
int main(){
scanf("%d%d", &n, &m); int x;
for (int i=0; i<=n; ++i) scanf("%lf", &A[i].x);
for (int i=0; i<=m; ++i) scanf("%lf", &B[i].x);
while (len<n+m) len<<=1; //idft需要至少l1+l2個點值
fdft(A, len, 1); fdft(B, len, 1);
for (int i=0; i<len; ++i) C[i]=A[i]*B[i];
fdft(C, len, -1); //idft
for (int i=0; i<=n+m; ++i){
x=C[i].x/len+0.5;
printf("%d ", x);
}
return 0;
}
題目是luogu的模板。注意給出的n和m都是多項式的最高次數，也就是說乘起來后的多項式最高次數為n+m，至少需要n+m個點。
迭代版FFT
遞歸版的太慢了，暗中觀察我們是如何處理序列的，可以發現：
把每個元素的編號二進制反轉一下，就是我們要求的序列編號！原因是原序列的最后1位決定了當前元素被分到前半區還是后半區，也就是轉換后元素編號的第1位。依次類推。
有一個O(n)推出n個數各自編號鏡像反轉的方法，大體思想是通過i<<1的反轉推出i的反轉。
由於各種原因，迭代版要比遞歸版快四倍左右～
#include <cmath>
#include <cctype>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
const int maxn=2e6+5;
const double pi=3.1415926535898;
int t, n, m, len=1, l, r[maxn*2];
struct Cpx{ //復數
double x, y;
Cpx (double t1=0, double t2=0){ x=t1, y=t2; }
}A[maxn*2], B[maxn*2], C[maxn*2];
Cpx operator +(Cpx a, Cpx b){ return Cpx(a.x+b.x, a.y+b.y); }
Cpx operator -(Cpx a, Cpx b){ return Cpx(a.x-b.x, a.y-b.y); }
Cpx operator *(Cpx a, Cpx b){ return Cpx(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x); }
void fdft(Cpx *a, int n, int flag){ //快速將當前多項式從系數表達轉換為點值表達
for (int i=0; i<n; ++i) if (i<r[i]) swap(a[i], a[r[i]]);
for (int mid=1; mid<n; mid<<=1){ //當前區間長度的一半
Cpx w1(cos(pi/mid), flag*sin(pi/mid)), x, y;
for (int j=0; j<n; j+=(mid<<1)){ //j:區間起始點
Cpx w(1, 0);
for (int k=0; k<mid; ++k, w=w*w1){ //系數轉點值
x=a[j+k], y=w*a[j+mid+k];
a[j+k]=x+y; a[j+mid+k]=x-y;
}
}
}
}
inline int getint(int &x){
char c; int flag=0;
for (c=getchar(); !isdigit(c); c=getchar())
if (c=='-') flag=1;
for (x=c-48; c=getchar(), isdigit(c);)
x=(x<<3)+(x<<1)+c-48;
return flag?x:-x;
}
int main(){
getint(n); getint(m); int x;
for (int i=0; i<=n; ++i) getint(x), A[i].x=x;
for (int i=0; i<=m; ++i) getint(x), B[i].x=x;
while (len<=n+m) len<<=1, ++l; //idft需要至少l1+l2個點值
for (int i=0; i<len; ++i) //編號的字節長度為l
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
fdft(A, len, 1); fdft(B, len, 1);
for (int i=0; i<len; ++i) C[i]=A[i]*B[i];
fdft(C, len, -1); //idft
for (int i=0; i<=n+m; ++i) printf("%d ", int(C[i].x/len+0.5));
return 0;
}
這樣可以做到1e6的數據最差也能跑進1s。我太菜了，並不會什么常數優化。
兩個月后的PS：注意n個點確定一個n-1次多項式。這是因為，對多項式求點值表達，相當於將一個范德蒙德矩陣乘上系數矩陣（前文有圖）。而范德蒙德矩陣是可逆的，所以在已知y的情況下，a也是唯一確定的。因此n個點一定唯一確定一個n-1次多項式。
五個月后的PS：qwq 借用了不少大佬的東西，侵刪。