FFT·快速傅立葉變換

---

title: 【學習筆記】從單位根到FFT
date: 2019-02-19 11:26:08
tags:
- 多項式基礎
top: 6009
categories:
- 學習筆記
- 多項式

---

青春的回憶啊…

Preface

這篇文章初寫於 $ 7/1/2018 $ ，是在陪同好友 $ yjk $ 與 $ wx $ 以及學長 $ rqy $ 一起去參加省隊集訓時寫的。今天突然來了興致，打算重新復習一遍 FFT 並且寫 MTT ，於是便有了這篇文章。

其實一開始我是不情願把這篇文章搬到這兒來的——這好像是一個時代的縮影，那個時代的orchidany特別喜歡扮演老師，每天仿佛來到奧賽室只是為了“為人師”的：心性浮躁，學習功利。但現在我則是想沉下心來，認真做學問。

但無論如何，我希望這篇原本冗雜繁長的文章可以更短、更新穎、從更高的角度審視一些問題。

Convolution

卷積 $ \boldsymbol{(Convolution)} $ ，准確來說是一種通過兩個函數 $ \boldsymbol f $ 和 $ \boldsymbol g $ 生成第三個函數的一種數學算子.

而廣義上其定義為：

\[\boldsymbol{h(x) = \int _{- \infty} ^{\infty}g(\tau) \cdot f(x - \tau)} \rm{d\tau} \]

我們稱 $ h(x) $ 是 $ g(x) $ 與 $ f(x) $ 的卷積。

而此處我們討論的是整次多項式域，那么就把其定義划歸進多項式，得到

\[A(x) \cdot B(x) =\sum\limits_{i = 0}^{n} \sum\limits_{j=0}^{i}{a_jb_{i-j}} \]

其中 $ A(x) $ 和 $ B(x) $ 均為 \(n-1\) 次多項式
比較顯然的是，在我們合並同類項的情況下，最終卷出來的結果應該有 $ 2n+1 $ 項。

Dot Method

我們知道，原本的多項式是系數表示法，現在我們將其轉化為**點值表示法 $ (\boldsymbol{dot~method} ) $ **。即我們可以把多項式 $ F(x) $ 轉化為多項式函數 $ f(x) $ ，那么這個 $ n $ 階函數就可以由 $ n+1 $ 個點唯一確定。即給定

\[f(x)\Longleftrightarrow{(x_0,y_0),(x_1,y_1),(x_2,y_2)....(x_n,y_n)}\\ \forall k,y_k = f(x_k) \]

這是很顯然的，並且這 $ n+1 $ 個點是隨意選取的——只要求它們相異即可。

Advanced Trick Point $\color{red}{1} $ Multiplication

假設我們有兩個關於 $ x $ 的 $ n+1 $ 次多項式 $ A(x) $ 和 $ B(x) $ ，我們要對它的點值表達式進行乘法操作。由於結果有 $ 2n+1 $ 項，我們考慮補上一堆項，並對

\[A(x) = {(x_0,y_0),(x_1,y_1)....(x_{2n},y_{2n})}\\ B(x) = {(x_0,y_0'),(x_1,y_1')....(x_{2n},y_{2n})} \]

做乘法可得

\[A(x)B(x) = {(x_0,y_0y_0'),(x_1,y_1y_1')\cdots(x_{2n},y_{2n}y_{2n}')} \]

我們觀察點乘法，它的時間復雜度達到了 $ \Theta(n) $ ，完全可以接受。那么不妨先看一下算法的大體思路：

對於每個因子多項式，選取 $ n+1 $ 個點，得出點值表達式(復雜度 $ \Theta(n^2) $)

$ \longrightarrow $ 點乘法(時間復雜度 $ \Theta(n) $ )

\(\longrightarrow\) 將得出來的 $ C(x) $ 的點值表達式再轉換成系數表達式(復雜度 $ \Theta(n^2) $ )。

這就是 FFT ​ 的大體流程。轉化之后怎么沒多快常數還大了

雖然其余部分的時間復雜度還是很麻煩的 $ O(n^2) $ ，但是都是可以優化成 $ O(n\log n) $ 的。

本質上的 $ FFT $ 包含 $ \mathbf{DFT} $ （離散傅立葉變換）和 $ \mathbf{IDFT} $ （逆離散傅立葉變換）實際上， $\rm DFT $ 對應着的就是把系數表達式映射到點值表達式的過程， $\rm IDFT $ 對應着的就是我們把點值表達式映射到系數表達式的過程。

Base of Optimization

因為實際上，我們的第一步——求值(系數轉點值)和我們的第三步(點值轉系數)都是可以做到 $ n\log n $ 的，那么總的時間復雜度，漸進意義下就是 $ O(n\log n) $ 的。

下面就讓我們來看看如何優化：

Advanced Trick Point $ \color{red}{2} $ Unit Complex Root

$ n $ 次單位復根是滿足 $ \omega^n = 1 $ 的復數 $ \omega $ ，其中我們可以由復數的運算法則(輻角相加，模相乘)很簡單地得出 $ n $ 次單位根有 $ n $ 個 這個結論——亦或者是用代數基本定理證，都可以。

而又因為復數 $ \omega^n $ 在復數平面上的模都是一，所以相乘之后還會是一，那么所有的 $ \omega_i,1 \leq i \leq n $ 就會均勻分布在單位圓上，類似當 $ n = 8 $ 時它是這樣的：

我們考慮歐拉公式：

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

我們取 $ x =2\pi $ ，可以得到如下關系式：

\[e^{2 \pi i} = 1 = \omega^n \Longleftrightarrow \omega = e^{\frac{2\pi i}{n}} \]

們把此時的單位根稱之為主次單位根，記作 $$ \omega_n = e^{\frac{2\pi i}{n}} $$ 。

那么對於其他的單位根，記作 $$ \omega_n^k=e^{\frac{2\pi ik}{n}},0 \leq k < n $$ 都是主次單位根的整次冪，也就是上圖中的一圈。

以下是單位根的性質：

\(\frak{Elimination~Lemma}\)

對任何整數 $ n \geq 0,k \geq 0,d >0 $ ,有

\[\omega_{dn}^{dk} = \omega_n^k \]

$ \mathcal{Proof.} $

這個好像很好證的樣子……代入定義可以獲得

\[\omega_{dn}^{dk} = \omega^{\frac{2\pi dk}{dn}} = e^{\frac{2\pi ik}{n}} = \omega_n^k \]

$ \frak{Binary~Lemma} $

引理：對於任何大於 $ 0 $ 的偶數 $ n $ ，都有 $ n $ 個 $ n $ 次單位復根的平方的集合，等於 $ \dfrac{n}{2} $ 個 $ \dfrac{n}{2} $ 次單位復根的集合。

$ \mathcal{Proof.} $

我們可以由消去引理得到

\[(\omega _n^k)^2 = \omega^{2k}_n=\omega_{\frac{n}{2}}^k \]

那么

\[(\omega_n^{k + \frac{n}{2}})^2 = \omega_n^{2k + n} \Longrightarrow \omega_n^{2k} \cdot \omega_n^n \Longrightarrow \omega_n^{2k} = (\omega_n^k)^2 \]

那么接下來，如果對所有的 $ n $ 次單位跟平方一下，我們會發現 $ \frac{n}{2} $ 次單位根每個都恰好出現了兩次。也就是說，在 $ n $ 個 $ n $ 次單位復數根平方的集合里，只有 $ \frac{n}{2} $ 個元素。我們參考上面那張圖，它的意義就在於方向相對的兩只向量，其平方相等。

那么把所有 $ n $ 單位根的平方畫到一個數列上就是這樣。

這個引理直接保證了我們求值的復雜度為 $ \Theta(n \log n) $ 。

而我們在代碼實現中，不能直接得到 $ e $ 或者虛數 $ i $ ，所以這個時候求單位根的任務就交給了我們上文中提到過的歐拉公式。

$ \frak{Sum~Lemma} $

引理：對於任意 $ n>0 $ 且 $ k $ 不能整除 $ n $ ，我們都有

\[\sum\limits_{j =0}^{n-1}{(\omega_n^k)^j} = 0 \]

由幾何級數的求和公式（等比數列求和公式)

\[\sum\limits_{j = 0}^{n}{x^j} = \frac{x^{j +1} -1}{x -1} \]

可得

\[\sum\limits_{j =0}^{n-1}{(\omega_n^k)^j} = \frac{(\omega_n^k)^n -1}{\omega_n^k -1} \Longrightarrow \frac{(\omega_n^n)^k -1}{\omega_n^k -1} = \frac{(1)^k -1}{\omega_n^k -1} \]

由於保證了 $ k $ 不可整除 $ n $ 所以分母一定不為 $ 0. $

DFT \(\to\) FFT

那么我們在了解完單位復數根之后，便可以正式地對 DFT 給出定義與操作方案了。

DFT

對於我們已知的一個多項式

\[A(x) = \sum\limits_{i =0}^{n - 1}{a_ix^i} \]

在

\[\omega_n^0,\omega_n^1,\omega_n^2 \cdots \omega_n^{n-1} \]

處的取值，我們可以假定 $ n $ 是 $ 2 $ 的冪，因為即使它本身不是 $ 2 $ 的冪，我們也可以通過向高次冪補值為 $ 0 $ 的項來解決這個問題。

那我們現在會有一個向量組 $ \vec{a} = {a_1, a_2, a_3 \cdots a_{n-1}} $ 。對於 $ k = 0, 1, 2, \cdots n -1 $ ，定義 $ y_k $ 如下：

\[y_k = A(\omega_n^k)=\sum\limits_{j =0}^{n -1}{a_j \cdot \omega_n^{kj}} \]

那么向量

\[\vec{y} = {y_0, y_1, y_2 \cdots y_{n-1}} \]

就稱作系數向量 $ \vec{a} = {a_1, a_2, a_3 \cdots a_{n-1}} $ 的離散型傅立葉變換(DFT)。

嗯，這個離散型我們可以由點乘法聯想意會一下：本來 $ A(x) $ 是一個優美的多項式，轉變成函數之后是一條優美的曲線，結果你突然把它拆成了一堆離散的點，把它用點值法表示，故稱之曰：“離散型” 。

FFT 優化 DFT

在上文中我們分析過，將系數表達式轉化為點值表達式需要的時間復雜度為 $ O(n^2) $ ，這是朴素算法。而我們只需要用一種被稱作快速傅立葉變換（ $\mathbf{Fast \ \ Fourier \ \ Transformation} $ ）的方式，就可以將其時間復雜度壓縮成 $ O(n\log n) $ 。而在這里我們就用到了剛才證明的引理——折半引理。

我們考慮將原來的多項式

\[A(x) = a_0+a_1x+ a_2x^2 \cdots +a_{n-1}x^{n-1} \]

重定義成兩個次數為 $ \frac{n}{2} $ 的小多項式 $ A^{[0]}(x) $ 和 $ A^{[1]}(x) $ ：

\[A^{[0]}(x) = a_0 + a_2x+a_4x^2 \cdots +a_{n-2}x^{\frac{n}{2} - 1} \\ A^{[1]}(x) = a_1 + a_3x+a_5x^2 \cdots +a_{n-1}x^{\frac{n}{2} - 1} \]

那么也就是說， $ A^{[0]}(x) $ 存儲的是所有偶數位（二進制位最后一位是 $ 0 $ ），而 $ A^{[1]}(x) $ 存儲的是所有的奇數位（二進制位最后一位是 $ 1 $ ），那么有下式

\[A(x) = A^{[0]}(x^2)+xA^{[1]}(x^2) \]

那我們求 $ A(x) $ 在

\[\omega_n^0,\omega_n^1,\omega_n^2 \cdots ,\omega_n^{n-1} \]

處的值，就變成了先求出 $ A^{[0]}(x^2) $ 和 $ A^{[1]}(x^2) $ 的值，然后根據上式進行合並即可。

而顯然的是，根據折半引理，我們根本不需要 $ O(n) $ 求，而是通過數據規模不斷減小使之成為 $ O(\log n) $ 。於是，我們成功通過 FFT​ 優化了求值的復雜度。

那么同時對於另一邊，我們可以根據

\[A(\omega_n^{k+\frac{n}{2}}) = A^{[0]}(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A^{[1]}(\omega_n^{2k+n}) \Longrightarrow A^{[0]}(\omega_n^{2k}\cdot \omega_n^{n})-\omega_n^{k}A^{[1]}(\omega_n^{2k}\cdot \omega_n^{n}) \]

得到

\[A(\omega_n^{k+\frac{n}{2}})=A^{[0]}(\omega_n^{2k})-\omega_n^{k}A^{[1]}(\omega_n^{2k}) \]

從而有偽代碼：

int Lim = 1, N, M ;
function FFT(int lenth, complex *A, int flag){
    IF (Lim == 1) return ;
    complex A0[lenth >> 1], A1[lenth >> 1] ;//分成兩部分
    for(int j : 0 to lenth by_grow 2) A0[j >> 1] = A[j], A1[j >> 1] = A[j + 1] ;
    FFT(lenth >> 1, A0, flag) ;
    FFT(lenth >> 1, A1, flag) ;
    complex Wn = unit(,) , w = (1, 0) ;
  //Wn是單位根，w用來枚舉冪，即我們令主次單位根不變，由於其余單位根都是其整次冪，所以可以用一個w來記錄到第幾次冪
  /*此處求單位根的時候會用到我們的參數flag……嗯沒錯就用這一次，並且flag的值域為(-1, 1)……是的，只會有兩個值*/
    for(int j : 0 to (lenth >> 1) by_grow 1 with w = w * Wn){
        A[i] = A0[i] + A1[i] * w ;//應用公式，下同 
        A[i + (lenth >> 1)] = A0[i] - A1[i] * w ; //順便求出另一半，由折半引理可顯然。 
    } 
} 
function Main{
    input(N), input(M) ;
    for(i : 0 to N by_grow 1) => input(A) ;
    for(i : 0 to M by_grow 1) => input(B) ; 
    while(Lim < N + M) Lim <<= 1 ;//Lim為結果多項式的長度（暫時），化為2的冪的便於分治（二分）
    FFT(Lim, A, 1) ;//兩遍FFT表示從系數化為點值 
    FFT(Lim, B, 1) ;
    for(i : 0 to Lim by_grow 2) => A[i] *= B[i] ;//點乘法，此處需要重定義乘號，因為每一項現在表示的是一個點，有x和y兩個屬性qwq 
}

以上是基於 $ pks $ 標准下的偽代碼你可以試試在c++標准下運行，其中 $ for $ 循環部分，$ grow $ 表示當前循環變量的單次增量，之后帶有 $ with $ 表示每次循環結束都會進行的運算（下同

嗯，這就是求值的方法，好像很 $ nice $ 地達到了 $ O(n \log n) $ 。

FFT 優化 IDFT

上文中我們曾經提及過的范德蒙德矩陣可以放到這兒用：

\[\left|\begin{array}{ccccc}1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_{n} & \omega_{n}^{2} & \cdots & w_{n}^{n-1} \\ 1 & \omega_{n}^{2} & \omega_{n}^{4} & \cdots & \omega_{n}^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_{n}^{n-1} & \omega_{n}^{2(n-1)} & \cdots & \omega_{n}^{(n-1)(n-1)}\end{array}\right| \quad\left|\begin{array}{c}a_{0} \\ a_{1} \\ a_{2} \\ \vdots \\ a_{n-1}\end{array}\right|=\left|\begin{array}{c}y_{0} \\ y_{1} \\ y_{2} \\ \vdots \\ y_{n-1}\end{array}\right| \]

那為了求出我們的 $ \vec{a} = {a_0, a_1 \cdots ,a_{n-1}} $ 我們應該讓剛剛求值算出的 $ \vec{y} $ 乘上 $ \vec{V}^{~-1} $ 即可。於是有

推論：對於 $ j,k = 0,1, 2 \cdots n-1,V_n^{-1} $ 的 $ (j, k) $ 處的值為 $\dfrac{\omega_n^{-kj}}{n} $ 。

我們考慮驗證這個斷言的正確性，已知 $ V_n' $ 是一個 $ (j,k) $ 處值為 $ \omega_n^{-kj}/n $ 的、與 $ V $ 形態相同的矩陣，那我們只需要證明 $ V' \cdot V = I_n $ 即可，其中 $ I_n $ 是 $ n $ 階單位矩陣，即主對角線都是 $ 1 $ ，其余位置上是 $ 0 $ 的矩陣。

那么我們考察 $ V' V $ 中的元素 $ (i, j) $ ，有如下的式子

\[V'V = \sum\limits^ {n-1}_{k=0}{(\dfrac{\omega_n^{-ki}}{n})} \cdot {\omega_n^{kj}} = \frac{1}{n} \sum\limits^ {n-1}_{k=0}{\omega_n^{k(j-i)}} \]

由求和引理當且僅當 $ i=j $ 時其值為一，其余的時刻均為零，所以有 $ V'V = I_n $ 。

那么我們把我們剛剛求出來的逆矩陣 $ V^{-1} $ 美化一下，提出每一項所除的 \(n\) :

\[\mathrm{IDFT}_n(y) = \frac{1}{n}\sum\limits_{k = 0}^{n-1}{y_k\omega_n^{-kj}},j\in [0,n-1] \]

誒，這個好像…跟我們求值時的公式差不多？沒錯，除了帶個負號，其余的都差不多。所以我們可以考慮打個標記：當 flag=1 時，他是正向 DFT ；當它等於 $ -1 $ 時，它是逆向的 IDFT 。這可以讓我們通過這一個函數解決兩個過程。我們只需要用 $ y $ 替換 $ a $ ，用 $ \omega_n^{-1} $ 替換 $ \omega_n $ ，其余的沒什么差別，於是復雜度還是 $ O(n \log n) $ 的。

void FFT(int Lim,complex *A,int flag){
    if(Lim == 1) return ;
    complex A0[Lim >> 1], A1[Lim >> 1] ;
    for(int i = 0; i <= Lim ; i += 2)
        A0[i >> 1] = A[i], A1[i >> 1] = A[i+1] ;
    FFT(Lim >> 1, A0, flag) ;
    FFT(Lim >> 1, A1, flag) ;
    complex unit = (complex){cos(2.0 * Pi / Lim) , flag * sin(2.0 * Pi / Lim)}, w = complex(1, 0) ;//歐拉公式 
    for(int i = 0;i < (Lim >> 1) ; i ++, w = w * unit) {
        A[i] = A0[i] + w * A1[i] ;
        A[i + (Lim>>1)] = A0[i] - w*A1[i];
    }
}

int main(){
......................
FFT(A, 1), FFT(B, 1) ;
for(i = 0; i <= Lim; i ++) A[i] = A[i] * B[i] ;
FFT(A, -1) ;
......................
}

好的，現在嘛……可以考慮撒花花啦！因為我們的 $ FFT $ 實際上已經結束了！ $ But $ ，這個遞歸版本的 $ FFT $ 由於牽扯到 $ sin/cos $ 的運算、 $ double $ 、遞歸時的入棧出棧（底層），所以常數特別的大 $ emmmmm $ ，那么——

Iterative Optimization

我們現在要引出的就是迭代版的 FFT 。

Advanced Trick Point $ \color{red}{3} $ The Butterfly Operation

$ emmm $ 先上一個不是特別卡常數的優化。我們觀察之前的代碼中，有這么一步：

   for(int i = 0;i < (Lim >> 1) ; i ++, w = w * unit) {
        a[i] = A0[i] + w * A1[i] ;
        a[i + (Lim>>1)] = A0[i] - w*A1[i];
}

我們會發現…… $ \omega \cdot A^{[1]}[i] $ 被執行了兩次，所以我們不妨用個變量記錄它：

   for(int i = 0;i < (Lim >> 1) ; i ++, w = w * unit) {
        int temp = w * A1[i] ;
        a[i] = A0[i] + t ;
        a[i + (Lim>>1)] = A0[i] - t ;
    }

嗯，這就是全部的優化啦！那么，FFT，完！

$ qwq $ 這分明是騙小孩子的啦……如果單單這一步就可以卡常數的話，那這個世界該多么美好 $ \mathcal{QAQ} $ 。

好吧，說這個的原因，只是為了引出我們關於蝴蝶操作的定義：

我們定義 $ \omega_n^k $ 為旋轉因子，那么每一次我們先將 $ y_k^{[1]} $ 與旋轉因子的乘積存儲在一個變量 $ t $ 里，並在 $ y_k^{[0]} $ 增加、減去 $ t $ 的操作稱為一次蝴蝶操作。

說白了，蝴蝶操作是一次 $ O(2) $ 的求出 $ A^{[0]}_k $ 與 $ A^{[1]}_k $ 的操作。

我們首先考慮按照遞歸的思路，將 FFT 的分治流程刻畫一下：

我們會發現，其實我們是可以對它進行反向迭代的。以上面的迭代樹為例，我們的步驟大體如下：

step 1 成對地取出兒子節點，用蝴蝶操作計算出其 DFT。
step 2 用這一步的 DFT 替換之前的。
step 3 直到我們迭代到根節點為止，否則返回 step 1。

而反向迭代似乎有規律可循。我們發現只要我們用迭代的過程模擬出回溯的過程即可。

那么思路便有了：三層 $ for $ ，先枚舉區間長度（1，2，4，8……），第二層枚舉長度為 $ j\times 2 $ 的每個區間的起點——意圖為同時枚舉兩個相鄰區間，便於相鄰區間之間 DFT 的合並，第三層負責遍歷每段區間，運用蝴蝶操作逐個合並：

for(j = 1; j < Lim; j <<= 1){//枚舉區間長度，從小區間到大區間依次合並。
        node T(cos(Pi / j), flag * sin(Pi / j)) ;
        for(k = 0; k < Lim; k += (j << 1) ){//兩段區間兩段區間的枚舉，用於合並
            node t(1, 0) ;
            for(l = 0 ; l < j; l ++, t = t * T){//枚舉k所枚舉的兩個區間內的值，並進行蝴蝶操作。
                node Nx = J[k + l], Ny = t * J[k + j + l] ;
                J[k + l] = Nx + Ny ; J[k + j + l] = Nx - Ny ;//一次蝴蝶操作
            }
        }
    }

嗯，好像…海星？哈，思維不嚴謹的我們漏了一個地方：我們在 DFT 的時候，為了保證時間復雜度是 $ \Theta(\log n) $ ，我們曾經進行過一次 $ A(x) = A^{[0]}(x^2)+xA^{[1]}(x^2) $ 的操作，所以我們需要自動調整順序。通俗一點，就是我們原來的序列順序是 $ 0,1,2,3,4,5,6,7 $ ，但是迭代版的 FFT 卻需要的順序應該跟葉子結點的順序吻合，即 $ 0, 4, 2, 6, 1, 5,3,7 $ 。所以 ——

Trick Point $\color{red}{4} $ The Butterfly Law

這個嘛……我們可以選擇打個表觀察：

原來的序號 $ 0 \ \ 1 \ \ 2 \ \ 3 \ \ 4 \ \ 5 \ \ 6 \ \ 7 $
現在的序號 $ 0 \ \ 4 \ \ 2 \ \ 6 \ \ 1 \ \ 5 \ \ 3 \ \ 7 $
原來的二進制表示 $ 000 \ \ 001 \ \ 010 \ \ 011 \ \ 100 \ \ 101 \ \ 110 \ \ 111 $
現在的二進制表示 $ 000 \ \ 100 \ \ 010 \ \ 110 \ \ 100 \ \ 101 \ \ 011 \ \ 111 $

誒，二進制好像是反序的嗷~這便是我們的最后一個 trick ，蝴蝶定理。而因為我們觀察到的蝴蝶定理是滿足一一對應性的，所以我們在 FFT 之前 swap 一遍即可。

嗯，然后我們可以將這個反序存在一個數組里面。類似這樣求出來：

for(i = 0; i < Lim; i ++ ) R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1)) ;

這大概是一個類似遞推的過程。

那么我們可以看到，這就簡化了很多冗余的步驟，並讓我們脫離遞歸的大常數。真開森啊

最后附迭代版的代碼(我寫的常數好像有點兒大 $ QAQ $ )

#include <cmath>
#include <cstdio>
#include <iostream>
#define il inline

using namespace std ;
int N, M, K ;
const int MAXN = 3000100 ;
const double Pi = acos(-1.0) ;
int i, j, k, l, Lim = 1, L, R[MAXN] ; 
struct node{
    double x, y ;
    node (double xx = 0, double yy = 0){
        x = xx, y = yy ;
    }
}A[MAXN], B[MAXN] ;
node operator * (node J, node Q){
    return node(J.x * Q.x - J.y * Q.y , J.x * Q.y + J.y * Q.x);
}
node operator + (node J, node Q){
    return node(J.x + Q.x , J.y + Q.y);
}
node operator - (node J, node Q){
    return node(J.x - Q.x , J.y - Q.y );
}

il int qr(){
    int k = 0, f = 1 ;
    char c = getchar() ;
    while(!isdigit(c)){if(c == '-') f = -1 ;c = getchar() ;}
    while(isdigit(c)) k = (k << 1) + (k << 3) + c - 48 ,c = getchar() ;
    return k * f ;
}
void FFT(node *J, int flag){
    for(i = 0; i < Lim; i ++)
        if(i < R[i]) swap(J[i], J[R[i]]) ;//前面的if保證只換一次
    for(j = 1; j < Lim; j <<= 1){
        node T(cos(Pi / j), flag * sin(Pi / j)) ;
        for(k = 0; k < Lim; k += (j << 1) ){
            node t(1, 0) ;
            for(l = 0 ; l < j; l ++, t = t * T){
                node Nx = J[k + l], Ny = t * J[k + j + l] ;
                J[k + l] = Nx + Ny ;
                J[k + j + l] = Nx - Ny ;
            }
        }
    }
}
int main(){
    N = qr(), M = qr() ; 
    for(i = 0; i <= N; i ++) A[i].x = qr() ;
    for(i = 0; i <= M; i ++) B[i].x = qr() ;
    while(Lim <= N + M) Lim <<= 1, L ++ ;
    for(i = 0; i < Lim; i ++ ) R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1)) ;
    FFT(A, 1), FFT(B, 1) ;
    for(i = 0; i <= Lim; i ++) A[i] = A[i] * B[i] ;
    FFT(A, -1) ;
    for(i = 0; i <= N + M; i ++)
        printf("%d ", (int)(A[i].x / Lim + 0.5)) ;//我們推過的公式里面有一個1/n這一項，最后輸出的時候添上即可 
    return 0 ;
}

啊……那就撒花花吧！！

Afterword

以下是原尾語，保留了下來：

嗯……怎么說呢，現在看這個算法，真是簡單的一匹啊……代碼這么短這么容易背過。但是當時理解起來卻花了很大心思呢！這篇博客我寫了整整三天 $ qwq $ ，由於要培訓和考試，所以拖拖沓沓地寫了三天，一邊寫一邊感嘆自己理解的簡直太淺顯了。每一個證明、每一個引理、甚至每一個符號，我都需要去和其他 $ DALAO $ 比對審核、或者纏着 $ rqy $ 問個沒完；每次一講到原理，我都發現自己原來並不理解那些，於是不得不推倒重來。這篇博客會持續更新，補充語意不明、證明難以理解的地方。

以下是溫馨提示：

1、好多自己當初不理解的地方在代碼里就只有半行qaq

2、三個引理中，只有消去引理跟算法的實現沒有關系——消去引理主要是用來證明其他引理的

真 · 結束語：

其實沒什么好說的，今天重新復習了一遍，發現自己以前有好多內容雖然如原尾語所言，看上去現在看這個算法，真是簡單的一匹啊，但實際上忽略了好多東西。我想大概只有一遍一遍地鑽研才能了解完全一件事情吧。

$ \rm{Reference} $

• $ [1] $ : $ rvalue $ 的 $ blog $ $ ^{^{[\nearrow ]}} $
• $ [2] $ :算法導論 $ ^{^{[\nearrow]}} $ 提取碼: txs2
• [3]*：鳴謝rqy

$ \mathfrak{writter:pks} $