wiki - 算法小站

主题名称

是否为父主题

是否为其添加wiki文章

选择主题

修改主题名称

修改父主题

修改wiki文章

算法基础搜索数据结构字符串数学图论动态规划计算几何杂项

数学

组合计数

二项式定理

生成函数（一）

FFT

[[ item.c ]]

FFT（快速傅立叶变换）

## 背景
关于多项式计数问题，很多时候我们只需要取前$$n+1$$项
```math
A(x) = a_0 + a_1x + a_2x^2 + \cdots + a_nx^n \Rightarrow A(x) \bmod x^{n+1}
```
FFT/NTT 可以在$$O(n\log n)$$的复杂度内，实现两个$$n$$次多项式的乘积$$A(x)B(x)$$

## 多项式的表示形式
假设$$f(x)$$是一个$$n$$次多项式，则$$f(x)$$的系数表示是
```math
\displaystyle
f(x) = a_0 + a_1x + a_2x^2 + \cdots + a_{n-1}x^{n-1} + a_nx^n
```
除此之外，$$f(x)$$还有**点值**表示，记为
```math
(x_0, f(x_0)), (x_1, f(x_1)), \cdots, (x_n, f(x_n))
```
这$$n+1$$个点值就可以表示一个$$n$$次多项式
在点值表示下，$$n$$次多项式乘法的复杂度是$$O(n)$$

```math
h(x) = f(x)g(x) \\
\quad \\
(x_0, f(x_0)) \cdots (x_n, f(x_n)) \\
(x_0, g(x_0)) \cdots (x_n, g(x_n)) \\

\Rightarrow (x_0, h(x_0)) \cdots (x_n, h(x_n)), \quad h(x_i) = f(x_i)g(x_i)
```
那么怎么实现**从点值表示到系数表示的相互转化呢？**，FFT 就可以解决这个问题

## 复数和单位根
复数具有指数形式$$a+bi = re^{i\theta}$$，其中$$r = \sqrt{a^2 + b^2}, \tan(\theta) = \dfrac{b}{a}$$

**单位根**，$$x^n = 1$$在复数域上的根称为$$n$$次单位根，$$n$$次单位根有$$n$$个
具有形式$$\displaystyle \omega_n^k = e^{i\frac{2k \pi}{n}}$$

> 证明
是因为欧拉公式$$e^{ix} = \cos x + i\sin x$$，比如复平面$$(1, 1)$$这个点表示$$1 + i$$
实际上，用极坐标形式可以写成$$\sqrt{2}e^{i \frac{\pi}{4}}$$
将$$x = \frac{\pi}{4}$$带入欧拉公式，$$\sqrt{2}e^{i \frac{\pi}{4}} = 1 + i$$

> 同理
取$$(\omega_n)^n = (e^{i \frac{2\pi}{n}})^n = e^{i \cdot 2\pi}$$
相当于长度为`1`的线段，绕远点转一圈，回到了实轴上，恰好表示`1`

单位根的性质
```math
\displaystyle
\omega_n^k = \omega_{2n}^{2k} \\
\omega_{2n}^{k+n} = -\omega_{2n}^k
```
> 证明
(1)很简单，考虑$$e^{i \frac{2k \pi}{n}}$$，`k, n`同时$$\times 2$$，答案不变

> (2)实际上只需要证明$$\omega\_{2n}^n = -1$$，实际上$$=e^{i\pi} = -1$$
(相当于绕着原点转`180`度，这里逆时针为正方向)

## DFT
将多项式
```math
\displaystyle
A(x) = a_0 + a_1x + \cdots + a_{n-1}x^{n-1}
```
转化成为点值形式$$(\omega\_{n}^k, A(\omega\_n^k)), \quad k = 0, 1, \cdots, n-1$$

考虑$$n = 2^l$$，是`2 的幂`的情况，奇偶分开
```math
\displaystyle
A(x) = (a_0 + a_2x^2 + \cdots + a_{n-2}x^{n-2}) + (a_1x + a_3x^3 + a_5x^5 + \cdots + a_{n-1}x^{n-1}) \\
B(x) = a_0 + a_2x + \cdots + a_{n-2}x^{n/2-1} \\
C(x) = a_1 + a_3x + \cdots + a_{n-1}x^{n/2-1}
```
那么$$A(x) = B(x^2) + xC(x^2)$$

对于$$0 \leqslant k < n/2$$
```math
\displaystyle
A(\omega_n^k) = B(\omega_n^{2k}) + \omega_n^k C(\omega_n^{2k}) = B(\omega_{n/2}^k) + \omega_n^k C(\omega_{n/2}^k)
```
对于$$k \geqslant n/2$$
```math
A(\omega_n^{k+n/2}) = B(\omega_n^{2k+n}) + \omega_n^{k+n/2} C(\omega_n^{2k+n}) \\
= B(\omega_n^{2k}) + \omega_n^{k+n/2} C(\omega_n^{2k}) = B(\omega_{n/2}^{k}) + \omega_n^{k+n/2} C(\omega_{n/2}^{k}) \\
= B(\omega_{n/2}^{k}) - \omega_n^k C(\omega_{n/2}^{k}) 
```
> 以$$n = 8$$为例
![](/media/content_img/dft01.png)

## DFT的蝴蝶变换
**bit-reversal permutation**
![](/media/content_img/dft02.png)
**注意保持`n`是`2 的幂`**，如果不满足，就尾项补`0`直到满足
```bash
x:         [  ......  ] | ?
rev(x):    ? | [ 高位的翻转 ]

把最低位放到最高位 | 剩余的位翻转后右移 1 位以插入最高位

std::vector<int> rev;
if (int(rev.size()) != n) {
    int k = __builtin_ctz(n) - 1;
    rev.resize(n);
    // coef index [a(0)....a(n-1)]
    for (int i = 0; i < n; i++) {
        rev[i] = (rev[i>>1] >> 1) | ((i & 1) << k);
    }
}
for (int i = 0; i < n; i++) {
    if (rev[i] < i) {
        // rev[i] < i 保证只交换一次，避免反复交换
        std::swap(a[i], a[rev[i]]);
    }
}
```
## DFT root 的构造
![](/media/content_img/dft03.png)
![](/media/content_img/dft04.png)
![](/media/content_img/dft05.png)

注意到因为每一层$$j \in [0, half)$$，所以构造$$j' = half + j \in [half, 2half)$$
可以保证这样的下标对应序列中**连续的一段**，`for`$$\forall i \in [2^{k-1}, 2^k)$$
```bash
root[2i] = root[i]
root[2i + 1] = root[i]
```
## DFT 迭代蝶形合并
```bash
for k = half, k <<= 1:
    for i = 起点 in [i, i+2k, ...], 考虑两段 [i...), [i+k...), then i += 2k:
        for j = offset in [0, half) = [0, k):
            auto u = a[i+j], v = [i+k+j] * root[k+j]
            迭代:
            a[i+j] = u + v;
            a[i+k+j] = u - v;
```

## IDFT
将多项式的点值表示$$(\omega_n^k, A(\omega_n^k)=b_k), k = \\{0, 1, \cdots, n-1\\}$$转化为其系数表示
```math
\displaystyle
A(x) = a_0 + a_1x + \cdots + a_{n-1}x^{n-1}
```
实际上
```math
\displaystyle
\begin{pmatrix}
(\omega_n^0)^0 & (\omega_n^0)^1 & \cdots & (\omega_n^0)^{n-1} \\
(\omega_n^1)^0 & (\omega_n^1)^1 & \cdots & (\omega_n^1)^{n-1} \\
\vdots & \vdots \\
(\omega_n^i)^0 & (\omega_n^i)^1 & \cdots & (\omega_n^i)^{n-1} \\
\vdots & \vdots \\
(\omega_n^{n-1})^0 & (\omega_n^{n-1})^1 & \cdots & (\omega_n^{n-1})^{n-1} 
\end{pmatrix}
\begin{pmatrix}
a_0 \\
a_1 \\
\vdots \\
a_i \\
\vdots \\
a_{n-1}
\end{pmatrix}
= \begin{pmatrix}
b_0 \\
b_1 \\
\vdots \\
b_i \\
\vdots \\
b_{n-1}
\end{pmatrix}
```
设$$n \times n$$的矩阵$$\Omega$$，其中$$\Omega\_{i, j} = \omega_n^{ij}$$

设向量$$\bold{a} = (a\_0, a\_1, \cdots, a\_{n-1}), \bold{b} = (b\_0, b\_1, \cdots, b\_{n-1})$$

那么 IDFT 相当于求解方程
```math
\displaystyle
\Omega\bold{a} = \bold{b}
```
> 问题转化成，求$$\Omega^{-1}$$

那么我们有
```math
\displaystyle
\begin{pmatrix}
a_0 \\
a_1 \\
\vdots \\
a_i \\
\vdots \\
a_{n-1}
\end{pmatrix}
= \Omega^{-1}
\begin{pmatrix}
b_0 \\
b_1 \\
\vdots \\
b_i \\
\vdots \\
b_{n-1}
\end{pmatrix}
```
考虑**共轭**
```math
\overline{\Omega} = 
\begin{pmatrix}
(\omega_n^{-0})^0 & (\omega_n^{-0})^1 & \cdots & (\omega_n^{-0})^{n-1} \\
(\omega_n^{-1})^0 & (\omega_n^{-1})^1 & \cdots & (\omega_n^{-1})^{n-1} \\
\vdots & \vdots \\
(\omega_n^{-i})^0 & (\omega_n^{-i})^1 & \cdots & (\omega_n^{-i})^{n-1} \\
\vdots & \vdots \\
(\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 & \cdots & (\omega_n^{-(n-1)})^{n-1} 
\end{pmatrix}
```
那么$$M = \overline{\Omega} \cdot \Omega$$，考虑$$\overline{\Omega}$$第`i`行乘上$$\Omega$$第`j`列，我们有
```math
\displaystyle
M_{i, j} = (\omega_n^{-i})^0 (\omega_n^0)^j + (\omega_n^{-i})^1 (\omega_n^1)^j + \cdots + (\omega_n^{-i})^j (\omega_n^j)^j \\ + \cdots + (\omega_n^{-i})^{n-1} (\omega_n^{n-1})^j  \\
\quad \\
M_{i, j} = (\omega_n^{j-i})^0 + (\omega_n^{j-i})^1 + \cdots + (\omega_n^{j-i})^{n-1} \\
\quad \\
M_{i, j} = \begin{cases} n & w_n^{j - i} = 1 \\
\dfrac{1-(w_n^{j-i})^n}{1- w_n^{j-i}} & w_n^{j - i} \neq 1
\end{cases}
```
注意到$$(w_n^{j-i})^n = (w_n^n)^{j - i} = 1$$，如果$$w_n^{j - i} \neq 1$$

有$$M\_{i, j} = 0$$，什么时候$$w\_n^{j - i} = 1$$呢？只有$$j = i$$

因为$$0 \leqslant i, j < n$$，只有$$\omega\_n^{kn} = 1$$，所以只能是$$j = i$$时候，$$M\_{i, j} = 1$$

那么$$\overline{\Omega} \cdot \Omega = n \bold{I}$$，也就是说$$\Omega^{-1} = \dfrac{1}{n} \overline{\Omega}$$

**IDFT**
令$$\overline{\Omega}$$满足$$\overline{\Omega}\_{i, j} = \omega\_n^{-ij}$$，我们有$$\overline{\Omega}\Omega = nI$$，因此
```math
\displaystyle
\bold{a} = \dfrac{1}{n} \overline{\Omega}\bold{b}
```
相当于给定
```math
B(x) = b_0 + b_1x + \cdots + b_{n-1}x^{n-1}
```
求点值$$B(\omega\_n^{-k}), 0 \leqslant k < n$$
```math
\displaystyle
\begin{pmatrix}
b_0(\omega_n^{-0})^0 + b_1(\omega_n^{-0})^1 +\cdots + b_{n-1}(\omega_n^{-0})^{n-1} \\
b_0(\omega_n^{-1})^0 + b_1(\omega_n^{-1})^1 + \cdots + b_{n-1}(\omega_n^{-1})^{n-1} \\
\vdots \\
b_0(\omega_n^{-i})^0 + b_1(\omega_n^{-i})^1 + \cdots + b_{n-1}(\omega_n^{-i})^{n-1} \\
\vdots \\
b_0(\omega_n^{-(n-1)})^0 + b_1(\omega_n^{-(n-1)})^1 + \cdots + b_{n-1}(\omega_n^{-(n-1)})^{n-1} 
\end{pmatrix} = 
\begin{pmatrix}
B(\omega_n^{-0}) \\
B(\omega_n^{-1}) \\
\vdots \\
B(\omega_n^{-i}) \\
\vdots \\
B(\omega_n^{-(n-1)})
\end{pmatrix}
```
第$$k$$行
```math
\displaystyle
b_0(\omega_n^{-k})^0 + b_1(\omega_n^{-k})^1 + \cdots + b_{n-1}(\omega_n^{-k})^{n-1} = a_k \\
\quad \\
b_j\omega_n^{-kj} = b_j \omega_n^{nk - kj} = b_j (\omega_n^{n - j})^k = b_j (\omega_n^k)^{n - j} \\
\quad \\
LHS(k)  = b_0 + \sum_{j = 1}^{n-1} b_j(\omega_n^k)^{n - j}
```
令$$j = n - j$$，上式就可以写成$$b\_{n-j}(\omega_n^k)^j$$
```math
\displaystyle
LHS(k) = b_0 + \sum_{j = 1}^{n-1}b_{n-j}(\omega_n^k)^j
```
我们原先执行`DFT`之后求出的系数`a[0...n-1]`
只要令$$b\_0 = a\_0, \quad b\_{j} = a\_{n - j}$$

**这实际上是在执行**
```bash
reverse(a.begin()+1, a.end())
```

## IDFT 的实现以及多项式卷积
**IDFT**
```bash
inline void idft(std::vector<cd> &a) {
    int n = a.size();
    if (n == 0) return;

reverse(a.begin()+1, a.end());
    dft(a);

double invn = 1.0 / n;
    for (int i = 0; i < n; i++) a[i] *= invn;
}
```
**多项式卷积**
```bash
struct Poly {
    std::vector<cd> a;
	
    friend Poly operator* (Poly A, Poly B) {
        if (A.size() == 0 or B.size() == 0) return Poly();
        if (A.size() < B.size()) std::swap(A, B);

if (B.size() < 128) {
            Poly C(A.size() + B.size() - 1);
            for (int i = 0; i < A.size(); i++) for (int j = 0; j < B.size(); j++) {
                C[i+j] += A[i] * B[j];
            }
            return C;
        }

int sz = 1, tot = A.size() + B.size() - 1;
        while (sz < tot) sz <<= 1;

A.a.resize(sz), B.a.resize(sz);

dft(A.a), dft(B.a);
        for (int i = 0; i < sz; i++) A.a[i] *= B.a[i];
        idft(A.a);

A.resize(tot);
        return A;
    }
};
```

看完文章有啥想法

算法小站

wiki 主题

数学

FFT（快速傅立叶变换）

目录