---
# System prepended metadata

title: Fourier Transform
tags: [講義]

---


## 多項式乘法

### 表示方法

- 係數表示法
- 
    $A(x) = \displaystyle\sum_{i = 0}^{n} a_i\,x^i = a_0 + a_1\,x + a_2\,x^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) = \displaystyle\sum_{i = 0}^{n} a_i\,x^i = a_0 + a_1\,x + a_2\,x^2 + \cdots + a_{n - 1}\,x^{n - 1} + a_n\,x^n$
$B(x) = \displaystyle\sum_{i = 0}^{n} b_i\,x^i = b_0 + b_1\,x + b_2\,x^2 + \cdots + b_{n - 1}\,x^{n - 1} + b_n\,x^n$

將這二個多項式相乘得到 $C(x) = \displaystyle\sum_{i = 0}^{2n} c_i\,x^i$，且 $c_i = \displaystyle\sum_{j = 0}^{i} a_j\,b_{i - j}\,x^i$

在係數表示下做乘法，時間複雜度為 $O(n^2)$；但若改為點值表示，我們可以挑出 $2n + 1$ 個點並且直接相乘，時間複雜度為 $O(n)$。
如果能找出一種快速對多項式插值的方法，我們就能更有效率地計算多項式乘法。
但高中所認識的牛頓、拉格朗日插值法，時間複雜度皆為 $O(n^2)$.

----

## 離散傅立葉變換 Discrete Fourier Transform, DFT

取樣時選擇夠好的點，讓取樣插值都超快!!!

假設現在有一個 $n - 1$ 次多項式
$A(x) = \displaystyle\sum_{i = 0}^{n - 1} a_i\,x^i$

為求方便，令 $n = 2^k, k \in \mathbb{Z}$，若高次項不足可補 $0$

我們選擇 $n$ 個 $n$ 次單位根 帶入 $A(x)$ 轉換成點值表示!

什麼是單位根? 它是 $x^n = 1$ 的解，共有 $n$  個複數解，分別為 $e^{2 \pi ki/n}, k = 0,1, \dots, n - 1$
令 $ω_n = e^{2 \pi i/n}$，則 $n$ 個單位根分別為 $ω_n^0, ω_n^1, \dots, ω_n^{n - 1}$

我們將代入的點值向量 $\vec{y} = (A(ω^0_n), A(ω^1_n), \dots, A(ω^{n - 1}_n))$ 稱作**係數向量**的**離散傅立葉變換**，記作 $\vec{y} = \text{DFT}_n(\vec{a})$

不過，這樣的變換還是需要 $O(n^2)$ 的時間複雜度。

----

## 快速傅立葉變換 Fast Fourier Transform, FFT

我們將點值向量的每一項按照次數奇偶分類

$\begin{split}
A(ω^k_n) &= \displaystyle\sum_{i = 0}^{n - 1} a_i\,ω^{ki}_n \\
&= \displaystyle\sum_{i = 0}^{\frac{n}{2} - 1} a_{2i}\,ω^{2ki}_n \,+\, ω^k_n\displaystyle\sum_{i = 0}^{\frac{n}{2} - 1}a_{2i + 1}\,ω^{2ki}_n
\end{split}$

配合單位根的性質 $ω^2_n = ω_{n/2}$、 $ω^{k + \frac{n}{2}}_n = ω^k_n \cdot ω^{n/2}_n= -\,ω^k_n$ 可發現 $k \lt \frac{n}{2}$ 的時候

$\begin{split}
A(ω^k_n) &= \displaystyle\sum_{i = 0}^{\frac{n}{2} - 1} a_{2i}\,ω^{ki}_{n / 2}\,+\,ω^k_n\displaystyle\sum_{i = 0}^{\frac{n}{2} - 1} a_{2i + 1}\,ω^{ki}_{n / 2} \\
&= A_1(ω_{n / 2}^k)\,+\,ω^k_n\,A_2\,(ω_{n / 2}^k)
\end{split}$

$\begin{split}
A(ω^{k + \frac{n}{2}}_n) &= \displaystyle\sum_{i = 0}^{\frac{n}{2} - 1} a_{2i}\,ω^{(k + \frac{n}{2}) i}_{n / 2}\,+\,ω^{k+\frac{n}{2}}_n\displaystyle\sum_{i = 0}^{\frac{n}{2} - 1} a_{2i + 1}\,ω^{(k + \frac{n}{2}) i}_{n / 2} \\
&= \displaystyle\sum_{i = 0}^{\frac{n}{2} - 1} a_{2i}\,ω^{ki}_{n / 2}\,-\,ω^k_n\displaystyle\sum_{i = 0}^{\frac{n}{2} - 1} a_{2i + 1}\,ω^{ki}_{n / 2} \\
&= A_1\,(ω_{n / 2}^k)\,-\,ω^k_n\,A_2\,(ω_{n / 2}^k)
\end{split}$

\* $\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{(\frac{n}{2})} + O(n) = O(n \log n)$

----
## 逆離散傅立葉變換 Inverse Discrete Fourier Transform

這個問題其實是相當於解一個線性方程組

$\begin{cases}
\quad{a_0\,(ω_n^0)^0} &+ \quad{a_1\,(ω_n^0)^1} &+ \quad{\cdots} &+ \quad{a_{n - 1}\,(ω_n^0)^{n - 1}} &= \quad{A(ω_n^0)} \\
\quad{a_0\,(ω_n^1)^0} &+ \quad{a_1\,(ω_n^1)^1} &+ \quad{\cdots} &+ \quad{a_{n - 1}\,(ω_n^1)^{n - 1}} &= \quad{A(ω_n^1)} \\
&\quad{\vdots} \\
\quad{a_0\,(ω_n^{n - 1})^0} &+ \quad{a_1\,(ω_n^{n - 1})^1} &+ \quad{\cdots} &+ \quad{a_{n - 1}\,(ω_n^{n - 1})^{n - 1}} &= \quad{A\,(ω_n^{n - 1})} \\
\end{cases}$

矩陣形式：

$\begin{equation}
    \begin{bmatrix}
        (ω_n^0)^0 &(ω_n^0)^1 &\cdots &(ω_n^0)^{n - 1} \\
        (ω_n^1)^0 &(ω_n^1)^1 &\cdots &(ω_n^1)^{n - 1} \\
        \vdots &\vdots &\ddots &\vdots \\
        (ω_n^{n - 1})^0 &(ω_n^{n - 1})^1 &\cdots &(ω_n^{n - 1})^{n - 1}
    \end{bmatrix}
    \begin{bmatrix} 
        a_0 \\
        a_1 \\
        \vdots \\
        a_{n - 1}
    \end{bmatrix}
\,=\,
    \begin{bmatrix} 
        A\,(ω_n^0) \\
        A\,(ω_n^1) \\
        \vdots \\
        A\,(ω_n^{n - 1}) 
    \end{bmatrix}
\end{equation}$

記上面的係數矩陣為 $V, \,v_{ij} = ω^{ij}_n$

構造矩陣 $D, \,d_{ij} = ω^{-ij}_n$

令他們相乘的結果為 $E = D \cdot V$

$\begin{split}
e_{ij} &= \displaystyle\sum_{k = 0}^{n - 1} d_{ik}\,v_{kj} \\
		&= \displaystyle\sum_{k = 0}^{n - 1} ω_n^{-ik}\,ω_n^{kj} \\
		&= \displaystyle\sum_{k = 0}^{n - 1} ω_n^{k(j - i)} \\
\end{split}$

當 $i = j$ ，$e_{ij} = n$

當 $i \ne j$ ，$e_{ij} = {\Large\frac{1 \,-\, (ω_n^{j - i})^n}{1 \,-\, ω_n^{j - i}}} = 0$

可得到 $I_n = \frac{1}{n}\,E, \,\frac{1}{n}\,D = V^{-1}$

所以

$\begin{equation}
    \begin{bmatrix}
        a_0 \\
        a_1 \\
        \vdots \\
        a_{n - 1}
    \end{bmatrix}
\,=\,
    {\Large\frac{1}{n}}
    \begin{bmatrix}
        (ω_n^{-0})^0 &(ω_n^{-0})^1 &\cdots &(ω_n^{-0})^{n - 1} \\
        (ω_n^{-1})^0 &(ω_n^{-1})^1 &\cdots &(ω_n^{-1})^{n - 1} \\
        \vdots &\vdots &\ddots &\vdots \\
        (ω_n^{-(n - 1)})^0 &(ω_n^{-(n - 1)})^1 &\cdots &(ω_n^{-(n - 1)})^{n - 1} \\
    \end{bmatrix}
    \begin{bmatrix}
    A\,(ω_n^0) \\
    A\,(ω_n^1) \\
    \vdots \\
    A\,(ω_n^{n - 1})
    \end{bmatrix}
\end{equation}$

IDFT 就相當於將 DFT 過程中 $ω_n^i$  換成 $ω_n^{-i}$，結果再除以 $n$

----
## 實作

### INIT

提前先將所有 $ω_n^k$ 計算出來

```cpp
#define N 100
#define complex_t complex&lt;double&gt;
#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 &lt; n; i++)
	    w[i] = w[i - 1] * delta;
}
```

----

### Recursion

`fft(n, l, step, flag)` 表示共 $n$ 個數要做 FFT，且這些數分別為 $a_l, a_{l + s}, \dots, a_{l + (n - 1)s}$

先遞迴處理 $a_l, a_{l + 2s}, \dots, a_{l + (n - 2)s}$ 和 $a_{l + s}, a_{l + 3s}, \dots, a_{l + (n - 1)s}$

再 $O(m)$ 列舉 $k$ 計算出 $a_{l + ks}, a_{l + (k + m)s}$

flag 用來標示做的是 DFT 還是 IDFT

$ω_n^k = ω_{nN/n}^{kN/n} = ω_N^{kN/n} = ω_N^{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 上去即可。

![](https://i.imgur.com/fkOOO7J.png)

實作 $\text{bit_reverse()}$ 的方式有很多種，下面的作法是 $i$ 從最低位 $+1$，$j$ 從最高位做 $+1$
```cpp
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 用到了單位根 $ω$ 的四個性質

1. $ω_n^n = 1$
2. $ω_n^0, ω_n^1, \dots, ω_n^{n - 1}$ 互不相同，這樣逆變換才存在
3. $ω_n^2 = ω_{n / 2}, ω_n^{k + \frac{n}{2}} = - ω_n^k$，這使得奇偶分類後能把代入的值減半，問題規模也減半
4. $\begin{equation}
	\displaystyle\sum_{k = 0}^{n - 1} (ω_n^{j - i})^k=
	\begin{cases}
		0, \quad i \ne j \\
		n, \quad i = j 
	\end{cases}
\end{equation}$，這點保證可以使用相同的方法做逆變換

FFT 可能出現精度問題，於是我們開始尋找是否有其他好的取樣點

----

### 原根

根據費馬小定理，對於一個質數 $p$

$\begin{equation}
    a^{p - 1} \equiv 1 \pmod{p}
\end{equation}$

我們定義 $p$ 的原根 $g$ 為滿足 $g^0, g^1, \dots, g^{p - 2} \pmod p$ 互不相同的數

如果我們取質數 $p = r \cdot 2^k + 1$，並找到它的原根 $g$，然後我們令 $g_n^k \equiv g^k \pmod p$，這樣就可以使得 $g_n^0, g_n^1, \dots, g_n^{n - 1} \pmod p$ 互不相同，並且 $g_n^n \equiv 1 \pmod p$，滿足性質一跟性質二。

由於 $p$ 是質數，這樣 $g_{\,n}^{n / 2} \mod p$ 必然是 $1$ 或 $-1$，根據 $g^k_n$ 互不相同，$g_{\,n}^{n / 2} \equiv -1 \pmod p$，滿足性質三。

性質四可自行驗證。

於是，原根可以在模這種特殊形式時直接替代單位根，只需要將複數運算改為剩餘系下的運算。

以下提供一些常用的 $p, g$

| p | r | k | g |
|---|---|---|---|
| 998244353 | 119 | 23 | 3 |
| 1004535809 | 479 | 21 | 3 |
| 2281701377 | 17 | 27 | 3 |
