https://hackmd.io/@peter12345678/BkZ8wIFvO --- # FFT --- ## 多項式乘法 ---- ### 表示方法 ---- ##### 係數表示法$$A(x) = \displaystyle\sum_{i=0}^na_i\,x^i = a_0 + a_1x + a_2x^2 + \cdots + a_{n-1}\,x^{n-1} + a_n\,x^n$$ ---- ##### 點值表示法 * 對於一個 $n$ 次多項式，我們可以隨機選擇 $n+1$ 個**不同的** $x_i$ ，並且對 $A$ 取值。 * 這種表示方法雖然不唯一，但**任何一種點值表示**都能藉由**插值法唯一**確定一個多項式。 $\{(x_i, A(x_i)) \mid 0 \le i \le n, i \in \mathbb Z \}$ --- #### 乘法 給定兩個多項式 $A(x), B(x)$ ##### $A(x) = \displaystyle\sum_{i=0}^na_i\,x^i = a_0 + a_1x + a_2x^2 + \cdots + a_{n-1}\,x^{n-1} + a_n\,x^n$ ##### $B(x) = \displaystyle\sum_{i=0}^nb_i\,x^i = b_0 + b_1x + b_2x^2 + \cdots + b_{n-1}\,x^{n-1} + b_n\,x^n$ ---- ##### 將$A(x), B(x)$這二個多項式相乘得到 $C(x) = \displaystyle\sum_{i=0}^{2n}c_i\,x^i$，且 $c_i = \displaystyle\sum_{j=0}^ia_j\,b_{i-j}\,x^i$ ---- 在係數表示下做乘法，時間複雜度為 $O(n^2)$；但若改為點值表示，我們可以挑出 $2n + 1$ 個點並且直接相乘，時間複雜度為 $O(n)$。 ---- 如果能找出一種快速對多項式插值的方法，我們就能更有效率地計算多項式乘法。 但顯然的，你高中所認識牛頓、拉格朗日插值法，時間複雜度皆為 $O(n^2)$ --- ## 快速傅立葉變換 ## Fast Fourier Transform ---- ### 離散傅立葉變換 ### Discrete Fourier Transform ---- 取樣時選擇**夠好**的點，讓取樣插值都超快!!! ---- 假設現在有一個 $n - 1$ 次多項式 $A(x) = \displaystyle\sum_{i=0}^{n-1} a_ix^i$ 為求方便，令 $n = 2 ^k, k \in Z$，若高次項不足可補 $0$ ---- 重點來了!!! 我們選擇 $n$ 個 $n$ 次單位根 帶入 $A(x)$ 轉換成點值表示! ---- 什麼是單位根? 它是 $x^n = 1$ 的解，共有 $n$ 個複數解，分別為 $e^{2\pi ki/n}, k = 0,1, \dots, n-1$ ---- 令 $\omega_n = e^{2 \pi i/n}$，則 $n$ 個單位根分別為 $\omega_n^0, \omega_n^1, \dots, \omega_n^{n-1}$ 我們將代入的**點值向量** $\vec y = (A(\omega^0_n), A(\omega^1_n), \dots, A(\omega^{n-1}_n))$ 稱作**係數向量**的**離散傅立葉變換**，記作 $\vec y = \text{DFT}_n(\vec a)$ ---- 不過，這樣的變換還是需要 $O(n^2)$ 的時間複雜度 ---- #### Cooley-Tukey Algorithm ---- 我們將點值向量的每一項按照次數奇偶分類 $A(\omega^k_n) = \displaystyle\sum_{i=0}^{n-1}a_i\omega^{ki}_n$ $\hphantom{A(\omega^k_n) }= \displaystyle\sum_{i=0}^{\frac n 2-1}a_{2i}\omega^{2ki}_n + \omega^k_n\displaystyle\sum_{i=0}^{\frac n 2-1}a_{2i+1}\omega^{2ki}_n$ 配合單位根的性質 $\omega^2_n = \omega_{\frac n 2}$、 $\omega^{k + \frac n 2}_n = \omega^k_n \cdot \omega^{\frac n 2 }_n= -\omega^k_n$ 可發現 $k < \frac n 2$ 的時候 $A(\omega^k_n) = \displaystyle\sum_{i=0}^{\frac n 2-1}a_{2i}\,\omega^{ki}_{n/2} + \omega^k_n\displaystyle\sum_{i=0}^{\frac n 2-1}a_{2i+1}\,\omega^{ki}_{n/2}$ $\hphantom{A(\omega^k_n)} = A_1(\omega_{n/2}^k) + \omega^k_n\,A_2(\omega_{n/2}^k)$ $A(\omega^{k+\frac n 2}_n)= \displaystyle\sum_{i=0}^{\frac n 2-1}a_{2i}\,\omega^{(k+\frac n 2)i}_{\frac n 2} + \omega^{k+\frac n 2}_n\displaystyle\sum_{i=0}^{\frac n 2-1}a_{2i+1}\,\omega^{(k+\frac n 2)i}_{\frac n 2}$ $\hphantom{A(\omega^{k+\frac n 2}_n)}= \displaystyle\sum_{i=0}^{\frac n 2-1}a_{2i}\,\omega^{ki}_{\frac n 2} - \omega^k_n\displaystyle\sum_{i=0}^{\frac n 2-1}a_{2i+1}\,\omega^{ki}_{\frac n 2}$ $\hphantom{A(\omega^{k+\frac n 2}_n)} = A_1(\omega_{\frac n 2}^k) - \omega^k_n\,A_2(\omega_{\frac n 2}^k)$ $\vec {a_1}$ 為 $(a_0, a_2, \dots, a_{n-2})$，$\vec {a_2}$ 為 $(a_1, a_3, \dots, a_{n-1})$ 我們可以藉由計算 $\vec {y_1} = \text{DFT}_{\frac n 2}(\vec {a_1})$、$\vec {y_2} = \text{DFT}_{\frac n 2}(\vec {a_2})$ $O(n)$ 合併出 $\vec y = \text{DFT}_n(\vec a)$ 的答案 問題就被分割成二個規模減半的子問題，時間複雜度為 $$T(n) = 2\,\,T{\small \left(\frac n 2\right)} + O(n) = O(n\log n)$$ 跟 merge sort 的複雜度一樣 #### 逆離散傅立葉變換 Inverse Discrete Fourier Transform 這個問題其實是相當於解一個線性方程組 $\begin{cases} \quad a_0\left(\omega_n^0\right)^0 &+ \quad a_1\left(\omega_n^0\right)^1 &+ \quad\cdots &+ \quad a_{n-1}\left(\omega_n^0\right)^{n-1} &= \quad A\left(\omega_n^0\right) \\ \quad a_0\left(\omega_n^1\right)^0 &+ \quad a_1\left(\omega_n^1\right)^1 &+ \quad\cdots &+ \quad a_{n-1}\left(\omega_n^1\right)^{n-1} &= \quad A\left(\omega_n^1\right) \\ &\quad\vdots \\ \quad a_0\left(\omega_n^{n-1}\right)^0 &+ \quad a_1\left(\omega_n^{n-1}\right)^1 &+ \quad\cdots &+ \quad a_{n-1}\left(\omega_n^{n-1}\right)^{n-1} &= \quad A\left(\omega_n^{n-1}\right) \\ \end{cases}$ 寫成矩陣形式的話 $\begin{bmatrix} \left(\omega_n^0\right)^0 & \left(\omega_n^0\right)^1 & \cdots & \left(\omega_n^0\right)^{n-1} \\ \left(\omega_n^1\right)^0 & \left(\omega_n^1\right)^1 & \cdots & \left(\omega_n^1\right)^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ \left(\omega_n^{n-1}\right)^0 & \left(\omega_n^{n-1}\right)^1 & \cdots & \left(\omega_n^{n-1}\right)^{n-1} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} = \begin{bmatrix} A\left(\omega_n^0\right) \\ A\left(\omega_n^1\right) \\ \vdots \\ A\left(\omega_n^{n-1}\right) \end{bmatrix}$ 記上面的係數矩陣為 $V,\ v_{ij} = \omega^{ij}_n$ 構造矩陣 $D, \ d_{ij} = \omega^{-ij}_n$ 令他們相乘的結果為 $E = D \cdot V$ $$e_{ij} = \displaystyle\sum_{k=0}^{n-1} d_{ik}\,v_{kj} \\ \ \ \ \ \ \ \ \ = \displaystyle\sum_{k=0}^{n-1} \omega_n^{-ik}\,\omega_n^{kj} \\ \ \ \ \ \ = \displaystyle\sum_{k=0}^{n-1} \omega_n^{k(j-i)} \\ $$ 當 $i = j$ ，$e_{ij} = n$ 當 $i \ne j$ ，$e_{ij} = {\Large\frac {1 - (\omega_n^{j-i})^n}{1 - \omega_n^{j-i}}} = 0$ 可得到 $I_n = \frac 1 n E, \ \frac 1 n D = V^{-1}$ 所以 $\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} = {\Large\frac{1}{n}} \begin{bmatrix} \left(\omega_n^{-0}\right)^0 & \left(\omega_n^{-0}\right)^1 & \cdots & \left(\omega_n^{-0}\right)^{n-1} \\ \left(\omega_n^{-1}\right)^0 & \left(\omega_n^{-1}\right)^1 & \cdots & \left(\omega_n^{-1}\right)^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ \left(\omega_n^{-(n-1)}\right)^0 & \left(\omega_n^{-(n-1)}\right)^1 & \cdots & \left(\omega_n^{-(n-1)}\right)^{n-1} \\ \end{bmatrix} \begin{bmatrix} A\left(\omega_n^0\right) \\ A\left(\omega_n^1\right) \\ \vdots \\ A\left(\omega_n^{n-1}\right) \end{bmatrix}$ IDFT 就相當於 DFT 過程中 $\omega_n^i$ 換成 $\omega_n^{-i}$，結果再除以 $n$ ### 實作 #### INIT 提前先將所有 $\omega_n^k$ 計算出來 ```cpp #define N 100 #define complex_t complex<double> #define pi 3.1415926535892 complex_t w[N]; void init(int n) { complex_t delta = complex_t(cos(2.0 * pi / n), sin(2.0 * pi / n)); w[0] = 1; for (int i = 1; i < n; i++) w[i] = w[i - 1] * delta; } ``` #### Recursion `fft(n, l, step, flag)` 表示共 $n$ 個數要做 FFT，且這些數分別為 $a_{\small l}, a_{\small l + s}, \dots, a_{\small l+(n-1)s}$ 先遞迴處理 $a_{\small l}, a_{\small l + 2s}, \dots, a_{\small l+(n-2)s}$ 和 $a_{\small l+s}, a_{\small l + 3s}, \dots, a_{\small l+(n-1)s}$ 再 $O(m)$ 列舉 $k$ 計算出 $a_{\small l+ks}, a_{\small l+(k+m)s}$ flag 用來標示做的是 DFT 還是 IDFT $\large \omega_{\normalsize n}^{\normalsize k} = \omega_{\normalsize nN/n}^{\normalsize kN/n} = \omega_{\normalsize N}^{\normalsize kN/n} = \omega_{\normalsize N}^{\normalsize ks}$ ```cpp complex_t a[N], b[N]; void fft(int n, int l, int s, bool flag = true) { if (n == 1) return; int m = n >> 1; fft(m, l, s << 1, flag); fft(m, l + s, s << 1, flag); for (int k = 0; k < n; k++) b[l + s * k] = a[l + s * k]; // buffer for (int k = 0; k < m; k++) { complex_t y = b[l + s * (k << 1)]; complex_t z = (flag ? w[k * s] : conj(w[k * s])) * b[l + s * ((k << 1) + 1)]; a[l + s * k] = y + z; a[l + s * (k + m)] = y - z; } } ``` #### Iterative 觀察遞迴的順序，不斷奇偶分割的情況下，數字會以項數的二進位後綴相同程度分在同一組。 將它們的二進位反轉，會發現每次分在同一組的數字都是連續的一段區間。 也就是說，最底層即是把 $a[i]$ 放在第 $\text{bit_reverse}(i)$ 項。 若事先將它們的順序安排好，那麼我們從下到上直接 merge 上去即可。  實作 $\text{bit_reverse}()$ 的方式有很多種，下面的作法是 $i$ 從最低位 $+1$，$j$ 從最高位做 $+1$ `for (int l = (n >> 1); (j ^= l) < l; l >>= 1);` 從最高位開始向右尋找第一個 $0$，把過程遇到的 $1$ 全部翻成 $0$ 這件事如果從最低位開始往左做，是不是跟 $+1$ 一模一樣呢? 為了只讓它們交換一次，當有一方大於(或小於)另一方時才進行 $\text{swap}$ 從寬度 $2, 4, \dots, n$，慢慢向上迭代，一些神奇的 index 就讓你們自己想怎麼填XD 而改值的順序可以省掉一個 buffer 的陣列! ```cpp void fft(int n, complex_t a[], bool inv=false) { init(n); // bit_reverse() for (int i = 1, j = 0; i < n; i++) { for (int l = (n >> 1); (j ^= l) < l; l >>= 1); if (i > j) swap(a[i], a[j]); } for (int i = 2; i <= n; i <<= 1) { int m = i >> 1; for (int j = 0; j < n; j += i) { for (int k = 0; k < m; k++) { complex_t z = (inv ? conj(w[n / i * k]) : w[n / i * k]) * a[j + k + m]; a[j + k + m] = a[j + k] - z; a[j + k] += z; } } } if (inv) for (int i = 0; i < n; i++) a[i] /= n; } ``` ### 大數乘法實作 ```cpp void str2cpx(string a, c_t b[], int n) { int len = a.size(), j = 0; for (int i = len - 1; i >= 0; i--) b[j++] = a[i] - '0'; while (j < n) b[j++] = 0; } string cpx2str(int n, complex_t a[]) { vector <int> tmp; tmp.resize(n); for (int i = 0; i < n - 1; i++) { int x = round(a[i].real()); a[i + 1] += x / 10; tmp[i] = x % 10; } tmp[n - 1] = round(a[n - 1].real()); while (!tmp.back()) tmp.pop_back(); string s = ""; while (!tmp.empty()) s += (char)(tmp.back() + '0'), tmp.pop_back(); return s; } string mul(string a_, string b_) { int n = 1; while (n < a_.size() + b_.size() - 1) n <<= 1; str2cpx(a_, a, n); str2cpx(b_, b, n); fft(n, a); fft(n, b); for (int i = 0; i < n; i++) a[i] *= b[i]; fft(n, a, 1); return cpx2str(n, a); } ``` ### 數論變換 Number Theory Transform FFT 用到了單位根 $\omega$ 的四個性質 1. $\omega_n^n = 1$ 2. $\omega_n^0, \omega_n^1, \dots, \omega_n^{n-1}$ 互不相同，這樣逆變換才存在 3. $\omega_n^2 = \omega_{(n/2)}, \omega_n^{(n/2) + k} = - \omega_n^k$，這使得奇偶分類後能把代入的值減半，問題規模也減半 4. $\displaystyle\sum_{k = 0}^{n - 1}\left(\omega_n^{j - i}\right)^k = \begin{cases} 0, \quad i \ne j \\ n, \quad i = j \end{cases}$，這點保證可以使用相同的方法做逆變換 FFT 可能出現精度問題，於是我們開始尋找是否有其他好的取樣點