炫酷数学魔术

壹：FFT, NTT

1. 域

(1) 定义

一个有两个运算的群，两个运算都有交换律，且一个运算对另一个运算有分配律。

(2) 域上多项式

在域 \(F\) 中形如 \(a_0 + a_1x + a_2x^2 + ... + a_{n - 1}x^{n -1}\) 的多项式，其中 \(a_0, a_1, ..., a_{n - 1} \in F\)。

2. FFT

(1) 引入

对于 \(n\) 次多项式 \(f\) 和 \(m\) 次多项式 \(g\)，如果我们要计算 \(f(x)g(x)\)，常规高精度是 \(\mathcal{O}(nm)\) 的。

首先，假设我们已经知道了对于一个 \(k\) 次多项式 \(h\) 在 \(k\) 个不同的 \(x\) 的 \(h(x)\)，那么我们就可以确定这个 \(h\)（高斯消元解矩阵）。

不过这样时间复杂度依然不变。我们的 \(x\) 可以取一些特殊值，这里我们先把 \(x\) 取为 \(\omega^0, \omega^1, ..., \omega^{k - 1}\)，其中 \(\omega\) 是 \(F\) 下的 \(k\) 次本原单位根。设矩阵 \(P\) 是 \(x\) 组成的 \(k \times k\) 的矩阵，\(Q\) 是长度为 \(k\) 的 \(h(x)\) 的列向量，则根据定义，\(P\cdot h = Q\) 如下：

\[\def\w{\omega} \begin{bmatrix} 1 & 1 & 1 &... & 1 \\ 1 & \w & \w ^ 2 & ... & \w ^ {k - 1} \\ 1 & \w ^ 2 & \w ^ 4 & ... & \w ^ {2(k - 1)} \\ 1 & \w ^ 3 & \w ^ 6 & ... & \w ^ {3(k - 1)} \\ ... \\ 1 & \w ^ {(k-1)} & \w ^ {2(k - 1)} & ... & \w ^ {(k - 1)(k - 1)} \\ \end{bmatrix} \begin{bmatrix} h_0 \\ h_1 \\ h_2 \\ h_3 \\ ... \\ h_{k - 1} \\ \end{bmatrix} = \begin{bmatrix} Q_0 \\ Q_1 \\ Q_2 \\ Q_3 \\ ... \\ Q_{k - 1} \\ \end{bmatrix} \]

显然，此时 \(h = Q \cdot P ^ {-1}\)，这时，\(P ^ {-1}\) 刚好等于：

\[\def\w{\omega} \begin{bmatrix} 1 & 1 & 1 &... & 1 \\ 1 & \w^{-1} & \w ^ {-2} & ... & \w ^ {-(k - 1)} \\ 1 & \w ^ {-2} & \w ^ {-4} & ... & \w ^ {-2(k - 1)} \\ 1 & \w ^ {-3} & \w ^ {-6} & ... & \w ^ {-3(k - 1)} \\ ... \\ 1 & \w ^ {-(k-1)} & \w ^ {-2(k - 1)} & ... & \w ^ {-(k - 1)(k - 1)} \\ \end{bmatrix} \]

这就是我们要选 \( \def\w{\omega} \w\) 为 \(x\) 的第一个原因。

(2) 整理

刚才已经强调过了，直接这么做时间复杂度依然不变。

准备了这么久，回到正题上，我们该怎么利用刚才的结论求出 \(f(x)g(x)\)？

顺序是这样的：首先，我们分别求出 \(Q_f\) 和 \(Q_g\)，然后乘起来就是 \(Q_{fg}\)（注意要带 \(n + m\) 个值才能解出来 \(fg\)）。接下来，我们再把求出来的 \(Q_{fg}\) 乘上 \(P^{-1}\)，就得到了 \(fg\)。

但现在的时间复杂度瓶颈就是如何通过 \(P, h\) 求出 \(Q\)（\(h\) 指系数矩阵）。直接矩乘时间复杂度肯定会炸，因此要找到更优的方法来求，FFT 就是优化了这一部分，使它能在 \(\mathcal{O}(n \log n)\) 内求出。

(3) 探究

我们现在唯一有的性质就是 \(P\) 它是和 \( \def\w{\omega} \w\) 相关的矩阵。

显然 \(Q_i = \sum\limits_{j = 0} ^ {k - 1} h_j\cdot x_i^j = h_0 + h_1x_i + h_2x_i^2 + ... + h_{k - 1}x_i^{k - 1}\)。

考虑奇偶分类，则有：

\[\begin{aligned} Q_i &= h_0 + h_2x_i^2 + h_4x_i^4 + ... + h_1x_i + h_3x_i^3 + h_5x^5 +...\\ &= h_0 + h_2x_i^2 + h_4x_i^4 + ... + x_i(h_1 + h_3x_i^2 + h_5x_i^4 + ...) \end{aligned} \]

观察到左右两边很相似，只有系数不一样。设 \(F_1 = h_0 + h_2x_i + h_4x_i^2 + ..., F_2 = h_1 + h_3x_i + h_5x_i^2 + ...\)，则 \(Q_i = F_1(x_i ^ 2) + x_iF_2(x_i ^ 2)\)。

\(F_1\), \(F_2\) 和 \(Q_i\) 一模一样。因此分成了两个子问题，现在我们利用分治可以在 \(\mathcal{O}(2n)\) 的时间复杂度求出一个 \(Q_i\) 了……？好像时间复杂度更差了。

先别急，我们要利用 \(x = \omega_k\) 这个性质来联系起不同的 \(Q_i\)。为了保证思维连贯性，这里直接说结论，推式子可以得到 \(Q_i = F_1(x_i ^ 2) + x_iF_2(x_i ^ 2)\)，则 \(Q_{i + \frac{k}{2}} = F_1(x_i ^ 2) - x_iF_2(x_i ^ 2)\)。因此我们只要知道左半边的 \(Q_i\)，就可以知道右半边的 \(Q_{i + \frac{n}{2}}\)。这就是选 \(x = \omega\) 的第二个原因。

也就是说，我们只需要求出 \(\frac{n}{2}\) 个取值的 \(F_1\) 和 \(F_2\)，就可以求出整个 \(Q\)。同理，对于 \(\frac{n}{2}\) 个取值的 \(F_1\) 来说，我们只需要求出 \(\frac{n}{4}\) 个取值的 \(F_{1_1}\) 和 \(F_{1_2}\) 就可以了（\(F_2\) 同理）。

因此我们的分治结构确定了，先把 \(h\) 奇偶分类，然后分治求左右半边 \(\frac{n}{2}\) 个取值，然后再 \(O(n)\) 合并。

总时间复杂度是 \(T(n) = 2T(\frac{n}{2}) + \mathcal{O}(n)\)。

还没结束，我们现在可以在 \(\mathcal{O}(n \log n)\) 内求出一个 \(Q\)，但是还不知道怎么把 \(Q\) 和 \(P^{-1}\) 乘到一起还原。其实答案已经出来了， \(Q \cdot P^{-1}\) 本质上和 \(Q \cdot P\) 类似，只不过是由 \(x = \omega\) 变成了 \(x = \omega^{-1}\)，因此也可以用 FFT 来求。

至此，FFT 已经结束了。不过还有一个问题没解决，就是怎么求 \(\omega\)。直接给公式：\(\omega = \cos(\frac{2\pi}{n}) + i\sin(\frac{2\pi}{n})\)。因此 FFT 需要复数运算。

(4) 回顾

总体来说 FFT 就是把 \(fg\) 从数值表示法到点值表示法再到数值表示法，利用分治来优化原本需要矩乘的东西。一个细节是之前的结论需要建立在 \(n\) 为 \(2\) 的整数幂，因此要提前把数组补 \(0\)。

此外 FFT 还有几个优化。

1. 迭代实现

我们只需要在 FFT 之前把每个元素移动到正确的位置，就可以直接从底到顶递推合并去求。还是直接给结论，假设 \(n = 2 ^ l\)，每个元素最终移到的位置 \(r_k = \lfloor \frac{r_{\lfloor \frac{k}{2} \rfloor}}{2} \rfloor + (k \bmod 2) \times 2^{l -1}\)。

2. 三次变两次

把 \(g\) 放到 \(f\) 的虚部，做完一遍 FFT 后乘法变成平方，最后答案就是虚部除以二。正确性易证明：\((a + bi) ^ 2 = (a ^ 2 - b ^ 2) + 2abi\)

(5) NTT

NTT 就是 \(\mathbb{F_p}\) 下的 FFT。因此不需要用虚数来运算，整数运算就可以了。\(998244353\) 的原根是 \(3\)。

贰：乘法逆元

1. 定义

若 \(aa^{-1} \equiv 1\pmod p\)，则称 \(a^{-1}\) 是 \(a\) 的逆元。
注：若 \(p | a\)，则不存在 \(a ^ {-1}\)。

2. 求解方法

(1) 快速幂

快速幂适用于 \(p\) 是质数的情况。根据费马小定理，\(a^p \equiv a \pmod p\)，则 \(a \cdot a ^ {p - 2} \equiv 1 \pmod p\)。因此此时 \(a\) 的逆元为 \(a ^ {p - 2}\)。

(2) exgcd

exgcd 适用于 \(p\) 和 \(a\) 互质的情况，根据定义 \(a \cdot a ^ {-1} \equiv 1\pmod p\)，则 \(a \cdot a ^ {-1} - py = 1\)。显然有解的情况是 \(\gcd(a, p) = 1\)。在这种情况下就可以解出 \(a ^ {-1}\) 和 \(y\)。

(3) 线性求解 \(1 \sim n\) 的逆元

适用于 \(p\) 是质数的情况。设 \(p = ki + j\)，其中 \(k = \lfloor \frac{p}{i} \rfloor, j = p \bmod i\)。则有 \(ki + j \equiv 0 \pmod p\)，同时乘 \(i^{-1}j^{-1}\)，变为 \(kj^{-1} + i^{-1} \equiv 0 \pmod p\)，即 \(i^{-1} \equiv -\lfloor \frac{p}{i} \rfloor (p \bmod i) ^ {-1}\)。

实际代码可写成：

inv[1] = 1;
for (int i = 2; i <= n; ++i) {
  inv[i] = (long long)(p - p / i) * inv[p % i] % p;
}

(4) 线性求任意 \(n\) 个数的逆元

\(s_i\) 表示前 \(i\) 个数的前缀积，先求解出 \(s_n\) 的逆元 \(sv_n\)，则 \(sv_n \cdot a_n = sv_{n - 1}\)。依次求解出 \(sv_{1 \sim n}\)。则 \(a_{i}^{-1} = sv_{i} \cdot s_{i - 1}\)。