==== €€£ WARNING ====
這篇博文內容相對偏少, 已經在后續博文中擴充.
大家可以看我的最新博文 [學習筆記&教程] 信號, 集合, 多項式, 以及各種卷積性變換 (FFT,NTT,FWT,FMT)
==== ====
引入
可能有不少OIer都知道FFT這個神奇的算法, 通過一系列玄學的變化就可以在 $O(nlog(n))$ 的總時間復雜度內計算出兩個向量的卷積(或者多項式乘法/高精度乘法), 而代碼量卻非常小. 博主一年半前曾經因COGS的一道叫做"神秘的常數 $\pi$"的題目而去學習過FFT, 但是基本就是照着板子打打完並不知道自己在寫些什么鬼畜的東西OwO
不過...博主這幾天突然照着算法導論自己看了一遍發現自己似乎突然意識到了什么OwO然后就打了一道板子題還1A了OwO再加上午考試差點AK以及日更頻率即將不保於是就有了這篇博文233
博主在寫這篇文字的過程中也發現了不少自己之前忽略的FFT細節, 感覺對FFT似乎有了更深刻的理解, 希望想學習FFT的讀者能認真看完這篇文章並仔細體味FFT的優雅性質.
本文可能並不會介紹關於使用FFT進行信號處理的相關信息, 主要介紹在OI中的應用即快速卷積.
Rush了好久終於寫出來了QAQ
那么現在, 就讓我們變玄學為科學了解 $FFT$ 背后的工作原理與它的優雅性質
前置知識: 多項式的表示
多項式的系數表達
相信大家都知道對於一個次數界為 $n$ 的多項式 $A(x)$ 可以表示為如下形式:
\[A(x) = \sum _{i=0} ^{n-1} a_i x^i\]
這時, 多項式 $A(x)$ 的系數所組成的向量 $\vec{a}=(a_0,a_1,...a_{n-1})$ 是這個多項式的系數表達. (實際上是列向量不是行向量OwO)
使用系數表達來表示多項式可以進行一些方便的運算, 比如對其進行求值, 這時我們可以采用秦九昭算法來在 $O(n)$ 時間復雜度內對多項式進行求值.
\[A(x_0) = a_0+x_0(a_1+x_0(a_2+...+x_0(a_{n-2}+x_0 a_{n-1})...))\]
但是當我們想把兩個采用系數表達的多項式乘起來的時候, 恭喜你, 問題來了: 一個多項式中的每個系數都要與另一個多項式中的每個系數相乘. 然后時間復雜度就變成了鬼畜的 $\Theta(n^2)$ ...
也就是對於兩個多項式的系數表示 $\vec{a}$ 和 $\vec{b}$ 來說, 兩個多項式乘積的系數表達 $\vec{c}$ 中的值的求值公式如下:
\[c_i=\sum_{j=0}^i a_jb_{i-j}\]
然而多項式乘法和卷積在很多地方都需要快速的實現, 而對於多項式乘法來說, 另一種表示方法似乎更為合適:
多項式的點值表達
首先我們對點值表達進行定義:
一個次數界為 $n$ 的多項式 $A(x)$ 的點值表達為一個由 $n$ 個點值對所組成的集合:
\[\{(x_0,y_0),(x_1,y_1),(x_2,y_2),...,(x_{n-1},y_{n-1})\}\]
使得對於 $k=0,1,2,...,n-1$ , 所有的 $x_k$ 各不相同
(我們可以先暫且把 $A(x)$ 看作一個函數)
其中 $y_k=A(x_k)$ , 也就是說選取 $n$ 個不同的值分別代入求值之后產生 $n$ 個值, 就會得到 $n$ 個未知數的值與多項式值的點對. 這 $n$ 個點對就是多項式的一個點值表達.
不難看出點值表達並不像多項式表達一樣具有唯一性, 因為 $x_k$ 是沒有限制條件的.
對於一個系數表達的多項式, 顯然求出它的一個點值表達非常方便, 只需取 $n$ 個不同的 $x$ 值代入並求出對應的點值即可. 計算出這些點值需要 $O(n^2)$ 的時間.
....
等等...? 怎么還是 $O(n^2)$ ? 別急, 很快我們就會發現, 如果我們像某些喪病出題人一樣精心選擇數據, 我們就可以在 $O(nlog(n))$ 的時間復雜度內完成這個運算.
而從多項式的點值表達轉換為唯一的系數表達就沒有那么顯然了. 首先我們介紹兩個概念:
求值與插值
從一個多項式的系數表達確定其點值表達的過程稱為求值(似乎很顯然的樣子? 畢竟我們求點值表達的過程就是取了 $n$ 個 $x$ 的值然后扔進了多項式求了 $n$ 個值出來233), 而求值運算的逆運算(也就是從一個多項式的點值表達確定多項式的系數表達)被稱為插值. 而插值多項式的唯一性定理告訴我們只有多項式的次數界等於已知的點值對的個數, 插值過程才是明確的(能得到一個確定的多項式表達) , 聯想之前的矩陣與線性方程組和增廣矩陣可以得到如下證明:
定理(插值多項式的唯一性): 對於任意 $n$ 個點值對組成的集合 $\{(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1}\}$ ,且其中所有 $x_k$ 不同, 則存在唯一的次數界為 $n$ 的多項式 $A(x)$ , 滿足 $y_k=A(x_k),k=0,1,...,n-1$ .
證明 證明主要是根據某個矩陣存在逆矩陣. 多項式系數表達等價於下面的矩陣方程:
\[
\begin{bmatrix}
1 & x_0 & x_0^2 & \dots & x_0^{n-1} \\
1 & x_1 & x_1^2 & \dots & x_1^{n-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & x_{n-1} & x_{n-1}^2 & \dots & x_{n-1}^{n-1}
\end{bmatrix}
\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix}
=
\begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_{n-1} \end{bmatrix}
\]左邊的矩陣表示為 $V(x_0,x_1,...,x_{n-1})$ , 稱為范德蒙矩陣. 而這個矩陣的行列式為:
\[\prod _{0\leq j < k \leq n-1} (x_k-x_j)\]
所以根據定理 "$n\times n$ 矩陣是奇異的, 當且僅當該矩陣的行列式為 $0$", 如果 $x_k$ 皆不相同, 則這個矩陣是可逆的. 因此給定點值表達 $\vec y$ , 則能確定唯一的系數表達 $\vec{a}$ 使:
\[\vec{a}=V(x_0,x_1,...,x_{n-1})^{-1}\vec{y}\]
$\blacksquare$
上面的證明過程實際上也給出了插值的一種算法, 即求出范德蒙矩陣的逆矩陣. 我們可以在 $O(n^3)$ 時間復雜度內求出逆矩陣所以相應地可以在 $O(n^3)$ 時間復雜度內計算出點值表達所對應的系數表達.
然而這特么比求值還慢...=_=
所幸我們還有一種基於 $n$ 個點的插值算法, 基於拉格朗日公式:
\[A(x)=\sum_{k=0}^{n-1} y_k \frac{\prod _{j\neq k} (x-x_j)} {\prod_{j\neq k} (x_k-x_j)}\]
基於這個公式我們可以在 $O(n^2)$ 時間復雜度內計算出點值表達所對應的系數表達.
然后插值勉強追上了求值的時間復雜度.
那么問題來了, 計算這個點值表達有啥用?
點值表達的運算
點值表達在很多多項式相關操作上比系數表達要方便很多, 比如對於加法, 如果 $C(x)=A(x)+B(x)$ ,則對於任意點值 $x_k$ , 滿足 $C(x_k)=A(x_k)+b(x_k)$ . 而正確性是顯而易見的.所以對於兩個點值表達的多項式進行加法操作, 時間復雜度為 $O(n)$.
與加法類似, 多項式乘法在點值表達的情況下也變得更加方便: 如果 $C(x)=A(x)B(x)$, 則對於任意點 $x_k$ , $C(x_k)=A(x_k)B(x_k)$ . 也就是說點值表達常常可以大大簡化多項式間的運算!
然而問題又來了: 多項式 $C$ 的次數界顯然是 $A$ 的次數界與 $B$ 的次數界的和. 所以如果我們要插值並獲得唯一的系數表達式就需要更多的點值對. 這時我們就有必要對 $A$ 與 $B$ 的點值表達進行"擴展", 至少擴展到 $C$ 的次數界.
但是我們似乎發現了一個嚴重的問題: 雖然點值表達中計算乘法速度比系數表達快很多, 但是到目前為止我們在兩種表達方式之間轉換的算法都是 $O(n^2)$ 的, 也就是說轉換之后時間復雜度相較於朴素的多項式乘法實現不但沒有提升, 甚至還增加了不少運算量.
對於這個問題, 我們可以采取一些方法: 因為我們在計算點值表達時可以采用任何點作為求值點, 通過精心挑選求值點就可以做到 $O(nlog(n))$ 的時間復雜度進行兩種表示法之間的轉換! 也就是說我們可以采取如下策略來加速卷積:
1.加倍次數界: 在多項式 $A$ 和多項式 $B$ 中加入若干值為 $0$ 的高階系數使其達到多項式 $C$ 所需的次數界, 這一過程為 $O(n)$ 時間復雜度.
2.求值: 在 $O(nlog(n))$ 時間復雜度內計算出多項式 $A$ 和多項式 $B$ 的點值表達.
3.逐點相乘: 在 $O(n)$ 時間復雜度內計算出多項式 $C$ 的點值表達.
4.插值: 在$O(nlog(n))$ 時間復雜度內計算出多項式 $C$ 的系數表達
也就是說我們可以在 $O(nlog(n))$ 的總時間復雜度內計算出多項式 $A$ 與多項式 $B$ 的乘積.
而這個關鍵的求值與插值算法又是什么呢? 沒錯, 就是快速傅里葉變換(Fast Fourier Transform, FFT)
DFT與FFT
我們剛剛說到, 當我們仔細選擇求值點的話就可以加速求值過程, 而我們所選擇的值就是一種具有特殊性質的數:n次單位復數根.
n次單位復數根
n次單位復數根是滿足 $\omega ^n=1$的復數 $\omega$ . 根據復數乘法運算的幾何意義(幅角相加, 模相乘)和同余, 我們可以得出 $n$ 次單位復數根恰好有 $n$ 個且均勻分布在以原點為圓心, 一單位長度為半徑的圓周上. 如下兩圖所示:
所以 . 根據復數在指數上的定義 $e^{iu}=cos(u) + i\:sin(u)$ , 我們可以得到對於 $k=0,1,...,n-1$ , $n$ 次單位根是 $e^{2\pi i k/n}$
其中, $\omega _n=e^{2\pi i /n}$ 為主 $n$ 次單位根. 因為所有其它單位根實際上都是它的整次冪.
實際上 $n$ 次單位根形成了一個乘法群(感興趣請自行查閱群論相關資料, 在這里寫的話估計又要扯很久).
我們注意到對於一個 $n$ 次單位根 $\omega_n^k$, 它在極坐標系中的幅角就是 $\frac{k}{n}$ 個周角. 因為它們均勻分布在單位圓上. 而對於分數, 我們是可以進行約分的, 這也就引出了對后續FFT十分重要的結論:
消去引理, 折半引理與求和引理
直接引用算法導論中給出的證明:
引理(消去引理): 對任意整數 $n \geq 0, k\geq 0 $, 以及 $d\geq 0$, 有:
\[\omega_{dn}^{dk}=\omega_n^k\]
證明: 由單位復數根的定義式可直接推出引理, 因為
\[\omega_{dn}^{dk}=(e^{2\pi i / dn})^{dk}=(e^{2\pi i /n})^k= \omega_n^k\]
$\blacksquare$
十分優雅而精簡的證明, 正是運用了單位復數根的特殊性質. 證明過程中使用了冪的冪中指數相乘的運算法則, 將 $d$ 乘進括號內, 與括號內在指數分母上的另一個 $d$ 相抵消.
根據消去引理, 我們還可以推出另一個引理:
引理(折半引理): 如果 $n>0$ 為偶數, 那么 $n$ 個 $n$ 次單位復數根的平方的集合就是 $n/2$ 個$n/2$ 次單位復數根的集合.
證明: 根據消去引理, 對任意非負整數 $k$ , 我們有 $(\omega_n^k)^2=\omega_{n/2}^k$ . 注意, 如果對所有 $n$ 次單位復數根進行平方, 那么獲得每個 $n/2$ 次單位根正好兩次. 因為:
\[(\omega_n^{k+n/2})^2=\omega_n^{2k+n}=\omega_n^{2k}\omega_n^n=\omega_n^{2k}=(\omega_n^k)^2\]
因此, \(\omega_n^{k+n/2}\) 與 \(\omega_n^k\) 平方相同.
$\blacksquare$
從消去引理推出的折半引理是后面FFT中使用的分治策略的正確性的關鍵, 因為折半引理保證了遞歸子問題的規模只是遞歸調用前的一半, 從而保證 $O(nlog(n))$ 時間復雜度.
另一個重要的引理是求和引理, 是后文中使用FFT快速插值的基礎.
引理(求和引理): 對於任意整數 $n\geq 1$ 與不能被 $n$ 整除的非負整數 $k$ , 有:
\[\sum_{j=0}^{n-1}(\omega_n^k)^j=0\]
證明: 根據幾何級數的求和公式:
\[\sum_{k=0}^nx^k=\frac{x^{n+1}-1}{x-1}\]
我們可以得到如下推導:
\[\sum_{j=0}^{n-1}(\omega_n^k)^j=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}=\frac{(\omega_n^n)^k-1}{\omega_n^k-1}=\frac{(1)^k-1}{\omega_n^k-1}=0\]
定理中要求 $k$ 不可被 $n$ 整除, 而僅當 $k$ 被 $n$ 整除時有 $\omega_n^k=1$, 所以保證分母不為 $0$ .
$\blacksquare$
介紹了這么多前置知識與理論基礎, 我們該開始今天的重頭戲了:
離散傅里葉變換 DFT
首先我們還是給出一個DFT的定義:
還記得我們需要FFT去做什么嗎? 我們想要快速從多項式的系數表達計算出它的點值表達, 而上一節我們說過, 我們要用 $n$ 個 $n$ 次單位復數根處進行這一求值過程.
對於 $A$ 的系數表達 $\vec{a}= (a_0,a_1,...,a_{n-1})$ 和 $k=0,1,...,n-1$ , 定義結果 $y_k$:
\[y_k=A(\omega_n^k)=\sum_{i=0}^{n-1}a_i\omega^{ki}_n\]
則 $\vec{y}=(y_0,y_1,...,y_{n-1})$ 就是系數向量 $\vec{a}$ 的離散傅里葉變換(DFT). 簡記為 $\vec{y}=DFT_n(\vec{a})$
然而問題一個接一個: 從這個定義式來看, 我們計算出一個多項式 $A(x)$ 的離散傅里葉變換依然需要 $O(n^2)$ 的時間復雜度, 依然沒有任何提升的樣子?
好了, 真正的主角該登場了:
快速傅里葉變換 FFT
如果我們使用FFT利用單位復數根的特殊性質,我們可以在 $O(nlog(n))$ 時間復雜度內計算出 $DFT_n(\vec a)$ . 如果直接按定義式計算的話依然是和求值相同的時間復雜度.
在下文中為了分治方便, 我們將假設 $n$ 是 $2$ 的整次冪, 雖然對於非整次冪的 $n$ 的算法已經出現, 但是平常 OI 中使用的時候一般采取補一大坨值為 $0$ 的高次項系數強行補到 $2$ 的整次冪23333
FFT使用的是分治策略. 根據系數下標的奇偶來拆分子問題, 將這些系數分配到兩個次數界為 $n/2$ 的多項式 $A^{[0]}(x)$ 和 $A^{[1]}(x)$ :
\[A^{[0]}(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{n/2-1}\]
\[A^{[1]}(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{n/2-1}\]
可以發現 $A^{[0]}(x)$ 中包含了 $A$ 中下標為偶數的系數, 而 $A^{[1]}(x)$ 則包含了 $A$ 中下標為奇數的系數. 所以實際上我們根據下標可以得到如下結論:
\[A(x)=A^{[0]}(x^2)+A^{[1]}(x^2)x\]
因為次數界折半了所以代入的是 $x^2$ , 然后對於奇數下標系數還要乘上一個 $x$ 作為補償.
於是我們就把求 $A(x)$ 在 $\omega_n^0,\omega_n^1,...,\omega_n^{n-1}$ 處的值轉化成了這樣:
1.求多項式 $A^{[0]}(x)$ 和 多項式 $A^{[1]}(x)$ 在下列點上的取值:
\[(\omega_n^0)^2,(\omega_n^1)^2,...(\omega_n^{n-1})^2\]
2.根據剛剛推導出的將兩個按奇偶拆開的多項式的值合並的公式合並結果.
但是根據折半引理, 這 $n$ 個要代入子多項式中求值的點顯然並不是 $n$ 個不同值, 而是 $n/2$ 個 $n/2$ 次單位復數根組成. 所以我們遞歸來對次數界為 $n/2$ 的多項式 $A^{[0]}(x)$ $A^{[1]}(x)$ 在 $n/2$ 個 $n/2$ 次單位復數根處求值. 然我們成功地縮小了了子問題的規模. 也就是說我們將 $DFT_n$ 的計算划分為兩個 $DFT_{n/2}$ 來計算. 這樣我們就可以計算出 $n$ 個元素組成的向量 $\vec a$ 的DFT了.
然后讓我們結合FFT的C++實現來理解一下:
1 void FFT(Complex* a,int len){ 2 if(len==1) 3 return; 4 Complex* a0=new Complex[len/2]; 5 Complex* a1=new Complex[len/2]; 6 for(int i=0;i<len;i+=2){ 7 a0[i/2]=a[i]; 8 a1[i/2]=a[i+1]; 9 } 10 FFT(a0,len/2); 11 FFT(a1,len/2); 12 Complex wn(cos(2*Pi/len),sin(2*Pi/len)); 13 Complex w(1,0); 14 for(int i=0;i<(len/2);i++){ 15 a[i]=a0[i]+w*a1[i]; 16 a[i+len/2]=a0[i]-w*a1[i]; 17 w=w*wn; 18 } 19 delete[] a0; 20 delete[] a1; 21 }
讓我們逐行解讀一下這個板子:
第2,3行是遞歸邊界, 顯然一個元素的 $DFT$ 就是它本身, 因為 $\omega_1^0=1$, 所以
\[y_0=a_0\omega_1^0=a_0\]
第4~9行將多項式的系數向量按奇偶拆成兩部分.
第12,13和第17行結合起來用來更新 $\omega$ 的值, 12行計算主 $n$ 次單位復數根, 13行將變量 $w$ 初始化為 $\omega_n^0$ 也就是 $1$ , 17行使在組合兩個子多項式的答案時可以避免重新計算 $\omega_n$ 的冪, 而是采用遞推方式最大程度利用計算中間結果(OI中的常用技巧, 對於主 $n$ 次單位復數根還可以打個表來節約計算三角函數時的巨大耗時).
第10,11行將拆出的兩個多項式遞歸處理. 我們設多項式 $A^{[0]}(x)$ 的 $DFT$ 結果為 $\vec{y^{[0]}}$ , $A^{[1]}(x)$ 的 $DFT$ 結果為 $\vec{y^{[1]}}$, 則根據定義, 對於 $k=0,1,...,n/2-1$, 有:
\[y_k^{[0]}=A^{[0]}(\omega_{n/2}^k)\]
\[y_k^{[1]}=A^{[1]}(\omega_{n/2}^k)\]
然后根據消去引理, 我們可以得到 $\omega_{n/2}^k=\omega_n^{2k}$, 所以
\[y_k^{[0]}=A^{[0]}(\omega_n^{2k})\]
\[y_k^{[1]}=A^{[1]}(\omega_n^{2k})\]
設整個多項式的 $DFT$ 結果為 $\vec y$, 則第15,16行用如下推導來將遞歸計算兩個多項式的計算結果綜合為一個多項式的結果:
第15行, 對於 $y_0,y_1,...,y_{n/2-1}$, 我們設 $k=0,1,...,n/2-1$ :
\begin{aligned} y_k&=y^{[0]}_k+\omega_n^ky_k^{[1]} \\ &=A^{[0]}(\omega_n^{2k})+\omega_n^kA^{[1]}(\omega_n^{2k}) \\ &= A(\omega_n^k) \end{aligned}
上面的推導中, 第一行是我們在代碼中進行的運算, 第二行將 $y_k^{[0]}$ 和 $y_k^{[1]}$ 按定義代入, 第三行基於我們剛剛推導出的組合兩個子多項式 $DFT$ 結果的式子.
第16行, 對於 $y_{n/2}, y_{n/2+1}, ..., y_{n-1}$ , 我們同樣設 $k=0,1,...,n/2-1$ :
\begin{aligned} y_{k+(n/2)} &= y^{[0]}_{k+(n/2)} - \omega_n^ky_{k+(n/2)}^{[1]} \\ &= y_k^{[0]} + \omega_n^{k+(n/2)} y_k^{[1]} \\ &= A^{[0]}(\omega_n^{2k}) +\omega_n^{k+(n/2)} A^{[1]}(\omega_n^{2k}) \\ &= A^{[0]}(\omega_n^{2k+n}) + \omega_n^{k+(n/2)} A^{[1]}(\omega_n^{2k+n}) \\ &= A(\omega_n^{k+(n/2)}) \end{aligned}
上面的推導中, 第一行同樣是我們在代碼中進行的運算, 第二行從 $\omega_n^{k+(n/2)}=-\omega_n^k$ 推導出, 第三行將 $y_k^{[0]}$ 和 $y_k^{[1]}$ 按定義代入, 第四行從 $\omega_n^n=1$ 推導出, 最后一行基於推導出的組合遞歸結果用的等式.
根據這兩段推導, 我們在代碼中所作的計算能夠得到正確的$\vec y$.
第19,20行用於釋放我們開出的臨時內存空間.
上面的幾段推導都印證了 $n$ 次單位復數根與DFT的優美對稱性, 所以我們得以減少重復計算來得到更優的時間復雜度.
其中我們把 $\omega_n^k$ 稱為旋轉因子.
上述的算法實際上對應於下面的遞歸式:
\[T(n)=2T(n/2)+\Theta(n)=\Theta(nlog(n))\]
然而還有一個問題
我們到現在為止一直在討論快速求值, 並且成功地通過FFT將時間復雜度降到了 $O(nlog(n))$, 然而我們還需要插值來將點值表達轉化為系數表達. 這時我們就得...
通過 FFT 在單位復數根處插值
根據我們的四步快速卷積計算方案(倍次/求值/點積/插值), 我們還剩下最后一步就可以完成這一切了.
首先, $DFT$ 的計算過程按照定義可以表示成一個矩陣, 根據上文關於范德蒙德矩陣的等式 , 我們可以將 $DFT$ 寫成矩陣乘積 $\vec y=V_n \vec a$, 其中 $V_n$ 是用 $\omega_n$ 的一定冪次填充的矩陣, 這個矩陣大概長這樣:
\[
\begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{bmatrix}=
\begin{bmatrix}
1 & 1 & 1 & 1 & \dots & 1 \\
1 & \omega_n & \omega_n^2 & \omega_n^ 3 & \dots & \omega_n^{n-1} \\
1 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \dots & \omega_n^{2(n-1)} \\
1 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \dots & \omega_n^{3(n-1)} \\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\
1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \dots & \omega_n^{(n-1)(n-1)}
\end{bmatrix}
\begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{bmatrix}
\]
觀察矩陣中 $\omega_n$ 的冪次, 我們似乎發現 $\omega_n$ 的指數似乎組成了一張....乘法表???
因為我們將 $DFT$ 表示為 $\vec y=V_n \vec a$ , 則對於它的逆運算 $\vec a=DFT^{-1}(\vec y)$ , 我們可以選擇直接求出 $V_n$ 的逆矩陣 $V^{-1}_n$ 並乘上 $\vec y$ 來計算.
然后就會有一個玄學的定理來告訴我們逆矩陣的形式, 算導上的證明如下:
定理: 對於 $j,k=0,1,...,n-1$, $V_n^{-1}$ 的 $(j,k)$ 處元素為 $\omega_n^{-kj}/n$.
證明: 我們證明 $V_n^{-1}V_n=I_n$, 其中 $I_n$ 為一個 $n\times n$ 的單位矩陣. 考慮 $V_n^{-1}V_n$ 中 $(j,j')$ 處的元素:
\[[V_n^{-1}V_n]_{jj'}=\sum_{k=0}^{n-1}(\omega_n^{-kj}/n)(\omega_n^{kj'})=\sum_{k=0}^{n-1}\omega_n^{k(j'-j)}/n \]
當 $j'=j$ 時, 此值為 $1$ , 否則根據求和引理, 此和為 $0$.注意, 只有 $-(n-1)\leq j'-j \leq n-1$ , 從而使得 $j'-j$ 不能被 $n$ 整除,才能應用求和引理.
$\blacksquare$
上面證明的思路是, 求出 $V_n^{-1}V_n$ 的各元素的值, 並和單位矩陣 $I_n$ 比較, 發現求出的值與單位矩陣相符, 原命題得證.
得到這個定理之后我們就可以推導出 $DFT^{-1}(\vec y)$ 了, 而我們發現它是長這樣的:
\[a_j=\frac{1}{n} \sum_{k=0}^{n-1} y_k\omega_n^{-kj}\]
我們發現它和 $DFT$ 的定義極其相似! 唯一的區別僅僅在於前面的多出的 $\frac{1}{n}$ 和 $\omega_n$ 指數上多出來的負號. 也就是說, 我們可以通過稍微修改一下 $FFT$ 就可以用來計算 $DFT^{-1}$ 了!
最終我們可以得到像這樣的代碼, 通過第三個參數指定計算 $DFT$ 還是 $DFT^{-1}$ . 參數為 $1$ 則計算 $DFT$ , 參數為 $-1$ 則計算 $DFT^{-1}$. 除以 $n$ 的過程在函數執行后在函數外進行即可.
1 void FFT(Complex* a,int len,int opt){ 2 if(len==1) 3 return; 4 Complex* a0=new Complex[len/2]; 5 Complex* a1=new Complex[len/2]; 6 for(int i=0;i<len;i+=2){ 7 a0[i/2]=a[i]; 8 a1[i/2]=a[i+1]; 9 } 10 FFT(a0,len/2); 11 FFT(a1,len/2); 12 Complex wn(cos(2*Pi/len),opt*sin(2*Pi/len)); 13 Complex w(1,0); 14 for(int i=0;i<(len/2);i++){ 15 a[i]=a0[i]+w*a1[i]; 16 a[i+len/2]=a0[i]-w*a1[i]; 17 w=w*wn; 18 } 19 delete[] a0; 20 delete[] a1; 21 }
FFT 的速度優化與迭代實現方式
看到上面的實現, 根據OIer的常識, 我們顯然可以意識到: 遞歸常數比較大. 而且上述實現中拆分多項式的時候還要重新分配一段空間, 造成了時間與空間上的多余開銷. 那么, 我們是否可以通過模擬遞歸時對系數向量的重排與操作來取得更好的執行效率呢?
答案是肯定的.
首先我們可以注意到, 在上面的實現代碼中計算了兩次 $\omega_n^k y^{[1]}_k$, 我們可以將它只計算一次並將結果放在一個臨時變量中. 然后我們的循環大概會變成這樣:
1 for(int i=0;i<(len/2);i++){ 2 int t=w*a1[i]; 3 a[i]=a0[i]+t; 4 a[i+len/2]=a0[i]-t; 5 w=w*wn; 6 }
我們定義將旋轉因子與 $y^{[1]}_k$ 相乘並存入臨時變量 $t$ 中, 並從 $y^{[0]}_k$ 中增加與減去 $t$ 的操作為一個蝴蝶操作.
接着我們考慮如何將遞歸調用轉化為迭代. 首先我們觀察遞歸時的調用樹:
我們可以發現如果我們自底向上觀察這棵樹, 如果將初始向量按照葉子的位置預先重排好的話, 完全可以自底向上一步一步合並結果. 具體的過程大概像這樣:
首先成對取出元素, 對於每對元素進行 $1$ 次蝴蝶操作計算出它們的 $DFT$ 並用它們的 $DFT$ 替換原來的兩個元素, 這樣 $\vec a$ 中就會存儲有 $n/2$ 個二元素 $DFT$ ;
接着繼續成對取出這 $n/2$ 個 $DFT$ , 對於每對 $DFT$ 進行 $2$ 次蝴蝶操作計算出包含 $4$ 個元素的 $DFT$ , 並用計算出的 $DFT$ 替換掉原來的值.
重復上面的步驟, 直到計算出 $2$ 個長度為 $n/2$ 的 $DFT$ , 最后使用 $n/2$ 次蝴蝶操作即可計算出整個向量的 $DFT$.
實際實現時, 我們發現計算該輪要進行蝴蝶操作的兩個元素的下標是比較方便的, 最外層循環維護當前已經計算的 $DFT$ 長度, 每次循環完成后翻倍; 次級循環維護當前所在的 $DFT$ 的開始位置的下標, 每次循環完成后加上 $DFT$ 長度值; 最內層循環枚舉 $DFT$ 內的偏移量. 如果三個循環的循環變量分別為 $i,j,k$, 則 $j+k$ 指向當前要操作的第一個值, $j+k+i$ 指向下一個 $DFT$ 中需要計算的值.
然后關鍵的問題就在於如何快速重排得到遞歸計算時的葉子順序.
我們以長度為 $8$ 的 $DFT$ 為例, 讓我們打表找規律看看能不能發現什么:
| 原位置 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
| 原位置的二進制表示 | 000 | 001 | 010 | 011 | 100 | 101 | 110 | 111 |
| 重排后下標的二進制表示 | 000 | 100 | 010 | 110 | 001 | 101 | 011 | 111 |
| 重排結果 | 0 | 4 | 2 | 6 | 1 | 5 | 3 | 7 |
看起來似乎...是把二進制表示給翻過來了?
這個過程在算導上似乎叫做位逆序置換, 這個過程可以直接遍歷一遍並交換二進制表示互為逆序的各個下標所對應的值, 也可以預處理出來保存在一個數組中來在多次定長 $FFT$ 中減少重復計算.
代碼實現類似這樣(我在這里將位逆序置換預處理掉存在 $rev$ 數組里了)
1 void FFT(Complex* a,int len,int opt){ 2 for(int i=0;i<len;i++) 3 if(i<rev[i]) 4 std::swap(a[i],a[rev[i]]); 5 for(int i=1;i<len;i<<=1){ 6 Complex wn=Complex(cos(PI/i),opt*sin(PI/i)); 7 int step=i<<1; 8 for(int j=0;j<len;j+=step){ 9 Complex w=Complex(1,0); 10 for(int k=0;k<i;k++,w=w*wn){ 11 Complex x=a[j+k]; 12 Complex y=w*a[j+k+i]; 13 a[j+k]=x+y; 14 a[j+k+i]=x-y; 15 } 16 } 17 } 18 }
預處理部分大概這樣:
1 for(int i=0;i<bln;i++){ 2 rev[i]=(rev[i>>1]>>1)|((i&1)<<(bct-1)); 3 }
然后我們就成功地給 $FFT$ 加了個速OwO
總結
$FFT$ 作為快速計算卷積的算法, 在現在OI中較少考察高精度計算的情況下也很有用, 常常用於加速轉移方程為卷積形式的動態規划問題.
而 $FFT$ 還有一個數論版本 $NTT$ , 利用的是模意義下的原根而不是單位復根. 我可能會在后續的博文中介紹關於 $NTT$ 的部分.
參考資料:
Introduction To Algorithms , 3rd Edition , Episode 30 , Polynomials and the FFT
(繼續厚臉皮求推薦)
如果你覺得文章中有錯誤也歡迎在評論區指正OwO