第1章 引言
傅立葉變換（Fourier Transform）是由數學家傅立葉提出的一套對函式進行變換的方法，其主要分為連續傅立葉變換（Continuous Fourier Transform，CFT）和離散傅立葉變換（Discrete Fourier Transform，DFT）兩種，在本文中，我們只研究離散傅立葉變換。
離散傅立葉變換雖然在數學層面很有用，但其演演算法的時間複雜度較高，在演演算法層面並不實用。繼而，後續研究者又提出了快速傅立葉變換（Fast Fourier Transform，FFT）演演算法，這才徹底解決了問題。
那麼，離散傅立葉變換到底有什麼用呢？它的用途十分直白：用於計算多項式乘法。
多項式乘法早在中學數學中就已經學過，例如：
相信讀者一定有這樣的感覺：多項式的加減法很好算，只需要做幾次合併同類項就行了，這是因為多項式的加減法是時間複雜度為\(\varTheta(N)\)的演演算法。但多項式的乘法可就很難算了，上式中，想要計算一次3*3的多項式乘法，則需要進行9次乘法計算，和少量的合併同類項計算。也就是說，如果想要計算一次N*N的多項式乘法，則需要進行\(N^2\)次乘法計算，和少量的合併同類項計算，即：這種計算多項式乘法的演演算法的時間複雜度為\(\varTheta(N^2)\)。
快速傅立葉變換正是這樣一個演演算法：其能夠突破上述演演算法的\(\varTheta(N^2)\)時間複雜度，將多項式乘法的時間複雜度最佳化至\(\varTheta(NlogN)\)。
第2章 多項式的係數表達與點值表達
想要研究多項式，就需要先把多項式寫出來。在本章中，我們研究多項式的兩種表達方式：係數表達與點值表達。
2.1 係數表達
多項式的係數表達我們都非常熟悉：指的就是上面的\(\left( 1+2x+3x^2 \right)\)和\(\left( 4+5x+6x^2 \right)\)這種形式，透過寫出多項式中每一項的係數，從而表達出一個多項式是什麼樣的。
這種寫法還可以再省略一些：由於每個係數後面的\(x^n\)寫不寫出來都一樣，所以可以只寫出每一項的係數，並構成一個向量：
此時，多項式乘法就有了一種新的表示：
這是一種全新的向量間的乘法運算，稱為：卷積（Convolution），用符號\(\otimes\)表示。
2.2 點值表達
多項式的係數表達很好理解，那麼，點值表達又是什麼呢？
我們都知道：兩點確定一條直線，而直線是一個包含常數項和一次項的多項式；所以，我們也可以說：兩點確定一個一次多項式。那麼，三點能不能確定一個二次多項式呢？四點又能不能確定一個三次多項式呢？更多的點呢？
答案是肯定的，請看如下引理：
引理1（點值表達的唯一性）：對於任意的n個點構成的點集：\(\left\{ \left( x_0, y_0 \right) , \left( x_1, y_1 \right) , \left( x_2, y_2 \right) , ..., \left( x_{n-1}, y_{n-1} \right) \right\}\)，如果\(x_0\ne x_1\ne x_2\ne ...\ne x_{n-1}\)，則該點集能夠唯一確定一個\(n-1\)次多項式。
證明：
也就是說，我們不僅可以用n個係數唯一的表達一個n-1次多項式，還可以用這個多項式上的n個點唯一表達。
現在，需要考慮三個問題：
- 已知一個多項式的係數表達，怎麼得到其點值表達？
- 已知一個多項式的點值表達，怎麼得到其係數表達？
- 點值表達存在的意義是什麼？
第一個問題很簡單，我們已經知道，想要得到一個多項式的點值表達，需要滿足兩個條件：
- 點的數量要夠，需要找n個點
- 每個點的x互不相同
這兩個條件可太好滿足了，專門挑幾個最好算的x帶入計算，點不就來了嗎？例如：對於多項式\(1+2x+3x^2\)，需要找3個點，那我們就選3個最好算的點：0、1、-1，帶進去，點就有了：\(\left\{ \left( 0, 1 \right) , \left( 1, 6 \right) , \left( -1, 2 \right) \right\}\)，這就是多項式\(1+2x+3x^2\)的點值表達。當然，還可以挑別的3個點，算出來的點集都能唯一確定\(1+2x+3x^2\)這個多項式。
第二個問題也不難解決，不是有很多點嗎，那就用這些點解上面證明過程中的那個線性方程組，解出來的\(a_0, a_1, a_2, ..., a_{n-1}\)就是係數表達了。
第三個問題很重要。我們之所以要研究點值表達，顯然是為瞭解決多項式的乘法問題的。說到這，點值表達存在的意義就出現了：點值表達下的多項式四則運算，全都是時間複雜度為\(\varTheta(N)\)的演演算法。這句話可能需要想一下才能明白，這是因為：兩個多項式做四則運算，其實就是在這兩個多項式的每一對y上都做一遍四則運算。那麼，如果我們選擇相同的一組x帶入兩個多項式，得到兩組點值表達，多項式乘法就可以變成這兩組點之間做一遍乘法了。舉個例子：比如我們對\(\left( 1+2x+3x^2 \right)\)和\(\left( 4+5x+6x^2 \right)\)這兩個多項式都選0、1、-1帶入，就能得到兩組x相同的點值表達：\(\left\{ \left( 0, 1 \right) , \left( 1, 6 \right) , \left( -1, 2 \right) \right\}\)和\(\left\{ \left( 0, 4 \right) , \left( 1, 15 \right) , \left( -1, 5 \right) \right\}\)；此時，多項式乘法就可以用這三個點的y進行計算，結果是：\(\left\{ \left( 0, 4 \right) , \left\{ 1, 90 \right\} , \left( -1, 10 \right) \right\}\)。
事實上，我們在這裡犯了一個錯誤：兩個二次多項式做乘法，結果應當是一個四次多項式，而四次多項式需要用五個點才能唯一表示，而上面只有三個點，顯然是不夠的。這個問題提醒了我們：在計算多項式乘法的時候，不能只看現在這個多項式是幾次的，還應當看乘法的結果多項式是幾次的，後者的次數才能決定一開始要取幾個點。
三個問題都討論完了，但是除了第三個問題的結論很有吸引力以外，前兩個問題的結論實在是平平無奇。
先看從係數表達轉點值表達的過程，想做這件事，就要選n個點，每個點依次帶入多項式算一遍，而多項式裡面全是各種高次冪，時間複雜度早已不可接受了。如果使用霍納法則，將計算一個多項式的時間複雜度最佳化到\(\varTheta(N)\)，那最終也是一個時間複雜度為\(\varTheta(N^2)\)的演演算法。
而從點值表達轉系數表達就更難算了，由於需要解矩陣方程，考慮線性代數中的LU分解演演算法，其時間複雜度是\(O(N^3)\)。拉格朗日插值法則更快一些，其基於以下公式：
不難看出，這個公式的時間複雜度為\(\varTheta(N^2)\)。
看來，我們遇到了木桶效應：雖然在點值表達下，多項式的乘法變得非常好算，但兩種表達方式間的來回轉換，全都是時間複雜度至少也要\(\varTheta(N^2)\)的演演算法。還不如用最普通的那個演演算法呢。
所以，現在的目標是：找到一對演演算法，能夠更快的在多項式的係數表達和點值表達之間進行轉換。
第3章 尤拉公式
由於傅立葉變換需要使用複數的相關知識，所以這一章中，我們研究尤拉公式及其推論。
尤拉公式是一個非常有名的公式，其將複數域上的指數函式和三角函式聯絡在了一起。讓我們從\(e^{i\theta}\)的麥克勞林級數展開開始：
這就是著名的尤拉公式：
引理2（尤拉公式）：
將\(\theta=\pi\)帶入尤拉公式，就能得到著名的尤拉恆等式：
引理3（尤拉恆等式）：
尤拉公式和尤拉恆等式在後續的研究中十分重要，讀者應熟練掌握。
第4章 n次單位複數根
這一章中，我們研究n次單位複數根。什麼是n次單位複數根呢？其指的是以下n次方程的根：
在實數域中，不管n是多少，這個方程都最多隻有兩個根：1和-1。但是，由代數學基本定理：任何一個n次方程都有且僅有n個根。那麼，對於\(n>2\)的這些高次方程來說，剩下的根去哪了呢？顯然，這些根都是虛數。
根據尤拉公式，我們可以湊出n次單位複數根的一般形式：
引理4（n次單位複數根）：
證明：
由複數在複平面上的極座標表示可知：各個n次單位複數根是複平面上單位圓的各個n等分點，故其具有一些特殊的性質。請看：
引理5（折半引理）：\(\omega _{n}^{k}=\omega _{n/2}^{k/2}\)
證明：
引理6：\(\omega _{n}^{n}=1\)
證明：
引理7：\(\omega _{n}^{n/2}=-1\)
證明：
引理8：\(\omega _{n}^{-k}=\overline{\omega _{n}^{k}}\)
證明：
n次單位複數根及其若干引理的用途將在後續章節展開。
第5章 離散傅立葉變換與離散傅立葉逆變換
上一章中，我們研究了n次單位複數根及其性質。那麼，這些根有什麼用呢？傅立葉變換現在正式登場：
使用所有的n次單位複數根將一個n-1次多項式從係數表達轉為點值表達的過程，就稱為離散傅立葉變換（Discrete Fourier Transform，DFT）。
原來，上一章中研究的這些n次單位複數根是用來帶入的。如果選用這一組根作為x帶入多項式，從而得到其點值表達，這個過程就稱為離散傅立葉變換。
那麼，離散傅立葉逆變換又是什麼呢？顧名思義，其是一個變換回來的過程：
將離散傅立葉變換得到的點值表達轉為係數表達的過程，就稱為離散傅立葉逆變換（Inverse Discrete Fourier Transform，IDFT）。
我們已經知道，想要將一個多項式從係數表達轉為點值表達，隨便選一組點都是可以的。那為什麼選n次單位複數根這組點進行轉換，就有了傅立葉變換和傅立葉逆變換這兩個專業術語呢？我們不禁猜想：n次單位複數根可能是一組非常特殊的點，將這組點帶入多項式時，是可以簡化計算的；此外，從這些點得到的點值表達向係數表達轉換時，也是可以簡化計算的。這樣一來，雙向的轉換過程就都能得到最佳化了。
事實上，這兩個猜想都是正確的。其演演算法分別被稱為：快速傅立葉變換（Fast Fourier Transform，FFT）和快速傅立葉逆變換（Inverse Fast Fourier Transform，IFFT）。在下面的兩章中，我們分別研究這兩種演演算法。
第6章 快速傅立葉變換
6.1 快速傅立葉變換的數學原理
上一章中，我們已經知道，如果將n次單位複數根帶入多項式，進行係數表達到點值表達的轉換，是可以簡化計算的，這樣的演演算法被稱為快速傅立葉變換。這一章中，我們具體研究這一演演算法。
既然是一種最佳化演演算法，那就不能一上來就帶入點進行計算，我們需要先做一些準備：
- 因為快速傅立葉變換是一個嚴格二分的演演算法，所以需要將多項式的項數補齊至2的整數次冪。什麼意思呢？比如想對一個二次多項式進行快速傅立葉變換，由於二次多項式有隻有三項，所以需要用係數0再補一項，將其補至四項；例如：\(1+2x+3x^2\)就需要補成\(1+2x+3x^2+0x^3\)。同理，具有5、6、7項的多項式都需要用係數0補齊至8項；以此類推
- 將多項式按奇偶順序重排並整理：
- 將多項式分解為兩個子多項式：
至此，準備工作就完成了。
現在，我們將n次單位複數根，即\(\omega _{n}^{k}\)帶入上式：
此時，\(\omega _{n}^{2k}\)可以使用折半引理，變成\(\omega _{n/2}^{k}\)。但需要小心：有一半的k是滿足\(k\geqslant n/2\)的，這些比較大的k如果直接使用折半引理的話，會變得不好處理，所以，需要分兩種情況討論：
對比這兩種情況得到的兩個結果，可以發現：這兩個結果唯一的差別就是一個正負號。也就是說，每當在\(k<n/2\)這半邊算出一個點，就可以立即在\(k+n/2\)處再算出一個點。這樣一來，計算量直接減少一半。這還沒完，在計算\(A^{\left[ 0 \right]}\left( \omega _{n/2}^{k} \right)\)和\(A^{\left[ 1 \right]}\left( \omega _{n/2}^{k} \right)\)的時候，這兩個多項式不僅長度減半，而且可以繼續用這個性質，計算量繼續減少一半。這個過程什麼時候停下來呢？不難想象：當一個只有兩項的多項式被拆成兩個只有一項的多項式時，這個過程就停了，因為只有一項的多項式就是一個常數，其值不僅不需要計算，而且與自變數的取值無關。
上述演演算法的時間複雜度可由主定理證明：為\(\varTheta(NlogN)\)。這就是快速傅立葉變換的數學原理。
6.2 一個手工計算快速傅立葉變換的例項
上一節中，我們研究了快速傅立葉變換的數學原理，但如果只看公式的話，讀者可能仍然不理解這個演演算法是怎麼進行的。所以，這一節中，我們就透過一個例子，來真正計算一次快速傅立葉變換。
例：計算多項式\(A\left( x \right) =1+2x+3x^2+4x^3\)的快速傅立葉變換。
按照快速傅立葉變換的演演算法流程，我們首先要做的是：將這個多項式的項數補齊至2的整數次冪。由於這個多項式的項數是4，其已經是2的整數次冪了，所以，不需要補齊。
接下來，需要按奇偶順序將多項式重排並整理：
此時，我們已經將一個四項的多項式分解為兩個兩項的多項式了。但這還不夠，我們還需要將這兩個多項式繼續分解為四個只有一項的多項式：
此時，\(A^{\left[ 0 \right] \left[ 0 \right]}\left( x^2 \right) \text{，}A^{\left[ 0 \right] \left[ 1 \right]}\left( x^2 \right) \text{，}A^{\left[ 1 \right] \left[ 0 \right]}\left( x^2 \right) \text{，}A^{\left[ 1 \right] \left[ 1 \right]}\left( x^2 \right)\)這四個多項式都是隻包含常數項的多項式，其值均與輸入無關，而分別恆等於其常數項的值。
至此，多項式已經完成分解。我們開始將n次單位複數根帶入，並向前倒推多項式的值：
同理：
接下來，讓我們回到最開始的目標：
至此，我們就完成了一次快速傅立葉變換的計算。讀者可以將4個自變數依次帶入多項式，來驗證結果的正確性。
6.3 過程更簡略的手工計算例項
上一節中，雖然我們花了大量的篇幅來演示一次快速傅立葉變換是怎麼計算的，但實際上讀者可以發現：其中的大多數過程都只是為了便於讀者理解而寫的，當熟練掌握後，完全可以只保留真正需要計算的部分，而這部分的計算量是非常少的。從這裡就能看出，快速傅立葉變換對多項式求值帶來的最佳化。
這一節中，我們再重新算一次上面這個多項式的快速傅立葉變換；但這一次，我們省去一切不必要的過程，只保留真正需要計算的部分，看看是什麼體驗。
首先，我們將多項式裡面的x全部省去，只留下係數：
1 2 3 4
這裡的1 2 3 4表示的就是上面的\(A\left( x \right) =1+2x+3x^2+4x^3\)
接下來，考慮對項數補齊。由於這個多項式的項數是4，剛好是2的整數次冪，所以不需要補齊。
接下來，將係數按奇偶順序分組：
1 2 3 4
1 3 | 2 4
1 | 3 | 2 | 4
這裡的1 | 3 | 2 | 4依次對應著上面的\(A^{\left[ 0 \right] \left[ 0 \right]}\left( x^2 \right) \text{，}A^{\left[ 0 \right] \left[ 1 \right]}\left( x^2 \right) \text{，}A^{\left[ 1 \right] \left[ 0 \right]}\left( x^2 \right) \text{，}A^{\left[ 1 \right] \left[ 1 \right]}\left( x^2 \right)\)，由於這四個多項式的值與輸入無關，所以其值分別就是1、3、2、4。
接下來，我們從第二行的係數開始，向上倒推多項式的值。第一次倒推以連續的兩個係數為一組，每一組中，相鄰的兩個係數之間做一對計算，需要用到係數：\(\omega _{2}^{0}=1\)：
1 3 | 2 4
1 + 1 * 3 = 4
1 - 1 * 3 = -2
2 + 1 * 4 = 6
2 - 1 * 4 = -2
4 -2 6 -2
第二次倒推以連續的四個係數為一組，每一組中，第n個係數和第n+2個係數之間做一對計算，需要用到係數：\(\omega _{4}^{0}=1\text{；}\omega _{4}^{1}=i\)：
4 -2 6 -2
4 + 1 * 6 = 10
4 - 1 * 6 = -2
-2 + i * -2 = -2-2i
-2 - i * -2 = -2+2i
10 -2-2i -2 -2+2i
至此，快速傅立葉變換就計算完成了（n次單位複數根的值不需要算出）。由此可見，一旦讀者熟練掌握了快速傅立葉變換的計算原理，就可以使用這種非常簡潔的計算過程進行快速傅立葉變換的計算了。
第7章 快速傅立葉逆變換
上一章中，我們研究瞭如何在\(\varTheta(NlogN)\)的時間複雜度下進行快速傅立葉變換，而快速傅立葉變換一旦完成，就可以進行一次時間複雜度為\(\varTheta(N)\)的多項式乘法計算，再之後，我們就需要將多項式的點值表達轉為係數表達了，即進行快速傅立葉逆變換。
7.1 快速傅立葉逆變換的數學原理
對於多項式：
我們現在已經知道其全部的快速傅立葉變換結果：
此時，我們需要構造一個新的多項式，這個多項式的係數由\(y_i\)給出，並將單位複數根全部取倒數後帶入：
這樣做有什麼用呢？接下來，將\(y_i\)展開：
此時，需要分兩種情況討論：
透過這一系列操作，我們找到了一種計算\(a_k\)的演演算法，基於這個演演算法，就能從一個多項式的點值表達轉為係數表達了。所以，現在只剩下最後一個問題：\(z_k\)表示的是將單位複數根的倒數（而不是單位複數根）帶入一個多項式，這並不是快速傅立葉變換的標準做法。那麼此時，快速傅立葉變換還能使用嗎？
為了研究這個問題，讓我們回到快速傅立葉變換推導過程的起點：
回顧一下：在快速傅立葉變換的推導過程中，我們是將一半的\(\omega _{n}^{k}\)以及另一半的\(\omega _{n}^{k^{\prime}+n/2}\)帶入到\(A\left( x \right)\)中的，得到的結論是：
即：只需要計算\(A^{\left[ 0 \right]}\left( \omega _{n/2}^{k} \right)\)和\(\omega _{n}^{k}A^{\left[ 1 \right]}\left( \omega _{n/2}^{k} \right)\)就行了。
那麼，如果將單位複數根的倒數帶入，又會怎麼樣呢？還是分兩種情況討論：離0比較近的一半，和離0比較遠的一半：
可以發現：如果將單位複數根的倒數帶入，那麼在進行快速傅立葉變換的過程中，只有一個係數的差別（\(A^{\left[ 0 \right]}\)和\(A^{\left[ 1 \right]}\)括號裡面的差別可以忽略，因為分解後的多項式最終將與輸入無關）。也就是說，只需要對快速傅立葉變換的計算過程做兩處微小的改動，我們就可以再一次利用快速傅立葉變換去計算所有\(z_k\)的值，從而計算出所有\(a_k\)的值了。這兩處改動分別為：
- 將\(\omega _{n}^{k}\)換成\(\overline{\omega _{n}^{k}}\)
- 將每個快速傅立葉變換的結果除以n
這就是快速傅立葉逆變換的數學原理。其時間複雜度與快速傅立葉變換一致，也是\(\varTheta(NlogN)\)。
7.2 一個手工計算快速傅立葉逆變換的例項
上一節中，我們研究了快速傅立葉逆變換的數學原理。這一節中，我們就使用上一章得到的快速傅立葉變換的計算結果10 -2-2i -2 -2+2i，進行一次快速傅立葉逆變換的手工計算，如果計算結果能夠回到1 2 3 4，就說明我們的計算是正確的。
首先，考慮是否需要對項數補齊。由於項數是4，所以不需要補齊。
接下來，對係數進行分組：
10 -2-2i -2 -2+2i
10 -2 | -2-2i -2+2i
接下來，進行第一次倒推，需要用到係數：\(\overline{\omega _{2}^{0}}=1\)：
10 -2 | -2-2i -2+2i
10 + 1 * -2 = 8
10 - 1 * -2 = 12
-2-2i + 1 * -2+2i = -4
-2-2i - 1 * -2+2i = -4i
8 12 -4 -4i
接下來，進行第二次倒推，需要用到係數：\(\overline{\omega _{4}^{0}}=1\text{；}\overline{\omega _{4}^{1}}=-i\)：
8 12 -4 -4i
8 + 1 * -4 = 4
8 - 1 * -4 = 12
12 + -i * -4i = 8
12 - -i * -4i = 16
4 8 12 16
上面得到的4 8 12 16，分別就是\(z_0, z_1, z_2, z_3\)。最後，我們將其都除以\(n=4\)，就能得到\(a_0, a_1, a_2, a_3\)了：
4 / 4 = 1
8 / 4 = 2
12 / 4 = 3
16 / 4 = 4
可見，計算結果完全符合預期。
第8章 快速傅立葉變換的實現
這一章中，我們研究快速傅立葉變換以及快速傅立葉逆變換的實現，有了這兩種變換，就能實現出時間複雜度為\(\varTheta(NlogN)\)的卷積演演算法。
8.1 對齊至2的整數次冪的演演算法
快速傅立葉變換的第一步是將多項式的項數對齊至2的整數次冪，所以，我們需要根據輸入的項數，來找到需要對齊到的項數。這一需求的樸素演演算法是使用一個迴圈，並使用一個從1開始，不斷自乘2的數字和輸入項數作比較，直至這個數字已經大於等於輸入項數時，演演算法終止。這個演演算法很簡單，讀者可以自行嘗試。
這裡給出一種更為高效的演演算法：
unsigned __nextPow2(unsigned N)
{
N--;
N |= N >> 1;
N |= N >> 2;
N |= N >> 4;
N |= N >> 8;
N |= N >> 16;
return N + 1;
}
這個演演算法不難理解，其要點在於：
- 如果N不是2的整數次冪，那麼，就將N從最高位的1開始，到最低位之間的所有位都變成1；然後，將這個全是1的數字再加1，這些1就都會變成0，並且一個新的1將出現在原最高位的更高一位上，這正是我們需要的數字。例如
0b101，我們希望將其變成0b111，再加1，就得到了0b1000，這就是我們需要的數字。那麼，具體要怎麼操作，才能將0b101變成0b111呢？我們可以從N的最高位的那個1開始，將其右移1位後與N位或，此時，N的最高兩位就一定都是1了；接下來，將N右移兩位後與N位或，使N的最高4位都變成1；以此類推：接下來使N的最高8、16、32位都變成1。讀者在理解這段話時要清楚：這裡所說的最高n位，都是在N足夠大，確實有這麼多位的前提下才成立，否則，N就會因為過多的右移而位或到一個0 - 如果N已經是2的整數次冪，那麼，演演算法直接返回N就行了。但是，判斷一個數字是不是2的整數次冪需要額外的代價，且很明顯，這個判定的失敗率是很高的，因為絕大多數的整數都不是2的整數次冪。所以，乾脆就不要判定這件事了，而是將N減去1。如果N是2的整數次冪，減去1後就會丟失其最高位的1，並變成一個全是1的數字，在經過多次（無用的）位或運算後，又被加上1，回到了原值。例如
0b1000，減去1後會變成0b111，其丟失了原數字最高位的1，最終，0b111又會因為加1而回到0b1000並返回。另一方面，如果N並不是2的整數次冪，那就說明N除了最高位的1以外，在低位還有1，這樣一來，減去1就不會使N丟失最高位的1，所以，其結果不受影響
8.2 位逆序演演算法
快速傅立葉變換的第二步是對係數進行分組，分組操作的樸素實現和前面的手工計算過程是一致的，讀者可以自行嘗試。這裡給出的是一種更為高效的演演算法。
仔細觀察分組前後，各個係數的索引值的二進位製表示，這裡以8個係數為例：
000 001 010 011 100 101 110 111 // 分組前
000 100 010 110 001 101 011 111 // 分組後
不難發現：分組前後的每一對索引值都是位逆序的。
這就意味著，對於輸入的每一個係數，我們都可以立即知道這個係數在分組後被放在哪裡了：只需要將係數的索引值進行位逆序即可。
那麼，怎麼實現位逆序呢？樸素的演演算法是：透過一個迴圈，將待轉換的數字不斷右移1位，同時將轉換後的數字不斷左移1位，並將兩個數字的最低位對接即可。讀者可以自行嘗試。
這裡給出的是一種更為高效的演演算法：
unsigned __bitReverse(unsigned N, unsigned bitWidth)
{
N = ((0xaaaaaaaa & N) >> 1) | ((0x55555555 & N) << 1);
N = ((0xcccccccc & N) >> 2) | ((0x33333333 & N) << 2);
N = ((0xf0f0f0f0 & N) >> 4) | ((0x0f0f0f0f & N) << 4);
N = ((0xff00ff00 & N) >> 8) | ((0x00ff00ff & N) << 8);
N = ((N >> 16) | (N << 16)) >> (32 - bitWidth);
return N;
}
這個演演算法不難理解，其要點在於：
0xaaaaaaaa是形如0b1010...的位掩碼，這個位掩碼會保留N的所有奇數位；而0x55555555是形如0b0101的位掩碼，這個位掩碼會保留N的所有偶數位；將二者的掩碼結果一個左移，一個右移，最後再位或到一起，就能使N的所有相鄰位發生交換- 類似的，
0xcccccccc是形如0b1100...的位掩碼；而0x33333333是形如0b0011的位掩碼；在這兩個位掩碼，以及後續的左右移位和位或的作用下，N會以每兩位為一組發生交換。以此類推，N又會以每4、8位為一組進行交換 - 最後，N需要以每16位為一組，完成最後一次交換，此時就不需要位掩碼了，直接交換即可。至此，N的所有位完成了逆序
- 上述演演算法完成的是32位無符號整數的位逆序，而實際輸入的數字很可能並沒有這麼多位。例如：
0b001的位逆序應該是0b100，而按照上述演演算法，最終的結果是：0b100...(後面還有29個0)，這不是我們需要的。所以，最後還需要做一次右移，將多餘的0去掉
8.3 快速傅立葉變換與快速傅立葉逆變換的實現
在前面的章節中我們已經知道，快速傅立葉變換和快速傅立葉逆變換的計算過程只有兩處微小的不同：
- 快速傅立葉變換使用的係數是\(\pm \omega _{n}^{k}\)，而快速傅立葉逆變換使用的係數是\(\pm \overline{\omega _{n}^{k}}\)
- 快速傅立葉逆變換需要在最後對所有的結果除以n
我們可以使用一個要麼是1，要麼是-1的數字來同時區別這兩處不同。請看：
vector<complex<double>> __FFT(const vector<complex<double>> &coefList, double conjNum)
{
vector<complex<double>> FFTList(coefList.size());
unsigned bitWidth = __builtin_ctz(FFTList.size());
for (unsigned idx = 0; idx < FFTList.size(); idx++)
{
FFTList[idx] = coefList[__bitReverse(idx, bitWidth)];
}
for (unsigned N = 2; N <= FFTList.size(); N *= 2)
{
for (unsigned startIdx = 0; startIdx < FFTList.size(); startIdx += N)
{
complex<double> curOmega(1.);
complex<double> mulOmega(cos(2 * M_PI / N), conjNum * sin(2 * M_PI / N));
for (unsigned leftIdx = startIdx, rightIdx = startIdx + N / 2; leftIdx < startIdx + N / 2; leftIdx++, rightIdx++)
{
auto leftNum = FFTList[leftIdx] + curOmega * FFTList[rightIdx];
auto rightNum = FFTList[leftIdx] - curOmega * FFTList[rightIdx];
FFTList[leftIdx] = leftNum;
FFTList[rightIdx] = rightNum;
curOmega *= mulOmega;
}
}
}
if (conjNum == -1.)
{
for (auto &FFTNum: FFTList)
{
FFTNum /= FFTList.size();
}
}
return FFTList;
}
__FFT函式用於計算快速傅立葉變換以及快速傅立葉逆變換。當conjNum = 1.時，其處於快速傅立葉變換模式；而當conjNum = -1.時，其處於快速傅立葉逆變換模式（由於這個函式不作為對外介面，所以沒有對conjNum使用布林值或列舉變數等程式設計手段限制其他錯誤的值，讀者如果對此感到介意，可以自行實現一個更嚴謹的介面）。
形參方面，coefList為係數列表，所有的係數已經由主調函式從double型別轉為了complex<double>型別，並已經進行了對齊處理；conjNum已在上文中說明，其只會傳入1.或-1.。
函式中，首先進行的是係數的分組操作。在進行這一操作之前，我們需要知道位逆序所需要的位寬，這是由GCC內建函式__builtin_ctz完成的，其返回輸入數字從最低位到第一個1之間的0的數量。分組操作透過迴圈進行，其將coefList中的係數重排至FFTList列表中。
接下來的程式碼是一個三重迴圈。
第一重迴圈用於遍歷N的取值，N從2開始，以不斷自乘2的方式遞增，直至與多項式的項數一致時終止。
第二重迴圈用於遍歷分組，startIdx存放的是當前分組的起始索引值；而分組的長度（決定了startIdx的迴圈增量）是恆等於N的。比如，第一次倒推時，以兩個數字為一組；第二次倒推時，以四個數字為一組；以此類推。
當確定了N以及當前分組後，就可以對分組內的每一對係數進行計算了。在計算過程中，\(\omega _{n}^{k}\)或\(\overline{\omega _{n}^{k}}\)需要伴隨迴圈而變化，具體來說，每計算一對係數，當前的\(\omega _{n}^{k}\)或\(\overline{\omega _{n}^{k}}\)就需要再乘一次\(\omega _{n}^{1}\)或\(\overline{\omega _{n}^{1}}\)；而\(\omega _{n}^{k}\)或\(\overline{\omega _{n}^{k}}\)的初始值為\(\omega _{n}^{0}\)或\(\overline{\omega _{n}^{0}}\)，均為1。程式碼方面，curOmega變數用於存放\(\omega _{n}^{k}\)或\(\overline{\omega _{n}^{k}}\)的當前值，其被初始化為1；而mulOmega變數用於存放\(\omega _{n}^{1}\)或\(\overline{\omega _{n}^{1}}\)，其用於在每一對係數計算完成後，將curOmega自乘一次mulOmega。
mulOmega的值基於尤拉公式：
conjNum用於控制上式中的正負號。
第三重迴圈用於計算當前分組內的每一對係數。分組的長度為N，將其分為左右兩半，左半邊的索引值由leftIdx維護，初始化為startIdx，即當前分組的起始索引值；右半邊的索引值由rightIdx維護，初始化為startIdx + N / 2，即當前分組右半邊的第一個索引值；這兩個索引值同步向前遞增，從而訪問到分組內的每一對係數。在迴圈體中，我們同時計算並更新一對係數，然後更新curOmega。
在這個函式的最後，實現的是快速傅立葉逆變換所需的額外操作：將每個係數都除以n。
8.4 卷積的實現
卷積的實現是前面所有準備工作的彙總。讓我們先梳理一下實現思路：
- 形參方面，卷積的輸入是兩個不保證等長的
vector<double> - 在進行快速傅立葉變換之前，我們需要先準備好足夠多的點來表示結果多項式，即：需要將輸入的兩個係數列表都用0擴充到足夠的長度。那麼，需要多少個係數呢？這裡需要做一個簡單的計算：一個長度為\(n\)的係數列表，表示的是一個\(n-1\)次多項式；而另一個長度為\(m\)的係數列表，表示的是一個\(m-1\)次多項式；這兩個多項式相乘的結果是一個\(n+m-2\)次多項式；而這樣的多項式一共有\(n+m-1\)項。此外，根據快速傅立葉變換的要求，係數列表的長度必須是2的整數次冪，所以，我們還需要將\(n+m-1\)這個數字對齊到2的整數次冪，作為兩個係數列表擴充後的長度
- 快速傅立葉變換需要的係數列表是
vector<complex<double>>型別的，而輸入的係數列表是vector<double>型別的，需要進行轉換 - 當兩個係數列表都準備好後，進行兩次快速傅立葉變換，將兩個多項式從係數表達轉為點值表達；然後，透過一個迴圈進行點值表達下的多項式乘法；最後，再進行一次快速傅立葉逆變換，將點值表達轉為係數表達
- 快速傅立葉逆變換的輸出是
vector<complex<double>>型別的係數列表，而我們最終需要的是vector<double>型別的係數列表，需要進行轉換。此外，還需要捨去由於對齊到2產生的係數擴充
下面請看實現：
vector<double> calcVectorConvolution(const vector<double> &leftCoefList, const vector<double> &rightCoefList)
{
unsigned coefSize = leftCoefList.size() + rightCoefList.size() - 1;
unsigned alignSize = __nextPow2(coefSize);
vector<complex<double>> leftAlignCoefList(alignSize);
vector<complex<double>> rightAlignCoefList(alignSize);
for (unsigned idx = 0; idx < leftCoefList.size(); idx++)
{
leftAlignCoefList[idx].real(leftCoefList[idx]);
}
for (unsigned idx = 0; idx < rightCoefList.size(); idx++)
{
rightAlignCoefList[idx].real(rightCoefList[idx]);
}
auto leftFFTList = __FFT(leftAlignCoefList, 1.);
auto rightFFTList = __FFT(rightAlignCoefList, 1.);
for (unsigned idx = 0; idx < leftFFTList.size(); idx++)
{
leftFFTList[idx] *= rightFFTList[idx];
}
auto resFFTList = __FFT(leftFFTList, -1.);
vector<double> resCoefList(coefSize);
for (unsigned idx = 0; idx < resCoefList.size(); idx++)
{
resCoefList[idx] = resFFTList[idx].real();
}
return resCoefList;
}
形參方面，leftCoefList和rightCoefList是兩個多項式的係數表達。
coefSize用於存放結果多項式的項數，其計算公式已經由上文討論過；alignSize用於存放將coefSize對齊到2的整數次冪後的多項式的項數，其決定了快速傅立葉變換需要的列表長度。
接下來，使用alignSize作為長度生成leftAlignCoefList和rightAlignCoefList，這兩個係數列表用於快速傅立葉變換；並將leftCoefList和rightCoefList中的係數分別放入這兩個列表的前面部分。
接下來，進行兩次快速傅立葉變換，將leftAlignCoefList和rightAlignCoefList從係數表達轉為點值表達，存放在leftFFTList和rightFFTList中；然後，使用一個迴圈進行點值表達下的多項式乘法，將rightFFTList乘入leftFFTList中；最後，再進行一次快速傅立葉逆變換，將leftFFTList從點值表達轉為係數表達，存放在resFFTList中。
現在，resFFTList中存放的是alignSize個complex<double>型別的係數，而我們需要的是這個列表中前coefSize個double型別的係數，所以，函式的最後一段用於提取這部分系數並返回。讀者可以自行驗證：在resFFTList中，除了我們提取出的部分外，其餘部分無論是實部還是虛部，都是0。
第9章 討論
本文中，對時間複雜度的描述多次使用了\(\varTheta\)符號而非\(O\)符號，讀者應予以關注。
第三章中研究的尤拉公式的推導過程不能作為其證明過程。這是因為，麥克勞林級數展開需要求出函式的高階導數，而複數域下的指數函式以及三角函式的導數均依賴於尤拉公式，從而造成迴圈論證。尤拉公式的嚴格證明超出了本文的範圍。
快速傅立葉變換的樸素實現基於：\(A\left( \omega _{n}^{k} \right) =A^{\left[ 0 \right]}\left( \omega _{n/2}^{k} \right) \pm \omega _{n}^{k}A^{\left[ 1 \right]}\left( \omega _{n/2}^{k} \right)\)這一結論；這是一個遞迴版本的演演算法，讀者可以自行嘗試。
已經存在不要求n為2的整數次冪的快速傅立葉變換演演算法，但其超出了本文的範圍。
本文中使用的"每次計算一對係數"的操作，在相關書籍和文獻中被稱為蝴蝶操作（Butterfly Operation；《演演算法導論》中將此術語拼寫為Bufferfly Operation，似為勘誤）；且\(\pm \omega _{n}^{k}\)和\(\pm \overline{\omega _{n}^{k}}\)被稱為旋轉因子（Twiddle Factor）。本文作者認為這兩個術語不夠生動形象，故未在正文中引入。
快速傅立葉變換在自然科學，電腦科學等諸多領域都有著廣泛的應用，希望本文能夠為讀者提供幫助。
櫻雨樓
2023.2