快速傅里葉變換(FFT)學習筆記(其一)

再探快速傅里葉變換(FFT)學習筆記(其一)

目錄

• 再探快速傅里葉變換(FFT)學習筆記(其一)
  • 寫在前面
    • 為什么寫這篇博客
    • 一些約定
  • 前置知識
    • 多項式卷積
    • 多項式的系數表達式和點值表達式
    • 單位根及其性質
  • DFT和IDFT
    • DFT的過程
    • IDFT的過程
  • FFT
    • FFT的數學證明及時間復雜度分析
    • FFT的遞歸實現
    • FFT的非遞歸實現
    • FFT的局限性
  • 例題

寫在前面

為什么寫這篇博客

筆者去年暑假剛剛學習過FFT,NTT的一些基礎應用。但當時對FFT和NTT的理解還不夠深入。本博客參考2016年國家集訓隊論文中雅禮中學毛嘯的《再探快速傅立葉變換》,對之前學習時的不足之處做了補充。

為了不使篇幅過長，預計將把學習筆記分為四部分:

1. DFT,IDFT,FFT的定義,實現與證明
2. NTT的實現與證明,任意模數NTT
3. FFT的優化技巧
4. 多項式相關操作

注意:本文為了嚴謹性可能會犧牲一些可讀性

陰晴朝夕態倏忽兮，紛變換乎秋春窈窕繚繞，別有天地兮，莽不知漢與秦之日月。
—— 〔明〕張萱 黃山歌為吳孝父作

一些約定

1. \([p(x)]=\begin{cases}1,p(x)為真 \\ 0,p(x)為假 \end{cases}\)

2. 本文中序列的下標從0開始

3. 若\(s\)是一個序列,\(|s|\)表示\(s\)的長度

4. 若大寫字母如\(F(x)\)表示一個多項式,那么對應的小寫字母如\(f\)表示多項式的每一項系數，即\(F(x)=\sum_{i=0}^{n-1} f_ix^i\)

5. 用\(DFT(a)\)表示a做DFT后的結果,IDFT同理

前置知識

多項式卷積

多項式\(A(x)=\sum_{i=0}^{n-1} a_ix^i\)(注意最高次項的指數為\(n-1\),至於為什么要這樣定義見下文)

定義1.1 對於兩個多項式\(A(x)=\sum_{i=0}^{n-1} a_ix_i,B(x)=\sum_{i=0}^{m-1} b_ix_i\),我們定義\(A(x)\)和\(B(x)\)的卷積為:

\[A(x) \cdot B(x)=\sum_{i=0}^{n+m-2} x^i\sum_{k=0}^i a_kb_{i-k} \]

實際上我們只要考慮系數，也就是給出兩個序列\(a,b\),求序列\(c\)使得$$c_i=\sum_{k=0}^i a_k b_{i-k} \tag{1.1}$$

直接計算卷積顯然是\(O(n^2)\)的,FFT可以用來優化這個過程

多項式的系數表達式和點值表達式

系數表達式就是多項式中每一項的系數\(a_0,a_1 \dots a_{n-1}\)

我們知道,n個點可以確定一個n-1次多項式.那么，一個n-1次多項式\(F(x)\)也可以表示成這樣的一個集合\(\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)) \dots (x_{n-1},f(x_{n-1}))\}\).比如多項式\(f(x)=x^2+x+1\)就可以表示成\(\{(1,3),(2,7),(-1,1)\}\)

這樣有什么用呢?如果已知\(A(x)\)和\(B(x)\)的點值表達式,可以\(O(n)\)計算\(C(x)=A(x)B(x)\)的點值表達式

比如\(A(x)=\{(x_0,A(x_0)),(x_1,A(x_1)),(x_2,A(x_2)),....,(x_n,A(x_n))\}\)

\(B(x)=\{(x_0,B(x_0)),(x_1,B(x_1)),(x_2,B(x_2)),....,(x_n,B(x_n))\}\)

\(C(x)=A(x)\times B(x)=\{(x_0,A(x_0)\times B(x_0)),(x_1,A(x_1)\times B(x_1)),(x_2,A(x_2)\times B(x_2)),....,(x_n,A(x_n)\times B(x_n))\}\)

所以:我們可以在\(O(n)\)的時間內計算出兩個點值表達式的乘積,只需要把對應點的\(y\)坐標相乘即可。但是,把系數表達式和點值表達式互相轉化仍然需要\(O(n^2)\)時間.(每代入一個\(x\)都需要\(O(n)\)時間計算\(A(x)\))

單位根及其性質

定義1.2:\(n\)次單位根為滿足:\(n\)是最小的滿足\(\omega^n=1\)的復數\(\omega\)

容易發現,\(n\)次單位根有\(n\)個,他們分布在復平面上的單位圓上。(這里我們只討論\(n\)為偶數的情況)

把單位圓\(n\)等分,圓上的每個分割點都對應一個單位根.從點\((1,0)\)開始,逆時針把這\(n\)個點從\(0\)開始編號,第\(k\)個點對應的復數記作\(w_n^k\),上圖呈現的就是\(\omega_{16}^0 \text{~} \omega_{16}^{15}\)

容易發現\(\omega_n^k\)的幅角為\(\frac{2k\pi}{n}\).(一整個圓為\(2\pi \ \text{rad}\),分為\(n\)份,再取\(k\)份)，因此，有:

\[\omega_n^k=\cos \frac{2k\pi}{n}+\text{i} \sin \frac{2k\pi}{n} \tag{1.2.1} \]

單位根有如下的性質:

定理1.2.1 \(\forall p \in \mathbb{Z}, \omega_n^k=\omega_{pn}^{pk}\)

代數證明:

\(\omega_{pn}^{pk}=\cos \frac{2pk \pi}{pn}+ \text{i}\sin \frac{2pk\pi}{pn}=\cos \frac{2k\pi}{n}+\text{i} \sin \frac{2k\pi}{n}=\omega_n^k\)

幾何證明:單位圓分成\(n\)份，數\(k\)次，和分成\(pn\)份，數\(pk\)次所代表復數是一樣的

定理1.2.2 \(\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^k\)

代數證明:

\[\begin{aligned} \omega_{n}^{k+\frac{n}{2}} &= \cos \frac{2(k+\frac{n}{2}\pi)}{n}+ \text{i} \sin \frac{2(k+\frac{n}{2}\pi)}{n}\\ &=\cos(\frac{2k \pi}{n}+\pi)+\text{i}\sin(\frac{2k \pi}{n}+\pi) \\&= -\cos \frac{2k\pi}{n}-\text{i} \sin \frac{2k\pi}{n}=-\omega_n^k\end{aligned} \]

倒數第二步用到了三角函數的誘導公式

幾何證明:

\(\frac{n}{2}\)代表半圈，所以這個式子的意思是從編號為\(k\)的復數開始走半圈，代表的復數和之前的正好相反（實部和虛部都相反）

定理1.2.3 \(\omega_n^p \times \omega_n^q=\omega_n^{p+q}\)

令\(\alpha=\frac{2p\pi}{n},\beta=\frac{2p\pi}{n}\)

則\(\omega_n^p=\cos \alpha+ \text{i} \sin \alpha,\omega_n^q=\cos \beta+ \text{i} \sin \beta\)

\[\begin{aligned} \omega_{n}^{p} \times \omega_{n}^{q}&=(\cos \alpha \cos \beta-\sin \alpha \sin \beta)+i(\sin \alpha \cos \beta+\cos \alpha \sin \beta)\\ &=\cos (\alpha+\beta)+i \sin (\alpha+\beta)\\&=\omega_{n}^{p+q} \end{aligned}\]

證明用到了復數乘法和三角函數的和角公式.注意到定理1.2.2實際上是定理1.2.3的推論，證明留給讀者。這個公式極為重要，在FFT的公式推導中會用到。

定理1.2.3的另一個推論，\(\omega_n^k=(\omega_n^1)^k\)

證明顯然

定理1.2.4

\[\forall v \in \mathbb{N},\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{v k}=[v \bmod n=0] \]

證明:

當\(v \bmod n =0\)時, \(\omega_{n}^{vk}=\cos\frac{2vk \pi}{n}+\text{i} \sin \frac{2vk \pi}{n}\)

​ 那么\(\frac{vk}{n}\)一定是整數,所以\(\cos\frac{2vk \pi}{n}=1,\sin \frac{2vk \pi}{n}=0,\omega_n^{vk}=1\)

​ 所以\(\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{v k}=\frac{1}{n} \times n =1\)

當\(v \mod n \neq 0\)時,容易發現\(\omega_n^v \neq 1\),

​ 那么上式實際上是一個等比數列,根據等比數列求和公式有:

​ $$\begin{aligned} \frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{v k}&=\frac{1-\omega_n^{{vn}}{n(1-\omega_n}v)} \ &= \frac{1-\omega_1^{v}{n(1-\omega_n}v)}\end{aligned}$$

​ 而\(\omega_1^v=\cos(2v\pi)+i\sin(2v \pi)=1\)

​ 因此\(\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{v k}=0\)

證畢.

這個定理會在IDFT的正確性證明中用到

DFT和IDFT

在上文中我們提到了多項式的系數表達式和點值表達式，並發現點值表達式可以簡化卷積運算。

離散傅里葉變換(Discrete Fourier Transform,DFT):是把系數表達式轉化為點值表達式的一種算法

逆離散傅里葉變換(Inverse Discrete Fourier Transform,IDFT):是把點值表達式轉化為系數表達式的一種算法。

那么求\(A\)和\(B\)的卷積就可以先對\(A\)和\(B\)做DFT,再將\(A\)和\(B\)的點值表達式相乘得到\(C\),最后再對\(C\)做IDFT

DFT的過程

對於多項式\(A,B\),我們找到最小的大於\(|A|+|B|\)的2的整數冪\(n\),至於為什么是2的整數冪后面會提到。給定序列\(a\),對於\(k \notin [0,|a|)\),我們定義\(a_k=0\)

要想得到一個系數表達式的點值表達式，我們可以任意選\(n\)個值代入。但是傅里葉想到,單位根有着特殊的性質，可以把\(n\)次單位根代入多項式。形式化的說,設多項式\(A(x)=\sum_{i=0}^{n-1}a_ix^i\),則\(A(x)\)的DFT變換結果為\(a'_k=A(\omega_n^k)=\sum_{i=0}^{n-1}a_i \omega_n^{ki}\)

我們后面提到的多項式與系數序列的對應關系都類似\(A(x)\)與\(a\)的對應關系，對\(A(x)\)做DFT就是對序列\(a\)做DFT;說一個東西做DFT得到了\(A(x)\),就是得到了序列\(a\)

容易發現DFT的時間復雜度為\(O(n^2)\)

IDFT的過程

代入\(n\)次單位根的作用在IDFT中顯現出來.

定理2.1 :把多項式\(A(x)\)的離散傅里葉變換結果\(f(\omega_n^0),f(\omega_n^1) \dots f(\omega_n^{n-1})\)作為另一個多項式\(B(x)\)的系數，取單位根的共軛復數即\(\omega_n^0,\omega_n^{-1},\omega_n^{-2} \dots \omega_n^{-(n-1)}\)作為\(x\)代入\(B(x)\)，得到的每個數再除以\(n\)，得到的是\(A(x)\)的各項系數

也就是說,把點值表達式中點的\(y\)坐標作為另一個多項式的系數后,IDFT的過程和DFT幾乎一樣,只是代入的是單位根的共軛復數。

我們暫時不證明定理2.1,而是證明一個更強的定理:

定理2.2:把多項式\(A(x)\)的DFT結果\(a\),和\(B(x)\)的DFT結果\(b\)對應相乘,得到序列\(t_i=a_ib_i\);再對\(t_i\)做IDFT,得到的序列就是\(C(x)=A(x)\cdot B(x)\)的各項系數\(c_i\)

證明:

由卷積的定義式\((1.1)\)我們有:

\[c_{r}=\sum_{p, q}[(p+q) \bmod n=r] a_{p} b_{q} \tag{2.1} \]

看到取模想到定理1.1.4

\[\forall v \in \mathbb{N},\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{v k}=[v \bmod n=0] \tag{2.2} \]

故:

\[\begin{aligned} &[(p+q) \bmod n=r] \\ =&[(p+q-r) \bmod n=0] \\ =& \frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{(p+q-r) k} \\ =& \frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{-r k} \omega_n^{p k} \omega_n^{q k} \end{aligned}\]

代入\((2.1)\)

\[\begin{align*} c_{r} &=\sum_{p, q}[(p+q) \bmod n=r] a_{p} b_{q} \\ &=\sum_{p, q} \frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{-r k} \omega_n^{p k} \omega_n^{q k} a_{p} b_{q} \\ &=\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{-r k} \sum_{p, q} \omega_n^{p k} a_{p} \omega_n^{q k} b_{q} \\ &=\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{-r k} \sum_{p} \omega_n^{p k} a_{p} \sum_{q} \omega_n^{q k} b_{q}\tag{2.3} \\&=\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{-r k} a'_k b'_k \\&=\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{-r k} t_k\end{align*} \]

容易發現定理2.1是定理2.2 的推論,即\(B(x)\)的系數都為1時的情況.式\((2.3)\)實際上指出了DFT和IDFT求卷積的過程。

和DFT類似,IDFT的時間復雜度為\(O(n^2)\)

FFT

上一部分我們提到了DFT的公式，但是這樣直接做還是太慢了。因此我們要利用單位根的性質來加速DFT,這個算法被稱為快速傅立葉變換(Fast Fourier Transform,FFT). 同樣,我們也可以加速IDFT,即IFFT.

這幾個算法的關系如下圖所示:

其中FFT,IFFT的時間復雜度均為\(O(n \log n)\)

FFT的數學證明及時間復雜度分析

我們對多項式\(A(x) = a_0 + a_1x + a_2x^2 +...+a_{n-1}x^{n-1}\)進行FFT,仍然把\(\omega_n^k\)代入

把\(a\)按照下標奇偶性分為兩部分:

\(A_1(x) = a_0 + a_2x + ... + a_{n - 2}x^{\frac{n}{2} - 1}\)

\(A_2(x) = a_1 + a_3x + ... + a_{n - 1}x^{\frac{n}{2} - 1}\)

那么我們有:

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

將\(x=\omega_n^k\)代入,

對於\(k < \frac{n}{2}\)(根據先前的定義,n為2的次冪,所以沒有取整問題)

\[\begin{align*} A(\omega_n^k) &= A_1(\omega_n^{2k}) + \omega_n^kA_2(\omega_n^{2k}) \\ &= A_1(\omega_{\frac{n}{2}}^{k}) + \omega_n^kA_2(\omega_{\frac{n}{2}}^{k}) \tag{3.1}\end{align*} \]

​ 證明用到了定理1.2.1:\(\forall p \in \mathbb{Z}, \omega_n^k=\omega_{pn}^{pk}\)

對於剩下的部分,即\(\omega_{n}^{k+\frac{n}{2}}(k<\frac{n}{2})\)

\[\begin{align*} A(\omega_n^{k + \frac{n}{2}}) &= A_1(\omega_n^{2k + n}) + \omega_n^{k + \frac{n}{2}}A_2(\omega_n^{2k + n}) \\ &= A_1(\omega_{\frac{n}{2}}^{k} \times \omega_n^n) + \omega_n^{k + \frac{n}{2}} A_2(\omega_{\frac{n}{2}}^{k} \times \omega_n^n) \\ &= A_1(\omega_{\frac{n}{2}}^{k}) - \omega_n^kA_2(\omega_{\frac{n}{2}}^{k}) \tag{3.2} \end{align*} \]

證明用到了定理1.2.1和定理1.2.3:\(\omega_n^p \times \omega_n^q=\omega_n^{p+q}\)

式\((3.1),(3.2)\)也被稱為“蝴蝶操作”, 通過蝴蝶操作,我們只要知道兩個多項式\(A_1(x),A_2(x)\)在\((\omega_{\frac{n}{2}}^{0}, \omega_{\frac{n}{2}}^{1}, \omega_{\frac{n}{2}}^{2}, ... , \omega_{\frac{n}{2}}^{\frac{n}{2} - 1})\)處的點值表達式,就可以求出\(A(x)\)在\((\omega_n^0, \omega_n^1, \omega_n^2, ..., \omega_n^{n-1})\)處的點值表達式了了。每次遞歸下去規模縮小一半,\(n=1\)時終止.由於長度\(n\)每次除2,且要求每次的\(n\)均為偶數，所以我們需要\(n\)是2的正整數次冪

設給長度為\(n\)的序列做\(FFT\)的時間為\(T(n)\),那么我們有

\[T(n)=2T(\frac{n}{2})+n \]

根據主定理,\(T(n)=\Theta(n \log n)\)

FFT的遞歸實現

FFT遞歸實現很好理解.注意可以使用C++自帶的復數類<complex>,也可以自己手寫復數類。

<complex>重載了+,-,*,/,=,==等運算符，在使用上很方便，幾乎和操作整型變量沒有區別。兩個成員函數.real(),.imag()分別返回復數的實部和虛部。更多它的用法請參考這個鏈接

#include<complex>//c++自帶復數類
const double pi=acos(-1.0);
typedef complex<double> com;
com a[maxn+5],b[maxn+5],c[maxn+5];
void fft(com *x,int n,int type){//orz Fourier 
	if(n==1) return;//遞歸邊界
	com l[n>>1],r[n>>1];
	for(int i=0;i<n;i++){//按奇偶分類
		if(!(i&1)) l[i>>1]=x[i]; 
		else r[i>>1]=x[i];
	}
	fft(l,n>>1,type);//遞歸對n/2的長度做FFT
	fft(r,n>>1,type);
	
	com wn1=com(cos(2*pi/n),type*sin(2*pi/n))/*w_n^0*/,wnk=com(1,0);	
	for(int i=0;i<(n>>1);i++){
		x[i]=l[i]+wnk*r[i];//式(3.1)
		x[i+(n>>1)]=l[i]-wnk*r[i];//式(3.2)
		wnk*=wn1;//遞推算w_n^k
	}
}
int main(){
     //....
     int m=n*2;//這個數根據題目算要求來定,這里是長度為n的兩個序列卷積,那就取2n
	 n=1;
	 while(n<m) n*=2;//計算出最小的2的次冪
	 fft(a,n,1);
	 fft(b,n,1);
	 for(int i=0;i<n;i++) c[i]=a[i]*b[i];
	 fft(c,n,-1);
     for(int i=0;i<m;i++) printf("%lf ",c[i].real()/n);//記得除n
    //....
}

FFT的非遞歸實現

很遺憾,FFT的遞歸實現常數極大(雖然非遞歸實現常數也很大),幾乎無法通過任何OI中FFT的題目。

我們觀察遞歸過程中根據奇偶改變系數位置的方式,例如長度為8的情況:

\[\begin{aligned} &\left(a_{0}, a_{1}, a_{2}, a_{3}, a_{4}, a_{5}, a_{6}, a_{7}\right)\\ &\left(a_{0}, a_{2}, a_{4}, a_{6}\right)\left(a_{1}, a_{3}, a_{5}, a_{7}\right)\\ &\left(a_{0}, a_{4}\right)\left(a_{2}, a_{6}\right)\left(a_{1}, a_{5}\right)\left(a_{3}, a_{7}\right)\\ &\left(a_{0}\right)\left(a_{4}\right)\left(a_{2}\right)\left(a_{6}\right)\left(a_{1}\right)\left(a_{5}\right)\left(a_{3}\right)\left(a_{7}\right) \end{aligned}\]

觀察\(0,4,2,6,1,5,3,7\)的二進制表示\(000,100,010,110,001,101,011,111\)，

比較\(0,1,2,3,4,5,6,7\)的二進制表示\(000,001,010,011,100,101,110,111\)

我們發現,初始時位置\(x\)上的數會被交換到x二進制位翻轉后的位置,記為\(rev(x)\)。注意這里的反轉指的是反過來寫,而不是每一位交換01。這種操作也被稱為“位逆序置換”。那么只要處理出這個序列,每次FFT的都是一段連續的區間,可以非遞歸求解。

\(rev(x)\)可以根據如下代碼遞推求出

for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));

證明: 對於數\(i\),我們已經求出\(i>>1\)的置換結果,如\(k=5,i=(10111)_2,i>>1=(1011)_2\);那么我們先反轉從左到右前\(k-1\)位,即\(rev[i>>1]=rev[(01011)_2]=(11010)_2\),右移一位去掉前導0(翻轉后到了末尾)就得到了前\(k-1\)位反轉的結果; 然后再加上最后一位,最后一位的值是\((i\text{&}1)\),要左移到最前面.

int rev[maxn+5];//數x位逆序置換后變成rev[x]
int n,k;//n=2^k
void fft(com *x,int n,int type){//type=1做FFT,type=-1做IFFT
    for(int i=0;i<n;i++) if(i<rev[i]) swap(x[i],x[rev[i]]);//做置換,i<rev[i]防止一個數被換兩次換回原位
    for(int len=1;len<n;len*=2){
        int sz=len*2;//當前處理的長度為len*2,len=1的情況不用計算,就是x[i]
        com wn1=com(cos(2*pi/sz),sin(2*pi/sz)*type);//計算n次單位根w_n^1,type=-1時得到共軛復數
        for(int l=0;l<n;l+=sz){//枚舉一對長度為len的區間，對這兩個區間做蝴蝶操作
            int r=l+len-1;//[l,l+len-1]為前半段,[l+len,r]為后半段
            com wnk=1;
            for(int i=l;i<=r;i++){
                com tmp=x[i+len];//用tmp臨時存儲值
                x[i+len]=x[i]-wnk*tmp;//式(3.1)
                x[i]=x[i]+wnk*tmp;//式(3.2)
                wnk*=wn1;//遞推計算w_n^k
            }
        }
    }
    if(type==-1) for(int i=0;i<n;i++) x[i]/=n;//IFFT記得除n
}
int main(){
    n=1,k=0;
	while(n<m){
        n*=2;
        k++;
    }//預處理n,k
    for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));//預處理rev
    fft(a,n,1);
    fft(b,n,1);
    for(int i=0;i<n;i++) ans[i]=a[i]*b[i];
    fft(ans,n,-1);
}

FFT的局限性

顯然FFT的過程中會出現浮點數的精度誤差。比如在遞推單位根\(w_n^k=w_n^{k-1}\times w_n^1\)時的精度誤差比較大， 但一般可以接受。如果追求更高精度，可以對於每個\(k\)重新計算\(\omega_n^k=\cos \frac{2k\pi}{n}+\text{i} \sin \frac{2k\pi}{n}\),但時間開銷較大。

如果有取模運算，可以使用第二節中提到的NTT算法

另外，在第三節中，我們會介紹減少FFT常數的幾種方法，以及對任意長度的序列做DFT的方法。

例題

[1.BZOJ 3527] [ZJOI2014]力(FFT)

2.[Codeforces 580D]Fizzy Search(FFT)