「模板」多项式全家桶
泰勒展开
对于多项式复合函数 \(G(f)\)，若已知其在 \(f_0\) 处的取值 \(G(f_0)\)， 那么其在任意 \(f\) 处的取值 \(G(f)\) 可以表示为
牛顿迭代
可以用于求解多项式方程。
给定多项式复合函数 \(G(x)\)，求多项式 \(F(x)\operatorname{mod} x^N\) 使得 \(G(F(x))\equiv 0(\operatorname{mod} x^N)\)。
考虑倍增法，假设已经求出 \(f_0\equiv f(\operatorname{mod} x^n)\)。在 \(f=f_0\) 处进行 \(\text{Taylor}\) 展开
因为 \(f_0\equiv f(\operatorname{mod} x^n)\)，所以 \((f_0-f)^2\equiv 0(\operatorname{mod} x^{2n})\)，由此将上式化简为
要求解的 \(f\) 满足 \(G(f)\equiv 0(\operatorname{mod} x^{2n})\)，因此
如此迭代即可求得多项式 \(f\)，称之为牛顿迭代法。
下面尝试用牛顿迭代解决一些常规的例子。
多项式求逆
已知 \(A(x)\)，求 \(F(x)\) 满足 \(A(x)F(x)\equiv 1(\operatorname{mod} x^n)\)。
构造 \(G(f)=\frac{1}{f}-A(x)\)，\(G'(f)=-\frac{1}{f^2}\)。
代入牛顿迭代式：
可以 \(O(n\log n)\) 求出。
多项式对数函数
已知 \(F(x)\)，求 \(G(x)\equiv \ln F(x)(\operatorname{mod} x^n)\)。
\(\ln f=\int \operatorname{d}\ln f=\int f'f^{-1} \text{d}x\)。
将 \(f\) 的导数和 \(f\) 的逆卷积，再进行积分即可。
时间复杂度 \(O(n\log n)\)。
多项式指数函数
已知 \(A(x)\)，求 \(F(x)\equiv e^{A(x)}(\operatorname{mod} x^n)\)。
两边同时取 \(\ln\)，\(\ln F(x)\equiv A(x)(\operatorname{mod} x^n)\)。
构造 \(G(f)=\ln f-A(x)\)，\(G'(f)=\frac{1}{f}\)。
时间复杂度 \(O(n\log n)\)，但常数巨大。
多项式开根
已知 \(A(x)\)，求 \(F(x)\equiv A(x)^{\frac{1}{2}}(\operatorname{mod} x^n)\)。
根据 \(F(x)^2\equiv A(x)(\operatorname{mod} x^n)\)，构造 \(G(f)=f^2-A(x)\)，\(G'(f)=2f\)。
时间复杂度 \(O(n\log n)\)。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef vector<int> Poly;
const int N = 1 << 20, mod = 998244353, inv2 = mod + 1 >> 1;
int r[N], Inv[N];
inline int add(int x, int y) { return x + y >= mod ? x + y - mod : x + y; }
inline int dec(int x, int y) { return x - y < 0 ? x - y + mod : x - y; }
inline int mul(int x, int y) { return (ll)x * y % mod; }
inline int qpow(int a, int b) {
int ans = 1;
for (; b; b >>= 1, a = mul(a, a)) if (b & 1) ans = mul(ans, a);
return ans;
}
inline void NTT(Poly &A, int len, int type){
for (int i = 0; i < len; ++ i) if(i < r[i]) swap(A[i], A[r[i]]);
for (int mid = 1; mid < len; mid <<= 1) {
int Wn = qpow(type == 1 ? 3 : (mod + 1) / 3, (mod - 1) / (mid << 1));
for (int j = 0; j < len; j += (mid << 1))
for (int k = 0, w = 1; k < mid; ++ k, w = mul(w, Wn)) {
int x = A[j + k], y = mul(w, A[j + k + mid]);
A[j + k] = add(x, y), A[j + k + mid] = dec(x, y);
}
}
if (type > 0) return;
for (int i = 0; i < len; ++ i) A[i] = mul(A[i], Inv[len]);
}
inline void init_rev(int len) {
for (int i = 0; i < len; ++ i) r[i] = r[i >> 1] >> 1 | ((i & 1) * (len >> 1));
}
inline Poly operator + (const Poly&a, const Poly&b) {
Poly c = a; c.resize(max(a.size(), b.size()));
for (int i = 0; i < b.size(); ++ i) c[i] = add(c[i], b[i]);
return c;
}
inline Poly operator - (const Poly&a, const Poly&b) {
Poly c = a; c.resize(max(a.size(), b.size()));
for (int i = 0; i < b.size(); ++ i) c[i] = dec(c[i], b[i]);
return c;
}
inline Poly operator * (Poly a, Poly b) {
int n = a.size(), m = b.size(), l = 1;
while (l < n + m - 1) l <<= 1;
init_rev(l);
a.resize(l), NTT(a, l, 1);
b.resize(l), NTT(b, l, 1);
for (int i = 0; i < l; ++ i) a[i] = mul(a[i], b[i]);
NTT(a, l, -1);
a.resize(n + m - 1);
return a;
}
inline Poly operator * (Poly a, int b) {
for (int i = 0; i < a.size(); ++ i) a[i] = mul(a[i], b);
return a;
}
inline Poly Deriv(Poly a) {
for (int i = 0; i + 1 < a.size(); ++ i) a[i] = mul(a[i + 1], i + 1);
return a.pop_back(), a;
}
inline Poly Integ(Poly a) {
a.emplace_back(0);
for (int i = a.size() - 1; i; -- i) a[i] = mul(a[i - 1], Inv[i]);
return a[0] = 0, a;
}
inline Poly Polyinv(const Poly&a, int lim) {
Poly c, b(1, qpow(a[0], mod - 2));
for (int l = 4; (l >> 2) < lim; l <<= 1) {
init_rev(l);
c = a, c.resize(l >> 1);
c.resize(l), NTT(c, l, 1);
b.resize(l), NTT(b, l, 1);
for (int i = 0; i < l; ++ i) b[i] = mul(b[i], dec(2, mul(b[i], c[i])));
NTT(b, l, -1), b.resize(l >> 1);
}
return b.resize(lim), b;
}
inline Poly Polyinv(const Poly&a) { return Polyinv(a, a.size()); }
inline Poly Polyln(Poly a, int lim) {
a = Integ(Deriv(a) * Polyinv(a, lim));
return a.resize(lim), a;
}
inline Poly Polyln(const Poly&a) { return Polyln(a, a.size()); }
inline Poly Polyexp(const Poly&a, int lim) {
Poly c, b(1, 1); int n = a.size();
for (int i = 2; (i >> 1) < lim; i <<= 1) {
c = Polyln(b, i);
for (int j = 0; j < i; ++ j) c[j] = dec(j < n ? a[j] : 0, c[j]);
c[0] = add(c[0], 1);
b = b * c, b.resize(i);
}
return b.resize(lim), b;
}
inline Poly Polyexp(const Poly&a) { return Polyexp(a, a.size()); }
inline Poly Polysqrt(const Poly&a, int lim) {
Poly c, d, b(1, 1);
for (int l = 4; (l >> 2) < lim; l <<= 1) {
init_rev(l);
c = a, c.resize(l >> 1);
d = Polyinv(b, l >> 1);
c.resize(l), NTT(c, l, 1);
d.resize(l), NTT(d, l, 1);
for (int j = 0; j < l; ++ j) c[j] = mul(c[j], d[j]);
NTT(c, l, -1), b.resize(l >> 1);
for (int j = 0; j < (l >> 1); ++ j) b[j] = mul(add(b[j], c[j]), inv2);
}
return b.resize(lim), b;
}
inline Poly Polysqrt(const Poly&a) { return Polysqrt(a, a.size()); }
inline Poly Polypow(Poly a, int k, int lim) {
a = Polyexp(Polyln(a) * k);
return a.resize(lim), a;
}
inline Poly Polypow(const Poly&a, int k) { return Polypow(a, k, a.size()); }
inline void precalc() {
Inv[0] = Inv[1] = 1;
for (int i = 2; i < N; ++ i) Inv[i] = mul(mod - mod / i, Inv[mod % i]);
}
int main() {
precalc();
}
