快速傅立葉變換 學習筆記
簡介：
卷積是形如 \(C_k=\sum_{i\oplus j==k} A_i*B_j\) 的式子，其中 \(\oplus\) 為表示某種運算。
而快速傅立葉變換（FFT）用於在 \(O(n\log n)\) 的時間複雜度內求加法卷積。
對於一個 \(k\) 次多項式，如果我們知道了它在 \(k+1\) 個點上的值，那麼我們可以求出這個多項式。而兩個多項式在同一個點上的值相乘的話，就可以得到 這兩個多項式相乘後的多項式 在這個點上的值。FFT 的主要思路就是先將兩個多項式在若干個點上的取值，然後相乘，再透過點值求出兩個多項式相乘的結果。
單位根
先考慮求點值。我們肯定不能暴力選擇點去代入求值，這樣的時間複雜度是 \(O(n^2)\) 的，要找一些有特殊性質的點去代入。
於是我們在複數域找到了單位根，它是在複平面上的單位圓上的點，表示為 \(\omega_n^k\) ，它的意義是將單位圓分割成 \(n\) 份，按順指標方向的第 \(k\) 個點所表示的值。兩個複數相乘時，它們的輻角相加（從實軸的正半軸開始，按逆時針旋轉），模長相乘（該點到原點的距離）。於是我們可以得到單位根的以下性質：
對於最後一條，我們可以將 \(n/2\) 看做旋轉半個圓周，所以就是取相反數。
對於一個 \(n\) 次多項式，我們看看將單位根代入會有什麼性質：
先將多項式 \(F[x]\) 拆成兩個部分： \(FL[x]=x^0+x^2+x^4...,FR[x]=x^1+x^3+x^5...\) 。然後就可以將多項式 \(F[x]\) 這樣表示：
再分別將 \(\omega_n^k\) 和 \(\omega_n^{k+n/2}\) 代入：
發現上下兩式只有後面一項的符號不同，所以如果我們求出了 \(FL[x]\) 和 \(FR[x]\) 在 \(\omega_{n/2}^{0..n/2-1}\) 上的值，我們就可以求出 \(F[x]\) 在 \(\omega_n^{0...n-1}\) 上的值了，這樣做一次的時間複雜度是 \(O(n)\) 的。同時我們可以發現，對於 \(FL[x]\) 和 \(FR[x]\) ，其實可以看做一個相同的子任務，再次遞迴下去求解。一共需要遞迴 \(log\ n\) 次。
但是此時我們又有一個問題：如果某個 \(n\) 不能被整除了怎麼辦？我們可以給它在後面補位，係數設為 0 就行。再觀察一下我們遞迴到最下面一層時，係數的編號與第一層相比有什麼特點。可以發現它是二進位制位的反轉，所以我們可以用數位DP的方法求出反轉後的編號，這樣就不用一遍一遍遞迴下去求了。
tr[i]=(tr[i>>1]>>1)|(i&1)?limit>>1:0
求多項式
透過上面的步驟，我們已經將 \(n\) 個點的值求出來了，接下來就要透過這 \(n\) 個點的值求多項式。
設點值構成的序列為 \(G[x]\) 。然後我們可以列出以下式子：
然後可以得到以下結論：
這個式子的證明就是將上面的式子代入。此時我們發現，這個式子和我們求點值的式子十分相似，只是在單位根上有點區別。我們可以發現，在單位圓上， \(\omega_n^{-1}\) 所在的點和 \(\omega_n^1\) 所在的點，它們關於實軸對稱。於是我們也可以用同樣的方法來求上面的式子。最後記得除以 \(n\) 就是我們最終的多項式。
Code
#include<cstdio>
#include<cmath>
using namespace std;
const int N=3e6+5;const double PI=6.28318530717958646;
struct Date{
double x,y;
Date operator*(Date v) {
Date res;res.x=x*v.x-y*v.y;res.y=x*v.y+y*v.x;return res;
};
Date operator+(Date v) {
Date res;res.x=x+v.x;res.y=y+v.y;return res;
};
Date operator-(Date v) {
Date res;res.x=x-v.x;res.y=y-v.y;return res;
};
};Date f[N],g[N],w[N];
int tr[N],n,m;
void swap(double &x,double &y)
{
double t=x;x=y;y=t;
}
void FFT(int flag)
{
int i,j,len;Date t;
w[0].x=1;w[0].y=0;
w[1].x=cos(PI/n);w[1].y=sin(PI/n)*flag;
for(i=1;i<n;i++)
w[i]=w[i-1]*w[1];
for(i=0;i<n;i++)
{
if(i<tr[i])
{
swap(f[i].x,f[tr[i]].x);swap(f[i].y,f[tr[i]].y);
}
}
for(len=1;len<n;len=len<<1)
{
for(i=0;i<n;i+=2*len)
{
for(j=0;j<len;j++)
{
t=f[i+j+len]*w[n/len/2*j];
f[i+j+len]=f[i+j]-t;f[i+j]=f[i+j]+t;
}
}
}
}
int main()
{
int i,h,len;
scanf("%d%d",&n,&m);
for(i=0;i<=n;i++)
scanf("%lf",&f[i].x);
for(i=0;i<=m;i++)
scanf("%lf",&g[i].x);
h=0;len=n+m;
while((1<<h)<n+m+1)
h++;
n=1<<h;
for(i=0;i<n;i++)
tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
FFT(1);
for(i=0;i<n;i++)
{
swap(f[i].x,g[i].x);swap(f[i].y,g[i].y);
}
FFT(1);
for(i=0;i<n;i++)
f[i]=f[i]*g[i];
FFT(-1);
for(i=0;i<=len;i++)
printf("%d ",int(f[i].x/n+0.49));
return 0;
}