淺談FFT(快速傅立葉變換)

發表於2022-02-07

前言

啊摸魚真爽哈哈哈哈哈哈 這個假期努力多更幾篇(

理解本演算法需對一些《 常 用 》數學概念比較清楚,如複數、虛數、三角函式等(不會的自己查去(其實就是懶得寫了(¬︿̫̿¬☆)
整理了一點點資料(確信 本文僅為作者的總結與完善和本人的理解與觀點,有任何誤導性錯誤請多多指出
[WARNING⚠]文筆極差,文章極度囉嗦且可能有些迷惑hhh,盡力了_(:з)∠)_


概述(可略過

離散傅立葉變換(Discrete Fourier Transform,縮寫為 DFT),是傅立葉變換在時域和頻域上都呈離散的形式,將訊號的時域取樣變換為其 DTFT 的頻域取樣。

FFT 是一種高效實現 DFT 的演算法,稱為 快速傅立葉變換(Fast Fourier Transform,FFT) 。它對傅立葉變換的理論並沒有新的發現,但是對於在計算機系統或者說數字系統中應用離散傅立葉變換,可以說是進了一大步。快速數論變換 (NTT) 是快速傅立葉變換(FFT)在數論基礎上的實現。[1]

The discovery of the Fast Fourier transformation (FFT) is attributed to Cooley and Tukey, who published an algorithm in 1965. But in fact the FFT has been discovered repeatedly before, but the importance of it was not understood before the inventions of modern computers. Some researchers attribute the discovery of the FFT to Runge and König in 1924. But actually Gauss developed such a method already in 1805, but never published it.
快速傅立葉變換(FFT)的發現歸功於Cooley和Tukey,他們在1965年發表了一種演算法。但事實上,FFT以前已經被反覆發現,但在現代計算機發明之前,人們並不瞭解它的重要性。一些研究人員將FFT的發現歸因於1924年的倫格和科尼格。但實際上,高斯早在1805年就提出了這種方法,但從未發表過。

Notice, that the FFT algorithm presented here runs in \(\Theta(n\log_2 n)\) time, but it doesn't work for multiplying arbitrary big polynomials with arbitrary large coefficients or for multiplying arbitrary big integers. It can easily handle polynomials of size \(10^5\) with small coefficients, or multiplying two numbers of size \(10^6\) , but at some point the range and the precision of the used floating point numbers will not no longer be enough to give accurate results. That is usually enough for solving competitive programming problems, but there are also more complex variations that can perform arbitrary large polynomial/integer multiplications. E.g. in 1971 Schönhage and Strasser developed a variation for multiplying arbitrary large numbers that applies the FFT recursively in rings structures running in \(\Theta(n)\) . And recently (in 2019) Harvey and van der Hoeven published an algorithm that runs in true \(\Theta(n\log_2 n)\).
請注意,這裡介紹的FFT演算法在\(\Theta(n\log_2 n)\)時間內執行,但它不適用於將任意大多項式與任意大系數相乘或將任意大整數相乘。它可以輕鬆地處理大小為\(10^5\)的多項式和小系數,或者乘以大小為\(10^6\)的兩個數字,但在某些情況下,所用浮點數的範圍和精度將不再足以給出準確的結果。這通常足以解決競爭性程式設計問題,但也有更復雜的變體可以執行任意大的多項式/整數乘法。例如,1971年,Schönhage和Strasser開發了一種將任意大數相乘的變體,該變體在\(\Theta(n)\)執行的環結構中遞回應用FFT。最近(2019年),Harvey和van der Hoeven釋出了一個演算法,該演算法以true\(\Theta(n\log_2 n)\)執行。[2]

另外,你可能會聽說過億些逼格很高的縮寫:

縮寫 全稱 作用 時間複雜度
DFT(Discrete Fourier Transform) 離散傅立葉變換 時頻域轉換 \(O(n^2)\)
FFT(Fast Fourier Transform) 快速傅立葉變換 時頻域轉換 ( 有精度誤差 ) \(O({\small\texttt{大常數}}+n\log_2n)\)
NTT/FNTT 快速數論變換 模意義下的時頻域轉換 \(O({\small\texttt{小常數}}+n\log_2n)\)
MTT 任意模數的NTT 任意模意義下的時頻域轉換 \(O(n\log_2n)\)
FWT(Fast Wavelet Transform) 快速沃爾什變換 快速集合卷積 \(O({\small\texttt{不定}})\)
FMT(Fast Möbius Transform) 快速莫比烏斯變換 逆莫比烏斯反演? \(O({\small\texttt{不定}})\)

今天主要來講前兩個
FFT(Fast Fourier Transform),譯為快速傅立葉變換,主要用於加速多項式乘法(高精乘)。樸素演算法,即DFT(Discrete Fourier Transform),譯為離散傅立葉變換,時間複雜度為\(\Theta(n^2)\),而FFT可加將其速到\(\Theta(n\log_2 n)\)


前置知識

有些基礎知識需要說一下-O-

多項式係數與點值表示方法

先上一段百度百科對FFT的定義:

快速傅立葉變換 (fast Fourier transform), 即利用計算機計算離散傅立葉變換(DFT)的高效、快速計算方法的統稱,簡稱FFT。快速傅立葉變換是1965年由J.W.庫利和T.W.圖基提出的。採用這種演算法能使計算機計算離散傅立葉變換所需要的乘法次數大為減少,特別是被變換的抽樣點數N越多,FFT演算法計算量的節省就越顯著。
FFT(Fast Fourier Transformation)是離散傅氏變換(DFT)的快速演算法。即為快速傅氏變換。它是根據離散傅氏變換的奇、偶、虛、實等特性,對離散傅立葉變換的演算法進行改進獲得的。
FFT是一種DFT的高效演算法,稱為快速傅立葉變換(fast Fourier transform)。傅立葉變換是時域一頻域變換分析中最基本的方法之一。在數字處理領域應用的離散傅立葉變換(DFT:Discrete Fourier Transform)是許多數字訊號處理方法的基礎 。

我們可以發現,FFT是在\(\Theta(n\log_2 n)\)的時間內將一個用係數表示的多項式轉化成這個多項式的點值表示演算法法。

特別地,我們通常將兩個多項式相乘稱作求卷積雖然沒啥用(:≡

係數表示

\(設一個一元n-1次n項式f(x)=\sum_{i=0}^{n-1}a_ix^i\)
係數表示法便是用每一項的係數來表示多項式\(f(x)\),即表示為$$f(x)={a_0,a_1,a_2,…,a_n-1}$$
可發現,運算過程是兩個多項式之間每個係數的乘積,複雜度為\(\Theta(n^2)\)

點值表示

\(同樣地,設一個一元n-1次n項式f(x)=\sum_{i=0}^{n-1}a_ix^i\)
可將多項式看為在平面直角座標系上的\(n\)次函式,而將\(n個互不相同的x帶入多項式當中,得到n個不同的y值\)

\[f(x_0)=y_0=\sum_{i=0}^{n-1}a_ix_0^i \\ f(x_1)=y_1=\sum_{i=0}^{n-1}a_ix_1^i \\ … \\ f(x_{n-1})=y_{n-1}=\sum_{i=0}^{n-1}a_ix_{n-1}^i\]

那麼這\(n\)個點唯一確定該多項式,此時有且僅有唯一一個多項式滿足\(\forall k(k\in [0,n)),使得y_k=f(x_k)\),即表示為

\[f(x)=\{(x_0,y_0),(x_1,y_1),…,(x_{n-1},y_{n-1})\} \]

Why?
不難想到高斯消元法,是想兩點確定直線。多一個點,可確定直線中另一引數,那麼也就是說\(n\)個點能確定\(n-1\)個引數(不考慮倍數點之類無意的點)。

假設在計算乘積的兩個多項式所選取的\(x\)序列相同,可得

\[f(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),…,(x_{n-1},f(x_{n-1}))\} \]

\[g(x)=\{(x_0,g(x_0)),(x_1,g(x_1)),…,(x_{n-1},g(x_{n-1}))\} \]

不妨設\(F(x)=f(x)\cdot g(x)\),可得

\[F(x)=\{(x_0,f(x_0)g(x_0)),(x_1,f(x_1)g(x_1)),…,(x_{n-1},f(x_{n-1})g(x_{n-1}))\} \]

可發現,上述過程便是把多項式從係數表示轉化為點值表示,將點值相乘之後,還原成係數表示,便仍是\(\Theta(n^2)\)的複雜度(瘋狂暗示(⊙x⊙;)

複數

玄學的一個東西
來自百度百科的解釋:

我們把形如a+bi(a,b均為實數)的數稱為複數,其中a稱為實部,b稱為虛部,i稱為虛數單位。當虛部等於零時,這個複數可以視為實數;當z的虛部不等於零時,實部等於零時,常稱z為純虛數。複數域是實數域的代數閉包,也即任何復係數多項式在複數域中總有根。 複數是由義大利米蘭學者卡當在十六世紀首次引入,經過達朗貝爾、棣莫弗、尤拉、高斯等人的工作,此概念逐漸為數學家所接受。

複數的定義

複數集 \(\mathbb{C}\) 中,定義 \(i=\sqrt{−1}\),一個複數 \(z\) 表示為 \(z=a+bi(a,b\in\mathbb{R})\),其中 \(a\) 稱為實部\(b\) 稱為虛部\(i\) 稱為虛數單位。特別地,當且僅當 \(a=0,b\not=0\) 時,我們稱 \(z\)純虛數,虛數和實數統稱為複數

複數的模(長) 對於虛數 \(z=a+b\text{i}\) ,定義 \(z\) 的模(長)為 \(|z|=\sqrt{a^2+b^2}\)。​其幾何意義是對應向量的模長,也是對應複平面上的點到原點的距離。

複數的輻角 假設以逆時針為正方向,從\(x\)軸正半軸到已知向量的轉角的有向角叫做幅角,
任意一個不為零的複數\(z=a+bi\)的輻角有無限多個值,且這些值相差 \(2\pi\) 的整數倍。特別地,我們把適合於 \(\pi\leq\theta<\pi\) 的輻角\(\theta\)的值,叫做輻角的主值,也稱主輻角記作\(\arg(z)\)。輻角的主值是唯一的。
若複數 \(z\) 所表示的向量與 \(x\) 軸正半軸的夾角為 \(\alpha\),則稱 \(\theta=\alpha+2k\pi,k\in \mathbb{Z}\) 為複數 \(z\) 的輻角

複數的表示

其實,可以將 \(i\) 理解成虛數單位(就像實數單位1一樣),在笛卡爾座標系中,可將x軸定義為實數軸,y軸定義為虛數軸,這樣的座標系叫做複平面。則可將複數表示在座標系上的一個點,如複數\(3i+2\)可表示成:

向量 是同時具有大小和方向的量。如上圖中的向量\(OB\)記作 \(\overrightarrow{OB}\)

複數的運算

設兩個複數分別為\(z_1=a_1+b_1i,z_2=a_2+b_2i\),則:

進行加法運算,實部與虛部分別相加,

\[z_1+z_2=(a_1+a_2)+(b_1+b_2)i \]

特別地,複數的加法滿足平行四邊形法則,即一組臨邊之和等於兩邊相夾的對角線。


表述為\(AB+AD=AC\)

進行乘法運算,拆括號後進行實部與虛部分別的合併(注:\(i^2=-1\)),

\[z_1z_2=(a_1a_2-b_1b_2)+(a_1b_2+a_2b_1)i \]

根據極座標系的相關知識,在確定複數 \(z\) 的輻角(主值) \(\theta\) 與模長 \(r=|z|\) 後複數有唯一解[3]

\[z=r(\cos\theta+\text{i}\sin\theta) \]

這便是複數的三角表示,顯然有

\[a=r\cos \theta,b=rsin\theta \]

儘管複數的三角形式在加減法時顯得十分笨拙,但複數的三角形式在乘除法方面上極具優勢,一般地,對於兩個複數

\[x=r_1(\cos\theta_1+\text{i}\sin\theta_1),y=r_2(\cos\theta_2+\text{i}\sin\theta_2) \]

根據三角恆等變換的知識,我們有

\[\begin{aligned} x×y&=r_1r_2(\cos\theta_1+i\sinθ_1) (\cosθ_2+i\sinθ_2) \\ &=r_1r_2[(\cosθ_1\cosθ_2−\sinθ_1\sinθ_2)+i(\sinθ_1\cosθ_2+\cosθ_1\sinθ_2)] \\ &=r_1r_2[\cos(θ_1+θ_2)+i\sin(θ_1+θ_2)]\end{aligned}\]

這說明兩個複數相乘得到的複數模長等於兩個複數模長之積,得到的複數輻角等於兩個複數的輻角之和。特別地,複數 \(z^n\) 的輻角為 \(n\arg(z)\),模長為 \(|z|^n\)
除法也同理,得到的複數主輻角相減,模長相除。

\[\begin{aligned} \frac{x}{y}&=\frac{r_1}{r_2}\frac{(\cosθ_1+i\sinθ_1)}{(\cosθ_2+i\sinθ_2)} \\ &=\frac{r_1}{r_2}\frac{(\cosθ_1+i\sinθ_1)(\cosθ_2−i\sinθ_2)}{(\cosθ_2+i\sinθ_2)(\cosθ_2−i\sinθ_2)} \\ &=\frac{r_1}{r_2}[\frac{(\cosθ_1\cosθ_2+\sinθ_1\sinθ_2)+i(\sinθ_1\cosθ_2−\cosθ_1\sinθ_2)}{\cos^2θ_2+\sin^2θ_2}] \\ &=\frac{r_1}{r_2}[\cos(θ_1−θ_2)+i\sin(θ_1−θ_2)]\end{aligned}\]

這便是複數相乘時的性質:模長相乘,極角相加

尤拉公式

便是這個(。・∀・)ノ

\[e^{i\theta}=\cos\theta+i\sin\theta \]

\(\sin x,\cos x\)泰勒展開

\[\begin{aligned}\sin(x)&=x-\frac{x^3}{3!}+\frac{x^5}{5!}- \frac{x^7}{7!}+\frac{x^9}{9!}\cdots \\ \cos(x)&=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\frac{x^8}{8!}\cdots\end{aligned} \]

\(e^x\)の展開:

\[e^x=1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!}+\frac{x^5}{5!}\cdots \]

因此,將 \(\cos\theta+i\sin\theta\) 用次公式展開,即

\[1+i\theta-\frac{\theta^2}{2!}-i\frac{\theta^3}{3!}+\frac{\theta^4}{4!}+i\frac{\theta^5}{5!}\cdots \]

可發現與\(e^x\)的展開式十分相似,令 \(x=i\theta\)並帶入,即

\[\begin{aligned}e^{i\theta}&=1+i\theta+\frac{i^2\theta^2}{2!}+\frac{i^3\theta^3}{3!}+\frac{i^4\theta^4}{4!}+\frac{i^5 \theta^5}{5!} \\ &=1+i\theta-\frac{\theta^2}{2!}-i\frac{\theta^3}{3!}+\frac{\theta^4}{4!}+i\frac{\theta^5}{5!}\cdots\end{aligned} \]

可發現 \(e^{i\theta}\)\(\cos\theta+i\sin\theta\) 的展開式完全相同,於是尤拉公式得證。
進一步地,尤拉公式可給出一推論:\(e^{\pi i}=-1\)

單位根

代數(學)基本定理

任何復係數一元 \(n\) 次多項式 方程在複數域上至少有一根\((n\geqslant 1)\)。對這個定理的證明目前無法繞開高等數學,由這個定理,我們可以發現任何實係數一元 \(n\) 次多項式方程在複數域上必有 \(n\) 個複數根,只需根據代數基本定理分解因式分解出所有因式即可。

一般地,關於 \(\omega\) 的方程 \(\omega^n=1,n\in \Z\) 在實數域內有且僅有一個根 \(1\),由代數基本定理,在複數域內,這個方程有 \(n\) 個根,它們統稱為 \(n\) 次單位根,本文用 \(\omega_n^i\)​ 來表示最小非負輻角第 \(i+1\) 大的 \(n\) 次單位根,可知 \(\omega_n^0=1\)
事實上,這個方程在複數域內的 \(n\) 個根為 \(\omega_n^i=\cos(i\cdot\frac {2\pi}{n})+\text{i}\sin(i\cdot\frac{2\pi}{n})\) ,因為其 \(n\) 次乘方對應的輻角為 \(n\cdot i\cdot\frac {2\pi}{n}=2i\pi\),故 \(\arg\left((\omega_n^1)^n\right)=0\)\(|(\omega_n^1)^n|=1\),因而 \((\omega_n^1)^n=1\),滿足方程。
把它們畫在複平面上可以發現,將單位圓 \(n\) 等分,做 \(n\) 個向量,那麼這 \(n\) 個向量對應複數即為 \(n\) 次單位根,以 \(1\) 為起點,逆時針分別為 \(\omega_n^i\)​。
例如當 \(n=8\) 時,在複平面上表示為(圓為單位圓)


規定這樣的次序還可以帶來一些便利,例如 \(\omega_n^i=(\omega_n^1)^i\) ,不妨再規定 \(\omega_n^1=\omega_n\) ​,那麼便可以使得 \(\omega_n^i\) ​ 既是第 \(i\)\(n\) 次單位根,又表示 \(\omega_n\) ​ 的 \(i\) 次方,同時將 \(\omega_n\)​ 的上標由 \([0,n-1]\)中的整數擴充到了 \(\mathbb{Z}\)

單位根的性質

性質1  \(w^0_n=w_n^n=1\)

性質2  \(w_n^0, w_n^1, \cdots,w_n^{n-1}\) 互不相同。

性質3  \(w^{2k}_{2n}=w^k_n\)
 代入定義式即可,也可直接觀察去掉單位圓上的奇數次單位根。

\[\begin{aligned} w^{2k}_{2n}&=\cos 2k×\frac{2\pi}{2n}+i\sin 2k×\frac{2π}{2n} \\ &=\cos k×\frac{2π}{n}+i\sin k×\frac{2π}{n} \\ &=w^k_n\end{aligned}\]

性質4  \(w_n^{k+\frac{n}{2}}=-w_n^k\)
 代入定義式即可,也可看出旋轉半圈前後關於原點對稱。

\[\begin{aligned}w_n^{k+\frac{n}{2}}&=w^k_n×w^{\frac{n}{2}}_{n} \\ &=w^k_n×[\cos(\frac{n}{2}×\frac{2π}{n})+i\sin(\frac{n}{2}×\frac{2π}{n})] \\ &=w^k_n×(\cos\pi+i\sin\pi) \\ &=w^k_n×(−1+0) \\ &=−w^k_n\end{aligned}\]

性質5  \(\sum\limits_{i=0}^{n−1}(ω^{j−k}_n)^i=\left\{ \begin{array}{rcl}0, & k≠j \\ n, & k=j\end{array}\right.\)
 1. 當 \(k\ne j\) 時,根據等比數列的求和公式,可得:

\[\sum_{i=0}^{n−1}(ω^{j−k}_n)^i=\frac{(w^{j−k}_n)^n−1}{w^{j−k}_n−1}=\frac{(w^n_n)^{j−k}−1}{w^{j−k}_n−1}=\frac{1−1}{w^{j−k}_n−1}=0 \]

2.當 \(k=j\) 時,可得:

\[\sum_{i=0}^{n−1}(ω^{j−k}_n)^i=\sum_{i=0}^{n−1} 1=n \]


正文

DFT(離散傅立葉變換)

對於任意係數的多項式來說,將其轉為點值表示法,暴力取點並計算便是 \(\Theta(n^2)\) 的複雜度,而DFT便是考慮帶入單位圓上的單位根,進而求出每個 \(\omega\) 的值

\[\omega^i_n=\cos\frac{i}{n}2\pi+i\sin\frac{i}{n}2\pi \]

可理解為

所推出的單位根便是所需要帶入的那些值

FFT(快速傅立葉變換)

\(\large\texttt{Cooley-Tukey} 演算法\)
\(FFT\) 最常見的演算法是 \(\texttt{Cooley-Tukey}\) 演算法,它的基本思路在 \(1965\) 年由 \(\texttt{J.W.Cooley}\)\(\texttt{J.W.Tukey}\) 提出的,它是一個基於DFT而採取分治策略的演算法。
上文中提到的點值表示提到過 \(n\) 次多項式可被 \(n+1\) 個點唯一確定,取點過程中是取複數點的,即 \(n\) 次單位根的 \(0\)\(n-1\) 次冪。

設多項式(函式) \(A(x)=\sum\limits_{i=0}^{n-1}a_ix^i=\) ,按下標 \(i\) 進行奇偶性分類,即

\[\begin{aligned}A(x)&=(a_0+a_2x^2+a_4x^4+...+a_{n-2}x^{n-2})+(a_1x+a_3x^3+a_5x^5+...+a_{n-1}x^{n-1}) \\ &=(a_0+a_2x^2+a_4x^4+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+a_5x^4+...+a_{n-1}x^{n-2})\end{aligned}\]

上述等式有著明顯的規律所在,很容易便想到設兩個多項式:

\[A_1(x)=a_0+a_2x+a_4x^2+...+a_{n−2}x^{\frac{n}{2}−1} \\ A_2(x)=a_1+a_3x+a_5x^2+...+a_{n−1}x^{\frac{n}{2}−1}\]

便有 \(A(x)=A_1(x^2)+xA_2(x^2)\)

ohhhhhhhhh,單位根來惹(。・ω・。)

接著,令 \(x=w_n^k (k<\frac{n}{2})\) 並將其代入函式 \(A(x)\),即

\[\begin{aligned}A(w^k_n)&=A_1(w^{2k}_n)+w^k_nA_2(w^{2k}_n)=A_1(w^{2k}_n)+w^k_nA_2(w^{2k}_n) \\ &=A_1(w^k_{\frac{n}{2}})+w^k_nA_2(w^k_{\frac{n}{2}})\end{aligned}\]

那麼,令 \(x=w_n^{k+\frac{n}{2}}\) 代入致函式 \(A(x)\),即

\[\begin{aligned}A(w^{k+\frac{n}{2}}_n)&=A_1(w^{2k+n}_n)+w^{k+\frac{n}{2}}_nA_2(w^{2k+n}_n)=A_1(w^{2k}_n\cdot w^n_n)−w^k_nA_2(w^{2k}_n\cdot w^n_n) \\ &=A_1(w^{2k}_n)−w^k_nA_2(w^{2k}_n)=A_1(w^k_{\frac{n}{2}})−w^k_nA_2(w^k_{\frac{n}{2}})\end{aligned}\]

神奇的是,上述的兩個函式 \(A(w^k_n)\)\(A(w^{k+\frac{n}{2}}_n)\) 所推出的多項式只相差一個加減符號
那麼就是說,這兩個函式的值只需在得到 $A_1​(w^​k_n​)和 \(A_2(w^k_{n\over 2})\) 便能求出其值,同時在遍歷求得其中一個函式的值時,便可用 \(\Theta(1)\) 的複雜度直接求得另一函式的值,且第一個式子的 \(k\) 在取遍了 \([0,\frac n2-1]\) 的時候,第二個式子取遍了 \([\frac n2,n-1]\) ,即原問題的規模縮小了一半,便可以用遞迴分治來解決FFT。
此時時間複雜度為 \(\Theta({\small \texttt{大常數}}+nlog_2n)\) 大常數便是複數運算得來的。因為複數乘法速度較慢,所以在資料規模不大時,應儘量避免使用FFT演算法。

時間複雜度 \(\Theta(n)\)

完了嗎?想?呢,當然還要逆回去_〆(´Д` )

IFFT(快速傅立葉逆變換)

我們上述的討論都是基於點值表示法而進行的,所以最後還需再次轉化為係數表示法,此過程便是所謂的逆變換

\((y_0,y_1,...,y_{n-1})\)\((a_0,a_1,a_2,\dots,a_{n-1})\) 的傅立葉變換(即點值表示),即 \((y_0,y_1,...,y_{n-1})\)\((a_0,a_1,a_2,\dots,a_{n-1})\)\((w_n^0,w_n^1,\dots,w_n^{n-1})\) 處的值。則有

\[\begin{aligned}y_i&=a_0+a_1w^1_n+a_2(w^1_n)^2+⋯+a_n(w^1_n)^n \\ &=a_0+a_1w^1_n+a_2w^2_n+⋯+a_nw^n_n=\sum_{j=0}^{n−1}a_j(w^i_n)^j\end{aligned}\]

設 向量\((c_0,c_1,c_2,\dots,c_{n-1})\)\((y_0,y_1,...,y_{n-1})\)\((w_n^0,w_n^{-1},\dots,w_n^{-(n-1)})\) 處的取值。
求C序列同樣用 \(\texttt{Cooley-Tukey}\) 演算法實現,只需要將只需要單位根變為 \(w_n^{-1}=w_n^{n-1}\) ,相當於順時針旋轉。可以發現,在複平面上, \(w_n^{n-1}\)\(w_n^1\)​ 的 \(x\) 座標相同,\(y\) 座標互為相反數,即為共軛複數

\((c_0,c_1,c_2,\dots,c_{n-1})\) 滿足

\[\begin{aligned}c_k&=\sum_{i=0}^{n−1}y_i(w^{−k}_n)^i=\sum_{i=0}^{n−1}(\sum_{j=0}^{n−1}a_j(w^i_n)^j)(w^{−k}_n)^i \\ &=\sum_{i=0}^{n−1}(\sum_{j=0}^{n−1}a_j(w^j_n)^i)(w^{−k}_n)^i=\sum_{i=0}^{n−1}(\sum_{j=0}^{n−1}a_j(w^j_n)^i(w^{−k}_n)^i) \\ &=\sum_{i=0}^{n−1}\sum_{j=0}^{n−1}a_j(w^j_n)^i(w^{−k}_n)^i=\sum_{i=0}^{n−1}\sum_{j=0}^{n−1}a_j(w^{j−k}_n)^i \\ &=\sum_{j=0}^{n−1}a_j(\sum_{i=0}^{n−1}(w^{j−k}_n)^i)=\sum_{j=0}^{n−1}a_j(\sum_{i=0}^{n−1}(w^{j−k}_n)^i)\end{aligned}\]

根據性質5,便有

  1. \(j\ne k\) 時,即 \(j-k\ne 0\) ,此時 \(\sum\limits_{i=0}^{n-1}(w_n^{j-k})^i=0\) ,共有 \(n-1\) 種情況
  2. \(j=k\) 時,即 \(j-k=0\) ,此時 \(\sum\limits_{i=0}^{n-1}(w_n^{j-k})^i=n\) ,共有 \(1\) 種情況。
    因此,$$c_k=na_k,\dfrac{c_k}n=a_k$$

如上,我們得到了點值表示法與係數表示法之間的關聯式,便可進行逆運算
總結性地來說,可概括為函式分治時所乘單位根的共軛複數,分治後每一項除以\(n\)即為原多項式的每一項係數,因此FFT與IFFT可同時完成。

特別地,C++中帶有複數模板庫#include <complex>


大概就是complex<double> cp定義,.real()實部.imag()虛部等等 [4]

同時,FFT只能處理 \(n\)\(2\) 的整數次冪的多項式(函式)

啊但是,樸素演算法的遞迴常數較大,需考慮優化(σ`д′)σ

FFT優化

非遞迴

通過觀察,我們可得出我們所求的序列便是原序列下表二進位制的翻轉後得到的結果,因此樸素演算法過程中可省略下標的奇偶性分類。此時,秩序 \(\Theta(n)\) 的複雜度來得到我們所求的序列,接著便是向上合併的過程了

蝴蝶操作

啊回頭來補嘛沒學會(逃(¬︿̫̿¬☆)

題表(坑

總結


後記

引用

  1. https://cp-algorithms.com/algebra/fft.html
  2. https://oi-wiki.org/math/poly/fft/
  3. https://zhuanlan.zhihu.com/p/143276207
  4. https://blog.csdn.net/linjiayang2016/article/details/80341958
  5. https://blog.csdn.net/enjoy_pascal/article/details/81478582
  6. https://blog.csdn.net/ggn_2015/article/details/68922404
  7. https://blog.csdn.net/leo_h1104/article/details/51615710
  8. https://www.cnblogs.com/RabbitHu/p/FFT.html
  9. https://www.cnblogs.com/zwfymqz/p/8244902.html
  10. https://baike.baidu.com/item/快速傅立葉變換/214957
  11. https://baike.baidu.com/item/FFT原理/8966333
  12. https://baike.baidu.com/item/離散傅立葉變換/6379901
  13. https://baike.baidu.com/item/複數/254365
    嗯,就挺多的(逃

天哪這算是彙總全了吧,碼字碼去世了?一定要留個贊啊啊QAQ緩兩天繼續更
-END-

相關文章