算法学习笔记：多项式
多项式基础
为了与幂次相对应，下文中系数均采用 \(\text{0-based}\) 下标计数。
多项式的表示方法
我们习惯于将一个 \(n\) 次多项式 \(f(x)\) 表示为 \(f(x)=\sum_{i=0}^{n-1}a_ix^i\)。那么有没有其他表示方法呢？
我们不妨考虑代入 \(n\) 个互不相同的 \(x_i\) 得到的 \(n\) 个点值 \((x_i,y_i)\)，很自然地猜测它们可以唯一确定这个 \(n\) 次多项式。可以发现 \(n\) 个点值对应一个关于系数的 \(n\) 元一次方程，显然有唯一解，因此这 \(n\) 个点值可以唯一确定这个 \(n\) 次多项式。我们称这种表示方法为点值表示法。而我们常用的 \(f(x)=\sum_{i=0}^{n-1}a_ix^i\) 的表示方法则被称作系数表示法。
既然这两种表示法都可以唯一确定一个多项式，那么必然可以相互转化。具体地，系数表示法转化为点值表示法的过程就是求值，点值表示法转化为系数表示法的过程就是插值。
离散卷积
在这里较为不严谨地给出离散卷积的概念。
给定函数 \(f,g\)，形如
的式子就是 \(f\) 与 \(g\) 的卷积，其中 \(\otimes\) 是某种运算。可以将其写作 \(h=f*g\)。
考虑多项式乘法，给出多项式 \(F(x),G(x)\)，令 \(f(n)=[x^n]F(x)\)，\(g(n)=[x^n]G(x)\) 则
因此多项式乘法本质上就是其系数序列的加法卷积。
拉格朗日插值
我们简单研究一下插值。一个显然的方法就是 \(O(n^3)\) 高斯消元，而下面介绍的拉格朗日插值可以做到 \(O(n^2)\) 插值。
现在给出 \(n\) 个点值 \((x_i,y_i)\)，保证 \(x_i\) 互不相同，我们想要找到对应的多项式 \(f(x)\) 的系数表示法。
拉格朗日插值的思想就是构造 \(n\) 个 \(n\) 次多项式 \(g_i(x)\)，满足 \(g_i(x_i)=y_i\)，\(g_i(x_j)=0(j\neq i)\)，那么所求 \(f(x)\) 即为 \(\sum_{i=0}^{n-1}g_i(x)\)。考虑如何构造 \(g_i(x)\)。显然 \(x_j(j\neq i)\) 给出了 \(g_i(x)\) 的 \(n-1\) 个根，因此它必然包含 \(\prod_{j\neq i}(x-x_j)\) 的因式，而我们还要求 \(g_i(x_i)=y_i\)，调整系数，乘上 \(\frac{y_i}{g_i(x_i)}\) 即可。综上
因此
实现时可以 \(O(n^2)\) 求出 \(\prod_{j=0}^{n-1}(x-x_j)\) 的各项系数，对于每个 \(i\) 直接 \(O(n)\) 除掉 \(x-x_i\)，再调整系数，最后对所有 \(i\) 系数相加即可。时间复杂度为 \(O(n^2)\)。
Luogu P4781 【模板】拉格朗日插值 只要求单点插值，直接带入拉格朗日插值公式计算即可，依然是 \(O(n^2)\) 的。
横坐标是连续整数的快速单点插值
若 \(x_i\) 为连续整数，则我们可以给出 \(O(n)\) 的单点插值的方法。
设 \(x_i=k+i\)，则拉格朗日插值公式变为
注意到分子是经典的去除一个单点的形式，考虑预处理前后缀积，即令 \(pre_i=\prod_{j=0}^i(x-j-k)\)，\(suf_i=\prod_{j=i}^{n-1}(x-j-k)\)。而对于分母，\(j<i\) 部分的乘积为 \(i!\)，\(j>i\) 部分的乘积为 \((-1)^{n-i-1}(n-i-1)!\)。容易得到
预处理阶乘逆元即可 \(O(n)\) 计算。
快速傅里叶变换
快速傅里叶变换（FFT）是多项式的基石。这里默认读者具有良好的线性代数与复数知识储备。
离散傅里叶变换
离散傅里叶变换（DFT）是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其 DTFT 的频域采样。在 OI 中，其恰好能用来将多项式的系数表示法转化为点值表示法。
序列 \(\{x_i\}_{i=0}^{n-1}\) 的 DFT 为
若将 \(\{x_i\}_{i=0}^{n-1}\) 视为多项式 \(f(x)\) 的系数序列，则我们有
通常以符号 \(\mathcal{F}\) 表示该变换，即 \(\hat{x}=\mathcal{F}x\)。而它的逆离散傅里叶变换为
也记作 \(x=\mathcal{F}^{-1}\hat{x}\)。
容易发现，离散傅里叶变换可以看作一个 \(n-1\) 次多项式在 \(n\) 个单位根处求值。
分治实现 FFT
FFT 是一种高效实现 DFT 的算法，不过其中的一些细节又与 DFT 略有区别，后文将会提及。
我们的目标就是，在可接受的时间复杂度下，快速求出一个 \(n-1\) 次多项式的 \(n\) 个点值。
FFT 的核心思想就是分治。我们考虑 \(f(x)\) 拆成一个偶函数 \(f_e(x)\) 和一个奇函数 \(f_o(x)\) 之和，也就是令 \(f_e(x)=a_0x^0+a_2x^2+\cdots\)，\(f_o(x)=a_1x^1+a_3x^3+\cdots\)，此时
这样并不能得到很好的复杂度，因为 \(f_e(x)\) 和 \(f_o(x)\) 依然是 \(n-1\) 次的多项式。但注意到项数减半，考虑换元，令 \(t=x^2\)，则 \(f_e(x)=a_0t^0+a_2t^1+\cdots\)，\(f_o(x)=x(a_1t^0+a_3t^1+\cdots)\)，换句话说
此时 \(f_e(x)\) 和 \(f_o(x)\) 降成了 \(\frac{n}{2}\) 次多项式，我们做到了规模减半。
接下来，由偶函数和奇函数的性质，考虑代入 \(\frac{n}{2}\) 对互为相反数的 \((x_i,-x_i)\)，因为这样会使得它们的 \(f_e(x^2)\) 和 \(f_o(x^2)\) 值是相同的。而我们需要保证子问题也能找到这样成对的相反数，也即要找到 \(n\) 个互不相同的 \(x_i\)，它们的 \(n\) 次幂相等，显然所有 \(n\) 次单位根就满足这个性质，因此我们代入 \(\omega_n^k\) 和 \(-\omega_n^k=\omega_n^{k+\frac{n}{2}}\)：
再利用 \(\omega_{n}^{2k}=\omega_{\frac{n}{2}}^k\)：
我们以 \(n=1\) 作为递归边界，在每个子问题中对系数 \(a_i\) 按 \(i\) 的奇偶性左右分组，再递归分治处理，最后统计答案，就可以分治实现 FFT 了！
\(T(n)=T(\frac{n}{2})+O(n)\)，由主定理，时间复杂度为 \(O(n\log{n})\)。
注意上述过程都是在 \(n=2^k(k\in\mathbb{N})\) 的前提下进行的。因此具体实现时，要把 \(n\) 向上补成 \(2^k\)。
我们还需要支持复数运算，可以手写结构体，也可以使用 STL complex。
分治实现的代码：
void FFT(complex<double> *a, int n) {
if (n == 1) return;
int mid = n >> 1;
for (int i = 0; i < n; i += 2)
tmp[i >> 1] = a[i], tmp[(i >> 1) + mid] = a[i + 1];
for (int i = 0; i < n; ++i) a[i] = tmp[i];
FFT(a, mid, tp), FFT(a + mid, mid, tp);
w[0] = { 1, 0 }, w[1] = { cos(PI * 2 / n), sin(PI * 2 / n) };
for (int i = 2; i < mid; ++i) w[i] = w[i - 1] * w[1];
for (int i = 0; i < mid; ++i)
tmp[i] = a[i] + w[i] * a[i + mid],
tmp[i + mid] = a[i] - w[i] * a[i + mid];
for (int i = 0; i < n; ++i) a[i] = tmp[i];
}
倍增实现 FFT/蝴蝶变换
分治实现 FFT 需要递归处理，效率较慢，我们希望找到非递归的实现方式。
考察 FFT 对系数位置的变换。对于一个规模为 \(n\) 的问题，位置 \(p\) 上的系数会在 \(p\) 为偶数时被分到左侧，为奇数时分到右侧。进一步地，位置 \(p\) 在第 \(i\) 次递归的方向决定了它最终下标二进制表示下第 \(\log_2n-i\) 位的值，这等价于 \(p\) 最终下标二进制表示下第 \(\log_2n-i\) 位的取值，就是 \(p\) 二进制表示下第 \(i\) 位的取值。因此，系数 \(a_p\) 在 \(\log_2n\) 次递归后的下标就是 \(p\) 的二进制表示反转的结果，不妨记作 \(rev_p\)。
显然 \(rev_p\) 可以 \(O(n)\) 递推：
我们在一开始就做一次位逆序置换，即令 \(a_i\leftarrow a_{rev_i}\)。然后考虑如何自底向上合并。假设我们已经完成了所有规模为 \(n\) 的变换，而想要转而完成所有规模为 \(2n\) 的变换，回到先前的分治式子
此时，\(f_e(\omega_n^k)\) 和 \(f_o(\omega_n^k)\) 的值分别存在下标为 \(k\) 和 \(k+n\) 的位置，而我们想求出的 \(f(\omega_{2n}^k)\) 和 \(f(\omega_{2n}^{k+n})\) 恰好将分别存在下标为 \(k\) 和 \(k+n\) 的位置，所以我们直接在原数组上覆写即可。具体来说，令 \(x=a_k\)，\(y=\omega_{2n}^ka_{k+n}\)，然后令 \(a_k\leftarrow x+y\)，\(a_{k+n}\leftarrow x-y\) 就完成了一对值的求解。这个过程就称为蝴蝶变换。
我们理一下算法的过程：先对系数数组做位逆序置换，然后枚举问题规模 \(2k(1\leq k<n)\)，枚举每个子问题 \(2ki(0\leq 2ki<n)\)，枚举子问题中的对应位置 \(x,y(x=2ki+j,y=2ki+j+k)\)，记 \(tx=a_x\)，\(ty=\omega_{2k}^ja_y\)，令 \(a_x\leftarrow tx+ty\)，\(a_y\leftarrow tx-ty\) 即可。
由 \(rev\) 的对称性，做位逆序置换时只需在 \(i<rev_i\) 时 \(\operatorname{swap}(a_i,a_{rev_i})\)。
时间复杂度还是 \(O(n\log{n})\) 的，但相比分治实现效率高了不少。
蝴蝶变换的代码：
void FFT(complex<double> *a, int n) {
int mid = n >> 1;
for (int i = 1; i < n; ++i) {
rev[i] = (rev[i >> 1] >> 1) + (i & 1 ? mid : 0);
if (i < rev[i]) swap(a[i], a[rev[i]]);
}
for (int k = 1; k < n; k <<= 1) {
w[0] = { 1, 0 }, w[1] = { cos(PI / k), sin(PI / k) };
for (int i = 2; i < k; ++i) w[i] = w[i - 1] * w[1];
for (int i = 0; i < n; i += k << 1) for (int j = 0; j < k; ++j) {
complex<double> x = a[i | j], y = w[j] * a[i | j | k];
a[i | j] = x + y, a[i | j | k] = x - y;
}
}
}
多项式求值和线性代数
在讲解快速傅里叶逆变换之前，我们先讲一下多项式求值和线性代数之间的关系。
其实多项式求值很显然可以用矩阵乘法描述。例如，把多项式 \(f(x)=\sum_{i=0}^{n-1}a_ix^i\) 对 \(x_0,x_1,\cdots,x_{n-1}\) 求值，就可以表示为
左侧的矩阵就是范德蒙德矩阵。
那么 FFT 也可以用这种方式表示：
快速傅里叶逆变换
设 FFT 中的范德蒙德矩阵为 \(T\)，即
由前面的矩阵乘法表示，容易想到逆变换就是左乘上这个矩阵的逆矩阵 \(T^{-1}\)。
考虑利用拉格朗日插值公式求逆。观察拉格朗日插值公式，容易得出
分子和分母都是
的形式，考虑研究其性质。由代数基本定理，\(\prod_{i=0}^{n-1}(x-\omega_n^i)=x^n-1\)，因此化为
模拟短除法，可以得到
因此 \([x^i]\prod_{k\neq j}(x-\omega_n^k)=(\omega_n^j)^{n-1-i}\)，\(g(\omega_n^j) =n(\omega_n^j)^{n-1}\)，于是
即我们所求的逆矩阵
所以 IFFT 的过程和 FFT 没有很大的区别。我们只需要将 FFT 中代入的单位根变为原来的共轭，并在最后对每一项系数除以 \(n\) 即可。
为了方便，我们可以把 FFT 与 IFFT 放在一个函数里，用 \(tp=1\) 表示是 FFT，用 \(tp=-1\) 表示是 IFFT，那么这个 \(tp\) 就恰好可以乘到单位根的虚部上来转成共轭。
FFT/IFFT 的代码：
void FFT(complex<double> *a, int n, int tp = 1) {
int mid = n >> 1;
for (int i = 1; i < n; ++i) {
rev[i] = (rev[i >> 1] >> 1) + (i & 1 ? mid : 0);
if (i < rev[i]) swap(a[i], a[rev[i]]);
}
for (int k = 1; k < n; k <<= 1) {
w[0] = { 1, 0 }, w[1] = { cos(PI / k), tp * sin(PI / k) };
for (int i = 2; i < k; ++i) w[i] = w[i - 1] * w[1];
for (int i = 0; i < n; i += k << 1) for (int j = 0; j < k; ++j) {
complex<double> x = a[i | j], y = w[j] * a[i | j | k];
a[i | j] = x + y, a[i | j | k] = x - y;
}
}
if (tp == -1) for (int i = 0; i < n; ++i) a[i] /= n;
}
DFT 和 FFT
回顾前面的内容，可以总结出 DFT 和 FFT 之间的几点差别：
- FFT 要求 \(n\) 是 \(2\) 的整数次幂，DFT 则不做要求。
- DFT 与 FFT 代入单位根的顺序相反。
快速数论变换
数论变换（NTT）是 DFT 在数论基础上的实现，快速数论变换（FNTT）则是 FFT 在数论基础上的实现。更具体地，FNTT 是在模 \(p\) 意义下的整数域进行的 FFT。
对次数为 \(n-1\) 的多项式进行变换，我们需要找到 \(n\) 次单位根 \(a\) 满足 \(\delta_p(a)=n\)，这要求模数 \(p\) 为一个质数，且 \(2^k\mid p-1\)，其中 \(2^k\leq\) 最大的可能 NTT 长度。
根据原根的性质，\(\delta_p(g)=p-1\)，即 \(g\) 的 \(0\sim p-2\) 次幂互不相同，因此 \(g^k\) 和 \(\omega_{p-1}^k(0\leq k<p-1)\) 形成的域是同构的。进一步地，\(\omega_{n}\) 等价于 \(g^{\frac{p-1}{n}}\)，其 \(0\sim n-1\) 次幂互不相同，满足我们做 FFT 时所需的性质。
因此我们只需要将 FFT 时的 \(\omega_n\) 换成 \(g^{\frac{p-1}{n}}\)，\(\overline{\omega_n}\) 换成 \((g^{-1})^\frac{p-1}{n}\) 就变成了 NTT 啦！由于没有了浮点数运算，NTT 的效率比 FFT 高不少。
下面列出一些 NTT 常见大模数：
- \(998244353=7\times 17\times 2^{23}+1\)，有原根 \(3\)。
- \(1004535809=479\times 2^{21}+1\)，有原根 \(3\)。
- \(469762049=7\times 2^{26}+1\)，有原根 \(3\)。
NTT 的代码（\(p=998244353\)）：
const int g1 = 3, g2 = (MOD + 1) / g1;
void NTT(int *a, int n, int tp = 1) {
int mid = n >> 1;
for (int i = 1; i < n; ++i) {
rev[i] = (rev[i >> 1] >> 1) + (i & 1 ? mid : 0);
if (i < rev[i]) swap(a[i], a[rev[i]]);
}
for (int k = 1; k < n; k <<= 1) {
w[0] = 1, w[1] = qpow(tp == 1 ? g1 : g2, (MOD - 1) / (k << 1));
for (int i = 2; i < k; ++i) w[i] = 1ll * w[i - 1] * w[1] % MOD;
for (int i = 0; i < n; i += k << 1) for (int j = 0; j < k; ++j) {
int x = a[i | j], y = 1ll * w[j] * a[i | j | k] % MOD;
a[i | j] = (x + y) % MOD, a[i | j | k] = (x - y + MOD) % MOD;
}
}
if (tp == -1) {
int iv = qpow(n, MOD - 2);
for (int i = 0; i < n; ++i) a[i] = 1ll * a[i] * iv % MOD;
}
}
FFT/NTT 的应用
多项式乘法
题意：给出一个次数为 \(n\) 的多项式 \(f(x)\) 和一个次数为 \(m\) 的多项式 \(g(x)\)，求它们的卷积。\(1\leq n,m\leq 10^6\)，\(0\leq [x^i]f(x),[x^i]g(x)\leq 9\)。
下面视 \(n,m\) 同阶。按照卷积的定义暴力做显然是 \(O(n^2)\) 的。不妨从点值表示法的角度考虑，容易发现两个多项式的卷积的点值表示即两个多项式的点值对应相乘。所以我们用 FFT 把 \(f(x),g(x)\) 转成点值表示法，\(O(n)\) 点值相乘，再 IFFT 转回系数表示法即可。时间复杂度为 \(O(n\log{n})\)。
进一步地，我们发现卷积后的系数不超过 \(nV^2\)，其中 \(V\) 为多项式系数的值域，所以我们可以把 FFT 改成对 \(998244353\) 取模的 NTT，从而优化常数。
代码中的函数calc_len 和 clr 在后面的代码中会用到：
int calc_len(int n) { return 1 << (int)ceil(log2(n)); }
void clr(int *a, int n) { memset(a, 0, n << 2); }
cin >> n >> m; ++n, ++m;
for (int i = 0; i < n; ++i) cin >> f[i];
for (int i = 0; i < m; ++i) cin >> g[i];
int len = calc_len(n + m - 1);
NTT(f, len), NTT(g, len);
for (int i = 0; i < len; ++i) f[i] = 1ll * f[i] * g[i] % MOD;
NTT(f, len, -1);
for (int i = 0; i < n + m - 1; ++i) cout << f[i] << " \n"[i == n + m - 2];
分治 FFT/NTT
半在线卷积
题意：给定长度为 \(n-1\) 的 序列 \(g\)，求出一个长度为 \(n\) 的序列 \(f\)，满足
答案对 \(998244353\) 取模。\(2\leq n\leq 10^5\)。
对于一类形如 \(f_i=\sum_{j=1}^if_{i-j}g_j\) 的递推式，\(f_i\) 的值依赖于 \(f[0,i)\) 的值，我们称其为半在线卷积问题。
这是一个在线问题，考虑用 CDQ 分治将其转化为离线问题。具体来说，假设我们求出了 \(f[l,mid]\)，现在想要计算其对 \(f[mid+1,r]\) 的贡献。对于 \(f_i(i\in[mid+1,r])\)，容易计算得到 \(f[l,mid]\) 对它的贡献为 \(\sum_{j=l}^{mid}f_{j}g_{i-j}\)，这就变成了很简单的卷积问题了，我们对 \(f[l,mid]\) 和 \(g[0,r-l]\) 做卷积即可得到贡献序列。
代码：
void solve(int l, int r) {
if (l == r) return;
static int a[N], b[N];
int mid = (l + r) >> 1;
solve(l, mid);
copy(f + l, f + mid + 1, a), copy(g, g + r - l + 1, b);
int len = calc_len(r - l + 1);
NTT(a, len), NTT(b, len);
for (int i = 0; i < len; ++i) a[i] = 1ll * a[i] * b[i] % MOD;
NTT(a, len, -1);
for (int i = mid + 1; i <= r; ++i) f[i] = (f[i] + a[i - l]) % MOD;
clr(a, len), clr(b, len);
solve(mid + 1, r);
}
单项式乘积
给出 \(n\) 个一次多项式 \(f_i(x)\)，我们希望计算它们的乘积多项式 \(\prod_{i=1}^n{f_i(x)}\)。
如果我们从左往右 NTT 合并，复杂度会是 \(O(n^2\log{n})\) 的。
考虑分治计算，令 \(g_{l,r}(x)=\prod_{i=l}^rf_i(x)\)，于是显然 \(g_{l,r}(x)=g_{l,mid}(x)g_{mid+1,r}(x)\)，递归合并即可。\(T(n)=2T(\frac{n}{2})+O(n\log{n})\)，由主定理，时间复杂度为 \(O(n\log^2{n})\)。
乍一看这似乎不太合理，为什么把从左到右合并改成分治合并就能得到更优的时间复杂度呢？究其原因，从左到右合并就是不断求一个 \(k\) 次多项式和一个单项式的乘积的过程，NTT 很浪费时间，甚至不如根据定义暴力合并。而分治则平衡了左右两侧的多项式的规模，使得复杂度得到了优化。
多项式牛顿迭代
泰勒级数和麦克劳林级数
这里简单介绍学习牛顿迭代的一些基础知识。读者需要有简单的微积分基础。
对于一个光滑的函数 \(f(x)\)，我们有
上式称为 \(f(x)\) 在 \(x=a\) 处的泰勒展开式。特别地，在 \(x=0\) 处的泰勒展开式被称为麦克劳林展开式。
这里列出一些常见的麦克劳林展开式，后文会用到一部分：
一般的牛顿迭代
一般的牛顿迭代按以下方式求出方程 \(f(x)=0\) 的根：
- 设我们得到的近似根为 \(x_0\)。
- 将 \(f(x)\) 在 \(x=x_0\) 处做泰勒展开。
- 取前两项的系数来近似成 \(f(x)\)，即 \(f(x)=f(x_0)+f'(x_0)(x-x_0)=0\)。
- 解得 \(x=x_0-\frac{f(x_0)}{f'(x_0)}\)，将其作为新的近似根继续迭代下去。
容易发现，一般的牛顿迭代的本质上就是不断求切线的过程。
多项式牛顿迭代
我们可以用牛顿迭代的方式解决这样一类问题：给出关于多项式 \(f(x)\) 的函数 \(g(f(x))\)，求一个多项式 \(f(x)\) 使得 \(g(f(x))\equiv 0\pmod{x^n}\) 成立。
与一般的牛顿迭代不同，由于多项式的一些特性，用牛顿迭代求解这类问题并无精度问题。
假设我们求出了 \(f_0(x)\) 使得 \(g(f_0(x))\equiv 0\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\)，我们把 \(g(f(x))\) 在 \(f(x)=f_0(x)\) 处做泰勒展开
注意到，\(f_0(x)\) 和 \(f(x)\) 的前 \(\left\lceil\frac{n}{2}\right\rceil\) 项相同，即
所以对于 \(n\geq 2\)，\(\frac{g^{(n)}(f_0(x))}{n!}(f(x)-f_0(x))^n\equiv 0\pmod{x^n}\)，于是我们的泰勒级数依然只需保留前两项：
类似解得
我们只需要在初始时给出 \(\bmod{x^1}\) 的解就能进行多项式牛顿迭代了。
有一些细节需要注意：
- 由于 \(g(f_0(x))\equiv 0\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\)，\(g'(f_0(x))\) 中次数 \(\geq\left\lceil\frac{n}{2}\right\rceil\) 的项会被 \(\bmod{x^n}\) 截断，所以 \(g'(f_0(x))\) 只需要做到 \(\bmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\) 的精度。
- \(g'(f_0(x))\) 是将 \(f_0(x)\) 视作主元，而非 \(x\)，即分母要求的是 \(\frac{\text{d}g(f_0(x))}{\text{d}f_0(x)}\)。
显然多项式牛顿迭代的精度成平方增长，这通常允许我们使用倍增求解。读者可以结合接下来的一些例子体会。
多项式乘法逆
题意：给定一个 \(n-1\) 次多项式 \(f(x)\)，求出 \(g(x)\) 使得 \(f(x)g(x)\equiv 1\pmod{x^n}\)，系数对 \(998244353\) 取模。\(1\leq n\leq 10^5\)。
这里给出两种求解方法：
方法一：平方法
显然 \(\bmod{x^1}\) 的解就是 \(a_0^{-1}\bmod{998244353}\)。接下来假设我们求出了 \(g_0(x)\) 使得 \(g_0(x)f(x)\equiv 1\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}\)，来快乐地推一下式子：
那么倍增求解到最小的 \(k\) 使得 \(2^k\geq n\) 即可。\(T(n)=T(\frac{n}{2})+O(n\log{n})\)，由主定理，时间复杂度为 \(O(n\log{n})\)。
方法二：牛顿迭代
我们相当于要让 \(h(g(x))\equiv g(x)f(x)-1\equiv 0\)。直接代入多项式牛顿迭代的公式求解：
这和前面推的式子是一样的，同样 \(O(n\log{n})\) 倍增求解即可。
回顾多项式牛顿迭代的推导过程，读者可以体会平方法本质上与牛顿迭代的相似性。
倍增法的细节
倍增法有一些很重要的细节，不加注意就会出现问题：
- 注意考虑好 NTT/FFT 的变换长度。
- 记得清空！！！记得清空！！！记得清空！！！
代码：
void poly_inv(const int *a, int *res, int n) {
static int ta[N], tr[N], b[N];
b[0] = qpow(a[0], MOD - 2);
for (int k = 1; k >> 1 < n; k <<= 1) {
clr(ta, k << 1), clr(tr, k << 1);
copy(a, a + min(n, k), ta), copy(b, b + min(n, k), tr);
NTT(ta, k << 1), NTT(tr, k << 1);
for (int i = 0; i < k << 1; ++i)
b[i] = (2 - 1ll * ta[i] * tr[i] % MOD + MOD) % MOD * tr[i] % MOD;
NTT(b, k << 1, -1);
clr(b + k, k);
}
copy(b, b + n, res);
}
多项式开根
Luogu P5205 【模板】多项式开根 | Luogu P5277 【模板】多项式开根（加强版）
题意：给定一个 \(n-1\) 次多项式 \(f(x)\)，求出 \(g(x)\) 使得 \(g(x)^2\equiv f(x)\pmod{x^n}\)，系数对 \(998244353\) 取模。非加强版保证 \(a_0=1\)，加强版则只保证 \(a_0\) 是模 \(998244353\) 意义下的二次剩余。\(1\leq n\leq 10^5\)。
多项式开根同样可以平方法做，读者可以自行推导。这里直接给出牛顿迭代做法。
构造函数 \(h(g(x))=g(x)^2-f(x)\)，代入公式计算：
对分母求逆即可计算上式。对于初值，若保证 \(a_0=1\)，直接取 \(g_0(x)\equiv 1\pmod{x^1}\) 即可，否则使用 Cipolla 求出模意义下开根的最小解即可。时间复杂度 \(O(n\log{n})\)。
加强版代码：
void poly_sqrt(const int *a, int *res, int n) {
int k;
static tmpa[N], tmpr[N], inv[N];
res[0] = cipolla(a[0]);
chk_min(res[0], MOD - res[0]);
for (k = 1; k >> 1 < n; k <<= 1) {
clr(tmpa, k << 1), clr(tmpr, k << 1), clr(inv, k << 1);
for (int i = 0; i < k << 1; ++i) tmpa[i] = tmpr[i] = inv[i] = 0;
copy(a, a + k, tmpa), copy(res, res + k, tmpr);
poly_inv(tmpr, inv, k), NTT(inv, k << 1);
NTT(tmpa, k << 1);
for (int i = 0; i < k << 1; ++i) res[i] = 1ll * tmpa[i] * inv[i] % MOD;
NTT(res, k << 1, -1);
for (int i = 0; i < k; ++i) res[i] = 1ll * (res[i] + tmpr[i]) % MOD * iv2 % MOD;
for (int i = k; i < k << 1; ++i) res[i] = 0;
}
for (int i = n; i < k; ++i) res[i] = 0;
}
多项式除法
题意：给定一个 \(n\) 次多项式 \(f(x)\) 和一个 \(m\) 次多项式 \(g(x)\)，求出一个 \(n-m\) 次多项式 \(q(x)\) 和一个 \(m-1\) 次多项式 \(r(x)\)，使得 \(f(x)\equiv g(x)q(x)+r(x)\pmod{x^n}\)。\(1\leq m<n\leq 10^5\)。
这题的做法还是比较牛的。
若保证能整除，直接求逆即可，因此考虑去除 \(f(x)\) 的余多项式。前面的题中，我们经常用 \(\bmod{x^k}\) 来消去无用的项，这种做法可以消掉高次项。但余多项式显然在低次项的位置，考虑反转系数。令 \(f_{R}(x)=x^nf(\frac{1}{x})\)，我们尝试推一下式子：
那么我们就可以多项式求逆求出 \(q_R(x)\)，反转系数后得到 \(q(x)\)，最后 \(r(x)\equiv f(x)-g(x)q(x)\) 即可。时间复杂度 \(O(n\log{n})\)。
实现细节还是比较多的。代码：
void poly_div(const int *a, const int *b, int *q, int *r, int n, int m) {
static int ar[N], br[N], inv[N], tq[N];
reverse_copy(a, a + n, ar), reverse_copy(b, b + m, br);
for (int i = n - m + 1; i < m; ++i) br[i] = 0;
poly_inv(br, inv, n - m + 1);
int len = calc_len(n * 2 - m);
NTT(ar, len), NTT(inv, len);
for (int i = 0; i < len; ++i) ar[i] = 1ll * ar[i] * inv[i] % MOD;
NTT(ar, len, -1);
reverse_copy(ar, ar + n - m + 1, q), reverse_copy(ar, ar + n - m + 1, tq);
clr(ar, len), clr(br, len), clr(inv, len);
len = calc_len(n);
copy(b, b + m, br);
NTT(br, len), NTT(tq, len);
for (int i = 0; i < len; ++i) r[i] = 1ll * br[i] * tq[i] % MOD;
NTT(r, len, -1);
for (int i = 0; i < m - 1; ++i) r[i] = (a[i] - r[i] + MOD) % MOD;
clr(br, len), clr(tq, len);
}
多项式对数函数
Luogu P4725 【模板】多项式对数函数（多项式 ln）
题意：给定一个 \(n-1\) 次多项式 \(f(x)\)，求出一个多项式 \(g(x)\) 使得 \(g(x)\equiv \ln(f(x))\pmod{x^n}\)，系数对 \(998244353\) 取模。保证 \(a_0=1\)，\(1\leq n\leq 10^5\)。
\(\ln\) 这个函数比较棘手，但注意到其导数具有很好的性质，两边同时求导再积分就做完了：
读者可能会想到求导/积分会带来常数项的损失，不过题目钦定了 \(a_0=1\)，所以我们无需担心这一点。
那这又导出了另一个问题：如果不保证 \(a_0=1\) 呢？我们声称
\(f(x)\) 的对数多项式存在当且仅当 \(a_0=1\)。
考虑用麦克劳林级数对 \(\ln(f(x))\) 展开：
考察其常数项，即
显然该级数收敛当且仅当 \(a_0=1\)。得证。
代码：
void poly_diff(const int *a, int *res, int n) {
for (int i = 1; i < n; ++i) res[i - 1] = 1ll * a[i] * i % MOD;
res[n - 1] = 0;
}
void poly_intg(const int *a, int *res, int n) {
for (int i = 1; i < n; ++i) res[i] = 1ll * a[i - 1] * qpow(i, MOD - 2) % MOD;
res[0] = 0;
}
void poly_ln(const int *a, int *res, int n) {
static int df[N], inv[N];
poly_inv(a, inv, n), poly_diff(a, df, n);
int len = calc_len(n << 1);
NTT(inv, len), NTT(df, len);
for (int i = 0; i < len; ++i) df[i] = 1ll * df[i] * inv[i] % MOD;
NTT(df, len, -1), poly_intg(df, res, n);
clr(inv, len), clr(df, len);
}
多项式指数函数
Luogu P4726 【模板】多项式指数函数（多项式 exp）
题意：给定一个 \(n-1\) 次多项式 \(f(x)\)，求出一个多项式 \(g(x)\) 使得 \(g(x)\equiv e^{f(x)}\pmod{x^n}\)，系数对 \(998244353\) 取模。保证 \(a_0=0\)，\(1\leq n\leq 10^5\)。
\(\exp\) 函数怎么求导/积分都是本身，怎么办？注意到其与 \(\ln\) 互为反函数，对等式两边同时取 \(\ln\)，得到
这个形式直接上牛顿迭代：
倍增 \(O(n\log{n})\) 做即可。
和多项式对数函数类似，考虑 \(\exp(f(x))\) 的麦克劳林级数展开即可证明
\(f(x)\) 的指数多项式存在当且仅当 \(a_0=0\)。
代码：
void poly_exp(const int *a, int *res, int n) {
static int ta[N], tr[N], b[N], lnr[N];
b[0] = 1;
for (int k = 1; k >> 1 < n; k <<= 1) {
clr(ta, k << 1), clr(tr, k << 1);
copy(a, a + min(n, k), ta), copy(b, b + min(n, k), tr);
poly_ln(tr, lnr, k);
for (int i = 0; i < k << 1; ++i) ta[i] = (ta[i] - lnr[i] + MOD) % MOD;
ta[0] = (ta[0] + 1) % MOD;
NTT(ta, k << 1), NTT(tr, k << 1);
for (int i = 0; i < k << 1; ++i) b[i] = 1ll * ta[i] * tr[i] % MOD;
NTT(b, k << 1, -1);
clr(b + k, k);
}
copy(b, b + n, res);
}
多项式幂函数
Luogu P5245 【模板】多项式快速幂 | Luogu P5273 【模板】多项式幂函数（加强版）
题意：给定一个 \(n-1\) 次多项式 \(f(x)\) 和 \(k\)，求出一个多项式 \(g(x)\) 使得 \(g(x)\equiv f(x)^k\pmod{x^n}\)，系数对 \(998244353\) 取模。\(1\leq n\leq 10^5\)，\(0\leq k\leq 10^{10^5}\)。非加强版保证 \(a_0=1\)，加强版则不做保证。
看到幂函数，一个思路就是两边同时取 \(\ln\) 再取 \(\exp\)，得到
不过 \(k\) 很大，如何处理？考虑右侧式子的麦克劳林级数展开：
因此我们可以令 \(k\leftarrow k\bmod{998244353}\)，而不影响答案正确性。
我们需要求 \(\ln(f(x))\)，而这要求 \(a_0=1\)，因此上述做法只能通过非加强版。
非加强版的代码：
void poly_pow(int *a, int *res, int n, int b) {
static int lna[N];
poly_ln(a, lna, n);
for (int i = 0; i < n; ++i) lna[i] = 1ll * lna[i] * b % MOD;
poly_exp(lna, res, n);
clr(lna, n);
}
cin >> n >> str + 1;
for (int i = 1; str[i]; ++i) k = (1ll * k * 10 + (str[i] ^ 48)) % MOD;
for (int i = 0; i < n; ++i) cin >> a[i];
poly_pow(a, ans, n, k);
for (int i = 0; i < n; ++i) cout << ans[i] << " \n"[i == n - 1];
既然上述做法要求 \(a_0=1\)，考虑对 \(f(x)\) 的每一项系数除掉 \(a_0\)，即
这就做完了吗？并没有。加强版还不保证 \(a_0\neq 0\)。那考虑给 \(f(x)\) 降次再升次，即找到最小的 \(t\) 使得 \(a_t\neq 0\)，然后转化成
也就是先把系数整体左移 \(t\) 个位置，沿用前面除掉常数项系数的做法，最后再整体右移 \(tk\) 个位置，这题就做完啦！
但还有很多细节。首先计算 \(a_0^k\) 需要利用欧拉定理，即
需要和多项式的幂次做区分，即要存两个 \(k\)，其中 \(k_1=k\bmod{998244353}\)，作为多项式的幂次，\(k_2=k\bmod{998244352}\)，作为计算 \(a_0^k\) 的幂次。
其次，我们需要特判偏移量 \(tk\geq n\) 的情况。由于 \(n\leq 10^5\)，所以直接把 \(k\) 在数位上截断低 \(6\) 位判断即可，即令 \(k_3=k\bmod{10^6}\)，判断此时是否有 \(k_3\geq n\)，然后再去判断是否有 \(tk_1\geq n\)。
代码：
void poly_pow(int *a, int *res, int n, int b1, int b2) {
int t = 0;
static int ta[N], tr[N], lna[N];
for (int i = 0; i < n && !a[i]; ++i, ++t);
if (1ll * t * b1 >= n) { clr(res, n); return; }
int a0 = a[t], iv = qpow(a0, MOD - 2);
for (int i = 0; i < n; ++i) ta[i] = 1ll * a[i + t] * iv % MOD;
poly_ln(ta, lna, n);
for (int i = 0; i < n; ++i) lna[i] = 1ll * lna[i] * b1 % MOD;
poly_exp(lna, tr, n);
t *= b1;
int pw = qpow(a0, b2);
clr(res, t);
for (int i = t; i < n; ++i) res[i] = 1ll * tr[i - t] * pw % MOD;
}
ios::sync_with_stdio(false), cin.tie(nullptr);
cin >> n >> str + 1;
for (int i = 1; str[i]; ++i) {
k1 = (1ll * k1 * 10 + (str[i] ^ 48)) % MOD;
k2 = (1ll * k2 * 10 + (str[i] ^ 48)) % (MOD - 1);
}
for (int i = 1; str[i] && i <= 6; ++i) k3 = k3 * 10 + (str[i] ^ 48);
for (int i = 0; i < n; ++i) cin >> a[i];
if (!a[0] && k3 >= n) {
for (int i = 0; i < n; ++i) cout << 0 << " \n"[i == n - 1];
return 0;
}
poly_pow(a, ans, n, k1, k2);
for (int i = 0; i < n; ++i) cout << ans[i] << " \n"[i == n - 1];
return 0;
多项式多点求值
题意：给定一个 \(n\) 次多项式 \(f(x)\) 和 \(m\) 个整数 \(p_i\)，\(\forall i\in[1,m]\) 求出 \(f(p_i)\)，答案对 \(998244353\) 取模。\(1\leq n,m\leq 6.4\times 10^4\)。
挺有趣的科技啊。下面默认 \(n,m\) 同阶。
方法一：多项式取模
构造多项式 \(g_i(x)=x-p_i\)，设 \(f(x)\bmod{g_i(x)}=r_i\)，可以发现此时恰好有 \(r_i=f(p_i)\)，因为考虑将 \(f(x)\) 对 \(g_i(x)\) 进行多项式除法，得到 \(f(x)=g_i(x)q_i(x)+r_i\)，代入 \(x=p_i\)，此时 \(g_i(x)=0\)，所以 \(r_i=f(p_i)\)。容易将其进一步推广：
令 \(g_{l,r}(x)=\prod_{i=l}^r(x-p_i)\)，\(f(x)\bmod{g_{l,r}(x)}=r_{l,r}(x)\)，则 \(\forall i\in[l,r],f(p_i)=r_{l,r}(p_i)\)。
于是我们得到了一个较为显然的做法。首先容易 \(O(n\log^2{n})\) 分治 NTT 计算出 \(g_{l,r}(x)\)，然后，我们依然分治解决问题，设当前递归到 \([l,r]\)，并得到了多项式 \(f_{l,r}(x)\)，我们只需将 \(f_{l,r}(x)\) 分别对 \(g_{l,mid}(x)\) 和 \(g_{mid+1,r}(x)\) 取模，得到 \(f_{l,mid}(x)\) 和 \(f_{mid+1,r}(x)\)，继续递归求解即可。当递归到边界 \([l,l]\) 时，由上述推论，显然此时 \(f_{l,l}(x)\) 仅有常数项且它就是 \(f(p_l)\) 的值。初始时只需令 \(f_{1,m}(x)=f(x)\bmod{g_{1,m}(x)}\)。取模保证了递归到 \([l,r]\) 时的多项式 \(f_{l,r}(x)\) 是一个 \(r-l\) 次多项式，复杂度得到了保证。由主定理，这部分复杂度还是 \(O(n\log^2{n})\) 的。
由于多项式取模的大常数，通常上述做法需要经过刻意卡常才能通过本题。
方法二：转置原理
咕咕咕。
代码：
void poly_mult(const int *a, const int *b, int *res, int n, int m) {
int len = calc_len(n);
static int ta[N], tb[N];
copy(a, a + n, ta), reverse_copy(b, b + m, tb);
NTT(ta, len), NTT(tb, len);
for (int i = 0; i < len; ++i) ta[i] = 1ll * ta[i] * tb[i] % MOD;
NTT(ta, len, -1);
copy(ta + m - 1, ta + n, res);
clr(ta, len), clr(tb, len);
}
namespace PolyEval {
#define ls(p) (p << 1)
#define rs(p) (p << 1 | 1)
int inv[N];
int len[N << 2], *g[N << 2], *q[N << 2];
void eval_build(int p, int l, int r, const int *a) {
if (l == r) {
len[p] = 1;
g[p] = new int[2], q[p] = new int[2];
g[p][0] = 1, g[p][1] = (-a[l] + MOD) % MOD;
return;
}
int mid = (l + r) >> 1;
eval_build(ls(p), l, mid, a), eval_build(rs(p), mid + 1, r, a);
len[p] = len[ls(p)] + len[rs(p)];
g[p] = new int[len[p] + 1], q[p] = new int[len[p] + 1];
poly_mul(g[ls(p)], g[rs(p)], g[p], len[ls(p)] + 1, len[rs(p)] + 1);
}
void eval_solve(int p, int l, int r, int *ans) {
if (l == r) { ans[l] = q[p][0]; return; }
int mid = (l + r) >> 1;
poly_mult(q[p], g[rs(p)], q[ls(p)], r - l + 1, len[rs(p)] + 1);
eval_solve(ls(p), l, mid, ans);
poly_mult(q[p], g[ls(p)], q[rs(p)], r - l + 1, len[ls(p)] + 1);
eval_solve(rs(p), mid + 1, r, ans);
}
void solve(const int *a, const int *b, int n, int m, int *ans) {
static int tmp[N];
eval_build(1, 1, m, b);
poly_inv(g[1], inv, m + 1), reverse(inv, inv + m + 1);
poly_mul(a, inv, tmp, n, m + 1), copy(tmp + n, tmp + n + m, q[1]);
eval_solve(1, 1, m, ans);
for (int i = 1; i <= m; ++i) ans[i] = (1ll * ans[i] * b[i] % MOD + a[0]) % MOD;
}
#undef ls
#undef rs
}
多项式快速插值
题意：给出 \(n\) 个点 \((x_i,y_i)\)，求出一个 \(n-1\) 次多项式 \(f(x)\) 使得 \(\forall i\in[1,n],f(x_i)\equiv y_i\pmod{998244353}\)，系数对 \(998244353\) 取模。\(1\leq n\leq 10^5\)。
回顾拉格朗日插值公式：
不妨将分母提出，得到
设 \(g(x)=\prod_{i=1}^n(x-x_i)\)，那么
而函数 \(\frac{g(x)}{x-x_i}\) 是连续可导的，考虑取极限并运用洛必达法则：
于是插值公式化为
这是一个比较典的可以分治计算的形式，具体来说，令
即 \(h_{l,r}(x)\) 表示 \((x_l,y_l),\cdots,(x_r,y_r)\) 插值得到的多项式。我们考虑缺口位于左半区间还是右半区间，推一下式子：
把 \(g'(x)\) 对 \(x_i(i\in[1,n])\) 进行多点求值，然后分治 \(O(n\log^2{n})\) 解决即可。
代码：
namespace PolyInter {
#define ls(p) (p << 1)
#define rs(p) (p << 1 | 1)
int val[N];
int len[N << 2], *g[N << 2], *h[N << 2];
void inter_build(int p, int l, int r, const int *a) {
if (l == r) {
len[p] = 1;
g[p] = new int[2], h[p] = new int[2];
g[p][0] = -a[l] + MOD, g[p][1] = 1;
return;
}
int mid = (l + r) >> 1;
inter_build(ls(p), l, mid, a), inter_build(rs(p), mid + 1, r, a);
len[p] = len[ls(p)] + len[rs(p)];
g[p] = new int[len[p] + 1], h[p] = new int[len[p] + 1];
poly_mul(g[ls(p)], g[rs(p)], g[p], len[ls(p)] + 1, len[rs(p)] + 1);
}
void inter_solve(int p, int l, int r, const int *y) {
if (l == r) { h[p][0] = 1ll * y[l] * qpow(val[l], MOD - 2) % MOD; return; }
int mid = (l + r) >> 1;
inter_solve(ls(p), l, mid, y), inter_solve(rs(p), mid + 1, r, y);
static int tmp[N];
poly_mul(h[ls(p)], g[rs(p)], tmp, len[ls(p)], len[rs(p)] + 1);
poly_mul(h[rs(p)], g[ls(p)], h[p], len[rs(p)], len[ls(p)] + 1);
for (int i = 0; i < len[p]; ++i) h[p][i] = (h[p][i] + tmp[i]) % MOD;
clr(tmp, len[p]);
}
void solve(const int *x, const int *y, int n) {
static int df[N];
inter_build(1, 1, n, x);
poly_diff(g[1], df, n + 1);
PolyEval::solve(df, x, n + 1, n, val);
inter_solve(1, 1, n, y);
clr(df, n);
}
#undef ls
#undef rs
}
例题
Luogu P5395 第二类斯特林数·行
题意：给定 \(n\)，对于所有整数 \(i\in[0,n]\)，求出 \({n\brace i}\)，答案对 \(167772161=5\times 2^{25}+1\) 取模。\(1\leq n\leq 2\times 10^5\)。
这里不加推导地给出第二类斯特林数的通项公式，其由二项式反演容易得到：
遇到类似 \(i\) 和 \(m-i\) 这类和为定值的形式，考虑构造成加法卷积的形式。对这题而言，
令 \(f_i=\frac{i^n}{i!},g_i=\frac{(-1)^i}{i!}\)，则 \({n\brace m}=\sum_{i=0}^mf_ig_{m-i}\)。这就是很显然的卷积形式了，构造多项式 \(f(x)\) 满足 \([x^i]f(x)=f_i\)，和 \(g(x)\) 满足 \([x^i]g(x)=g_i\)，NTT 计算它们的卷积即可。时间复杂度 \(O(n\log{n})\)。
AtCoder ABC196F Substring 2
题意：给定 \(0/1\) 串 \(s,t\)，求最少修改多少次 \(t\) 可以使其成为 \(s\) 的子串。\(1\leq |t|\leq|s|\leq 10^6\)。
令 \(n=|s|,m=|t|\)，枚举把 \(t\) 改成子串 \(s[k,k+m-1]\)，此时的修改次数为
操作数是 \(0/1\)，此时我们有经典结论 \(x\oplus y=(x-y)^2=x-2xy+y\)，代回上面的式子
\(\sum_{i=0}^{m-1}t_i\) 显然是一个定值，\(\sum_{i=0}^{m-1}s_{k+i}\) 可以前缀和处理，于是重点在于计算 \(\sum_{i=0}^{m-1}t_is_{k+i}\)。式子中出现了乘法，但不是我们所熟知的和为定值的形式，而是差为定值，对于这种形式，有经典套路：将其中一个多项式反转。对这题而言，令 \(s'_i=s_{n-1-i}\)，则 \(\sum_{i=0}^{m-1}t_is_{k+i}=\sum_{i=0}^{m-1}t_is_{n-1-k-i}'\)，此时和为定值 \(n-1-k\)，就可以 NTT 直接做了。
浙公网安备 33010602011771号