多项式学习笔记（Updating）
1. 复数，单位根，原根
1.1 复数的表示
- 代数表示：\(z = a+b_i\)，记 \(\overline z\) 的共轭为 \(\overline z = a-bi\)。
- 三角表示：记复数的模为 \(|z| = r = \sqrt{a^2+b^2}\)，复数的辐角为 \(\mathrm{Arg}\space z = \theta = \frac b a\)，\(\theta \in (-\pi, \pi]\) 的辐角称作辐角主值，记作 \(\arg z = \arctan \frac b a\)，那么 \(z = r(\cos\theta+i\sin\theta)\)。
- 用欧拉公式表示：\(z=re^{i\theta}\)。
1.2 复数的运算
1.3 单位根
\(r = 1\) 的复数在复平面的单位圆上，其 \(n\) 次幂为 \(\cos n\theta + i\sin n\theta\)，即在单位圆上从 \((1, 0)\) 开始以 \(\theta\) 的角度顺时针旋转了 \(n\) 次。
考虑 \(n\) 等分单位圆，记逆时针第 \(k\) 个（从 \(0\) 开始标号）等分点为 \(\omega_n^k = \cos(\frac{2k\pi} n)+i\sin(\frac{2k\pi}n)\)。所有的 \(\omega_n^k\) 记为 \(n\) 次单位根。\(\omega_n^1\) 简记为 \(\omega_n\)。显然 \((\omega_n^k)^n=1\)。若 \(\gcd(n, k) = 1\)，称 \(\omega_n^k\) 为本源单位根，所有本源单位根的 \(0 \sim n - 1\) 次幂是互不相同的。
另外，根据欧拉公式，我们有 \(\omega_n^k=e^{\frac{2k\pi}ni}\)。
性质：
- \(\omega_n^k = \omega_n^{k+tn}, t\in \mathbb Z\)。
- \(\omega_n^k=\omega_{tn}^{tk}, t \in \mathbb Z\)。由定义直接得出。
- \(n\) 为偶数时，\(-\omega_n^k=\omega_n^{k\pm\frac n2}\)。相当于在单位圆上顺时针或逆时针转半圈。
- \(\sum_{i=0}^{n-1}\omega_n^i=0\)，等比数列求和可得。
1.4 原根
原根可以看做是模 \(p\) 意义下的单位根，设 \(n = \varphi(p)\)，若 \(\delta_p(g)=n\)，则 \(g\) 为模 \(p\) 的原根，如果模 \(p\) 的原根存在，则 \(g^i \bmod p, i \in [0, n - 1]\) 两两不同，这和单位根的性质极其相似，模 \(p\) 意义下我们可以直接用 \(g\) 代替 \(n\) 次单位根，对于 \(d\mid n\)，也可以用 \(g^{\frac n d}\) 代替 \(d\) 次单位根（可以用 \(\omega_d=\omega_n^{\frac n d}\) 理解）。
2. 多项式
形如 \(\sum_{i = 0}^n a_ix^i\) 的和式称为多项式，记为 \(f(x) = \sum_{i=0}^n a_ix^i\)，其最高系数非零项记作 \(\deg f\)，若 \(a_n \neq 0\) 则 \(\deg f = n\)。所有 \(f(x) = 0\) 的 \(x\) 称为多项式的根，一个\(n\) 次多项式有 \(n\) 个复数根，这些根可能有重合。也可以用 \(f_i\) 表示 \(f\) 的 \(i\) 次项系数。
2.1 表示方法
系数表示
用 \(a_i\) 组成的序列来表示。
点值表示
我们有 \(n\) 个点值 \((x_0, y_0), (x_1, y_1), \dots, (x_{n - 1}, y_{n - 1})\)，若 \(x_i\) 互不相同，则这 \(n\) 个点值可以唯一确定一个最高次数为 \(n - 1\) 的多项式。
从系数表示法转为点值表示法称为求值。
从点值表示法转为系数表示法称为插值。
2.2 运算
设 \(f(x) = \sum_{i=0}^n a_ix^i, g(x) = \sum_{i=0}^m b_ix_i\)，则
\(\deg(f+g) = \max(\deg f, \deg g)\)。
\(\deg(f* g)=\deg f + \deg g\)。
记 \(f* g\) 或 \(fg\) 为加法卷积，而 \(f \cdot g\) 表示点乘，两者系数相乘。
计算两个多项式相加，在系数表示法下时间复杂度为 \(\mathcal{O}(n + m)\)，点值表示法下要 \(\max(n, m) + 1\) 个点值才能确定，时间复杂度 \(\mathcal{O}(n + m)\)。计算两个多项式相乘，在系数表示法下时间复杂度为 \(\mathcal{O}(nm)\)，点值表示法下要 \(n + m + 1\) 个点值才能确定，时间复杂度 \(\mathcal{O}(n + m)\)。
2.3 次数上界
我们用 \(\pmod {x^n}\) 来截断一个多项式，只保留 \(0 \sim n - 1\) 次项，\(F(x) \equiv G(x) \pmod {x^n}\) 则 \(G(x) = \sum_{i=0}^{n-1}F_ix^i\)。
3. 拉格朗日插值
3.1 算法介绍
我们有 \(n + 1\) 个点 \((x_0, y_0), (x_1, y_1), \dots, (x_n, y_n)\)，要求出一个 \(n\) 次多项式 \(f(x)\) 满足 \(f(x_i) = y_i\)。
3.1.1 常规形式
我们构造 \(n\) 个函数 \(f_1(x), \dots, f_n(x)\)，\(f_i(x)\) 的图象过 \(\begin{cases} (x_j, 0), & (j\neq i) \\ (x_j, y_j), & (j = i) \end{cases}\)，那么 \(f(x) = \sum_{i=1}^n f_i(x)\)。
我们设 \(f_i(x) = a\prod_{j\neq i}(x - x_j)\)，带入 \((x_i, y_i)\) 得 \(a = \frac{y_i}{\prod_{j\neq i} (x_i-x_j)}\)。即
时间复杂度 \(\mathcal{O}(n^2)\)。但是可以优化到 \(\mathcal{O}(n \log^2 n)\)，具体见 5.?.?，我还没写。
为了求出 \(f(x)\) 的系数，考虑预处理 \(M = \prod\limits_{i=1}^n (x - x_i)\) 的各项系数，\(px_i = \prod\limits_{j\neq i}(x_i - x_j)\)，然后就好求了。
洛谷模版题代码，0-indexed：
namespace Loop1st {
int n, k;
ll qpow(ll a, ll b) {
ll res = 1;
for (; b; b >>= 1, a = a * a % mod) if (b & 1) res = res * a % mod;
return res;
}
vector<ll> lagrange(const vector<ll> &x, const vector<ll> &y) {
vector<ll> M(n + 1), px(n, 1), f(n);
M[0] = 1;
for (int i = 0; i < n; i++) {
for (int j = i; j >= 0; j--) {
M[j + 1] = (M[j + 1] + M[j]) % mod;
M[j] = M[j] * (mod - x[i]) % mod;
}
}
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) if (i != j) {
px[i] = px[i] * (x[i] + mod - x[j]) % mod;
}
}
for (int i = 0; i < n; i++) {
ll t = y[i] * qpow(px[i], mod - 2) % mod, now = M[n];
for (int j = n - 1; j >= 0; j--) {
f[j] = (f[j] + now * t) % mod;
now = (now * x[i] + M[j]) % mod;
}
}
return f;
}
void Main(int tc) {
cin >> n >> k;
vector<ll>x(n), y(n);
for (int i = 0; i < n; i++) cin >> x[i] >> y[i];
const auto f = lagrange(x, y);
ll pw = 1, ans = 0;
for (int i = 0; i < n; i++) ans = (ans + pw * f[i]) % mod, pw = pw * k % mod;
cout << ans << '\n';
}
}
3.1.2 取值连续的形式
即给出 \(n\) 个点，但是 \(x_i\) 是从小到大连续的 \(n\) 个数。
则公式简化为
对于分子容易通过前后缀预处理得到，分母考虑对 \(i > j\) 和 \(i < j\) 分讨，\(i < j\) 的时候可能是负的，讨论一下得到
于是可以 \(\mathcal{O}(n)\) 计算单点的值。但是似乎不能 \(\mathcal{O}(n)\) 计算出所有系数。
代码见例题。
3.1.3 动态问题
我们插入一个点值时，容易 \(\mathcal{O}(n)\) 重新插值，可以 \(\mathcal{O}(n)\) 求任意点值，但是求系数似乎还是要 \(\mathcal{O}(n^2)\)。
3.2 例题
CF622F
考虑到 \(f_n = \sum\limits_{i=1}^n i^k\)，则 \(f_n\) 的差分为 \(n^k\)，是一个 \(k\) 次多项式，所以 \(f_n\) 是一个 \(k + 1\) 次多项式，所以找 \(k + 2\) 个连续的点插值，最后求 \(f_n\) 即可。时间复杂度 \(\mathcal{O}(k)\)。
namespace Loop1st {
int n, k, b[N], p[N], tot;
ll f[N], pre[N], suf[N], fac[N], ifac[N];
ll qpow(ll a, ll b) {
ll res = 1;
for (; b; b >>= 1, a = a * a % mod) if (b & 1) res = res * a % mod;
return res;
}
void Main(int tc) {
cin >> n >> k;
f[1] = 1;
for (int i = 2; i <= k + 2; i++) {
if (!b[i]) { p[++tot] = i; f[i] = qpow(i, k); }
for (int j = 1; j <= tot && i * p[j] <= k + 2; j++) {
b[i * p[j]] = 1;
f[i * p[j]] = f[i] * f[p[j]] % mod;
if (i % p[j] == 0) break;
}
f[i] = qpow(i, k);
}
for (int i = 2; i <= k + 2; i++) f[i] = (f[i] + f[i - 1]) % mod;
if (n <= k + 2) { cout << f[n] << '\n'; return ; }
pre[0] = suf[k + 3] = 1;
for (int i = 1; i <= k + 2; i++) pre[i] = pre[i - 1] * (n - i) % mod;
for (int i = k + 2; i; i--) suf[i] = suf[i + 1] * (n - i) % mod;
fac[0] = ifac[0] = 1;
for (int i = 1; i <= k + 2; i++) fac[i] = fac[i - 1] * i % mod;
ifac[k + 2] = qpow(fac[k + 2], mod - 2);
for (int i = k + 1; i; i--) ifac[i] = ifac[i + 1] * (i + 1) % mod;
ll ans = 0;
for (int i = 1; i <= k + 2; i++) {
ll sgn = (k + 2 - i) & 1 ? mod - 1 : 1;
ans = (ans + sgn * ifac[i - 1] % mod * ifac[k + 2 - i] % mod * pre[i - 1] % mod * suf[i + 1] % mod * f[i] % mod) % mod;
}
cout << ans << '\n';
}
}
calc
我们不考虑顺序，从小到大填数，设 \(f_{i, j}\) 为考虑到 \([1, i]\)，填了 \(j\) 个数的答案，则 \(f_{i, j} = f_{i - 1, j} + if_{i - 1, j - 1}\)，答案即 \(n!f_{n, k}\)，把 \(f_{i, j}\) 当成关于 \(i\) 的多项式，差分一下，\(f_j(i) - f_j(i - 1) = if_{j - 1}(i - 1)\)，所以 \(f_j(i)\) 的差分的次数是 \(f_{j-1}(i-1)\) 的次数 \(+1\)，\(f_j(i)\) 的次数是 \(f_{j-1}(i-1)\) 的次数 \(+2\)，而 \(f_0(i)\) 的次数是 \(0\)，所以 \(f_n(k)\) 是 \(2n\) 次多项式。带入 \(2n + 1\) 个点值即可，时间复杂度 \(\mathcal{O}(n^2)\)。
namespace Loop1st {
int n, m, mod, b[N], p[N], tot;
ll k, f[N][N], pre[N], suf[N], fac[N], ifac[N];
ll qpow(ll a, ll b) {
ll res = 1;
for (; b; b >>= 1, a = a * a % mod) if (b & 1) res = res * a % mod;
return res;
}
void Main(int tc) {
cin >> k >> n >> mod;
m = n << 1 | 1;
for (int i = 0; i <= m; i++) f[i][0] = 1;
for (int i = 1; i <= m; i++) {
for (int j = 1; j <= n; j++) {
f[i][j] = (f[i - 1][j] + f[i - 1][j - 1] * i) % mod;
}
}
fac[0] = ifac[0] = 1;
for (int i = 1; i <= m; i++) fac[i] = fac[i - 1] * i % mod;
ifac[m] = qpow(fac[m], mod - 2);
for (int i = m - 1; i; i--) ifac[i] = ifac[i + 1] * (i + 1) % mod;
if (k <= m) { cout << fac[n] * f[k][n] % mod << '\n'; return ; }
pre[0] = suf[m + 1] = 1;
for (int i = 1; i <= m; i++) pre[i] = pre[i - 1] * (k - i) % mod;
for (int i = m; i; i--) suf[i] = suf[i + 1] * (k - i) % mod;
ll ans = 0;
for (int i = 1; i <= m; i++) {
ll sgn = (m - i) & 1 ? mod - 1 : 1;
ans = (ans + sgn * ifac[i - 1] % mod * ifac[m - i] % mod * pre[i - 1] % mod * suf[i + 1] % mod * f[i][n] % mod) % mod;
}
cout << fac[n] * ans % mod << '\n';
}
}
4. FFT & NTT
这个部分的代码参考了 Alex_Wei。
4.1 点值、系数表示的转化
设 \(f(x) = \sum_{i=0}^{n-1} a_ix^i\)，我们考虑现在给定互不相同的 \(m - 1\) 个位置 \(x_0, x_1, x_2, \dots, x_{m - 1}\)，需要求 \(f(x_0), f(x_1), f(x_2), \dots, f(x_{m - 1})\) 的值，用矩阵乘法的形式描述问题：
左侧矩阵叫做范德蒙德矩阵（\(n = m\) 时称作范德蒙德方阵），\(m < n\) 时无法求逆，\(m \ge n\) 时可以选出 \(x_i\) 互不相同的 \(n\) 行求逆，通过点值乘上范德蒙德矩阵的逆还原系数。
朴素求值的复杂度为 \(\mathcal{O}(nm)\)。
4.2 DFT
DFT 就是把一个复数序列 \(x_0, x_1, \dots, x_{N - 1}\) 用如下公式转化为另一个复数序列 \(X_0, X_1, \dots, X_{N-1}\)：
设 \(f(x) = \sum_{i=0}^{N-1}a_ix^i\)，则
即可以看成 \(f(x)\) 在 \(N\) 个单位根处求值。
4.3 FFT
朴素地求出 \(n^2\) 个点值仍然是 \(\mathcal{O}(n^2)\) 的。我们现在的目标是给定一个 \(n - 1\) 次多项式 \(f(x) = \sum_{i=0}^{n-1}a_ix^i\)，求出其至少 \(n\) 个点值。但 FFT 所实现的变换又和 DFT 有着细微的差别，具体体现在两者的范德蒙德矩阵带入单位根的顺序是相反的，且 DFT 不要求长度为 \(2\) 的整幂。
4.3.1 朴素的 FFT
一个多项式函数可以表示为一个奇函数和一个偶函数的和。奇函数的偶次项系数都为 \(0\)，偶函数的奇次项系数都为 \(0\)。
不妨令 \(n\) 变为 \(2^{\left\lceil\log n\right\rceil}\)，\(a_i = 0(i \ge n)\)。这样让问题变得更加方便。
考虑将 \(f(x)\) 拆成偶函数 \(f_e(x)\) 和奇函数 \(f_o(x)\) 之和，则
我们现在选取 \(n\) 个位置 \(x_0, x_1, \dots, x_{n-1}\)，并且 \(x_{i+\frac n2}=-x_i(i\in [0, \frac n2-1])\)，求解 \(f_e(x_i)\) 和 \(f_o(x_i)(i \in [0, \frac n2-1])\) 就可以求出 \(f(x)\) 的至少 \(n\) 个点值。由于 \(f_e(x)\) 只有偶次项，\(f_o(x)\) 只有奇次项，我们可以把 \(f_e(x)\) 和 \(f_o(x)\) 变为 \(\frac n2\) 次项，即 $$g(x) = a_0+a_2x+a_4x2+\dots+a_{n-2}x=\sum_{i=0}^{\frac n2-1}a_{2i}x^i$$
那么 \(f(x) = g(x) + xh(x), f(-x)=g(x)-xh(x)\)。
这样我们就把规模减半了。
如果我们选取的 \(x_i\) 满足前 \(\frac n2\) 个数和后 \(\frac n2\) 个数互为相反数，并且前 \(\frac n2\) 个数平方后相反数仍然满足前 \(\frac n4\) 个数和后 \(\frac n4\) 个数互为相反数，依次类推，即前 \(\frac n{2^i}\) 个数前 \(\frac n{2^{i+1}}\) 个数的 \(2^i\) 次方和后 \(\frac n{2^{i+1}}\) 个数的 \(2^i\) 次方互为相反数，那么我们就可以层层递归下去，递归 \(\log n\) 层后求出所有点值。
我们发现 \(x_i = \omega_n^i\) 恰好满足这个需求，用 \(\omega_n^k=-\omega_n^{k+\frac n2}\) 这个性质即可证明。
然后我们再尝试化简一下，
然后递归下去求解即可，时间复杂度 \(\mathcal{O}(n \log n)\)。
递归的常数比较大，下面介绍非递归的写法。
4.3.2 蝴蝶变换
考虑到递归的最后一层，若当前所处理的位置的系数为 \(a_p\)，则这次递归处理的 \(f(x) = a_p\)，我们如果知道最后一层 \(a_p\) 的顺序，就可以递推地进行 FFT 了，比如，初始系数为 \((a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7)\)，递归第一次，左边 \(g(x)\) 的系数为 \((a_0, a_2, a_4, a_6)\)，右边 \(h(x)\) 的系数为 \((a_1, a_3, a_5, a_7)\)，递归第二次，变为 \((a_0, a_4), (a_2, a_6), (a_1, a_3), (a_5, a_7)\)。在第 \(k\) 层向下递归，\(a_i\) 被分到左边还是右边（\(g\) 还是 \(h\)），取决于 \(\left\lfloor \frac i{2^{k-1}} \right\rfloor\) 的奇偶性，即 \(i\) 从低到高第 \(k\) 位是 \(0\) 还是 \(1\)。即考虑令 \(r_i\) 表示 \(i\) 二进制表示的前 \(\log n\) 位（不足用 \(0\) 补齐）翻转后的数，如 \(r_{11}=(\overline{1011})_2=(1101)_2=13\)。那么最后 \(a_p\) 的顺序即 \(r_p\) 从小到大的顺序。\(r_i\) 容易 \(\mathcal{O}(n)\) 递推得到。
有了最后一层的顺序我们就可以自底向上递推了，先同时让 \(a_{r_i} \gets a_i\)，由于 \(r_{r_i}=i\)，也要让 \(a_i \gets a_{r_i}\)，这可以通过交换所有 \(a_i, a_{r_i}(i < r_i)\) 来实现。然后从小到大枚举问题规模 \(2d\)，每个子问题为 \([i, i + 2d)\)，其中 \(i = 2kd\)，\([i, i + d)\) 和 \([i + d, i + 2d)\) 都已在之前处理过，那么直接按照上述递归 FFT 的公式,枚举一个 \(j \in [0, d)\)，令 \(x = i + j, y = i + j + d\)，\(f'_k\) 表示当前规模下的 \(f(\omega_{2d}^k)\)，\(f_k\) 表示上一规模下的 \(f(\omega_d^k)\)，则
这个过程就是 FFT。
4.4 IFFT
如何将点值转化为系数？范德蒙德矩阵的本质是系数转化为点值，而范德蒙德矩阵的逆就是点值转化为系数，即插值。
设范德蒙德矩阵为 \(F\)，其逆为 \(F^{-1}\)。
两边左乘 \(F^{-1}\)，则 \(f(x_j)\) 对 \(a_i\) 的贡献为 \((F^{-1})_{i, j}\)。将 \(x_{i} = \omega_n^i\) 带入拉格朗日插值的式子：
所以 \((F^{-1})_{i, j} = [x^i]\prod\limits_{k \neq j}\frac{x-\omega_n^k}{\omega_n^j-\omega_n^k}\)。考虑令
把 \(F\) 和 \(F^{-1}\) 对照一下：
我们发现，\(F^{-1}\) 就是 \(F\) 的所有单位根次数取反，再将系数除以 \(n\)。
所以，IFFT 的过程就是 FFT 的过程上进行一些简单的修改。
4.5 多项式乘法
4.5.1 朴素实现
现在，给定一个 \(n - 1\) 次多项式和一个 \(m - 1\) 次多项式 \(a, b\)，要求出 \(c = a \times b\)。\(c\) 是 \(n + m - 2\) 次多项式，我们先将 \(n + m - 2\) 补全到 \(2\) 的幂次 \(L\)，然后对 \(a, b\) 分别做长度为 \(L\) 的 FFT，点值相乘得到 \(\hat{c}\)，再对 \(\hat{c}\) 做 IFFT 得到 \(c\)。时间复杂度 \(\mathcal{O}(n \log n)\)。下面给出 P3803 的代码。
namespace Loop1st {
int n, m, r[N];
struct Complex {
double a, b;
Complex(double _a = 0, double _b = 0): a(_a), b(_b) {}
Complex operator + (const Complex &A) const { return {a + A.a, b + A.b}; }
Complex operator - (const Complex &A) const { return {a - A.a, b - A.b}; }
Complex operator * (const Complex &A) const { return {a * A.a - b * A.b, a * A.b + b * A.a}; }
} f[N], g[N], h[N];
void FFT(int L, Complex *f, bool type) { // type = 1: FFT type = 0: IFFT
for (int i = 1; i < L; i++) {
r[i] = (r[i >> 1] >> 1) + (i & 1 ? L >> 1 : 0);
if (i < r[i]) swap(f[i], f[r[i]]);
}
for (int d = 1; d < L; d <<= 1) {
Complex w1(cos(pi / d), sin(pi / d));
if (!type) w1.b = -w1.b; // 指数取反相当于取倒数，对单位根来说相当于取共轭
static Complex w[N];
w[0] = {1, 0};
for (int i = 1; i < d; i++) w[i] = w[i - 1] * w1;
for (int i = 0; i < L; i += (d << 1)) {
for (int j = 0; j < d; j++) {
Complex x = f[i | j], y = w[j] * f[i | j | d];
f[i | j] = x + y, f[i | j | d] = x - y;
}
}
}
if (!type) for (int i = 0; i < L; i++) f[i].a /= L; // IFFT 最后只要输出实部就只除实部
}
void Main(int tc) {
cin >> n >> m;
for (int i = 0; i <= n; i++) cin >> f[i].a;
for (int i = 0; i <= m; i++) cin >> g[i].a;
int L = 1;
while (L <= n + m) L <<= 1;
FFT(L, f, 1); FFT(L, g, 1);
for (int i = 0; i < L; i++) h[i] = f[i] * g[i];
FFT(L, h, 0);
for (int i = 0; i <= n + m; i++) cout << (int)(h[i].a + 0.5) << " \n"[i == n + m];
}
}
4.5.2 三次变两次优化
对于 \(a, b\) 系数为实数的情况下，根据 \((a + bi)^2 = (a^2 - b^2) +2abi\)，我们将 \(a + bi\) 这个复系数多项式平方后取出虚部，除以 \(2\) 即可，进一步减小了常数。
namespace Loop1st {
int n, m, r[N];
struct Complex {
double a, b;
Complex(double _a = 0, double _b = 0): a(_a), b(_b) {}
Complex operator + (const Complex &A) const { return {a + A.a, b + A.b}; }
Complex operator - (const Complex &A) const { return {a - A.a, b - A.b}; }
Complex operator * (const Complex &A) const { return {a * A.a - b * A.b, a * A.b + b * A.a}; }
} f[N];
void FFT(int L, Complex *f, bool type) {
for (int i = 1; i < L; i++) {
r[i] = (r[i >> 1] >> 1) + (i & 1 ? L >> 1 : 0);
if (i < r[i]) swap(f[i], f[r[i]]);
}
for (int d = 1; d < L; d <<= 1) {
Complex w1(cos(pi / d), sin(pi / d));
if (!type) w1.b = -w1.b;
static Complex w[N];
w[0] = {1, 0};
for (int i = 1; i < d; i++) w[i] = w[i - 1] * w1;
for (int i = 0; i < L; i += (d << 1)) {
for (int j = 0; j < d; j++) {
Complex x = f[i | j], y = w[j] * f[i | j | d];
f[i | j] = x + y, f[i | j | d] = x - y;
}
}
}
if (!type) for (int i = 0; i < L; i++) f[i].b /= L;
}
void Main(int tc) {
cin >> n >> m;
for (int i = 0; i <= n; i++) cin >> f[i].a;
for (int i = 0; i <= m; i++) cin >> f[i].b;
int L = 1;
while (L <= n + m) L <<= 1;
FFT(L, f, 1);
for (int i = 0; i < L; i++) f[i] = f[i] * f[i];
FFT(L, f, 0);
for (int i = 0; i <= n + m; i++) cout << (int)(f[i].b / 2 + 0.5) << " \n"[i == n + m];
}
}
4.5.3 合并两次实系数 FFT/IFFT
鸽。想学的请见 Alex_Wei 老师的 blog。
4.5.4 拆系数 FFT
用于解决模 \(p\) 意义下的多项式乘法，且 \(p\) 不为 NTT 模数。
我们取 \(B = \lfloor\sqrt p\rfloor\)，将系数拆成 \(f = f_0B + f_1, g = g_0B + g_1\) 的形式，其中 \(f_1, g_1\) 的所有系数 \(<B\)，两两相乘得到 \(fg = (f_0g_0)B^2 + (f_0g_1+f_1g_0)B + f_1g_1\)，然后对 \(f_0, f_1, g_0, g_1\) 进行 FFT，\(\hat{f_0}\hat{g_0}, \hat{f_0}\hat{g_1}+\hat{f_1}\hat{g_0},\hat{f_1}\hat{g_1}\) 进行 IFFT，得到一个 \(7\) 次 FFT 的做法，用三次变两次的优化可以变成 \(5\) 次 FFT。系数的值域为 \(nB^2 = np\)，一般来说可以接受，但是卡精度要 long double。如果用 4.5.3 的技巧可以优化到 \(2\) 次 FFT 和 \(2\) 次 IFFT，效率会更优秀，我没有写代码。
4.6 NTT
4.6.1 朴素的 NTT
事实上，算法竞赛中的多项式题常常需要对一个固定模数 \(p\) 取模，这时我们就要引入 NTT 了。NTT 的本质是在模 \(p\) 意义下进行的 FFT。一般情况下我们也不写 FFT 只写 NTT。
设变换的长度为 \(n\)，我们需要存在一个用于代替单位根 \(\omega_n\) 的数 \(a\) 满足 \(\delta_p(a) = n\)，对于 \(p\) 为质数的情况，我们有 \(\delta_p(g) = p - 1\)，我们发现 \(g^k\) 和 \(\omega_{p-1}^k\) 是等价的。那么，\(q = g^{\frac{p-1}n}\) 等价于 \(\omega_{p-1}^{\frac{p-1}n}=\omega_n\)，换句话说，\(q\) 满足以下性质：
- \(q^0 \sim q^{n-1}\) 互不相同
- \(q^k=q^{k\bmod n}\)
- \((g^{\frac{p-1}{2n}})^2=g^{\frac{p-1}n}\)。
这也隐性要求了 NTT 的模数 \(p\) 满足 \(p - 1\) 是 \(n\) 的倍数，而 \(n\) 是 \(2\) 的整幂，所以 \(2^N \mid p - 1\)，其中 \(N \ge \log n\)。
常见的 NTT 模数：
- \(p=469762049 = 7\times 2^{26}+1\)
- \(\color{red}{p=998244353=7\times 17 \times 2^{23} + 1}\)
- \(p=1004535809=479\times 2^{21}+1\)
这三个模数的最小原根都是 \(3\)。
如果不是常见模数，但是 \(p\) 形如 \(q2^N+1\) 的形式，且 \(N \ge \log n\)，则也可以找一个原根然后 NTT。如果 \(p\) 是一般模数，请参见 4.5.4 的拆系数 FFT 或 4.6.2 三模 NTT。
当然其实如果存在 \(\delta_p(a) = n\) 的单位根 \(a\) 和 \(n\) 的逆元 \(n^{-1}\)，也是可以做 NTT 的。
下面给出多项式乘法的代码实现：
namespace Loop1st {
int n, m, r[N], f[N], g[N], h[N];
int qpow(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = (ll)a * a % mod) if (b & 1) res = (ll)res * a % mod;
return res;
}
void NTT(int L, int *a, bool type) {
static ull w[N], f[N]; // 用 ull 减小常数
w[0] = 1;
for (int i = 0; i < L; i++) {
r[i] = (r[i >> 1] >> 1) + (i & 1 ? L >> 1 : 0);
f[i] = a[r[i]];
}
for (int d = 1; d < L; d <<= 1) {
int w1 = qpow(type ? 3 : iv3, ((mod - 1) / d) >> 1);
for (int i = 1; i < d; i++) w[i] = w[i - 1] * w1 % mod;
for (int i = 0; i < L; i += (d << 1)) {
for (int j = 0; j < d; j++) {
int y = w[j] * f[i | j | d] % mod;
f[i | j | d] = f[i | j] + mod - y;
f[i | j] += y; // 注意顺序
}
}
if (d == (1 << 16)) for (int i = 0; i < L; i++) f[i] %= mod; // 16 次之后 y 最大是 18mod * mod, 还在 ull 范围内
}
int iv = qpow(L, mod - 2);
for (int i = 0; i < L; i++) a[i] = f[i] % mod * (type ? 1 : iv) % mod;
}
void Main(int tc) {
cin >> n >> m;
for (int i = 0; i <= n; i++) cin >> f[i];
for (int i = 0; i <= m; i++) cin >> g[i];
int L = 1;
while (L <= n + m) L <<= 1;
NTT(L, f, 1); NTT(L, g, 1);
for (int i = 0; i < L; i++) h[i] = (ll)f[i] * g[i] % mod;
NTT(L, h, 0);
for (int i = 0; i <= n + m; i++) cout << h[i] << " \n"[i == n + m];
}
}
4.6.2 任意模 NTT
这被称为 MTT。
4.6.2.1 三模 NTT
我们取 \(p_1 = 469762049, p_2 = 998244353, p_3 = 1004535809\)，先分别计算出结果模 \(p_1, p_2, p_3\) 的值，再用 exCRT 合并，我们答案模 \(p_1, p_2, p_3\) 的值分别为 \(h_1, h_2, h_3\)，先合并前两个得到答案模 \(p_1p_2\) 的值为 \(x\)，然后根据第三个求出答案为 \(x+yp_1p_2\)，注意此时我们可以边乘边取模，所以不需要 __int128，如果用普通 CRT 是需要的。注意此处模数互质，可以直接设 \(h_1 + kp_1 \equiv h_2 \pmod {p_2}\)，这样就可以直接移项然后求逆，而不需要 exgcd。
最后会用到 \(9\) 次 NTT，还是挺慢的。
下面给出 P4245 的 \(9\) 次 NTT 代码：
namespace Loop1st {
int n, m, p, r[N], f[N], g[N], h1[N], h2[N], h3[N];
int qpow(int a, int b, int mod) {
int res = 1;
for (; b; b >>= 1, a = (ll)a * a % mod) if (b & 1) res = (ll)res * a % mod;
return res;
}
void NTT(int L, int *a, bool type, int mod) {
int iv3 = qpow(3, mod - 2, mod);
static ull w[N], f[N]; // 用 ull 减小常数
w[0] = 1;
for (int i = 0; i < L; i++) {
r[i] = (r[i >> 1] >> 1) + (i & 1 ? L >> 1 : 0);
f[i] = a[r[i]] % mod;
}
for (int d = 1; d < L; d <<= 1) {
int w1 = qpow(type ? 3 : iv3, ((mod - 1) / d) >> 1, mod);
for (int i = 1; i < d; i++) w[i] = w[i - 1] * w1 % mod;
for (int i = 0; i < L; i += (d << 1)) {
for (int j = 0; j < d; j++) {
int y = w[j] * f[i | j | d] % mod;
f[i | j | d] = f[i | j] + mod - y;
f[i | j] += y; // 注意顺序
}
}
if (d == (1 << 16)) for (int i = 0; i < L; i++) f[i] %= mod; // 16 次之后 y 最大是 18mod * mod, 还在 ull 范围内
}
int iv = qpow(L, mod - 2, mod);
for (int i = 0; i < L; i++) a[i] = f[i] % mod * (type ? 1 : iv) % mod;
}
void mul(int L, int *f, int *g, int *h, int mod) {
static int a[N], b[N];
memcpy(a, f, sizeof a); memcpy(b, g, sizeof b);
NTT(L, a, 1, mod); NTT(L, b, 1, mod);
for (int i = 0; i < L; i++) h[i] = (ll)a[i] * b[i] % mod;
NTT(L, h, 0, mod);
}
void Main(int tc) {
cin >> n >> m >> p;
for (int i = 0; i <= n; i++) cin >> f[i];
for (int i = 0; i <= m; i++) cin >> g[i];
int L = 1;
while (L <= n + m) L <<= 1;
mul(L, f, g, h1, p1);
mul(L, f, g, h2, p2);
mul(L, f, g, h3, p3);
int iv1 = qpow(p1, p2 - 2, p2), iv12 = qpow((ll)p1 * p2 % p3, p3 - 2, p3);
for (int i = 0; i <= n + m; i++) {
ll x = (ll)(h2[i] - h1[i] % p2 + p2) * iv1 % p2 * p1 + h1[i];
ll y = (ll)(h3[i] - x % p3 + p3) * iv12 % p3;
cout << (y * p1 % p * p2 + x) % p << " \n"[i == n + m];
}
}
}
4.6.2.2 用单个大模数的 NTT
如果我们欣然接受 __int128，可以尝试找一个大一点的质数，这样就可以优化到 \(6\) 次甚至 \(3\) 次 NTT，比如我们用模数 \(25\times 2^{78}+1\)，原根是 \(3\)，不过这时取模是显然需要优化的，可以用 Montgomery 模乘。不过相信编译器的优化也是卡常的一环。代码以后补吧。
4.7 分治 NTT
不同于寻常 \(h_i = \sum_{j=0}^{i} f_jg_{i-j}\) 的卷积，现在我们需要求 \(f_i = \sum_{j=0}^{i-1}f_jg_{i-j}\)，即给定 \(f_0\) 和 \(g\)，求出 \(f\) 的所有系数。这被称为半在线卷积。
直接 NTT 是不可行的，因为 \(f_i\) 的值和 \(f_0, \dots, f_{i - 1}\) 有关。我们考虑（CDQ）分治，设当前分治区间为 \([l, r]\)，我们先递归计算 \(f_l, \dots, f_{mid}\) 的值，再计算 \(f_l, \dots, f_{mid}\) 对 \(f_{mid + 1}, \dots, f_r\) 的值，由于此时 \(f_l, \dots, f_{mid}\) 已知，我们可以直接用多项式乘法算出贡献，然后再递归处理 \(f_{mid + 1}, \dots, r\) 的值即可，时间复杂度 \(\mathcal{O}(n \log^2 n)\)。
下面给出 P4721 的代码。
namespace Loop1st {
int n, r[N];
int qpow(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = (ll)a * a % mod) if (b & 1) res = (ll)res * a % mod;
return res;
}
void NTT(int L, int *a, bool type) {
static ull w[N], f[N]; // 用 ull 减小常数
w[0] = 1;
for (int i = 0; i < L; i++) {
r[i] = (r[i >> 1] >> 1) + (i & 1 ? L >> 1 : 0);
f[i] = a[r[i]];
}
for (int d = 1; d < L; d <<= 1) {
int w1 = qpow(type ? 3 : iv3, ((mod - 1) / d) >> 1);
for (int i = 1; i < d; i++) w[i] = w[i - 1] * w1 % mod;
for (int i = 0; i < L; i += (d << 1)) {
for (int j = 0; j < d; j++) {
int y = w[j] * f[i | j | d] % mod;
f[i | j | d] = f[i | j] + mod - y;
f[i | j] += y; // 注意顺序
}
}
if (d == (1 << 16)) for (int i = 0; i < L; i++) f[i] %= mod; // 16 次之后 y 最大是 18mod * mod, 还在 ull 范围内
}
int iv = qpow(L, mod - 2);
for (int i = 0; i < L; i++) a[i] = f[i] % mod * (type ? 1 : iv) % mod;
}
int f[N], g[N];
void solve(int l, int r) { // a 需要用到 [l, mid], b 需要用到 [0, r - l]
if (l == r) return ;
int mid = (l + r) >> 1;
solve(l, mid);
int L = 1;
while (L <= r - l) L <<= 1;
static int a[N], b[N], c[N];
memset(a, 0, L << 2); memcpy(a, f + l, (mid - l + 1) << 2); memcpy(b, g, L << 2);
NTT(L, a, 1); NTT(L, b, 1);
for (int i = 0; i < L; i++) c[i] = (ll)a[i] * b[i] % mod;
NTT(L, c, 0);
for (int i = mid + 1; i <= r; i++) f[i] = (f[i] + c[i - l]) % mod;
solve(mid + 1, r);
}
void Main(int tc) {
cin >> n;
for (int i = 1; i < n; i++) cin >> g[i];
f[0] = 1;
solve(0, n - 1);
for (int i = 0; i < n; i++) cout << f[i] << " \n"[i == n - 1];
}
}
5. 多项式全家桶
这个部分我会将多项式板子封装成一个 namespace，默认使用 long long 类型，并放在这部分的最下方，参考了 lzyqwq 和 jiangly 的代码。
5.1 多项式加、减、乘法 P3803
不再赘述。
const ll mod = 998244353, iv3 = (mod + 1) / 3;
using poly = vector<ll>;
ll qpow(ll a, ll b) {
ll res = 1;
for (; b; b >>= 1, a = a * a % mod) if (b & 1) res = res * a % mod;
return res;
}
int deg(const poly &a) { return a.size() - 1; } // 如果不传 const & 就是 O(n) 的
poly rev(poly a) { reverse(a.begin(), a.end()); return a; }
void NTT(int L, poly &a, bool type = 0) {
vector<ull>f(L), w(L);
static int r[N];
w[0] = 1;
for (int i = 0; i < L; i++) {
r[i] = (r[i >> 1] >> 1) + (i & 1 ? L >> 1 : 0);
f[i] = a[r[i]];
}
for (int d = 1; d < L; d <<= 1) {
ll w1 = qpow(type ? 3 : iv3, ((mod - 1) / d) >> 1);
for (int i = 1; i < d; i++) w[i] = w[i - 1] * w1 % mod;
for (int i = 0; i < L; i += d << 1) {
for (int j = 0; j < d; j++) {
ll y = w[j] * f[i | j | d] % mod;
f[i | j | d] = f[i | j] + mod - y;
f[i | j] += y;
}
}
if (d == (1 << 16)) for (int i = 0; i < L; i++) f[i] %= mod;
}
ll iv = qpow(L, mod - 2);
for (int i = 0; i < L; i++) a[i] = f[i] % mod * (type ? 1 : iv) % mod;
}
poly modxk(poly a, int k) { a.resize(k); return a; }
poly operator + (poly a, poly b) {
int n = max(deg(a), deg(b)); a.resize(n + 1); b.resize(n + 1);
poly c = poly(n + 1);
for (int i = 0; i <= n; i++) c[i] = (a[i] + b[i]) % mod;
return c;
}
poly operator - (poly a, poly b) {
int n = max(deg(a), deg(b)); a.resize(n + 1); b.resize(n + 1);
poly c = poly(n + 1);
for (int i = 0; i <= n; i++) c[i] = (a[i] + mod - b[i]) % mod;
return c;
}
poly operator * (poly a, poly b) {
int n = deg(a), m = deg(b), L = 1;
while (L <= n + m) L <<= 1;
a.resize(L); b.resize(L);
NTT(L, a, 1); NTT(L, b, 1);
poly c = poly(L);
for (int i = 0; i < L; i++) c[i] = a[i] * b[i] % mod;
NTT(L, c, 0);
c.resize(n + m + 1);
return c;
}
5.2 多项式求逆 P4238
给定 \(f(x)\)，求 \(g(x)\) 使得 \(f(x)g(x)\equiv 1 \pmod {x^n}\)。
5.2.1 递推法
朴素的 \(\mathcal{O}(n^2)\) 递推求逆：
可以用半在线卷积优化至 \(\mathcal{O}(n \log^2 n)\)。
5.2.2 倍增法
我们设已经求出 \(f(x)\) 在模 \(x^{\lceil \frac n2\rceil}\) 意义下的逆元 \(f_0^{-1}(x)\)，。
由于 \(f(x)f_0^{-1}(x) \equiv 1 \pmod {x^{\lceil \frac n2\rceil}}, f(x)f^{-1}(x) \equiv 1 \pmod {x^{\lceil \frac n2\rceil}}\)，得到 \(f^{-1}(x)-f_0^{-1}(x) \equiv 0 \pmod{x^{\lceil \frac n2\rceil}}\)。
两边平方：
两边同乘 \(f(x)\)，移项：
递归或递推计算，时间复杂度 \(T(n) = T(\frac n2) + \mathcal{O}(n \log n) = \mathcal{O}(n \log n)\)。
poly modxk(poly a, int k) { a.resize(k); return a; }
poly inv(poly a) {
poly c = poly(1), d = poly(1);
c[0] = qpow(a[0], mod - 2); d[0] = 2;
for (int i = 1; (1 << (i - 1)) < a.size(); i++) {
c = modxk(c * (d - modxk(a, 1 << i) * c), 1 << i);
}
return c;
}
5.2.3 牛顿迭代法
以后再说
5.3 多项式除法 P4512
对于 \(f(x) = \sum_{i=0}^n f_ix^i\)，我们发现 \(x^nf(x^{-1})=\sum_{i=0}^nf_ix^{n-i}=\sum_{i=0}^n f_{n-i}x^i\)。记 \(f^R(x) = x^nf(x^{-1})\)。
我们有
两边同乘 \(x^n\)
考虑两边同乘 \(G^R(x)\) 在模 \(x^{n - m + 1}\) 意义下的逆 \(H(x)\)
而 \(\deg Q^R(x)=n-m<n-m+1\)，所以
然后根据 \(Q^R(x)\) 求出 \(Q(x)\)，再用 \(F(x) - G(x)Q(x)\) 得到 \(R(x)\)。
int deg(poly a) { return a.size() - 1; }
poly rev(poly a) { reverse(a.begin(), a.end()); return a; }
poly inv(poly a, int n) {
poly c = poly(1), d = poly(1);
c[0] = qpow(a[0], mod - 2); d[0] = 2;
for (int i = 1; (1 << (i - 1)) < n; i++) {
c = modxk(c * (d - modxk(a, 1 << i) * c), 1 << i);
}
return modxk(c, n);
}
poly operator / (poly a, poly b) {
int n = deg(a), m = deg(b);
return rev(modxk(rev(a) * inv(rev(b), n - m + 1), n - m + 1));
}
poly operator % (poly a, poly b) {
return modxk(a - a / b * b, deg(b));
}
5.4 多项式求导、积分
很显然。
poly direv(poly a) {
poly b = poly(deg(a));
for (int i = 0; i < deg(a); i++) b[i] = a[i + 1] * (i + 1) % mod;
return b;
}
poly integr(poly a) {
poly b = poly(deg(a) + 2);
for (int i = 1; i <= deg(a) + 1; i++) b[i] = a[i - 1] * qpow(i, mod - 2) % mod;
return b;
}
5.5 多项式 \(\ln\) P4725
只需要知道 \(B'(x) \bmod x^{n-1}\) 就可以积分得到 \(B(x) \bmod x^n\)，那么只需要知道 \(\dfrac{A'(x)}{A(x)} \bmod x^{n-1}\)，多项式求逆即可。
这里保证了 \(a_0 = 1\) 的深层原因，是因为如果 \(a_0 \ne 1\) 的话 \(\ln\) 没有意义。
poly ln(poly a, int n) {
a = modxk(a, n);
assert(a[0] == 1);
return integr(modxk(direv(a) * inv(a, n - 1), n - 1));
}
5.x 封装好的板子
持续更新中……
参考资料
快速傅里叶变换
[多项式 I：拉格朗日插值与快速傅里叶变换 Alex_Wei](多项式 I：拉格朗日插值与快速傅里叶变换 - qAlex_Weiq - 博客园)
【学习笔记】多项式 1：基础操作
NTT与多项式全家桶
浙公网安备 33010602011771号