多项式半家桶速成

叫半家桶是因为多项式太多了，很多板子都不会，但又不是完全不会（雾

主要以 vector 写多项式模板，之前使用数组写多项式太痛苦了。

令 \([x^n]F(x)\) 和 \(F[x]\) 都表示 \(F(x)\) 的第 \(n\) 次项系数。

后文会附代码，其中 Z 是一个自动取模的数（结构体），同时 deg(a) 指的是 \(a\) 的度数\(+1\)，因为我比较习惯于让多项式 \(\bmod~x^{deg(a)}\)。

using poly = vector < Z >;
#define deg(a) int(a.size())

基本运算

加

\[F(x) +G(x) = \sum_{i = 0} (F[i] + G[i])x^i \]

poly operator + (const poly &a, const poly &b) {
	poly c = a; c.resize(max(deg(a), deg(b)));
	rep(i, 0, deg(b) - 1) c[i] += b[i];
	return c;
}

减

\[F(x) - G(x) = \sum_{i = 0} (F[i] - G[i])x^i \]

poly operator - (const poly &a, const poly &b) {
	poly c = a; c.resize(max(deg(a), deg(b)));
	rep(i, 0, deg(b) - 1) c[i] -= b[i];
	return c;
}

乘

\[F(x) G(x) = \sum_{i = 0} x^i \sum_{j = 0}^i F[j] G[i - j] \]

需要 \(\rm NTT\) 做到 \(\mathcal O(n\log n)\)。

void init_NTT(int l) {
	for(len = 1; len <= l + 2; len <<= 1);
	rep(i, 1, len - 1) rev[i] = rev[i >> 1] >> 1 | (i & 1 ? len >> 1 : 0);
	rep(i, 0, len - 1) A[i] = B[i] = 0;
}

void NTT(Z *f, int fl = 1) {
	rep(i, 1, len - 1) if(i < rev[i]) swap(f[i], f[rev[i]]);
	for(int h = 2; h <= len; h <<= 1) {
		Z wn = qp(Z(G), (mod - 1) / h);
		for(int i = 0, j = h >> 1; i < len; i += h) {
			Z ww = 1;
			for(int k = i; k < i + j; k++) {
				Z u = f[k], v = f[k + j] * ww; ww *= wn;
				f[k] = u + v; f[k + j] = u - v;
			}
		}
	} if(fl == -1) {
		reverse(f + 1, f + len); Z t = len; t = t.inv();
		rep(i, 0, len - 1) f[i] *= t;
	}
}
poly operator * (const poly &a, const poly &b) {
	init_NTT(max(deg(a), deg(b)) * 2);
	rep(i, 0, deg(a) - 1) A[i] = a[i]; rep(i, 0, deg(b) - 1) B[i] = b[i];
	NTT(A), NTT(B); rep(i, 0, len - 1) A[i] = A[i] * B[i]; NTT(A, -1);
	poly c; rep(i, 0, len - 1) c.eb(A[i]); c.resize(deg(a) + deg(b) - 1); return c;
}

除

就是求多项式乘法逆。

递推

因为 \(F(x)G(x) \equiv 1 \pmod{x^n}\)，于是有 \(\sum_{i = 0}^{n} F[i]G[n - i] = [n = 0]\)，对于 \(n = 0\) 的时候显然就是直接求逆，对于 \(n \neq 0\) 的时候 \(F[n]G[0] + \sum_{i = 0}^{n - 1} F[i]G[n - i] = 0\)，于是有 \(F[n] = -\frac{1}{G[0]} \sum_{i = 0}^{n - 1} F[i]G[n - i]\)。然后可以直接 \(\mathcal O(n^2)\) 递推求出，但是速度太慢了。

倍增

考虑我们已经求出了 \(H(x)\) 满足 \(F(x) H(x) \equiv 1 \pmod {x^{n / 2}}\)，这里的 \(n / 2\) 是向上取整。

因为 \(F(x) G(x) \equiv 1 \pmod {x^n}\)，一定能够满足 \(F(x) G(x) \equiv 1 \pmod {x^{n / 2}}\)，于是有：

\[G(x) - H(x) \equiv 0 \pmod {x^{n / 2}} \]

考虑平方，然后乘上 \(F(x)\)，于是有：

\[G(x)^2 + H(x)^2 - 2G(x)H(x) \equiv 0 \pmod {x^n} \\ G(x) + H(x)^2F(x) - 2H(x) \equiv 0 \pmod {x^n}\\ G(x) \equiv H(x) (2 -H(x)F(x)) \pmod {x^n} \\ \]

然后可以倍增求出 \(G(x)\) 了。

---

后面会写牛顿迭代的解法。

之前写的是数组版本，于是没有什么问题，直接用第二种即可。

如果使用 vector，可以考虑在项数为奇数的时候递推，项数为偶数的时候倍增。（但是经过实测，开了 \(\rm O2\) 后感觉速度只是略快，于是这个只是卡常的时候再考虑，另外，如果真的要卡常，可以考虑将 \(2\) 次多项式乘法变为 \(1\) 次，NTT 后的点值可以直接操作）

void getinv(const poly &a, poly &b) {
	if(deg(a) == 1) b.eb(qp(a.front()));
	else {
		poly t = a; t.resize(deg(a) + 1 >> 1); getinv(t, b);
		b = b * 2 - b * b * a; b.resize(deg(a));
	}
}

poly inv(const poly &a) { poly b; getinv(a, b); return b; }

poly operator / (const poly &a, const poly &b) { return a * inv(b); }

求导

\[F'(x) = \sum_{i = 1} F[i]x^{i - 1} \]

poly der(const poly &a) { poly b; rep(i, 1, deg(a) - 1) b.eb(a[i] * i); return b; }

积分

\[\int F(x) = C + \sum_{i = 0} \frac{F[i]}{i+ 1} x^{i + 1} \]

poly inte(const poly &a) { poly b; b.eb(0); rep(i, 0, deg(a) - 1) b.eb(a[i] / (i + 1)); return b; }

多项式牛顿迭代

已知 \(G(x)\) ，并且满足 \(G(F(x)) \equiv 0 \pmod {x^n}\)，求 \(F(x) \pmod {x^n}\)。

当 \(n =1\) 的时候，需要特殊处理。

然后考虑我们已经求出来了在 \(\bmod~x^{n / 2}\) 意义下的解 \(H(x)\)，考虑将 \(F(x)\) 作为主元，将 \(G(F(x))\) 在 \(H(x)\) 处泰勒展开：

\[G(F(x)) = \sum_{i = 0} \frac{G^{(i)}(F(x))}{ i!} (F(x) - H(x))^i \]

因为 \(F(x) - H(x)\) 的最低次至少是 \(n / 2\)，于是当 \(i \ge 2\) 时候，都是 \(0\)。

于是有：

\[G(F(x)) = G(H(x)) + G'(H(x))(F(x) - H(x))\\ F(x) = H(x) - \frac{G(H(x))}{G'(H(x))} \]

只要记住上面的公式，就可以快速推一些多项式公式，比较省事。

之后的例子只是推导，具体实现以及细节会在后文补充。

求逆

我们假设给定的多项式是 \(T(x)\)，然后 \(F(x)T(x) - 1 \equiv 0 \pmod {x^n}\)，于是有 \(G(F(x)) = F(x)T(x) - 1 \equiv 0 \pmod {x^n}\)，因为我们实际上是以 \(F(x)\) 作的主元，于是有 \(G'(F(x)) = T(x)\)，假设我们求出在 \(\bmod ~ x^{n / 2}\) 意义下的解为 \(H(x)\)，那么有：

\[F(x) = H(x) - \frac{G(H(x))}{G'(H(x))} = H(x) - \frac{H(x) T(x) - 1}{T(x)} \]

注意到 \(\frac{1}{T(x)}\) 的精度只要到 \(x^{n / 2}\) 即可，于是让 \(\frac{1}{T(x)}\) 为 \(H(x)\) 即可，于是有：

\[F(x) = 2H(x) - H(x) ^2T(x) \]

开根

令给定多项式为 \(T(x)\)，然后 \(F(x)^2 - T(x) \equiv 0 \pmod { x^ n}\)，于是有 \(G(F(x)) = F(x)^ 2- T(x) \pmod {x^n}\)，有 \(G'(F(x)) = 2F(x)\)，假设已经求出在 \(\bmod~x^{n / 2}\) 意义下的解为 \(H(x)\)，有：

\[F(x) = H(x) - \frac{G(H(x))}{G'(H(x))} = H(x) - \frac{H(x)^2 - T(x)}{2H(x)} = \frac{H(x)}{2} + \frac{T(x)}{2H(x)} \]

Exp

令给定多项式为 \(T(x)\)，然后 \(\ln F(x) - T(x) \equiv 0 \pmod { x^ n}\)，于是有 \(G(F(x)) = \ln F(x)- T(x) \pmod {x^n}\)，有 \(G'(F(x)) = \frac{1}{F(x)}\)，假设已经求出在 \(\bmod~x^{n / 2}\) 意义下的解为 \(H(x)\)，有：

\[F(x) = H(x) - \frac{G(H(x))}{G'(H(x))} = H(x) - \frac{\ln H(x) - T(x)}{\frac{1}{H(x)}} = H(x)(1 - \ln H(x) + T(x)) \]

进阶操作

开根

假设我们要求的是 \(T(x)\) 的开根后的结果为 \(F(x)\)，假设我们已经求出了在 \(\bmod ~ n / 2\) 意义下的解 \(H(x)\)，有：

\[H^2(x) - T(x) \equiv 0 \pmod {x^{n / 2}}\\ (H^2(x) - T(x))^2 \equiv 0 \pmod {x^n}\\ (H^2(x) + T(x))^2 \equiv 4H(x)^2T(x) \pmod {x^n}\\ (H^2(x) + T(x))^2 \equiv 4H(x)^2T(x) \pmod {x^n}\\ (\frac{H^2(x) + T(x)}{2H(x)})^2 \equiv T(x) \pmod {x^n}\\ \frac{H^2(x) + T(x)}{2H(x)} \equiv F(x) \pmod {x^n}\\ \]

然后可以倍增推出。要注意的是，使用 vector 版的开根的 \(H(x)\) 的逆必须要在 \(\bmod~{x^n}\) 意义下。

---

嫌麻烦的话可以直接用多项式牛顿迭代。

void getsqrt(const poly &a, poly &b) {
	if(deg(a) == 1) assert(a.front().val() == 1), b.eb(1);
	else {
		poly t = a; t.resize(deg(a) + 1 >> 1); getsqrt(t, b); b.resize(deg(a));
		b = (b + a / b) * Z(2).inv(); b.resize(deg(a));
	}
}

poly sqrt(const poly &a) { poly b; getsqrt(a, b); return b; }

ln

定义

由泰勒展开 \(f(x) = \sum_{i = 0} \frac{f^{(i)}(x_0)}{i!} (x - x_0)^ i\)，我们将 \(\ln F(x)\) 在 \(x = 1\) 处展开（把 \(F(x)\) 当做主元）。

首先 \(\ln^{(i)}( F(x))\) 为 \(-(-1)^i(i - 1)! F(x)^{-i}\)，然后有：

\[\ln F(x) =- \sum_{i = 1} \frac{(-1)^i (F(x) - 1)^i}{i} = - \sum_{i = 1} \frac{(1 - F(x))^i}{i} \]

求法

首先，只有 \(F[0] = 1\) 的才可能有对数函数。

考虑对于 \(\ln F(x)\) 先求导再积分，于是有：

\[(\ln F(x))' \equiv \frac{F'(x)}{F(x)} \pmod {x^n}\\ \int (\ln F(x))' \equiv F(x) \equiv \int \frac{F'(x)}{F(x)} \pmod {x^n} \]

只要求逆，然后乘法，然后积分即可。

poly ln(const poly &a) { assert(a.front().val() == 1); return inte(der(a) / a); }

Exp

定义

\[\exp F(x) = \sum_{i = 0} \frac{F(x)^i}{i!} \]

求法

如果 \(\exp F(x)\) 存在，那么必须 \(F[0] = 0\)，否则常数项不收敛。

鉴于分治 NTT 的 \(\exp F(x)\) 是 \(\log^2\) 的，于是这里不做推导。其实是不想写，虽然分治 NTT 的常数比较小，但是常数不重要，好写才是王道，于是就不写了。

void getexp(const poly &a, poly &b) {
	if(deg(a) == 1) assert(a.front().val() == 0), b.eb(1);
	else {
		poly t = a; t.resize(deg(a) + 1 >> 1); getexp(t, b); b.resize(deg(a));
		b = b + b * (a - ln(b)); b.resize(deg(a));
	}
}

poly exp(const poly &a) { poly b; getexp(a, b); return b; }

多项式快速幂

也就是计算 \(F(x)^k\) 的 \(\mathcal O(n\log n)\) 做法：

• \(F[0] = 1\) 的时候，\(F(x)^k = \exp(k\ln F(x))\)。
• \(F[0] \neq 1\) 的时候，我们设最低项为第 \(d\) 项，那么 \(F(x)^k = F[d]^k x^{dk} \exp(k \ln \frac{F(x)}{F[d]x^d})\)。

多项式带余除法

给定一个 \(n\) 次多项式 \(F(x)\) 和一个 \(m\) 次多项式 \(G(x)\)，求 \(Q(x), R(x)\)，满足：

• \(Q(x)\) 次数为 \(n -m\)，\(R(x)\) 次数 \(<m\)。
• \(F(x) = Q(x)G(x)+R(x)\)。

只要我们能不考虑 \(R(x)\) 的存在，就可以直接求逆了，于是接下来考虑消除 \(R(x)\)。

定义操作 \(F^T(x)\) 表示将 \(n\) 次多项式 \(F(x)\) 函数系数翻转，\(F^T(x) = x^nF(x^{-1})\)，然后可以将 \(F(x) = Q(x)G(x)+R(x)\) 进行一次操作：

\[x^n F(x) = x^n Q(x)G(x) + x^n R(x)\\ F^T(x) = Q^T(x)G^T(x) + x^{n - m + 1}R^T(x) \]

对于整个等式 \(\bmod~x^{n - m + 1}\)，于是：

\[F^T(x) \equiv Q^T(x)G^T(x) \pmod {x^{n - m + 1}}\\ Q^T(x) \equiv F^T(x) (G^T(x))^{-1} \pmod {x^{n - m + 1}} \]

然后可以直接通过多项式求逆求出 \(Q^T(x)\) 了，再翻转一下就是 \(Q(x)\) 了。

poly div(poly a, poly b) { // a / b = c
	int t = deg(a) - deg(b) + 1; reverse(a.begin(), a.end()); reverse(b.begin(), b.end());
	a.resize(t); b.resize(t); poly c = a / b; c.resize(t); reverse(c.begin(), c.end()); return c;
}

最后根据式子算一算 \(R(x)\) 即可。

最后的板子

目前的目标只是能够快点打出板子，因此常数不重要！

#include <bits/stdc++.h>

#define eb emplace_back
#define ep emplace
#define fi first
#define se second
#define in read<int>()
#define lin read<ll>()
#define rep(i, x, y) for(int i = (x); i <= (y); i++)
#define per(i, x, y) for(int i = (x); i >= (y); i--)
#define deg(x) int(x.size())

using namespace std;

using ll = long long;
using db = double;
using pii = pair < int, int >;
using vec = vector < int >;
using veg = vector < pii >;

template < typename T > T read() {
	T x = 0; bool f = 0; char ch = getchar();
	while(!isdigit(ch)) f |= ch == '-', ch = getchar();
	while(isdigit(ch)) x = x * 10 + (ch ^ 48), ch = getchar();
	return f ? -x : x;
}

template < typename T > void chkmax(T &x, const T &y) { x = x > y ? x : y; }
template < typename T > void chkmin(T &x, const T &y) { x = x < y ? x : y; }

const int N = 1e6 + 10;
constexpr int mod = 998244353;
constexpr int G = 3;
//constexpr int mod = 1e9 + 7;

int reduce(int x) {
	x += x >> 31 & mod;
	return x;
}

template < typename T > T qp(T x, ll t = mod - 2) { T res = 1; for(; t; t >>= 1, x *= x) if(t & 1) res *= x; return res; }

struct Z { // modint
	int x;
	Z(int x = 0) : x(reduce(x)) {}
	Z(ll x) : x(reduce(x % mod)) {}
	Z operator -() const { return Z(reduce(mod - x)); }
	int val() const { return x; }
	Z inv() const { assert(x); return qp(*this); }
	Z &operator += (const Z &t) { x = reduce(x + t.x - mod); return *this; }
	Z &operator -= (const Z &t) { x = reduce(x - t.x); return *this; }
	Z &operator *= (const Z &t) { x = (ll)x * t.x % mod; return *this; }
	Z &operator /= (const Z &t) { return *this *= t.inv(); }
	friend Z operator + (const Z &a, const Z &b) { Z res = a; res += b; return res; }
	friend Z operator - (const Z &a, const Z &b) { Z res = a; res -= b; return res; }
	friend Z operator * (const Z &a, const Z &b) { Z res = a; res *= b; return res; }
	friend Z operator / (const Z &a, const Z &b) { Z res = a; res /= b; return res; }
};

Z fac[N], ifac[N];
Z C(int x, int y) { return x < 0 || y < 0 || x < y ? Z(0) : fac[x] * ifac[y] * ifac[x - y]; }
void init(int l) {
	fac[0] = 1; rep(i, 1, l) fac[i] = fac[i - 1] * Z(i); ifac[l] = fac[l].inv();
	per(i, l - 1, 0) ifac[i] = ifac[i + 1] * Z(i + 1);
}

using poly = vector < Z >;

int rev[N << 2], len;
Z A[N << 2], B[N << 2];

void init_NTT(int l) {
	for(len = 1; len <= l + 2; len <<= 1);
	rep(i, 1, len - 1) rev[i] = rev[i >> 1] >> 1 | (i & 1 ? len >> 1 : 0);
	rep(i, 0, len - 1) A[i] = B[i] = 0;
}

void NTT(Z *f, int fl = 1) {
	rep(i, 1, len - 1) if(i < rev[i]) swap(f[i], f[rev[i]]);
	for(int h = 2; h <= len; h <<= 1) {
		Z wn = qp(Z(G), (mod - 1) / h);
		for(int i = 0, j = h >> 1; i < len; i += h) {
			Z ww = 1;
			for(int k = i; k < i + j; k++) {
				Z u = f[k], v = f[k + j] * ww; ww *= wn;
				f[k] = u + v; f[k + j] = u - v;
			}
		}
	} if(fl == -1) {
		reverse(f + 1, f + len); Z t = len; t = t.inv();
		rep(i, 0, len - 1) f[i] *= t;
	}
}

poly operator + (const poly &a, const poly &b) {
	poly c = a; c.resize(max(deg(a), deg(b)));
	rep(i, 0, deg(b) - 1) c[i] += b[i];
	return c;
}

poly operator - (const poly &a, const poly &b) {
	poly c = a; c.resize(max(deg(a), deg(b)));
	rep(i, 0, deg(b) - 1) c[i] -= b[i];
	return c;
}

poly operator * (const poly &a, Z t) { poly c; for(const auto &v : a) c.eb(v * t); return c; }

poly operator * (const poly &a, const poly &b) {
	init_NTT(max(deg(a), deg(b)) * 2);
	rep(i, 0, deg(a) - 1) A[i] = a[i]; rep(i, 0, deg(b) - 1) B[i] = b[i];
	NTT(A), NTT(B); rep(i, 0, len - 1) A[i] = A[i] * B[i]; NTT(A, -1);
	poly c; rep(i, 0, len - 1) c.eb(A[i]); c.resize(deg(a) + deg(b) - 1); return c;
}

void getinv(const poly &a, poly &b) {
	if(deg(a) == 1) b.eb(qp(a.front()));
	else {
		poly t = a; t.resize(deg(a) + 1 >> 1); getinv(t, b);
		b = b * 2 - b * b * a; b.resize(deg(a));
	}
}

poly inv(const poly &a) { poly b; getinv(a, b); return b; }

poly operator / (const poly &a, const poly &b) { return a * inv(b); }

void getsqrt(const poly &a, poly &b) {
	if(deg(a) == 1) assert(a.front().val() == 1), b.eb(1);
	else {
		poly t = a; t.resize(deg(a) + 1 >> 1); getsqrt(t, b); b.resize(deg(a));
		b = (b + a / b) * Z(2).inv(); b.resize(deg(a));
	}
}

poly sqrt(const poly &a) { poly b; getsqrt(a, b); return b; }

poly der(const poly &a) { poly b; rep(i, 1, deg(a) - 1) b.eb(a[i] * i); return b; }

poly inte(const poly &a) { poly b; b.eb(0); rep(i, 0, deg(a) - 1) b.eb(a[i] / (i + 1)); return b; }

poly ln(const poly &a) { assert(a.front().val() == 1); return inte(der(a) / a); }

int n, m;

void getexp(const poly &a, poly &b) {
	if(deg(a) == 1) assert(a.front().val() == 0), b.eb(1);
	else {
		poly t = a; t.resize(deg(a) + 1 >> 1); getexp(t, b); b.resize(deg(a));
		b = b + b * (a - ln(b)); b.resize(deg(a));
	}
}

poly exp(const poly &a) { poly b; getexp(a, b); return b; }

poly div(poly a, poly b) { // a / b = c
	int t = deg(a) - deg(b) + 1; reverse(a.begin(), a.end()); reverse(b.begin(), b.end());
	a.resize(t); b.resize(t); poly c = a / b; c.resize(t); reverse(c.begin(), c.end()); return c;
}

int main() {
#ifndef ONLINE_JUDGE
	freopen("1.in", "r", stdin);
#endif

	return 0;
}