快速傅里葉變換(FFT)詳解

本文只討論FFT在信息學奧賽中的應用

文中內容均爲我的理解，若有錯誤請指出，不勝感激

前言

先解釋幾個比較容易混淆的縮寫吧

DFT：離散傅里葉變換—>$O(n^2)$計算多項式乘法

FFT：快速傅里葉變換—>$O(n*\log(n)$計算多項式乘法

FNTT/NTT：快速傅里葉變換的優化版—>優化常數及偏差

FWT：快速沃爾什變換—>利用相似FFT的東西解決一類卷積問題

MTT：毛爺爺的FFT—>很是nb/任意模數

FMT 快速莫比烏斯變化—>感謝stump提供

多項式

係數表示法

設$A(x)$表示一個$n-1$次多項式

則$A(x)=\sum_{i=0}^{n} a_i * x^i$

例如：$A(3)=2+3*x+x^2$

利用這種方法計算多項式乘法複雜度爲$O(n^2)$

（第一個多項式中每一個係數都須要與第二個多項式的每一個係數相乘）

點值表示法

將$n$互不相同的$x$帶入多項式，會獲得$n$個不一樣的取值$y$

則該多項式被這$n$個點$(x_1,y_1),(x_2,y_2),\dots,(x_n,y_n)$惟一肯定

其中$y_i=\sum_{j=0}^{n-1} a_j*x_i^j$

例如：上面的例子用點值表示法能夠爲$(0,2),(1,6),(2,12)$

利用這種方法計算多項式乘法的時間複雜度仍然爲$O(n^2)$

（選點$O(n)$，每次計算$O(n)$）

 

咱們能夠看到，兩種方法的時間複雜度都爲$O(n^2)$，咱們考慮對其進行優化

對於第一種方法，因爲每一個點的係數都是固定的，想要優化比較困難

對於第二種方法，貌似也沒有什麼好的優化方法，不過當你看完下面的知識，或許就不這麼想了

 

複數

在介紹複數以前，首先介紹一些可能會用到的東西

向量

同時具備大小和方向的量

在幾何中一般用帶有箭頭的線段表示

圓的弧度制

等於半徑長的圓弧所對的圓心角叫作1弧度的角，用符號rad表示，讀做弧度。用弧度做單位來度量角的制度叫作弧度制

公式:

$1^{\circ }=\dfrac{\pi}{180}rad$

$180^{\circ }=\pi rad$

平行四邊形定則

 

(好像畫的不是很標準。。)

平行四邊形定則：AB+AD=AC

 

複數

定義

設$a,b$爲實數，$i^2=-1$，形如$a+bi$的數叫複數，其中$i$被稱爲虛數單位，複數域是目前已知最大的域

在複平面中，$x$表明實數，$y$軸（除原點外的點）表明虛數，從原點$(0,0)$到$(a,b)$的向量表示複數$a+bi$

模長：從原點$(0,0)$到點$(a,b)$的距離，即$\sqrt{a^2+b^2}$

幅角：假設以逆時針爲正方向，從$x$軸正半軸到已知向量的轉角的有向角叫作幅角

運算法則

加法：

由於在複平面中，複數能夠被表示爲向量，所以複數的加法與向量的加法相同，都知足平行四邊形定則（就是上面那個）

乘法：

幾何定義：複數相乘，模長相乘，幅角相加

代數定義：$$(a+bi)*(c+di)$$

$$=ac+adi+bci+bdi^2$$

$$=ac+adi+bci-bd$$

$$=(ac-bd)+(bc+ad)i$$

 

單位根

下文中，默認$n$爲$2$的正整數次冪

在複平面上，以原點爲圓心，$1$爲半徑做圓，所得的圓叫單位圓。以圓點爲起點，圓的$n$等分點爲終點，作$n$個向量，設幅角爲正且最小的向量對應的複數爲$\omega_n$，稱爲$n$次單位根。

根據複數乘法的運算法則，其他$n-1$個複數爲$\omega_n^2,\omega_n^3,\ldots,\omega_n^n$

注意$\omega_n^0=\omega_n^n=1$（對應複平面上以$x$軸爲正方向的向量）

那麼如何計算它們的值呢？這個問題能夠由歐拉公式解決$$\omega_{n}^{k}=\cos\ k *\frac{2\pi}{n}+i\sin k*\frac{2\pi}{n}$$

 

 

例如

 

圖中向量$AB$表示的複數爲$8$次單位根

單位根的幅角爲周角的$\frac{1}{n}$

在代數中，若$z^n=1$，咱們把$z$稱爲$n$次單位根

單位根的性質

• $\omega _{n}^{k}=\cos k\dfrac{2\pi}{n}+i\sin k\dfrac {2\pi }{n}$（即上面的公式）

• $\omega _{2n}^{2k}=\omega _{n}^{k}$

證實：

$$\omega _{2n}^{2k}=\cos 2k*\frac{2\pi}{2n}+i\sin2k*\frac{2\pi}{2n}$$

$$=\omega _{n}^{k}$$

• $\omega _{n}^{k+\frac{n}{2}}=-\omega _{n}^{k}$

$$\omega _{n}^{\frac{n}{2}}=\cos\frac{n}{2}*\frac{2\pi}{n}+i\sin\frac{n}{2}*\frac{2\pi}{n}$$

$$=\cos \pi+i\sin\pi$$

$$=-1$$

• $\omega _{n}^{0}=\omega _{n}^{n}=1$

 講了這麼多，貌似跟咱們的正題沒啥關係啊。。

 OK！各位坐穩了，前方高能！

快速傅里葉變換

咱們前面提到過，一個$n$次多項式能夠被$n$個點惟一肯定。

那麼咱們能夠把單位根的$0$到$n-1$次冪帶入，這樣也能夠把這個多項式肯定出來。可是這樣仍然是$O(n^2)$的呀！

咱們設多項式$A(x)$的係數爲$(a_o,a_1,a_2,\ldots,a_{n-1})$

那麼$$A(x)=a_0+a_1*x+a_2*{x^2}+a_3*{x^3}+a_4*{x^4}+a_5*{x^5}+ \dots+a_{n-2}*x^{n-2}+a_{n-1}*x^{n-1}$$

將其下標按照奇偶性分類

$$A(x)=(a_0+a_2*{x^2}+a_4*{x^4}+\dots+a_{n-2}*x^{n-2})+(a_1*x+a_3*{x^3}+a_5*{x^5}+ \dots+a_{n-1}*x^{n-1})$$

設

 

$$A_1(x)=a_0+a_2*{x}+a_4*{x^2}+\dots+a_{n-2}*x^{\frac{n}{2}-1}$$

$$A_2(x)=a_1+a_3*{x}+a_5*{x^2}+ \dots+a_{n-1}*x^{\frac{n}{2}-1}$$

那麼不可貴到

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

咱們將$\omega_n^k (k<\frac{n}{2}) $代入得

$$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})$$

同理，將$\omega_n^{k+\frac{n}{2}}$代入得

$$A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}(\omega_n^{2k+n})$$

$$=A_1(\omega_n^{2k}*\omega_n^n)-\omega_n^kA_2(\omega_n^{2k}*\omega_n^n)$$

$$=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k})$$

 

你們有沒有發現什麼規律？

沒錯！這兩個式子只有一個常數項不一樣！

那麼當咱們在枚舉第一個式子的時候，咱們能夠$O(1)$的獲得第二個式子的值

又由於第一個式子的$k$在取遍$[0,\frac{n}{2}-1]$時，$k+\frac{n}{2}$取遍了$[\frac{n}{2},n-1]$

因此咱們將原來的問題縮小了一半！

而縮小後的問題仍然知足原問題的性質，因此咱們能夠遞歸的去搞這件事情！

直到多項式僅剩一個常數項，這時候咱們直接返回就好啦

 

時間複雜度：

不難看出FFT是相似於線段樹同樣的分治算法。

所以它的時間複雜度爲$O(nlogn)$

 

快速傅里葉逆變換

不要覺得FFT到這裏就結束了。

咱們上面的討論是基於點值表示法的。

可是在日常的學習和研究中不多用點值表示法來表示一個多項式。

因此咱們要考慮如何把點值表示法轉換爲係數表示法，這個過程叫作傅里葉逆變換

 

$(y_0,y_1,y_2,\dots,y_{n-1})$爲$(a_0,a_1,a_2,\dots,a_{n-1})$的傅里葉變換（即點值表示）

設有另外一個向量$(c_0,c_1,c_2,\dots,c_{n-1})$知足

$$c_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i$$

即多項式$B(x)=y_0,y_1x,y_2x^2,\dots,y_{n-1}x^{n-1}$在$\omega_n^{0},\omega_n^{-1},\omega_n^{-2},\dots,\omega_{n-1}^{-(n-1)}$處的點值表示

emmmm又到推公式時間啦

$(c_0,c_1,c_2,\dots,c_{n-1})$知足
$$c_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i$$

$$=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i$$

$$=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^j)^i)(\omega_n^{-k})^i$$

$$=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^j)^i(\omega_n^{-k})^i)$$

$$=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^j)^i(\omega_n^{-k})^i$$

$$=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^{j-k})^i$$

$$=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i)$$

 

設$S(x)=\sum_{i=0}^{n-1}x^i$

將$\omega_n^k$代入得

$$S(\omega_n^k)=1+(\omega_n^k)+(\omega_n^k)^2+\dots(\omega_n^k)^{n-1}$$

當$k!=0$時

等式兩邊同乘$\omega_n^k$得

$$\omega_n^kS(\omega_n^k)=\omega_n^k+(\omega_n^k)^2+(\omega_n^k)^3+\dots(\omega_n^k)^{n}$$

兩式相減得

$$\omega_n^kS(\omega_n^k)-S(\omega_n^k)=(\omega_n^k)^{n}-1$$

$$S(\omega_n^k)=\frac{(\omega_n^k)^{n}-1}{\omega_n^k-1}$$

$$S(\omega_n^k)=\frac{(\omega_n^n)^{k}-1}{\omega_n^k-1}$$

$$S(\omega_n^k)=\frac{1-1}{\omega_n^k-1}$$

觀察這個式子，不難看出它分母不爲0，可是分子爲0

所以，當$K!=0$時，$S(\omega^{k}_{n})=0$

那當$k=0$時呢？

很顯然，$S(\omega^{0}_{n})=n$

 

繼續考慮剛剛的式子

$$c_k=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i)$$
當$j \neq k$時，值爲$0$
當$j=k$時，值爲$n$
所以，
$$c_k=na_k$$
$$a_k=\frac{c_k}{n}$$

這樣咱們就獲得點值與係數之間的表示啦

 

理論總結

至此，FFT的基礎理論部分就結束了。

咱們來小結一下FFT是怎麼成功實現的

 

首先，人們在用係數表示法研究多項式的時候遇阻

因而開始考慮可否用點值表示法優化這個東西。

而後根據複數的兩條性質（這個思惟跨度比較大）獲得了一種分治算法。

最後又推了一波公式，找到了點值表示法與係數表示法之間轉換關係。

 

emmmm

其實FFT的實現思路大概就是

係數表示法—>點值表示法—>係數表示法

引用一下遠航之曲大佬的圖

 

 

固然，在實現的過程當中還有不少技巧

咱們根據代碼來理解一下

 

遞歸實現

遞歸實現的方法比較簡單。

就是按找咱們上面說的過程，不斷把要求的序列分紅兩部分，再進行合併

在c++的STL中提供了現成的complex類，可是我不建議你們用，畢竟手寫也就那麼幾行，並且萬一某個毒瘤卡STL那豈不是很ＧＧ？

 

// luogu-judger-enable-o2
// luogu-judger-enable-o2
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int MAXN = 4 * 1e6 + 10;
inline int read() {
    char c = getchar(); int x = 0, f = 1;
    while (c < '0' || c > '9') {if (c == '-')f = -1; c = getchar();}
    while (c >= '0' && c <= '9') {x = x * 10 + c - '0'; c = getchar();}
    return x * f;
}
const double Pi = acos(-1.0);
struct complex {
    double x, y;
    complex (double xx = 0, double yy = 0) {x = xx, y = yy;}
} a[MAXN], b[MAXN];
complex operator + (complex a, complex b) { return complex(a.x + b.x , a.y + b.y);}
complex operator - (complex a, complex b) { return complex(a.x - b.x , a.y - b.y);}
complex operator * (complex a, complex b) { return complex(a.x * b.x - a.y * b.y , a.x * b.y + a.y * b.x);} //不懂的看複數的運算那部分
void fast_fast_tle(int limit, complex *a, int type) {
    if (limit == 1) return ; //只有一個常數項
    complex a1[limit >> 1], a2[limit >> 1];
    for (int i = 0; i <= limit; i += 2) //根據下標的奇偶性分類
        a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];
    fast_fast_tle(limit >> 1, a1, type);
    fast_fast_tle(limit >> 1, a2, type);
    complex Wn = complex(cos(2.0 * Pi / limit) , type * sin(2.0 * Pi / limit)), w = complex(1, 0);
    //Wn爲單位根，w表示冪
    for (int i = 0; i < (limit >> 1); i++, w = w * Wn) //這裏的w至關於公式中的k
        a[i] = a1[i] + w * a2[i],
               a[i + (limit >> 1)] = a1[i] - w * a2[i]; //利用單位根的性質，O(1)獲得另外一部分
}
int main() {
    int N = read(), M = read();
    for (int i = 0; i <= N; i++) a[i].x = read();
    for (int i = 0; i <= M; i++) b[i].x = read();
    int limit = 1; while (limit <= N + M) limit <<= 1;
    fast_fast_tle(limit, a, 1);
    fast_fast_tle(limit, b, 1);
    //後面的1表示要進行的變換是什麼類型
    //1表示從係數變爲點值
    //-1表示從點值變爲係數
    //至於爲何這樣是對的，能夠參考一下c向量的推導過程，
    for (int i = 0; i <= limit; i++)
        a[i] = a[i] * b[i];
    fast_fast_tle(limit, a, -1);
    for (int i = 0; i <= N + M; i++) printf("%d ", (int)(a[i].x / limit + 0.5)); //按照咱們推倒的公式，這裏還要除以n
    return 0;
}

遞歸版

 

update:遞歸版我本地是能夠AC的，只要開大數組就能夠了，可是交到洛谷上會WA0

 

 

woc?  臉好疼。。。。。。

咳咳。

速度什麼的纔不是關鍵呢？

關鍵是咱們AC不了啊啊啊

表着急，AC不了不表明我們的算法不對，只能說這種實現方法太low了

下面介紹一種更高效的方法

 

迭代實現

再盜一下那位大佬的圖

 

觀察一下原序列和反轉後的序列？

聰明的你有沒有看出什麼顯而易見的性質？

沒錯！

咱們須要求的序列實際是原序列下標的二進制反轉！

所以咱們對序列按照下標的奇偶性分類的過程實際上是沒有必要的

 

這樣咱們能夠$O(n)$的利用某種操做獲得咱們要求的序列，而後不斷向上合併就行了

 

// luogu-judger-enable-o2
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int MAXN = 1e7 + 10;
inline int read() {
    char c = getchar(); int x = 0, f = 1;
    while (c < '0' || c > '9') {if (c == '-')f = -1; c = getchar();}
    while (c >= '0' && c <= '9') {x = x * 10 + c - '0'; c = getchar();}
    return x * f;
}
const double Pi = acos(-1.0);
struct complex {
    double x, y;
    complex (double xx = 0, double yy = 0) {x = xx, y = yy;}
} a[MAXN], b[MAXN];
complex operator + (complex a, complex b) { return complex(a.x + b.x , a.y + b.y);}
complex operator - (complex a, complex b) { return complex(a.x - b.x , a.y - b.y);}
complex operator * (complex a, complex b) { return complex(a.x * b.x - a.y * b.y , a.x * b.y + a.y * b.x);} //不懂的看複數的運算那部分
int N, M;
int l, r[MAXN];
int limit = 1;
void fast_fast_tle(complex *A, int type) {
    for (int i = 0; i < limit; i++)
        if (i < r[i]) swap(A[i], A[r[i]]); //求出要迭代的序列
    for (int mid = 1; mid < limit; mid <<= 1) { //待合併區間的長度的一半
        complex Wn( cos(Pi / mid) , type * sin(Pi / mid) ); //單位根
        for (int R = mid << 1, j = 0; j < limit; j += R) { //R是區間的長度，j表示前已經到哪一個位置了
            complex w(1, 0); //冪
            for (int k = 0; k < mid; k++, w = w * Wn) { //枚舉左半部分
                complex x = A[j + k], y = w * A[j + mid + k]; //蝴蝶效應
                A[j + k] = x + y;
                A[j + mid + k] = x - y;
            }
        }
    }
}
int main() {
    int N = read(), M = read();
    for (int i = 0; i <= N; i++) a[i].x = read();
    for (int i = 0; i <= M; i++) b[i].x = read();
    while (limit <= N + M) limit <<= 1, l++;
    for (int i = 0; i < limit; i++)
        r[i] = ( r[i >> 1] >> 1 ) | ( (i & 1) << (l - 1) ) ;
    // 在原序列中 i 與 i/2 的關係是 ： i能夠看作是i/2的二進制上的每一位左移一位得來
    // 那麼在反轉後的數組中就須要右移一位，同時特殊處理一下奇數
    fast_fast_tle(a, 1);
    fast_fast_tle(b, 1);
    for (int i = 0; i <= limit; i++) a[i] = a[i] * b[i];
    fast_fast_tle(a, -1);
    for (int i = 0; i <= N + M; i++)
        printf("%d ", (int)(a[i].x / limit + 0.5));
    return 0;
}