多项式算法

FFT

算法描述，求解对于 $A(x)B(x)$ 的卷积，对于某一个多项式 $A(x)$ 而言
插值多项式的唯一性
对于 $n$ 个点值对组成的集合 $\{(x_0, y_0), (x_1, y_1), \cdots, (x_{n-1}, y_{n-1})\}$
其中所有的 $x_k$ 均不同，那么可以唯一确定次数为 $n$ 的多项式 $A(x)$，满足 $y_k = A(x_k)$

带入多项式，等价于方程

$$\begin{pmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^{n-1} \\ 1 & x_1 & x_1^2 & \cdots & x_1^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n-1} & x_{n-1}^2 & \cdots & x_{n-1}^{n-1} \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{pmatrix} = \begin{pmatrix} y_0 \\ y_1 \\ \vdots \\ y_{n-1} \end{pmatrix}$$

对于 $a_k$ 有唯一解，其中左边的矩阵表示为范德蒙行列式，该行列式值为

$$\prod_{0 \leqslant j < k \leqslant n-1} (x_k - x_j)$$

矩阵可逆，求出唯一的 $a$ 为 $V(x_0, x_1, \cdots, x_{n-1})^{-1} y$

给定 $n$ 次多项式，取 $n+1$ 个点可以将行列式唯一表示出来

根据点表示法，对于 $C(x) = A(x)B(x)$，可以用 $A(x), B(x)$ 的点表示法，假设 $A, B$ 的次数分别为 $n, m$
取 $n+m+1$ 个点

$$\{(x_1, A(x_1)), (x_2, A(x_2)) \cdots, (x_{n+m+1}, A(x_{n+m+1})) \} \\ \{(x_1, B(x_1)), (x_2, B(x_2)) \cdots, (x_{n+m+1}, B(x_{n+m+1})) \} \\ \{(x_1, A(x_1)B(x_1)), (x_2, A(x_2)B(x_2)) \cdots, (x_{n+m+1}, A(x_{n+m+1})B(x_{n+m+1})) \}$$

这样可以在 $O(n+m)$ 的时间复杂度内求出多项式 $C(x)$

普通乘法和点值乘法快速转换

如果能将普通多项式乘法转换成为点值乘法，那么 FFT 是很容易实现的，为此取单位复数根
即满足 $\omega^{n} = 1$ 的复数 $\omega$，$n$ 次单位根恰好有 $n$ 个，对于 $k = 0, 1, \cdots, n-1$
这些根是 $e^{2\pi ik / n}$，利用复数的指数表达，可以写成

$$e^{iu} = \cos(u) + i\sin(u)$$

$n$ 个单位根均匀分布在圆周上，值 $\omega_{n} = e^{2\pi i / n}$ 为主 $n$ 次单位根
所有其他负数根均为 $\omega_{n}$ 的幂次，这些根构成了一个模 $n$ 意义下的乘法群，结构如下
$\omega_{n}^{j} \omega_{n}^{k} = \omega_{n}^{j+k} = \omega_{n}^{(j+k) \bmod n}$，类似地
$\omega_{n}^{-1} = \omega_{n}^{n-1}$

消去引理
对于任意 $n \geqslant 0, k \geqslant 0, d > 0$，我们有

$$\omega_{dn}^{dk} = \omega_{n}^{k}$$

这是显然的，因为 $(e^{2\pi i / dn})^{dk} = (e^{2\pi i / n})^{k} = \omega_{n}^{k}$

那么可以根据消去引理，推出折半引理
对于任意偶数 $n > 0$，有 $\omega_{n}^{n/2} = \omega_{2} = -1$
证明也很简单，$\omega_{n}^{n/2} = \omega_{n/2}^{n/4} = \cdots = \omega_{n/(n/2)}^{n/((n/2) \cdot 2)} = \omega_{2}$

根据折半引理

$$\omega_{n}^{k + n / 2} = - \omega_{n}^{k}$$

求和引理
对于任意整数 $n \geqslant 1$ 和不能被 $n$ 整除的非负整数 $k$，有

$$\sum_{j = 0}^{n-1} (\omega_{n}^{k})^{j} = 0$$

可以用等比数列求和

$$\sum_{j = 0}^{n-1} (\omega_{n}^{k})^{j} = \frac{(\omega_{n}^{k})^{n} - 1}{\omega_{n}^{k} - 1} = \frac{(\omega_{n}^{n})^{k} - 1}{\omega_{n}^{k} - 1} = \frac{(1)^{k} - 1}{\omega_{n}^{k} - 1} = 0$$

DFT

假设 $n$ 为 $2$ 的整数幂
根据 $A(x)$ 中偶数下标和奇数下标的系数，定义最高次数为 $n/2-1$ 的多项式 $A^{[0]}(x)$ 和 $A^{[1]}(x)$

$$A^{[0]}(x) = a_0 + a_2 x + a_4 x^2 + \cdots + a_{n-2}x^{n/2-1} \\ A^{[1]}(x) = a_1 + a_3 x + a_5 x^2 + \cdots + a_{n-1}x^{n/2-1}$$

那么有 $A(x) = A^{[0]}(x^2) + x\cdot A^{[1]}(x^2)$

那么可以采用分治策略，对于 $k \in [0, n/2-1]$，可以得到

$$A(\omega_{n}^{k}) = A_0(\omega_{n}^{2k}) + \omega_{n}^{k} A_1(\omega_{n}^{2k}) = A_0(\omega_{n/2}^{k}) + \omega_{n}^{k} A_1(\omega_{n/2}^{k})$$

而对于 $k \in [n/2, n-1]$，令 $k' \leftarrow k+n/2, k \in [0, n/2-1]$，问题可以转为上述情况，即求 $A(\omega_{n}^{k+n/2})$

$$A(\omega_{n}^{k + n/2}) = A_0(\omega_{n}^{2k+n}) - \omega_n^k A_1(\omega_n^{2k+n}) = A_0(\omega_{n/2}^{k}) - \omega_{n}^{k} A_1(\omega_{n/2}^{k})$$

这里用了折半引理

IDFT

接下来需要根据点表示法，逆变换求出卷积系数
已知点表示 $(\omega_{n}^{k}, A(\omega_{n}^{k}))$，令 $y_k \leftarrow A(\omega_{n}^{k})$
则 $A(x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_{n-1}x^{n-1}$，推导出

$$a_k = \frac{1}{n} \cdot \sum_{i = 0}^{n-1} y_i (\omega_{n}^{-k})^{i}$$

求出 $a_k$，实际上也是在求多项式，构造 $B(x) = y_0 + y_1 x + y_2 x^2 + \cdots y_{n-1} x^{n-1}$
则 $a_k = \displaystyle\frac{1}{n}B(\omega_{n}^{-k})$
只要求出 $B(x)$ 在 $(\omega_{n}^{0}, \omega_{n}^{-1}, \cdots, \omega_{n}^{-(n-1)})$ 上的值
这也可以通过 $\text{DFT}$ 求出

下面给出证明

$$\begin{aligned} c_k &= \sum_{i = 0}^{n-1} y_i (\omega_{n}^{-k})^{i} = \sum_{i = 0}^{n-1} \left(\sum_{j = 0}^{n-1}a_j(\omega_{n}^{i})^{j} \right) \cdot \left(\omega_{n}^{-k}\right)^{i} \\ &= \sum_{i = 0}^{n-1}\left(\sum_{j = 0}^{n-1}a_j (\omega_{n}^{j})^{i} (\omega_{n}^{-k})^{i} \right) = \sum_{i = 0}^{n-1}\left(\sum_{j = 0}^{n-1}a_j (\omega_{n}^{j-k})^{i}\right) \\ &= \sum_{j = 0}^{n-1} a_j \left( \sum_{i = 0}^{n-1} (\omega_{n}^{j-k})^{i} \right) \end{aligned}$$

注意到项 $S(\omega_{n}^{k'}) = \displaystyle \left( \sum_{i = 0}^{n-1} (\omega_{n}^{j-k})^{i} \right)$

其中 $\omega_{n}^{k'} \leftarrow \displaystyle \omega_{n}^{j-k}$，根据求和引理，当
$k' \neq 0$，$S(\omega_{n}^{k'}) = 0$
$k' = 0$ 时，$\omega_{n}^{0} = 1$，那么 $S(1) = n$

也就是对应 $c_k$ 表达式中，$j = k$ 时才有贡献，我们有 $c_k = na_k$，证毕

证明 2，可以考虑用矩阵来证明

根据 $y_k = A(\omega_n^k)$，系数 $a_k$ 对应的项就是 $\omega_n^k$

$$\begin{pmatrix} 1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1} \\ 1 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)} \\ 1 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)} \\ \vdots & \vdots & \ddots & \ddots & \vdots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \\ \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{pmatrix} = \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{pmatrix}$$

不妨设 $\text{DFT}$ 的 $\omega_n^k$ 对应的矩阵是 $V_n$，只需要证明
引理，对于 $j, k \in [0, 1, \cdots, n-1]$，$V_n^{-1}$ 的 $(j, k)$ 处元素为 $\omega_{n}^{-kj} / n$

考虑矩阵乘法，$V_n^{-1}V_n$，对应 $(j, k) \times (k, j') = (j, j')$，考虑 $(j, j')$ 的元素

只要证明，我们有 $(j, k)$ 元素为 $\omega_{n}^{-kj} / n$ 的矩阵 $V'$，能推出该矩阵 $V'$ 为 $V_n^{-1}$ 即可

$\displaystyle [V_n^{-1}V_n]_{jj'} = \sum_{k = 0}^{n-1} \left(\omega_n^{-kj} / n \right) \cdot \omega_{n}^{kj'} = \sum_{k = 0}^{n-1} \omega_n^{k(j - j')} / n$

注意此时 $-(n-1) \leqslant j'-j \leqslant n-1$，可以用求和引理
（说明，如果此时 $-(n-1) \leqslant j'-j < 0$，可以用 $\omega_n^{-p} = \omega_n^{n-p}$ 变换）
如果 $j = j'$，右边等于 $1$，否则根据求和引理，和为 $0$
所以 $V_n^{-1}V_n$ 为单位矩阵 $I_n$，证毕

FFT 实现

蝴蝶变换
输入向量下标为 $(0, 1, 2, 3, 4, 5, 6, 7)$，那么 DFT 递归树叶子结点的下标是 $(0, 4, 2, 6, 1, 5, 3, 7)$
这是根据之前的分析，奇数项写在前面，偶数项写在后面

恰好是二进制表示的位逆序置换，比如 $6 = 110$，位逆序置换成 $011 = 3$

这个很好证明，假设最终系数 $x$ 二进制表示是 $(\cdots \ 1)$，那么它是奇数项，一定在 $[n/2, n-1]$ 这一半的区间
这一半区间中，对应回原来未打乱顺序的数，是较大的一半，最高位都是 $1$
如果二进制表示最低位是 $0$，即 $(\cdots \ 0)$，那么是偶数项
对应回原来的前一半，此时最高位为 $0$

那么核心是位翻转的实现，这可以递归来做

$$r_i = \left(r_{i \gg 1} \gg 1 \right) \mid (i \ \& \ 1) \ll (\text{bit} - 1)$$

$\displaystyle \text{bit}$ 表示最高位的 $1$ 在第几位上

综上所述，$\text{FFT}$ 实现的过程如下

• 预处理出卷积多项式有 $tot$ 项，其中 $tot$ 满足 $2$ 的 $b$ 次幂，不足可以补 $0$
  比如 $A(x)$ 有 $n+1$ 项（最高次幂为 $n$），$B(x)$ 有 $m+1$ 项，那么必须满足 $(1 \ll b) \geqslant m+n+2$
  $tot = (1 \ll b)$，则卷积多项式的系数下标范围是 $[0, tot)$
• $[mid = 1, 2, 2^{2} \cdots, tot]$ 逐层归并计算，对于每一层，分治区间起始点是 $[i = 0, i + 2mid, i + 4mid, \cdots, tot)$
  接下来计算 $\text{DFT}(\bm{v_1})$ 和 $\text{DFT}(\bm{v_2})$，根据如上所述的 $\text{DFT}$ 公式
  用折半引理，$j$ 遍历长度为 $mid$ 的前后两半区间，其中前一半下标是 $i+j$，后一半下标是 $i+j+mid$
  $\bm{v1} = A_{i+j} + \omega_{2mid}^{k} \cdot A_{i+j+mid}$，$\bm{v2} = A_{i+j} - \omega_{2mid}^{k} \cdot A_{i+j+mid}$
  具体来说，在算法实现上
  每一次迭代循环 $\text{mid}$ 的时候，令 $\displaystyle \omega n = \omega_{2mid} = \text{e}^{2\pi i / 2 \text{mid}} = \text{e}^{\pi i / \text{mid}} = \cos(\pi / \text{mid}) - i \cdot \sin(\pi / \text{mid})$
  $\bm{x} = A(i+j), \ \bm{y} = \omega n^{k} \cdot A(i+j+mid)$
  然后令 $A(i+j) \leftarrow \bm{x} + \bm{y}, \ A(i+j+mid) \leftarrow \bm{x} - \bm{y}$，特别注意，这里 $\omega n^{k}$ 是可以实时计算的
  只要在归并的开始，即对 $j$ 循环 $\textbf{for} \ j \in [0, mid)$ 前，令 $\omega n^k = 1$，然后每一次 $j++$ 的时候，令 $\omega n^k \leftarrow \omega n^k \cdot \omega n$