## 多項式乘法
### 表示方法
- 係數表示法
-
$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<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_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 上去即可。
實作 $\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 |
Sign in with Wallet
Connect another wallet