Templates of Polynomial
最近在写一些多项式的模板题,代码总是显得十分混乱。
深感代码通用的重要性,因此开一篇小结用以记录比较通用的模板。
如果模板出了问题请联系我,万分感谢。
写在前面
为了体现简洁,在每一部分只会放关键代码。
关于完整代码,在每一部分的代码标题处都放了链接。
可以保证代码之间函数的调用是合法的,会给出参数说明。
因为放代码可能会比较长,可以点右边的小火箭回到目录。
如果 MathJax 加载不出来或加载有误,请您多刷新几次。
模板出了问题请联系我,万分感谢。
参考资料
Picks 的博客 Picks’s Blog
Miskcoo 的博客 Miskcoo’s Space
多项式乘法
问题描述
给定两个多项式 $A(x) $ 和 $B(x)$ :
求卷积 $C(x) =A(x) * B(x)$ ,满足
解决方法
这部分上面提到的神仙们讲解的都十分详细。
我的小结限于篇幅,只挂一下 pdf 版本的连接了:快速傅里叶变换初步
代码实现
使用 [ UOJ 34 ] 多项式乘法 作为测试题。
Fast Fourier Transform,使用的是 3 次变换的最基本写法,用时 362 ms
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15inline void FFT(Complex *f, int len, int o) {
for (int i = 0; i < len; ++i)
if (rev[i] > i) swap(f[i], f[rev[i]]);
for (int i = 1; i < len; i <<= 1) {
Complex wn = Complex(cos(PI / i) , o * sin(PI / i));
for (int j = 0; j < len; j += (i << 1)) {
Complex w = Complex(1, 0), x, y;
for (int k = 0; k < i; ++k, w = w * wn) {
x = f[j + k]; y = w * f[i + j + k];
f[i + j + k] = x - y; f[j + k] = x + y;
}
}
}
if (o == -1) for (int i = 0; i < len; ++i) f[i].x /= len;
}Fast Number-Theoretic Transform,模数为 998244353,用时 403 ms
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19inline void NTT(int *f, int len, int o) {
for (int i = 1; i < len; ++i)
if (i > rev[i]) swap(f[i], f[rev[i]]);
for (int i = 1; i < len; i <<= 1) {
int wn = qpow(3, (mod - 1) / (i << 1));
if (o == -1) wn = qpow(wn, mod - 2);
for (int j = 0; j < len; j += (i << 1)) {
int w = 1, x, y;
for (int k = 0; k < i; ++k, w = 1ll * w * wn % mod) {
x = f[j + k]; y = 1ll * w * f[i + j + k] % mod;
f[j + k] = mo(x + y); f[i + j + k] = mo(x - y + mod);
}
}
}
if (o == -1) {
int invl = qpow(len, mod - 2);
for (int i = 0; i < len; ++i) f[i] = 1ll * f[i] * invl % mod;
}
}
多项式求逆
问题描述
给定一个 $n$ 次多项式 $A(x)$ ,求出一个多项式 $B(x)$, 满足
系数对 998244353 取模。
解决方法
采用倍增的思想。
考虑只有常数项的时候,$A(x)\equiv c\pmod{x}$ ,那么 $A^{-1}(x)$ 即为 $c^{-1}$。
对于 $n>1$ 的时候,设 $B(x)=A^{-1}(x)$ ,有
因为此时模 $x^n$ 相当于只保留多项式前 $n$ 项,所以该同余式在模 $x^k,0\le k\le n$ 时都成立
假设在 $\pmod {x^{\lceil \frac n2\rceil}}$ 意义下 $A(x)$ 的逆元是 $B′(x)$ 并且我们已经求出,那么
两式相减,得
两边平方,得
模数平方的合法性在于,大于 $n$ 的系数卷积中每组乘法必然有一项下标小于 $n$ ,因此乘起来必然为 $0$ 。
两侧同称 $A(x)$ ,整理得
因此只需将 $B’(x)$ 和 $A(x)$ 在模 $x^n$ 意义下的插值求出,有
因此一遍 NTT 就可以由 $B’(x)$ 求出 $B(x)$ 了。
总的时间复杂度为
由此过程也可以得到一个结论:一个多项式有没有逆元完全取决于其常数项是否有逆元。
代码实现
使用 [ Luogu P4238 ] 多项式求逆 作为测试题。
递归版本,使用 O2 优化,用时 562 ms
1
2
3
4
5
6
7
8
9
10
11
12
13inline void Inv(int *a, int *b, int n) {
if (n == 1) {b[0] = qpow(a[0], mod - 2); return;}
Inv(a, b, (n + 1) >> 1);
int len = Rev(n << 1);
for (int i = 0; i < n; ++i) tmp[i] = a[i];
for (int i = n; i < len; ++i) b[i] = tmp[i] = 0;
NTT(b, len, 1); NTT(tmp, len, 1);
for (int i = 0; i < len; ++i)
b[i] = (2ll - 1ll * tmp[i] * b[i] % mod + mod) * b[i] % mod;
NTT(b, len, -1);
for (int i = 0; i < len; ++i) tmp[i] = 0;
for (int i = n; i < len; ++i) b[i] = 0;
}迭代版本,使用 O2 优化,用时 570 ms
1
2
3
4
5
6
7
8
9
10
11
12
13inline void Inv(int *a, int *b, int n) {
b[0] = qpow(a[0], mod - 2);
int len;
for (int l = 1; l < (n << 1); l <<= 1) {
len = Rev(l << 1);
for (int i = 0; i < l; ++i) tmp[i] = a[i];
NTT(b, len, 1); NTT(tmp, len, 1);
for (int i = 0; i < len; ++i)
b[i] = (2ll - 1ll * tmp[i] * b[i] % mod + mod) * b[i] % mod;
NTT(b, len, -1);
for (int i = l; i < len; ++i) b[i] = 0;
}
}
多项式开根
问题描述
给定一个 $n$ 次多项式 $A(x)$,求一个在 $\bmod{x^{n+1}}$ 意义下的多项式 $B(x)$,使得
多项式的系数在模 998244353 意义下进行运算,保证常数项 $a_0=1$。
解决方法
同样采用倍增的思想。
考虑只有常数项的时候,$A(x)\equiv c\pmod x$ ,那么 $\sqrt {A(x)}$ 即为 $\sqrt c \equiv 1\pmod{x}$ (二次剩余)。
对于 $n>1$ 的时候,同样根据上一题的结论,我们可以把问题范围缩小到 $\bmod {x^{\lceil \frac{n}2\rceil}}$ ,有
不妨设我们已经求出来了 $\bmod{x^{\lceil \frac{n}2\rceil}} $ 意义下的根 $D(x)$,即
因此 $B(x)$ 与 $D(x)$ 在模 $x^{\lceil \frac{n}2\rceil}$ 意义下同余,移项得
两侧平方,得
模数能平方的原因与上一题相同。
我们知道$\bmod{x^n}$ 时 $B^2(x)$ 即为 $A(x)$ ,因此
移项,得
因此倍增时进行多项式求逆即可,总的时间复杂度为
代码实现
使用 [ Luogu P5205 ] 多项式开根 作为测试题,多项式求逆部分均采用递归版本。
递归版本,使用 O2 优化,用时 3081 ms
1
2
3
4
5
6
7
8
9
10
11
12
13inline void Sqrt(int *a, int *b, int n) {
if (n == 1) {b[0] = 1; return;}
Sqrt(a, b, (n + 1) >> 1);
Inv(b, b0, n);
int len = Rev(n << 1);
for (int i = 0; i < n; ++i) a0[i] = a[i];
for (int i = n; i < len; ++i) a0[i] = 0;
NTT(a0, len, 1); NTT(b0, len, 1);
for (int i = 0; i < len; ++i) a0[i] = 1ll * a0[i] * b0[i] % mod;
NTT(a0, len, -1);
for (int i = 0; i < n; ++i) b[i] = 1ll * (b[i] + a0[i]) % mod * inv2 % mod;
for (int i = n; i < len; ++i) b[i] = 0;
}迭代版本,使用 O2 优化,用时 3104 ms
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15inline void Sqrt(int *a, int *b, int n) {
b[0] = 1;
int len;
for (int l = 1; l < (n << 1); l <<= 1) {
Inv(b, b0, l);
len = Rev(l << 1);
for (int i = 0; i < l; ++i) a0[i] = a[i];
for (int i = l; i < len; ++i) a0[i] = 0;
NTT(a0, len, 1); NTT(b0, len, 1);
for (int i = 0; i < len; ++i) a0[i] = 1ll * a0[i] * b0[i] % mod;
NTT(a0, len, -1);
for (int i = 0; i < l; ++i) b[i] = 1ll * (b[i] + a0[i]) % mod * inv2 % mod;
for (int i = l; i < len; ++i) b[i] = 0;
}
}
多项式除法和取模
问题描述
给定一个 $n$ 次多项式 $A(x)$ 和一个 $m$ 次多项式 $B(x)$,求出多项式 $D(x)$, $R(x)$,满足
$D(x)$ 次数为 $n-m$,$R(x)$ 次数小于 $m$ ,所有的运算在模 998244353 意义下进行。
解决方法
注意到带着 $R(x)$ 在这里很麻烦,前人们想到了一个神奇的解决办法。
设 $A^R(x)=x^nA(\frac{1}{x})$ ,我们将右侧展开:
所以 $A^R(x)$ 就是将 $A(x)$ 的系数反转 。
我们将所求的等式中 $x$ 全部换成 $\frac{1}{x}$ ,然后两侧同乘 $x^n$ :
注意到 $D^R(x)$ 最高次反转后不变,依然为 $n-m$。
而右侧的 $R^R(x)$ 因为前面有 $x^{n-m+1}$ 所以最低次为 $n-m+1$ 。
所以我们可以把多项式运算在 $\bmod{x^{n-m+1}}$ 意义下进行,这样 $R(x)$ 就消失了:
因此就可以得到 $D^R(x)$ 的解法:
一个多项式求逆就可以求出 $D^R(x)$了,再将 $D^R(x)$ 进行反转就得到了答案。
将求出的 $D(x)$ 回代,再进行一次减法即可求出 $R(x)$。
复杂度与多项式求逆同阶,为 $\mathcal O(n \log n)$ 。
代码实现
使用 [ Luogu P4512 ] 多项式除法 作为测试题,不再区分多项式求逆部分的实现方式。
多项式求逆递部分归版本,使用 O2 优化,用时 710 ms
注意求 $D(x)$ 部分的那次卷积是在 $\bmod{x^{n-m+1}}$意义下的,所以 $A^R(x)$ 和 $B^R(x)$ 中高于 $n-m+1$ 的项需要清空(因为 FFT 卷积过程中高次系数也会对低次系数造成影响)
1 | inline void Div(int *a, int *b, int n, int m) { |
分治 FFT
问题描述
给定长度为 $n-1$ 的数组 $g[1],\dots,g[n-1]$,求长度为 $n$ 的数组 $f[0],\dots,f[n-1]$,其中
边界为 $f[0]=1$,运算在模 998244353 下进行。
解决方法
分治求解
暴力做是 $\mathcal O(n^2)$ 的,考虑将相同的转移一起做以达到优化的目的。
考虑使用类似 CDQ 分治的思想,每次我们求出 $[L, mid]$ 范围内的 $f$ 数组之后,把这部分 $f$ 对 $[mid+1, R]$ 范围内 $f$ 的贡献一起做。
考虑对 $x\in[mid + 1, R]$ 的 $f[x]$ 的贡献 $w_x$,有
因此 $w$ 数组可以卷积求了,注意求 $w_x$ 时后半段的 $f$ 需要认为是 $0$,否则就存在右区间内部的贡献了。
总的时间复杂度为
多项式求逆
一阶分治 FFT 是可以看作卷积处理的。
不妨设将数组看成多项式,有
将两个多项式卷积,有
后一个等式成立的原因是,注意到后一个求和就是 $f[i]$ 的形式,所以只有 $f[0]$ 没有被计数
求 $f$ 数组可以看作是 $\bmod{x^n}$ 意义下进行的,因此有
于是一遍多项式求逆就可以求出来了,复杂度为 $\mathcal O(n\log n)$
代码实现
使用 [ Luogu P4721 ] 分治 FFT 作为测试题,多项式求逆采用递归版本。
分治版本,使用 O2 优化,用时 974 ms
多项式求逆版本,使用 O2 优化,用时 261 ms
1
2
3
4
5
6inline void solve(int *a, int n) {
a[0] = 1;
for (int i = 1; i < n; ++i) a[i] = mo(mod - a[i]);
Inv(a, b, n);
for (int i = 0; i < n; ++i) print(b[i], 0);
}
多项式求导和积分
问题描述
给定一个 $n$ 次多项式 $A(x)$ ,求一个 $n-1$ 次多项式 $B(x)$ ,和一个 $n+1$ 次多项式 $C(x)$,满足
解决方法
直接按照定义做即可,有
我们一般认为 $C(x)$ 的常数项为 $0$ ,复杂度显然为 $\mathcal O(n)$。
代码实现
1 | inline void Der(int *a, int n) { |
多项式 ln
问题描述
给出 $n$ 次多项式 $A(x)$,求一个 $\bmod{x^{n+1}}$ 下的多项式 $B(x)$,满足
所有运算在模 998244353 下进行。
解决方法
设 $F(x)=\ln x$,则 $B(x)=F(A(x))$ 。
对 $B(x)$ 求导,根据链式法则,有
因此对 $A(x)$ 分别进行求导和求逆,卷积即可求出 $B’(x)$ ,再对其进行积分即可。
复杂度与多项式求逆同阶,为 $\mathcal O(n \log n)$ 。
代码实现
使用 [ Luogu P4725 ] 多项式对数函数 作为测试题,不再区分多项式求逆部分的实现方式。
多项式求逆递部分归版本,使用 O2 优化,用时 682 ms
1 | inline void Ln(int *a, int *b, int n) { |