快速傅立葉變換(FFT)
FFT 是之前學的,現在過了比較久的時間,終於打算在回顧的時候系統地整理一篇筆記,有寫錯的部分請指出來啊 qwq。
卷積
卷積、旋積或褶積(英語:Convolution)是通過兩個函式 \(f\) 和 \(g\) 生成第三個函式的一種數學運算元。
定義
設 \(f,g\) 在 \(R1\) 上可積,那麼 \(h(x) = \int_{-∞}^∞f(\tau)g(x-\tau)d\tau\) 稱為 \(f\) 與 \(g\) 的卷積。
對於整係數多項式域,\(n-1\) 次多項式 \(A,B\) 的卷積 \(h(x) = A(x)B(x) = \sum_{i=0}^{2(n-1)}\sum_{j=0}^{\min(i,n-1)} a_{j}b_{i-j}x^i\)。
係數表示法
即用多項式各項係數來刻畫這個多項式,例如 \(n-1\) 次多項式就可以寫成這樣:\(A(x) = a_0 + a_1x+...+a_{n-1}x^{n-1}\)
點值表示法
我們知道,\(n\) 個不同的點可以確定一個 \(n-1\) 次的多項式,所以我們可以使用 \(n\) 個(不同)點來刻畫一個 \(n-1\) 次多項式。
這樣做會有什麼方便呢?
例如 \(f(x) = (x_0, f(x_0)),...(x_n,f(x_n)),g(x) = (x_0, g(x_0)),...(x_n, g(x_n))\),那麼它們的卷積 \(h(x) = (x_0, f(x_0)g(x_0)),...(x_n, f(x_n)g(x_n))\)。
這意味著在係數表示法中需要 \(O(n^2)\) 次的乘法運算在點值表示法中只需要 \(O(n)\) 次。
係數表示法轉點值表示法(DFT)
下面考慮如何將 \(n-1\) 次多項式從係數表示法轉為點值表示法。
因為用普通的方法選取 \(n\) 個點然後將係數表示法轉為點值表示法的複雜度為 \(O(n^2)\)(因為需要選 \(n\) 個點,然後對於每個點 \(x\) 需要計算共 \(n\) 項的結果),我們考慮如何優化這一步。
注意到滿足 \(w^n=1\) 的單位根 \(w\) 有 \(n\) 個,故從這裡入手。
我們記方程 \(w^n = 1\) 的第 \(k\) 個單位根為 \(w_n^k\)。
方便起見,設 \(n\) 為 \(2\) 的冪(就算不是也可以看作是高次項的係數為 \(0\))。
將 \(A(x)\) 按照次數的奇偶性分別分成兩組 \(F(x),G(x)\),並表示為 \(A(x) = F(x^2) + xG(x^2)\)
例如 \(A(x) = a_0+a_1x+a_2x^2+a_3x^3\),那麼 \(F(x) = a_0+a_2x,G(x)=a_1+a_3x\)。
將 \(x=w_n^k\) 代入 \(A(x)\),由複數的性質,
\(A(w_n^k) = F(w_{\frac{n}{2}}^k) + w_n^k G(w_{\frac{n}{2}}^k)\) ,類似地 \(A(w_{n}^{k+\frac{n}{2}}) = F(w_{\frac{n}{2}}^k) - w_n^k G(w_{\frac{n}{2}}^k)\)。
推導:
\(A(w_{n}^{k+\frac{n}{2}}) = F(w_n^{2k+n}) + w_n^{k+\frac{n}{2}}G(w_n^{2k+n}) \\= F(w_{\frac{n}{2}}^k) + w_n^{k+\frac{n}{2}} G(w_{\frac{n}{2}}^k) \\= F(w_{\frac{n}{2}}^k) - w_n^k G(w_{\frac{n}{2}}^k)\)
可以發現對於兩個相應的單位根 \(w_n^k,w_n^{\frac{n}{2}+k}\),可以用對應的 \(F,G\) 算出(可以遞迴地實現這個過程),而且計算的範圍折半,所以一共需要計算 \(O(logN)\) 層,每一層執行 \(O(n)\) 次運算,所以複雜度為 \(O(NlogN)\)。
點值表示法轉系數表示法(IDFT)
下面考慮如何將 \(n-1\) 次多項式從點值表示法轉為係數表示法。
因為對於每個點值 \(y_i = \sum_{k=0}^{n-1}w_n^{ki}\),其中 \(i\in[0, n-1]\),我們可以寫出等式:
現在我們已經有向量 \(y\) 了(就是右式),因此,如果要得到向量 \(a\),只需要兩邊乘上 \(w\) 矩陣的逆即可。
這裡的 \(w\) 矩陣正是著名的範德蒙矩陣,它的逆正好是每一項都取倒數,然後除以 \(n\)。
因此有 \(a_i = \frac{1}{n}\sum_{k=0}^{n-1}w_n^{-ki}\),其中 \(i\in[0, n-1]\)。
有沒有發現 \(a_i,y_i\) 的形式非常接近?據此,我們可以在實現的時候在同一個函式中寫出逆變換和正變換,然後在得到的結果 res
中除以 \(n\) 就可以了。(參照下面的程式碼)
至此,FFT 的基本原理講述完畢,下面是優化。
位逆序置換
按照上文的講述,如果不看下面的程式碼,那麼編寫出來的是遞迴版本,但是這個版本的常數太大了,因此執行起來的效果不好,故使用位逆序置換來降低常數。
我們看看遞迴過程是什麼樣的,以 \(n=8\) 為例:
這裡就有一個非常神奇的規律:在最後一行中,原下標所對應的二進位制數翻轉正好是在最後一行的序數。例如 \(x_6\) 的下標是 \(6 = 110_{(2)}\),那麼它的序數正好是 \(011_{(2)} = 3\)。
據此,可以處理出 rev
陣列,它記錄的正是最後一行所有元素對應的下標。
簡單地說,遞迴形式是自上而下地做 FFT,而利用位逆序置換我們可以自下而上地做 FFT,它們在實際執行中有著常數上的區別。
模板題及程式碼
https://www.luogu.com.cn/problem/P3803
給定一個 \(n\) 次多項式 \(F(x)\),和一個 \(m\) 次多項式 \(G(x)\)。
請求出 \(F(x)\) 和 \(G(x)\) 的卷積。
#include<bits/stdc++.h>
using namespace std;
const int N=3e5+5;
const double pi=acos(-1);
int n, m;
// 複數類
struct Complex{
double x, y;
Complex operator + (const Complex &o)const { return {x+o.x, y+o.y}; }
Complex operator - (const Complex &o)const { return {x-o.x, y-o.y}; }
Complex operator * (const Complex &o)const { return {x*o.x-y*o.y, x*o.y+y*o.x}; }
};
Complex a[N], b[N];
int res[N];
int rev[N], bit, tot;
void fft(Complex a[], int inv){ // inv 指示正變換、逆變換。
for(int i=0; i<tot; i++) if(i<rev[i]) swap(a[i], a[rev[i]]);
for(int mid=1; mid<tot; mid<<=1){
auto w1=Complex({cos(pi/mid), inv*sin(pi/mid)});
for(int i=0; i<tot; i+=mid*2){
auto wk=Complex({1, 0});
for(int j=0; j<mid; j++, wk=wk*w1){
auto x=a[i+j], y=wk*a[i+j+mid];
a[i+j]=x+y, a[i+j+mid]=x-y;
}
}
}
}
int main(){
cin>>n>>m;
for(int i=0; i<=n; i++) cin>>a[i].x;
for(int i=0; i<=m; i++) cin>>b[i].x;
while((1<<bit)<n+m+1) bit++; // 結果次數分佈在 [0, n+m] 內,一共有 n+m+1 位。
tot=1<<bit; // 得到上文所說的 n
for(int i=0; i<tot; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
fft(a, 1), fft(b, 1); // 正變換 DFT
for(int i=0; i<tot; i++) a[i]=a[i]*b[i];
fft(a, -1); // 逆變換 IDFT
for(int i=0; i<=n+m; i++) res[i]=(int)(a[i].x/tot+0.5), printf("%d ", res[i]);
return 0;
}