---
title: 快速傅立葉轉換
tag: 數論
disqus: hackmd
---

# 快速傅立葉轉換 (FFT)

## 序

本篇包含大量 $\LaTeX$ 數學表達式，如果爛掉請嘗試重新載入。

~~因為讀學測讀得很累，所以換換口味~~
最近突然腦子撞到，所以上網查了大數乘法的資料，不過似乎並不是特別的豐富(~~肯定是我不會查~~)，而且相關的文章對「為何這麼做的原因」的解釋都不太多，所以決定自己寫一份筆記，因為是筆記所以會寫的特別白話，希望對看這篇的人有幫助。
如果有疏漏歡迎留言指正。

這篇的目標是讓人一次就看懂 FFT，不過當然要一些國高中的先備知識。

#### 本篇會用到的數學與演算法：
- 數學
- 多項式的基本知識
- 簡單的三角函數
- 虛數（只要知道 $i^2 = -1$ 就可以了）
- [矩陣乘法](https://zh.wikipedia.org/zh-tw/矩陣乘法)
- 演算法
- 時間複雜度
- 分治

## 大數計算

如果我們需要計算的數字大於 C++ 裡 long long 的儲存範圍，那麼我們就只能用字串來存放數字，並利用模擬的方式來計算答案，不過感覺正常來說這種需求不是很多，所以主要當個切入點。

假設現在有兩個數 $A, B$，並稱兩個數字的位數是 $a, b$。
普通的加減計算很簡單，你平常怎麼算程式就怎麼寫，複雜度是 $O(a+b)$。
不過到了乘法，我們平常的計算方式所需的複雜度是 $O(ab)$。

\begin{array}{r}
111 \\
\underline{\times \ \color{red}2 \color{green}3 \color{darkorange}4} \\
\color{darkorange}{444} \\
\color{green}{333} \phantom{1} \\
\color{red}{222} \phantom{11} \\
\hline{25974}
\end{array}

這聽起來不太妙，所以我們如何優化它？
其中一個方式是把它轉化成多項式乘法。

我們知道一個整數都可以有以下的表示方法：
$$
1234 = 1 \times 10^3 + 2 \times 10^2 + 3 \times 10^1 + 4 \times 10^0
$$

這種寫法就和多項式一樣
\begin{split}
Num(x) &= 1x^3 + 2x^2 + 3x + 4 \\
Num(10) &= 1234
\end{split}

因此，我們可以**把大數乘法直接變成多項式乘法**，將多項式乘法做完之後，再把 10 代進結果中就可以得到答案了，即：
$$
\begin{split}
A(x) &= s_0 + s_1x + s_2x^2 + s_3x^3 \\
B(x) &= t_0 + t_1x + t_2x^2 + t_3x^3 \\
C(x) &= A(x) \times B(x) \\
C(10) &= A(10) \times B(10) = ans
\end{split}
$$

所以就像剛剛說的，大數乘法只是一個切入點，實際上我們在意的是**如何快速的計算出兩個多項式的乘積**。計算多項式的乘積在很多地方都有用途，而且多項式的項數可以非常多，所以計算的時間複雜度就是一個需要在意的事情了。

## 多項式乘法

在討論如何快速的多項式乘法前，我們先來看看如何表示一個多項式。

### 多項式表示法

1. 係數表示法
一個多項式 $f(x) = p_0 + p_1x + p_2x^2 + p_3x^3$，只要規定好書寫的方向（升冪或降冪），我們就可以只寫係數而不寫後面的 $x$，實際上我們在學校學綜合除法時就用過這種表示方式，讓算式變得比較乾淨，如：用 $[p_0, p_1, p_2, p_3]$ 表示上面的 $f(x)$。
這種寫法也可以稱為一個「**向量**」

2. 點值表示法
我們知道，在座標平面上，可以用兩個相異的點來決定一條唯一的直線（$x$ 的一次式）
實際上，我們可以利用 $n+1$ 個點來決定一個 $x$ 的 $n$ 次多項式，所以我們也可以利用 $n+1$ 個點來表示一個多項式。
而利用這 $n+1$ 個點反推回多項式的方法稱為「**插值法**」。
:::spoiler 插值法的唯一性證明
如果我們選取的每個點都對決定這個多項式有貢獻（我們選的點都滿足線性獨立），那麼，最少只需要 $n+1$ 個點就可以決定一個唯一的多項式，因為一個 $n$ 次多項式有 $n+1$ 個係數。

$$
\begin{split}
f(x) = p_0 + p_1x + p_2x^2 + \cdots + p_nx^n
\end{split}
$$

而我們可以將這 $n+1$ 個點帶入 $f(x)$ 而得到 $n+1$ 個方程式。

$$
\begin{split}
y_1 = p_0 + p_1x_1 + &p_2x^2_1 + \cdots + p_nx^n_1 \\
y_2 = p_0 + p_1x_2 + &p_2x^2_2 + \cdots + p_nx^n_2 \\
\vdots \\
y_{n+1} = p_0 + p_1x_{n+1} + &p_2x^2_{n+1} + \cdots + p_nx^n_{n+1} \\
\end{split}
$$

藉由解這 $n+1$ 個方程式就可以解出係數了。
:::

### 多項式乘法

既然用係數表示的多項式沒辦法高效的完成乘法，那我們來看看如果多項式使用點值表示法，做乘法的效率會如何？

一樣是兩個多項式 $A(x)$ 和 $B(x)$，不過這次他們的最高次都是 $n$
那我們可以利用 $n+1$ 個點來表示這兩個多項式：

$$
\begin{matrix}
(0, A(0)) & (1, A(1)) & \cdots & (n, A(n)) \\
(0, B(0)) & (1, B(1)) & \cdots & (n, B(n)) \\
\end{matrix}
$$

而我們要的 $C(x) = A(x) \times B(x)$ 可以寫成這樣：

$$
\begin{array}{l}
(0, C(0)) & (1, C(1)) & \cdots & (n, C(n)) \\
(0, A(0) \times B(0)) & (1, A(1) \times B(1)) & \cdots & (n, A(n) \times B(n)) \\
\end{array}
$$

但是我們知道 $A(x) \times B(x)$ 後，$C(x)$ 的最高次次數是 $2n$，所以我們應該準備 $2n+1$ 個點來確定出唯一的 $C(x)$，所以應該改寫成：

$$
\begin{array}{l}
(0, C(0)) & (1, C(1)) & \cdots & \color{red}{(2n, C(2n))} \\
(0, A(0) \times B(0)) & (1, A(1) \times B(1)) & \cdots & \color{red}{(2n, A(2n) \times B(2n))} \\
\end{array}
$$

如此這般，如果利用點值來表示 $A(x)$ 和 $B(x)$，那我們的多項式乘法就可以在 $O(n)$ 的時間內做完了，可喜可賀。

所以現在的問題就變成了「**如何快速的求出點值**」和「**如何快速的將點值轉換回多項式**」

## 多項式求值

對多項式求值這件事我們在寫數學題目時在熟悉不過了，但是仔細想想就會發現，我們使用的方法時間複雜度其實非常糟糕，以下是一些求值的方法。

1. 直接一個一個算
如果對 $f(x) = p_0 + p_1x + \cdots + p_nx^n$ 直接一個一個算，而且計算 $2n+1$ 個點的值，那麼時間複雜度將會是 $O(\frac{n\times(n+1)}{2} \times (2n+1)) = O(n^3)$，非常不妙。
2. 簡單的優化
如果對 $f(x)$ 做一些簡單的合併，就會得到這個：
$$
\begin{split}
f(x) = p_0 + x(p_1 + x(p_2 + x(\cdots + x(p_{n-1} + p_nx))))
\end{split}
$$
這樣一來，求全部值的複雜度就變成了 $O(n \times (2n+1)) = O(n^2)$，可以更好嗎？

3. 使用魔法
來到本篇的重點，如果使用快速傅立葉轉換，我們可以將這個複雜度在優化到 $O(n\log{n})$

所以終於來到了這篇的主題了，什麼是「**快速傅立葉轉換**」？為何他可以快速的求值？

## 快速傅立葉轉換

快速傅立葉轉換 Fast Fourier transform 簡稱 FFT

什麼是傅立葉轉換？最簡化的理解就是**對多項式求值**，不過用的地方不一樣所以名字比較酷。

所以快速傅立葉轉換就是「快速的對多項式求值」，也就是剛剛的多項式乘法中，時間複雜度的瓶頸之處。而快速傅立葉轉換利用函數的特性，將多項式求值的複雜度從 $O(n^2)$ 降到 $O(n\log{n})$。

接著就來看看他是如何做到的

### 函數值的觀察

下面是 $f(x) = x^2$ 的函數圖形


我們可以發現 $f(-x) = f(x)$，所以如果我們選的是這六個點，其實只需要算三次值
而對於 $x$ 項都是偶數次方的函數來說這件事都成立，這種函數叫做「**偶函數**」

而這是 $f(x) = x^3$ 的函數圖形


我們又可以發現 $f(-x) = -f(x)$，一樣可以減少計算次數
而對於 $x$ 項都是奇數次方的函數來說這件事都成立，這種函數叫做「**奇函數**」

### 原函數的修改

假設函數總共有 $n$ **項**（$n$ 項多項式最高次為 $n-1$）
我們原來的函數長這樣：
$$
f(x) = p_0 + p_1x + p_2x^2 + \cdots + p_{n-1}x^{n-1}
$$

先假設 $n$ 是一個二的倍數，於是我們可以將他分解成一個奇函數和一個偶函數的和：
\begin{split}
f_{even}(x) &= p_0 + p_2x^2 + p_4x^4 + \cdots + p_{n-2}x^{n-2} \\
f_{odd}(x) &= p_1x + p_3x^3 + p_5x^5 + \cdots + p_{n-1}x^{n-1} \\
f(x) &= f_{even}(x) + f_{odd}(x)
\end{split}

並利用剛剛的結論：
\begin{split}
f_{even}(-x) &= f_{even}(x) \\
f_{odd}(-x) &= -f_{odd}(x) \\
\end{split}

我們就可以得到這個等式：
\begin{array}{l}
f(x) &= f_{even}(x) + f_{odd}(x) \\
f(-x) = f_{even}(-x) + f_{odd}(-x) &= f_{even}(x) \color{red}- f_{odd}(x)
\end{array}

這樣我們就只需要計算 $n$ 個 $f_{even}(x)$ 和 $f_{odd}(x)$ 就可以得到 $2n$ 個點，雖然一乘一除看起來計算次數不變，不過這一步的意義在於我們將函數分解了。

### 新函數的計算

我們現在成功的將原函式分成了兩個新的多項式，那有沒有可能繼續將新函式分解成兩個，如果可以一直分到函式只剩一項，我們就可以直接回傳常數項了。
而且這樣最多分解成 $\log n$ 個新函式，複雜度變可以下降成 $O(n\log n)$。

讓我們先來觀察兩個新函數的不同：
\begin{split}
f_{even}(x) &= p_0x^\color{red}0 + p_2x^\color{red}2 + p_4x^\color{red}4 + \cdots + p_{n-2}x^\color{red}{n-2} \\
f_{odd}(x) &= p_1x^\color{red}1 + p_3x^\color{red}3 + p_5x^\color{red}5 + \cdots + p_{n-1}x^\color{red}{n-1} \\
\end{split}

除了係數不一樣外，如果將 $f_{odd}(x)$ 的 $x$ 提出來，兩個式子就長得很像了：
\begin{split}
f_{even}(x) &= \ \ \ \ p_0x^\color{red}0 + p_2x^\color{red}2 + p_4x^\color{red}4 + \cdots + p_{n-2}x^\color{red}{n-2} \\
f_{odd}(x) &= \color{red}x(p_1x^\color{red}0 + p_3x^\color{red}2 + p_5x^\color{red}4 + \cdots + p_{n-1}x^\color{red}{n-2}) \\
\end{split}

既然兩個式子長的差不多，那接下來就以 $f_{even}(x)$ 當作例子。

我們可以發現，如果將 $k = x^2$ 對新函數做變數代換：
\begin{split}
f_{even}(k) &= p_0 + p_2k^\color{red}1 + p_4k^\color{red}2 + \cdots + p_nk^\color{red}{\frac{n}{2}-1} \\
\end{split}

就可以得到一個與原函數幾乎一模一樣的新函數。
如此一來我們就可以用分治的方法繼續分下去了$\cdots$ 嗎？

### 繼續分下去？

我們分治的地方在於每次將函數分成奇偶函數，而這樣做的前提是我們需要將代入函數的值分成**正負的點對**，如果不能做到這一點，那我們基於奇偶函數特性所創造的算法就無法成立。

然而我們剛剛對新函數 $f_{even}(x)$ 用 $k$ 進行了代換變成了 $f_{even}(k)$ ，而 $k$ 是 $x^2$ 所以 $k$ 不可能是負的，也就無法將代入 $f_{even}(k)$ 的值分成「正負點對」，除非...

除非 $x$ 是一個平方過後還會出現負數的怪物

### 對 $x$ 使用虛數吧

虛數的定義就是 $\sqrt{-1}$ 並且以 $i$ 作為代號，所以根據定義 $i^2$ 就等於 $-1$。
而我們的 $x$ 值必然會出現虛數，因為這樣才能滿足我們的需求。

不過虛數世界千千萬萬，具體來說會出現哪些虛數呢？
從原函數來看我們不好決定要代哪些虛數回去，所以不妨從遞迴的中止條件來看看哪些數比較簡單。

我們的遞迴會不斷地將原函數分成奇偶兩個，所以理所當然遞迴的中止條件就是分到最後的函數只剩一項時，我們就可以停下了，而這個只有常數項的函數，其實不管 $x$ 代什麼結果都一樣，所以我們代個最簡單的數「$1$」就可以了。

而我們是透過不斷用 $x^2$ 代換來得到新的函數，所以中止條件所代入的 $1$，在回傳的途中就會被不斷的開根號，所以最後我們原函數需要代入的 $x$ 值就是 $\sqrt[2^k]{1}$，總共 $2^k$ 個數 $(k \in \mathbb{Z})$。
:::spoiler 根號的解的個數
$\sqrt[n]{1}$ 也可以寫成 $x^n = 1$ 的解，而一個 $n$ 次多項式有 $n$ 個解，所以 $\sqrt[n]{1}$ 會有 $n$ 個解。
:::

<br/>
然而，我們原函式的項數不一定是 $2^k$ 個，所以我們必須找到一個**比項數還要大的二的冪次**，而且為了避免在遞迴途中還要一直做這件事（找冪次），所以我們可以一開始就將所求的函數補成 $2^k$ 項，我們總是可以透過將係數設為 $0$ 來做到，所以剛剛才會假設多項式的項數 $n$ 是一個二的倍數。

說了這麼多，所以什麼是 $\sqrt[2^k]{1}$，它究竟長怎樣？

### $\sqrt[2^k]{1}$ 長怎樣

我們會利用複數平面來表示一個複數， $x$ 軸是實數而 $y$ 軸是複數，這樣就可以在直角坐標系中表示出任意一個複數，而特別的是，**$\sqrt[n]{z}$ 的 $n$ 個解在複數平面上是一個圓的內接正 $n$ 邊形**。

那如何可以表示出這正 $n$ 邊形每個頂點的數呢？

在平常的 $xy$ 平面上，我們會用 $(x, y)$ 來表示出一個點，而我們也可以用三角函數來表示這個點。
假設 $r = \sqrt{x^2 + y^2}$，也就是 $(x, y)$ 到原點的距離，那這個點也可以表示成 $(r\cos{\theta}, r\sin{\theta})$，其中 $\theta$ 是 $(x, y)$ 向量與 $x$ 軸正向的夾角。

而在複數平面上我們也可以這樣做，不過 $y$ 要多寫一個 $i$ 來表示複數軸。
這樣一來，複數 $z$ 就等於 $r(\cos{\theta} + i\sin{\theta})$。
好消息是，因為我們選的 $z$ 是 $1$，所以 $r = 1$，而且這 $n$ 個點是正多邊形，所以 $\theta$ 每次增加 $\frac{2\pi}{n}$。

這樣我們就得到了解出來的 $n$ 個點：
\begin{split}
(\cos{0} + i\sin{0}), (\cos{\frac{2\pi}{n}} + i\sin{\frac{2\pi}{n}}), (\cos{\frac{2\pi(2)}{n}} + i\sin{\frac{2\pi(2)}{n}}), \cdots, (\cos{\frac{2\pi(n-1)}{n}} + i\sin{\frac{2\pi(n-1)}{n}})
\end{split}

不過這樣子的表示法看著很不優雅，所以我們可以引入歐拉公式來簡化點的表示：
\begin{split}
e^{i\theta} = \cos{\theta} + i\sin{\theta}
\end{split}

原本的點就可以表示成：
\begin{split}
e^{i\theta(0)}, e^{i\theta(1)}, e^{i\theta(2)}, \cdots, e^{i\theta(n-1)}
\end{split}

如果我們令 $w$ 等於 $e^{i\theta}$，我們就可以得到這些點的簡單寫法：
\begin{split}
w^0, w^1, w^2, \cdots, w^{n-1}
\end{split}

這些就是我們所有需要的點值了。

### 實作快速傅立葉轉換

因為我們利用分治拆解多項式，所以多項式的項數需要是 2 的冪次，所以我們一開始假設項數 $n$ 是一個二的倍數並沒有問題，而且我們總是可以透過為零的係數來增加多項式的項數。

就此，我們就分解出了完成快速傅立葉轉換的所有步驟：


所以就來實作吧，這次使用 python 因為比較好讀

e  = 2.71828182845 # 自然常數
pi = 3.14159265358 # 圓周率

def FFT(f):
    n = len(f) # n 需為 2 的冪次
    if n == 1: # 如果只剩一項就直接回傳
        return f
    
    w = e ** (2*pi*1j / n) # j 為 python 的虛數單位 
    even, odd = f[::2], f[1::2]  # 分成奇偶函數

    even, odd = FFT(even), FFT(odd) # 遞迴計算

    y = [0] * n
    for i in range(n//2): # 合併
        y[i] = even[i] + w**i * odd[i]
        y[i + n//2] = even[i] - w**i * odd[i]

    return y # 回傳解答

到此，我們將計算點值的複雜度優化到了 $O(n\log{n})$，完成了多項式乘法的第一個部分「**如何快速的求出點值**」，在進到第二個「**如何快速的將點值轉換回多項式**」之前，我們先來看看函數求值的另一種表示法。

## 矩陣表示法

在剛剛的多項式表示法中我們提到，$[p_0, p_1, p_2, \cdots, p_n]$ 這種表示法稱為向量，也就是一個一維的矩陣，所以我們剛剛的函數求值可以寫成這樣的矩陣乘法：

\begin{equation}
\begin{bmatrix}
f(x_0) \\
f(x_1) \\
f(x_2) \\
\vdots \\
f(x_{n}) \\
\end{bmatrix} =
\underbrace{
\begin{bmatrix}
1 & x_0 & x_0^2 & \cdots & x_0^n \\
1 & x_1 & x_1^2 & \cdots & x_1^n \\
1 & x_2 & x_2^2 & \cdots & x_2^n \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & x_n & x_n^2 & \cdots & x_n^n \\
\end{bmatrix}}_{\Huge{M}}
\times
\begin{bmatrix}
p_0 \\
p_1 \\
p_2 \\
\vdots \\
p_n
\end{bmatrix}
\end{equation}

這樣就可以乾淨的表示出函數求值的過程，那我們能否把「**點值轉換回多項式**」也表示成矩陣乘法？

當然可以，下面是讓點值經過一個矩陣 $m$ 的作用後得到係數：
\begin{equation}
\begin{bmatrix}
p_0 \\
p_1 \\
p_2 \\
\vdots \\
p_n
\end{bmatrix} =
{\Huge{m}}
\times
\begin{bmatrix}
f(x_0) \\
f(x_1) \\
f(x_2) \\
\vdots \\
f(x_{n}) \\
\end{bmatrix}
\end{equation}

而這個 $m$ 就是 $M$ 的反矩陣 $M^{-1}$

:::spoiler 什麼是反矩陣
在平常的數學運算 $a = b\times c$ 我們可以用 $a / b$ 來求出 $c$，但是在矩陣中我們無法定義「除法」，而 $a / b$ 也可以寫成 $a \times b^{-1}$，所以我們用類似倒數的方法來將除法變成乘法（除以 $n$ 就是乘以 $\frac{1}{n}$ 嘛）

如果定義 $I$ 為單位矩陣（可以理解成矩陣的 $1$ ：$IV = V$），那 $M$ 和 $M^{-1}$ 的關係就是：
\begin{split}
MM^{-1} = I
\end{split}

**簡單理解 $M^{-1}$ 就是一個矩陣 $M$ 的倒數，與 $M$ 相乘後就會兩個一起消失**。
:::

<br/>

所以我們傅立葉轉換的矩陣乘法寫成這樣

\begin{equation}
\begin{bmatrix}
f(w^0) \\
f(w^1) \\
f(w^2) \\
\vdots \\
f(w^{n-1}) \\
\end{bmatrix} =
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & w & w^2 & \cdots & w^{n-1} \\
1 & w^2 & w^4 & \cdots & w^{2(n-1)} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & w^{n-1} & w^{2(n-1)} & \cdots & w^{(n-1)(n-1)} \\
\end{bmatrix}
\times
\begin{bmatrix}
p_0 \\
p_1 \\
p_2 \\
\vdots \\
p_{n-1}
\end{bmatrix}
\end{equation}

在我們對 $A(x)$ 和 $B(x)$ 都求出值並乘出 $C(x)$ 後，我們可以利用反矩陣求回多項式的係數

\begin{equation}
\begin{bmatrix}
p_0 \\
p_1 \\
p_2 \\
\vdots \\
p_{n-1}
\end{bmatrix} =
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & w & w^2 & \cdots & w^{n-1} \\
1 & w^2 & w^4 & \cdots & w^{2(n-1)} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & w^{n-1} & w^{2(n-1)} & \cdots & w^{(n-1)(n-1)} \\
\end{bmatrix}^{-1}
\times
\begin{bmatrix}
C(w^0) \\
C(w^1) \\
C(w^2) \\
\vdots \\
C(w^{n-1}) \\
\end{bmatrix}
\end{equation}

好消息是，對於這個每列都是等比數列的大矩陣有個特殊的名字叫「范德蒙矩陣」，而他的反矩陣可以用公式簡單的求出來，這裡就直接給出結果：

\begin{equation}
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & w & w^2 & \cdots & w^{n-1} \\
1 & w^2 & w^4 & \cdots & w^{2(n-1)} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & w^{n-1} & w^{2(n-1)} & \cdots & w^{(n-1)(n-1)} \\
\end{bmatrix}^{-1} =
\frac{1}{n}
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & w^{-1} & w^{-2} & \cdots & w^{-(n-1)} \\
1 & w^{-2} & w^{-4} & \cdots & w^{-2(n-1)} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & w^{-(n-1)} & w^{-2(n-1)} & \cdots & w^{-(n-1)(n-1)} \\
\end{bmatrix}
\end{equation}

所以可以得到：

\begin{equation}
\begin{bmatrix}
p_0 \\
p_1 \\
p_2 \\
\vdots \\
p_{n-1}
\end{bmatrix} =
\frac{1}{n}
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & w^{-1} & w^{-2} & \cdots & w^{-(n-1)} \\
1 & w^{-2} & w^{-4} & \cdots & w^{-2(n-1)} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & w^{-(n-1)} & w^{-2(n-1)} & \cdots & w^{-(n-1)(n-1)} \\
\end{bmatrix}
\times
\begin{bmatrix}
C(w^0) \\
C(w^1) \\
C(w^2) \\
\vdots \\
C(w^{n-1}) \\
\end{bmatrix}
\end{equation}

這個矩陣乘法看起來跟我們原本 FFT 的矩陣乘法一模一樣，只是 $w$ 變成了 $w^{-1}$，最後在每項乘以 $\frac{1}{n}$
所以我們完全可以直接修改剛剛 FFT 的程式中的 $w$ 值就可以做到「**快速的將點值轉換回多項式**」。

## Inverse FFT

因為是以 FFT 的反矩陣來構成的算法，所以叫做 Inverse FFT，也就是反向的 FFT，我們只需修改 $w$ 就可以了，並在最後將每項都乘以 $\frac{1}{n}$，其他程式碼都不變。
時間複雜度一樣是 $O(n\log{n})$。

e  = 2.71828182845 # 自然常數
pi = 3.14159265358 # 圓周率

def Inv_FFT(y, original = True):# 傳入函數值 y
    n = len(y)
    if n == 1:
        return y
    
    w = e ** (2*pi*1j / n)
    w = 1 / w # 新的 w
    # or you can write w = e ** (-2*pi*1j / n)
    even, odd = y[::2], y[1::2]

    even, odd = Inv_FFT(even, False), Inv_FFT(odd, False)

    f = [0] * n
    for i in range(n//2):
        f[i] = even[i] + w**i * odd[i]
        f[i + n//2] = even[i] - w**i * odd[i]

    if original: # 如果回傳最終答案需逐項除以 n
        f = [x/n for x in f]
    return f # 回傳多項式係數 f

就此，我們完成了多項式乘法的優化，將時間複雜度降到了 $O(n\log{n})$，可喜可賀。

:::spoiler 完整的多項式乘法 code

e  = 2.71828182845 # 自然常數
pi = 3.14159265358 # 圓周率

def FFT(f):
    n = len(f) # n 需為 2 的冪次
    if n == 1: # 如果只剩一項就直接回傳
        return f
    
    w = e ** (2*pi*1j / n) # j 為 python 的虛數單位 
    even, odd = f[::2], f[1::2]  # 分成奇偶函數

    even, odd = FFT(even), FFT(odd) # 遞迴計算

    y = [0] * n
    for i in range(n//2): # 合併
        y[i] = even[i] + w**i * odd[i]
        y[i + n//2] = even[i] - w**i * odd[i]

    return y # 回傳解答

def Inv_FFT(y, original = True):# 傳入函數值 y
    n = len(y)
    if n == 1:
        return y
    
    w = e ** (2*pi*1j / n)
    w = 1 / w # 新的 w
    # or you can write w = e ** (-2*pi*1j / n)
    even, odd = y[::2], y[1::2]

    even, odd = Inv_FFT(even, False), Inv_FFT(odd, False)

    f = [0] * n
    for i in range(n//2):
        f[i] = even[i] + w**i * odd[i]
        f[i + n//2] = even[i] - w**i * odd[i]

    if original: # 如果回傳最終答案需逐項除以 n
        f = [x/n for x in f]
    return f # 回傳多項式係數 f

# or you can combine FFT and Inv_FFT in a single function
# def FFT(f, inverse = False, original = True):
#     n = len(f)
#     if n == 1:
#         return f
    
#     w = e ** (2*pi*1j / n)
#     if inverse: w = 1 / w

#     even, odd = f[::2], f[1::2]
#     if inverse: even, odd = FFT(even, True, False), FFT(odd, True, False)
#     else: even, odd = FFT(even), FFT(odd)

#     y = [0] * n
#     for i in range(n//2):
#         y[i] = even[i] + w**i * odd[i]
#         y[i + n//2] = even[i] - w**i * odd[i]

#     if inverse and original:
#             y = [x/n for x in y]
#     return y # 回傳解答


if __name__ == "__main__":

    # 輸入多項式
    a = list(map(complex, input().split()))
    b = list(map(complex, input().split()))

    # 相乘後的項數
    n = len(a) + len(b) - 1
    
    # 將項數補成 2 的冪次個
    x = 1
    while x < n: x *= 2
    a += [complex(0)] * (x - len(a))
    b += [complex(0)] * (x - len(b))
    # print(len(a))

    # FFT 求值
    a = FFT(a)
    b = FFT(b)

    # 多項式乘法
    c = [a[i] * b[i] for i in range(x)]

    # Inverse FFT 求係數
    c = Inv_FFT(c)
    # c = FFT(c, inverse=True)

    # 輸出四捨五入後的整數係數
    for i in c:
        print(round(i.real), end=' ')

:::

## 大數計算

終於，在經過前面的煉獄後，我們終於可以來寫大數乘法的 code 了
既然剛剛實作 FFT 是用 Python，那我們也用 Python 來實作大數乘法吧

根據時間複雜度，我們假設輸入最多有 $10^6$ 位

import sys
sys.set_int_max_str_digits(2000000)
x = int(input())
print(x * x)

好ㄟ，實作完了，那同學們我們這節課就到這
~~python 你好噁（抱頭痛哭~~

## 參考資料

[The Fast Fourier Transform (FFT): Most Ingenious Algorithm Ever?](https://youtu.be/h7apO7q16V0?si=TVuPRhL3wRtNAi_9)
[大数乘法（快速傅立叶变换）上](https://blog.csdn.net/u013351484/article/details/48739415)
[大数乘法（快速傅立叶变换）下](https://blog.csdn.net/u013351484/article/details/48809943)