开头先扔板子：多项式板子们
定义
多项式（polynomial）是形如 \(P(x) = \sum \limits_{i = 0}^{n} a_i x ^ i\) 的代数表达式。其中 \(x\) 是一个不定元。
\(\partial(P(x))\) 称为这个多项式的次数。
多项式的基本运算
- 多项式的加减法
- 多项式的乘法
- 多项式除法
这里讨论多项式的带余除法。
可以证明，一定存在唯一的 \(q(x), r(x)\) 满足 \(A(x) = q(x) B(x) + r(x)\)，且 \(\partial(r(x)) < \partial(B(x))\)。
\(q(x)\) 称为 \(B(x)\) 除 \(A(x)\) 的商式，\(r(x)\) 称为 \(B(x)\) 除 \(A(x)\) 的余式。记作：
\[A(x) \equiv r(x) (\bmod \ B(x)) \]特别的，若 \(r(x) = 0\)，则称 \(B(x)\) 整除 \(A(x)\)，\(A(x)\) 是 \(B(x)\) 的一个倍式，\(B(x)\) 是 \(A(x)\) 的一个因式。记作 \(B(x) | A(x)\)。
有关多项式的引理
-
对于 \(n + 1\) 个点可以唯一确定一个 \(n\) 次多项式。
-
如果 \(f(x), g(x)\) 都是不超过 \(n\) 次的多项式，且它对于 \(n +1\) 个不同的数 \(\alpha_1, \alpha_2 \cdots \alpha_n\) 有相同的值，即 \(\forall i \in [1, n + 1], i \in \mathbb{Z}, f(\alpha_i) = g(\alpha_i)\)。则 \(f(x) = g(x)\)。
多项式的点值表示
如果选取 \(n + 1\) 个不同的数 \(x_0, x_1, \cdots x_n\) 对多项式进行求值，得到 \(A(x_0), A(x_1) \cdots A(x_n)\)，那么就称 \(\{(x_i, A(x_i)) \ | \ 0 \le i \le n, i \in \mathbb{Z} \}\) 为 \(A(x)\) 的点值表示。
快速傅里叶变换（FFT）
快速傅里叶变换是借助单位根进行求值和插值，从而快速进行多项式乘法的算法。
单位根
将复平面上的单位圆均分成 \(n\) 份，从 \(x\) 轴数，第 \(i\) 条线与单位圆的交点称为 \(i\) 次单位根，记作 \(\omega_{n}^{i}\)。
根据定义，可以得到：\(\omega_{n}^{1} = i \sin \alpha +\cos \alpha, \alpha = \frac{2 \pi}{n}\)。
根据欧拉恒等式，可以得到：\(\omega_{n}^{1} = e ^ {\frac{2 \pi i}{n}}\)。
由此那么可以得到下面的性质：
-
\(\omega_{n}^{i} \times \omega_{n}^{j} = \omega_{n}^{i + j}\)
-
\(\omega_{n}^{i + \frac{n}{2}} = - \omega_{n}^{i}\)
-
\(\omega_{n}^{i + n} = \omega_{n}^{i}\)
离散傅里叶变换（DFT）
离散傅里叶变换，是将 \(\omega_{n}^{k}, 0 \le k < n\) 代入 \(f(x)\) 和 \(g(x)\) 中求值的过程。
对于朴素的方法，每次代入一个单位根，然后用 \(O(n)\) 的复杂度计算函数值。时间复杂度 \(O(n ^ 2)\)。
离散傅里叶变换利用了单位根的性质巧妙优化了这个过程。离散傅里叶变换过程如下：
首先将 \(f(x)\) 根据次数的奇偶性拆成两部分，分别分为：
\[\begin{cases} e(x) = \sum \limits_{i = 0}^{\frac{n - 2}{2}} a_{2i} x ^ {2i} = a_0 + a_2 x^2 + a_4 x^4 \cdots a_{n - 2} x^{n - 2} \\ o(x) = \sum \limits_{i = 0}^{\frac{n - 2}{2}} a_{2i + 1} x ^ {2i + 1} = a_1 x + a_3 x^3 + a_5 x^5 \cdots a_{n - 1} x^{n - 1} \end{cases} \]设
\[\begin{cases} e'(x) = \sum \limits_{i = 0}^{\frac{n - 2}{2}} a_{2i} x ^ i = a_0 + a_2 x + a_4 x^2 \cdots a_{n - 2} x^{\frac{n - 2}{2}} \\ o'(x) = \sum \limits_{i = 0}^{\frac{n - 2}{2}} a_{2i + 1} x ^ {i} = a_1 + a_3 x^1 + a_5 x^2 \cdots a_{n - 1} x^{\frac{n - 2}{2}} \end{cases} \]则 \(f(x) = e'(x ^ 2) + x o'(x ^ 2)\)。
将 \(\omega_{n}^{k}\) 代入得到：
\[\begin{cases} f(\omega_{n}^{k}) = e'((\omega_{n}^{k}) ^ 2) + \omega_{n}^{k} o'((\omega_{n}^{k}) ^ 2) = e'(\omega_{n}^{2k}) + \omega_{n}^{k} o'(\omega_{n}^{2k}) \\\\ f(\omega_{n}^{k + \frac{n}{2}}) = e'((\omega_{n}^{k+ \frac{n}{2}}) ^ 2) + \omega_{n}^{k+ \frac{n}{2}} o'((\omega_{n}^{k+ \frac{n}{2}}) ^ 2) = e'(\omega_{n}^{2k}) - \omega_{n}^{k} o'(\omega_{n}^{2k}) \end{cases} \]然后你发现，\(f(\omega_{n}^{k})\) 和 \(f(\omega_{n}^{k + \frac{n}{2}})\) 仅仅差了一个符号！！！
所以只需要求出 \(e'(x)\) 和 \(o'(x)\) 对 \(\omega^{k}_{n}\)（\(0 \le k \le \frac{n}{2}\)）上的取值，就可以推出 \(f(x)\) 在两倍点数上的取值。
每次问题规模缩小一半，因此时间复杂度 \(O(n \log n)\)。
离散傅里叶逆变换（IDFT）
假设对于两个多项式都得到了 \(t = n + m - 1\) 个点值，设为 \(\{(x_i, A(x_i)) \ | \ 0 \le i < t, i \in \mathbb{Z} \}, \{(x_i, B(x_i)) \ | \ 0 \le i < t, i \in \mathbb{Z} \}\)。
那么可以知道，多项式 \(C(x) = A(x) \times B(x)\) 的点值表示一定为：
\[\{(x_i, A(x_i) B(x_i)) \ | \ 0 \le i < t, i \in \mathbb{Z} \} \]现在，只需要将这 \(t\) 个点插值回去，就可以得到 \(A(x)B(x)\) 了。
先设这 \(t\) 个点值分别是：\(\{(x_i, v_i) \ | \ 0 \le i < t, i \in \mathbb{Z} \}\)，设最后的多项式为 \(C(x) = \sum \limits_{i = 0}^{n + m - 2} c_i x^i\)，这里直接给出结论：
\[c_k = \dfrac{1}{n} \sum_{i = 0}^{t - 1} v_i \omega_{t}^{-ki} \]由此可见，IDFT 和 DFT 仅仅有一个负号的区别。只要将所有的单位根从 \(\omega_{n}^{k}\) 变成 $ - \omega_{n}^{k}$ 即可。
void FFT(cp a[], int n, int op) {
if (n == 1) return;
cp a1[n], a2[n];
rop(i, 0, n >> 1) a1[i] = a[i << 1], a2[i] = a[(i << 1) + 1];
FFT(a1, n >> 1, op), FFT(a2, n >> 1, op);
cp Wn = {cos(2 * pi / n), op * sin(2 * pi / n)};
cp Wk = {1, 0};
rop(i, 0, n >> 1) {
a[i] = a1[i] + Wk * a2[i];
a[i + (n >> 1)] = a1[i] - Wk * a2[i];
Wk = Wk * Wn;
}
}
void FFT(cp a[], cp b[], int n, int m) {
m = n + m; n = 1;
while (n <= m) n <<= 1;
FFT(a, n, 1); FFT(b, n, 1);
rop(i, 0, n) a[i] = a[i] * b[i];
FFT(a, n, -1);
rep(i, 0, m) a[i].x = a[i].x / n;
}
FFT 优化
- 三次变两次优化
原本的朴素 FFT，将 \(\{a\}, \{b\}\) 两个序列分别求值，乘起来再 IDFT 插值一下，一共跑了三次 FFT。这是不好的。
三次变两次优化是这样的：将原序列合并成一个复数：\(\{a_i + b_i \times i\}\)。做一遍 DFT 把求出的每个函数值平方。因为 \((a + bi) ^ 2 = (a ^ 2 - b ^ 2) + (2abi)\)。因此把虚部取出来以后除以 \(2\) 就是答案。
void FFT(cp a[], int n, int op) {
if (n == 1) return;
cp a1[n], a2[n];
rop(i, 0, n >> 1) a1[i] = a[i << 1], a2[i] = a[(i << 1) + 1];
FFT(a1, n >> 1, op), FFT(a2, n >> 1, op);
cp Wn = {cos(2 * pi / n), op * sin(2 * pi / n)};
cp Wk = {1, 0};
rop(i, 0, n >> 1) {
a[i] = a1[i] + a2[i] * Wk;
a[i + (n >> 1)] = a1[i] - a2[i] * Wk;
Wk = Wk * Wn;
}
}
int main() {
read(n, m);
rep(i, 0, n) scanf("%lf", &a[i].x);
rep(i, 0, m) scanf("%lf", &a[i].y);
m = n + m; n = 1;
while (n <= m) n <<= 1;
FFT(a, n, 1);
rop(i, 0, n) a[i] = a[i] * a[i];
FFT(a, n, -1);
rep(i, 0, m) printf("%d ", (int)(a[i].y / 2 / n + 0.5));
}
- 蝴蝶变换优化
后面再补吧。其实本人感觉这个优化不是那么必要，因为三次变两次实在太快了。
FFT 例题
可以设 \(x = 10\)，把 \(a\) 写成 \(A(x) = \sum \limits_{i = 0}^{n} a_i x^i\) 的形式（\(n = \log_{10} a\)）。同理可以把 \(b\) 转化为多项式 \(B(x)\)。
这样求两个数相乘就是求 \(A(x) \times B(x)\) 啊。
所以直接 \(O(n \log n)\) 做完了。
给出 \(n\) 个数 \(q_1,q_2, \dots q_n\)，定义
\[F_j~=~\sum_{i = 1}^{j - 1} \frac{q_i \times q_j}{(i - j)^2}~-~\sum_{i = j + 1}^{n} \frac{q_i \times q_j}{(i - j)^2} \]\[E_i~=~\frac{F_i}{q_i} \]对 \(1 \leq i \leq n\)，求 \(E_i\) 的值。
首先发现这个除以 \(q_i\) 就是没用。所以可以化简成：
\[E_j = \sum_{i = 1}^{j - 1} \frac{q_i}{(i - j)^2}~-~\sum_{i = j + 1}^{n} \frac{q_i}{(i - j)^2} \]先看前面这个式子。答案就是：
\[(j - 1) ^ 2 q_1 + (j - 2) ^ 2 q_2 + (j - 3) ^ 2 q_3 \cdots \]设 \(f(x) = \sum q_i x ^ i, g(x) = \dfrac{1}{i ^ 2} x ^ i\)。可以发现，\(E_j' = (f \times g)[j]\)
再看后面这一块的式子。我们把 \(f(x)\) 的系数翻转，变成 \(f'(x) = \sum q_{n - j + 1} x ^ j\)。那么可以发现 \(E_{j}'' = (f' \times g)[n - j + 1]\)。
跑两次 FFT 就完事了。
首先发现加减相对于两个手环是对称的。因此可以把对一个手环的减法转化成对另一个手环的加法。这样可以假设全是在第一个手环上执行的加减操作。
第一个手环执行了加 \(c\) 的操作，且旋转过之后的序列为 \([x_1, x_2 \cdots x_n]\)，第二个手环为 \([y_1, y_2 \cdots y_n]\)。计算差异值并化简，可以得到差异值是：
\(\sum x ^ 2 + \sum y ^ 2 + c ^ 2 n + 2c(\sum x - \sum y) - 2 \sum xy\)
可以发现，这个序列只有最后一项是不定的。
因此将 \(y\) 序列翻转后再复制一倍，与 \(x\) 卷积，答案就是卷积后序列的 \(n + 1 \sim 2n\) 项系数的 \(\max\)。
直接暴力枚举 \(c\)，加上前面依托就行了。
快速数论变换（NTT）
快速数论变换就是基于原根的快速傅里叶变换。
首先考虑快速傅里叶变换用到了单位根的什么性质。
-
\(\omega_{n}^{k}, 0 \le k < n\) 互不相同。
-
\(\omega_{n}^{k} = \omega_{n}^{k + \frac{n}{2}}\)。
-
\(\omega_{n}^{k} = \omega_{2n}^{2k}\)。
数论中，原根恰好满足这些性质。
对于一个素数的原根 \(g\)，设 \(g_n = g ^ {\frac{p - 1}{n}}\)。那么：
\[g_n ^ n = g ^ {p - 1} \equiv 1(\bmod \ p) \]\[g_n ^ {\frac{n}{2}} = g ^ {\frac{p - 1}{2}} \equiv - 1(\bmod ~ p) \]\[g ^ \alpha + g ^ \beta = g ^ {\alpha + \beta} \]\[g_{\alpha n}^{\alpha k} = g_{n}^{k} \]我们发现它满足 \(\omega_{n}^{k}\) 的全部性质！
因此，只需要用 \(g_{n}^k = g_{}^{\frac{p - 1}{n} k}\) 带替全部的 \(\omega_{n}^{k}\) 就可以了。
tips：对于一个质数，只有当它可以表示成 \(p = 2 ^ \alpha + 1\)，且需要满足多项式的项数 \(< \alpha\) 时才能使用 NTT。\(p\) 后面有个加一，是因为 \(g_n\) 指数的分子上出现了 \(-1\)；\(p - 1\) 需要时二的整数次幂，是因为每次都要除以 \(2\)。
bonus：常用质数及原根
\(p = 998244353 = 119 \times 2 ^ {23} + 1, g = 3\)
\(p = 1004535809 = 479 \times 2 ^ {21} + 1, g = 3\)
void NTT(int *a, int n, int op) {
if (n == 1) return;
int a1[n], a2[n];
rop(i, 0, n >> 1) a1[i] = a[i << 1], a2[i] = a[(i << 1) + 1];
NTT(a1, n >> 1, op), NTT(a2, n >> 1, op);
int gn = qpow((op == 1 ? g : invg), (mod - 1) / n), gk = 1;
rop(i, 0, n >> 1) {
a[i] = (a1[i] + 1ll * gk * a2[i]) % mod;
a[i + (n >> 1)] = (a1[i] - 1ll * gk * a2[i] % mod + mod) % mod;
gk = 1ll * gk * gn % mod;
}
}
NTT 优化：蝴蝶变换
盗大佬一张图
这是进行 NTT 的过程中数组下标的变化。
这样看似乎毫无规律。但是把他们写成二进制：
变换前：\(000 ~ 001 ~ 010 ~ 011 ~ 100 ~ 101 ~ 110 ~ 111\)
变换后：\(000 ~ 100 ~ 010 ~ 110 ~ 001 ~ 101 ~ 011 ~ 111\)
可以发现，就是对每个下标进行了二进制翻转。
因此可以先预处理出每个下标在递归底层对应的新下标。然后从底层往顶层迭代合并。
预处理每个下标的二进制翻转：
rop(i, 0, n) rev[i] = rev[i >> 1] >> 1 | (i & 1) << bit;
优化后的 NTT：
void NTT(int *a, int n, int op) {
rop(i, 0, n) if (i < rev[i]) swap(a[i], a[rev[i]]);
for (int mid = 1; (mid << 1) <= n; mid <<= 1) {
int gn = qpow((op == 1 ? g : invg), (mod - 1) / (mid << 1));
for (int i = 0, gk = 1; i < n; i += (mid << 1), gk = 1)
for (int j = 0; j < mid; j ++ , gk = 1ll * gk * gn % mod) {
int x = a[i + j], y = 1ll * a[i + j + mid] * gk % mod;
a[i + j] = Mod(x + y), a[i + j + mid] = Mod(x - y);
}
}
}
当然了，FFT 也可以用蝴蝶变换来优化。实践证明，速度变快了至少二分之一。
FFT 的迭代实现
void FFT(cp *a, int n, int op) {
rop(i, 0, n) if (i < rev[i]) swap(a[i], a[rev[i]]);
for (int mid = 1; (mid << 1) <= n; mid <<= 1) {
cp Wn = {cos(pi / mid), op * sin(pi / mid)};
for (int i = 0; i < n; i += (mid << 1)) {
cp Wk = {1, 0};
for (int j = 0; j < mid; j ++ , Wk = Wk * Wn) {
cp x = a[i + j], y = Wk * a[i + j + mid];
a[i + j] = x + y, a[i + j + mid] = x - y;
}
}
}
}
任意模数多项式乘法（MTT）
计算 \(f \times g(\bmod ~ p)\) 的结果（\(p\) 不一定能够表示成 \(2 ^ \alpha + 1\) 的形式）。
这个东西有两种做法，但是我只学会了三模 NTT。
首先，卷积之后每个系数最多达到 \(\max \{V\} ^ 2 \times n\) 的大小。拿模板题举例，这个东西就是 \(10 ^ {23}\)。因此只需要找三个模数 \(p_1, p_2, p_3\)，满足 \(p_1p_2p_3 > \max \{V\} ^ 2 \times n\)，然后用他们分别 NTT，最后再 CRT / exCRT 合并即可。
void NTT(int *a, int n, int op, int p) {
rop(i, 0, n) if (i < rev[i]) swap(a[i], a[rev[i]]);
for (int mid = 1; (mid << 1) <= n; mid <<= 1) {
int gn = qpow((op == 1 ? g : invg), (p - 1) / (mid << 1), p);
for (int i = 0, gk = 1; i < n; i += (mid << 1), gk = 1)
for (int j = 0; j < mid; j ++ , gk = gk * gn % p) {
int x = a[i + j], y = a[i + j + mid] * gk % p;
a[i + j] = Mod(x + y, p), a[i + j + mid] = Mod(x - y, p);
}
}
}
int CRT(int a, int b, int c) {
int k = Mod(b - a, p2) * inv1 % p2;
int x = k * p1 + a;
k = Mod(c - x, p3) * inv2 % p3;
x = (x + k * p1 % p * p2 % p) % p;
return x;
}
signed main() {
read(n, m, p);
rep(i, 0, n) read(a[i]);
rep(i, 0, m) read(b[i]);
memcpy(a1, a, sizeof a); memcpy(a2, a, sizeof a);
memcpy(a3, a, sizeof a); memcpy(b1, b, sizeof b);
memcpy(b2, b, sizeof b); memcpy(b3, b, sizeof b);
m = n + m; n = 1; int bit = 0;
while (n <= m) n <<= 1, bit ++ ; bit -- ;
rop(i, 0, n) rev[i] = rev[i >> 1] >> 1 | (i & 1) << bit;
invg = inv(g, p1); NTT(a1, n, 1, p1); NTT(b1, n, 1, p1);
invg = inv(g, p2); NTT(a2, n, 1, p2); NTT(b2, n, 1, p2);
invg = inv(g, p3); NTT(a3, n, 1, p3); NTT(b3, n, 1, p3);
rop(i, 0, n) a1[i] = 1ll * a1[i] * b1[i] % p1;
rop(i, 0, n) a2[i] = 1ll * a2[i] * b2[i] % p2;
rop(i, 0, n) a3[i] = 1ll * a3[i] * b3[i] % p3;
invg = qpow(g, p1 - 2, p1); NTT(a1, n, -1, p1);
invg = qpow(g, p2 - 2, p2); NTT(a2, n, -1, p2);
invg = qpow(g, p3 - 2, p3); NTT(a3, n, -1, p3);
inv1 = inv(n, p1); inv2 = inv(n, p2); inv3 = inv(n, p3);
rop(i, 0, n) a1[i] = a1[i] * inv1 % p1, a2[i] = a2[i] * inv2 % p2, a3[i] = a3[i] * inv3 % p3;
inv1 = inv(p1, p2); inv2 = inv(p1 * p2 % p3, p3);
rep(i, 0, m) printf("%lld ", CRT(a1[i], a2[i], a3[i]));
return 0;
}
这里选择的模数为：\(p_1 = 998244353, p_2 = 469762049, p_3 = 1004535809\)。他们的原根都为 \(g = 3\)。
多项式求逆
求多项式 \(f(x)\) 的逆元 \(f^{-1}(x)\)。\(f(x)\) 的逆元定义为满足 \(f(x) g(x) = 1(\bmod \ x ^ n)\) 的多项式 \(g(x)\)。
使用倍增法即可求出多项式的逆元。
首先假设已经求出了 \(f(x) g'(x) \equiv 1(\bmod \ x ^ {\lceil \frac{n}{2} \rceil})\)。再假设 \(f(x) \bmod \ x ^ n\) 意义下逆元为 \(g(x)\)。那么有：
\[g(x) \equiv g'(x) (\bmod \ x ^ {\lceil \frac{n}{2} \rceil}) \]\[(g(x) - g'(x)) \equiv 0 (\bmod \ x ^ {\lceil \frac{n}{2} \rceil}) \]\[(g(x) - g'(x)) ^ 2 \equiv 0 (\bmod \ x ^ n) \]\[g^2(x) - 2 g(x) g'(x) + g'^2(x) \equiv 0 (\bmod \ x ^ {n}) \]两边同时乘以 \(f(x)\)，可以得到：
\[g(x) - 2 g'(x) + f(x) g'^2(x) \equiv 0 (\bmod \ x ^ n) \]\[g(x) \equiv 2g'(x) - f(x) g'^2(x) (\bmod \ x ^ n) \]\[g(x) \equiv g'(x) (2 - f(x)g'(x)) (\bmod \ x ^ n) \]倍增求即可。
void Inv(int *f, int *g, int n) {
if (n == 1) {
g[0] = inv(f[0]); return;
} Inv(f, g, (n + 1) >> 1);
int m = 1, bit = 0;
while (m < (n << 1)) m <<= 1, bit ++ ; bit -- ;
rop(i, 0, n) tmp[i] = f[i];
rop(i, n, m) tmp[i] = 0;
rev = new int [m + 5]();
rop(i, 1, m) rev[i] = (rev[i >> 1] >> 1) | (i & 1) << bit;
NTT(tmp, m, 1); NTT(g, m, 1);
rop(i, 0, m) g[i] = (2ll - 1ll * g[i] * tmp[i] % p + p) % p * g[i] % p;
NTT(g, m, -1); int invn = inv(m);
rop(i, 0, m) g[i] = 1ll * g[i] * invn % p;
rop(i, n, m) g[i] = 0;
free(rev); rev = NULL;
}