[學習筆記] 多項式與快速傅里葉變換(FFT)基礎

時間 2019-11-13

標籤學習筆記多項式快速傅里葉變換 fft 基礎简体版

原文原文鏈接

==== €€£ WARNING ====html

這篇博文內容相對偏少, 已經在後續博文中擴充.
算法

你們能夠看個人最新博文 [學習筆記&教程] 信號, 集合, 多項式, 以及各類卷積性變換 (FFT,NTT,FWT,FMT)數組

==== ====函數

引入

可能有很多OIer都知道FFT這個神奇的算法, 經過一系列玄學的變化就能夠在 $O(nlog(n))$ 的總時間複雜度內計算出兩個向量的卷積(或者多項式乘法/高精度乘法), 而代碼量卻很是小. 博主一年半前曾經因COGS的一道叫作"神祕的常數 $\pi$"的題目而去學習過FFT, 可是基本就是照着板子打打完並不知道本身在寫些什麼鬼畜的東西OwO學習

不過...博主這幾天忽然照着算法導論本身看了一遍發現本身彷佛忽然意識到了什麼OwO而後就打了一道板子題還1A了OwO再加上午考試差點AK以及日更頻率即將不保因而就有了這篇博文233優化

博主在寫這篇文字的過程當中也發現了很多本身以前忽略的FFT細節, 感受對FFT彷佛有了更深入的理解, 但願想學習FFT的讀者能認真看完這篇文章並仔細體味FFT的優雅性質.spa

本文可能並不會介紹關於使用FFT進行信號處理的相關信息, 主要介紹在OI中的應用即快速卷積.code

Rush了很久終於寫出來了QAQorm

那麼如今, 就讓咱們變玄學爲科學瞭解 $FFT$ 背後的工做原理與它的優雅性質htm

前置知識: 多項式的表示

多項式的係數表達

相信你們都知道對於一個次數界爲 $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

相關標籤/搜索

每日一句

每一个你不满意的现在，都有一个你没有努力的曾经。