多项式进阶操作

多点求值

问题：给定一个 \(n-1\) 次多项式 \(f(x)\)，求在 \(a_0,a_2,...,a_{m-1}\) 处分别求得的点值。

\(n,m\le 10^5\)

首先我们先钦定 \(n=m\)，否则也可以适当补，下文中用 \(n\) 来代替 \(m\)。

设 \(F=[f_0,f_1,...,f_{n-1}]\)，\(A = \begin{bmatrix} a_0^0 & a_0^1 & ... & a_0^{n - 1} \\ a_1^0 & a_1^1 & ... & a_1^{n - 1} \\ ... & ... & ... & ... \\ a_{n - 1}^0 & a_{n - 1} ^ 1 & ... & a_{n - 1} ^ {n - 1}\end{bmatrix}\)。

那么我们要求的答案为 \(Y_i = \sum_{j = 0} ^ {n - 1} a_i^j f_j\)，即 \(Y = AF\)。

考虑转置原理，我们考虑求 \(Y'=A^T F\)，有 \(Y_i' = \sum_{j = 0} ^ {n - 1} a_j^i f_j\)。

考虑 \(Y_{0..n - 1}'\) 的生成函数：

\[\begin{aligned} \sum_{i = 0} x^i\sum_{j = 0} ^ {n - 1} a_j^i f_j &= \sum_{j=0} ^{n - 1} f_j\sum_{i = 0} x^ia_j^i \\ &= \sum_{j = 0} ^{n - 1} f_j\sum_{i = 0} (a_jx)^i \\ &= \sum_{j = 0} ^{n - 1} f_j \cdot \dfrac 1 {1 - a_jx} \\ \end{aligned}\]

我们很容易把这个式子写成 \(\dfrac {P(x)}{Q(x)}\) 的形式。

具体的，我们分治计算，对于当前区间 \([l,r]\)，合并左右 \(\dfrac{P_0(x)}{Q_0(x)}\) 和 \(\dfrac{P_1(x)}{Q_1(x)}\) 的结果为：

\[\dfrac {P_0(x)Q_1(x) + P_1(x)Q_0(x)}{Q_0(x)Q_1(x)} \]

即 \(P(x) = P_0(x)Q_1(x) + P_1(x)Q_0(x), \ \ Q(x) = Q_0(x)Q_1(x)\)。

对于叶子 \(i\)，\(P(x) = f_i, \ \ Q(x) = 1 - a_ix\)。

接下来我们需要把线性变换转置并倒着计算，就能得到答案。

注意到 \(Q(x)\) 和 \(F\) 无关，看作常量，我们无需关心。

对于 \(P(x)\) 我们把它看作 \(F\) 进行一些线性变换，以及一些相加后的结果。

不妨假设叶子 \(i\) 的线性变换为 \(t[i]\)，其中 \(t[i]_{0,i} = 1\)，其余位置为 \(0\)，容易发现叶子的 \(P(x)\) 为 \(t[i]F\)。

只有叶子有贡献，那么我们最终的 \(Y'\) 就是所有叶子到根路径的所有矩阵的乘积对 \(F\) 施加变换后得到的向量的求和，最后再用总的 \(\dfrac 1{Q_{rt}(x)}\) 变换。

这里，我们的变换矩阵都是 \(Q(x)\)，即为多项式乘法，转置后就变成了差卷积。

对于 \(t[i]\)，转置后的 \(t[i]^T\) 为：只有 \(t[i]_{i,0}=1\)，其余位置为 \(0\)。

至于倒着计算，我们先对 \(F\) 先与 \(\dfrac 1{Q_{rt}(x)}\) 进行差卷积。

然后分治，对于分治区间 \([l,r]\)，我们传递进一个多项式，然后对左边差卷积 \(Q_1(x)\)，对右边差卷积 \(Q_0(x)\)，递归计算。

最后每个叶子都有一个 \(t[i]^T\) 的线性变换，注意到只有 \(t[i]_{i,0}=1\)，意味这我们只能用该叶子传递进的多项式的 \(0\) 次项对 \(Y_i\) 做贡献。

换句话说：\(Y_i\) 就是该多项式的 \(0\) 次项。

至此，我们以 \(O(n\log^2n)\) 的复杂度完成了多项式多点求值操作。

点击查看代码

#include <bits/stdc++.h>
#define ll long long
#define ull unsigned ll
#define pir pair <ll, ll>
#define mkp make_pair
#define fi first
#define se second
#define pb push_back
using namespace std;
const ll maxn = 64010, mod = 998244353;
ll power(ll a, ll b = mod - 2) {
	ll s = 1;
	while(b) {
		if(b & 1) s = s * a %mod;
		a = a * a %mod; b >>= 1;
	} return s;
}
ll pls(const ll x, const ll y) {return (x + y >= mod? x + y - mod : x + y);}
ll mus(const ll x, const ll y) {return (x < y? x + mod - y : x - y);}
void add(ll &x, const ll y) {x = (x + y >= mod? x + y - mod : x + y);}
void sub(ll &x, const ll y) {x = (x < y? x + mod - y : x - y);}
#define Poly vector <ll>
namespace Polygon {
	ll rev[maxn << 2], tr;
	void Getrev(ll n) {
		if(tr == n) return;
		tr = n;
		for(ll i = 1; i < n; i++)
			rev[i] = (rev[i >> 1] >> 1) | (i & 1? n >> 1 : 0);
	}
	void ntt(ll *a, ll n) {
		Getrev(n);
		for(ll i = 1; i < n; i++)
			if(i < rev[i]) swap(a[i], a[rev[i]]);
		for(ll i = 1; i < n; i <<= 1) {
			ll g = power(3, (mod - 1) / (i << 1));
			for(ll j = 0; j < n; j += (i << 1)) {
				for(ll k = 0, t = 1; k < i; k++, t = t * g %mod) {
					ll x = a[j|k], y = a[i|j|k] * t %mod;
					a[j|k] = pls(x, y), a[i|j|k] = mus(x, y);
				}
			}
		}
	}
	ll a[maxn << 2], b[maxn << 2];
	void Mul(const Poly &A, const Poly &B, Poly &C) {
		ll n = A.size(), m = B.size(), k = n + m - 1;
		for(ll i = 0; i < n; i++) a[i] = A[i];
		for(ll i = 0; i < m; i++) b[i] = B[i];
		ll l = 1;
		while(l < k) l <<= 1;
		ntt(a, l), ntt(b, l);
		for(ll i = 0; i < l; i++) a[i] = a[i] * b[i] %mod;
		ntt(a, l);
		reverse(a + 1, a + l);
		ll Inv = power(l);
		C.resize(k);
		for(ll i = 0; i < l; i++) {
			if(i < k) C[i] = a[i] * Inv %mod;
			a[i] = b[i] = 0;
		}
	}
	Poly Mul(const Poly &A, const Poly &B) {
		Poly C;
		Mul(A, B, C);
		return C;
	}
	ll o[maxn], tt;
	Poly tmp, fw;
	void Getinv(const Poly &A, Poly &B) {
//		puts("need inv");
//		for(ll i = 0; i < A.size(); i++) printf("%lld ", A[i]);
//		puts("");
		ll n = A.size();
		while(n > 1) {
			o[++tt] = n;
			n = (n + 1) >> 1;
		}
		tmp.resize(1);
		tmp[0] = power(A[0]);
		o[tt + 1] = 1;
		for(ll i = tt; i; i--) {
			fw = tmp;
			Mul(fw, fw, fw);
			Poly _A; _A.resize(o[i]);
			for(ll j = 0; j < o[i]; j++) _A[j] = A[j];
			Mul(fw, _A, fw);
			fw.resize(o[i]);
			tmp.resize(o[i]);
			for(ll j = 0; j < o[i]; j++)
				tmp[j] = mus(pls(tmp[j], tmp[j]), fw[j]);
		}
		B = tmp;
//		puts("inv result");
//		for(ll i = 0; i < B.size(); i++)
//			printf("%lld ", B[i]); puts("");
	}
	Poly Getinv(const Poly &A) {
		Poly B;
		Getinv(A, B);
		return B;
	}
	void Mus_Mul(const Poly &A, Poly B, Poly &C) {
		reverse(B.begin(), B.end());
		Mul(A, B, C);
		ll n = A.size(), m = B.size();
		for(ll i = 0; i < n; i++)
			C[i] = C[i + m - 1];
		C.resize(A.size());
	}
	Poly Mus_Mul(const Poly &A, const Poly &B) {
		Poly C;
		Mus_Mul(A, B, C);
		return C;
	}
};
using Polygon::Mul;
using Polygon::Mus_Mul;
using Polygon::Getinv;
ll n, m, a[maxn];
Poly Q[maxn << 4], f;
void calc_Q(ll p, ll l, ll r) {
	if(l == r) {
		Q[p].resize(2);
		Q[p][0] = 1;
		Q[p][1] = mod - a[l];
		return;
	}
	ll mid = l + r >> 1;
	calc_Q(p << 1, l, mid);
	calc_Q(p << 1|1, mid + 1, r);
	Mul(Q[p << 1], Q[p << 1|1], Q[p]);
}
void solve(ll p, ll l, ll r, Poly ret) {
	ret.resize(r - l + 1);
	if(l == r) {
		if(l <= m) printf("%lld\n", ret[0]);
		return;
	}
	ll mid = l + r >> 1;
	solve(p << 1, l, mid, Mus_Mul(ret, Q[p << 1|1]));
	solve(p << 1|1, mid + 1, r, Mus_Mul(ret, Q[p << 1]));
}
int main() {
	scanf("%lld%lld", &n, &m);
	++n;
	f.resize(n);
	for(ll i = 0; i < n; i++) {
		scanf("%lld", &f[i]);
	}
	for(ll i = 1; i <= m; i++) {
		scanf("%lld", a + i);
	}
	if(n < m) f.resize(m), n = m;
	calc_Q(1, 1, n);
	solve(1, 1, n, Mus_Mul(f, Getinv(Q[1])));
    return 0;
}

快速插值

问题：给定一个 \(n-1\) 次多项式在 \(x_0,x_1,...,x_{n - 1}\) 的取值 \(y_0,y_1,...,y_{n - 1}\)，求这个多项式。

考虑拉格朗日插值，我们可以得到

\[ \begin{aligned} f(x) &= \sum_{i = 0} ^ {n - 1} y_i \prod_{j \not = i} \dfrac {x - x_j} {x_i - x_j} \\ &= \sum_{i = 0} ^ {n - 1} y_i \cdot \frac 1{\prod_{j \not = i} (x_i - x_j)} \cdot \prod_{j \not = i} (x - x_j) \end{aligned}\]

我们分为三部分，前面一部分的 \(y_i\) 是常量。

对于中间的 \(\prod_{j \not = i} (x_i - x_j)\)，我们考虑设 \(M(x) = \prod_{i = 0} ^ {n - 1} (x - x_i)\)。

那么 \(\prod_{j \not = i} (x_i - x_j) = \dfrac {\lim_{x\to x_i} M(x)}{\lim_{x\to x_i} x - x_i}\)。

显然分式上下两个式子极限取到 \(0\)，根据洛必达法则，可化为 \(M'(x_i)\)。

\(M'(x)\) 是已知的，我们利用多项式多点求值即可求出 \(M'(x_0),M'(x_1),...,M'(x_{n - 1})\)。

对于最后一部分 \(\prod_{j \not = i} (x - x_j)\)，我们考虑分治解决：

设 \(m = \lfloor \dfrac {n - 1} 2 \rfloor\)。

\[\begin{aligned} f(x) &= \sum_{i = 0} ^ m \frac {y_i} {M'(x_i)} \cdot \prod_{j \not = i} (x - x_j) + \sum_{i = m + 1} ^ {n - 1} \frac {y_i} {M'(x_i)} \cdot \prod_{j \not = i} (x - x_j) \\ &= \left (\prod_{j = m + 1} ^ {n - 1} (x - x_j) \right )\sum_{i = 0} ^ m \frac {y_i} {M'(x_i)} \cdot \prod_{0\le j \le m, j \not = i} (x - x_j) + \left (\prod_{j = 0} ^ m (x - x_j) \right )\sum_{i = m + 1} ^ {n - 1} \frac {y_i} {M'(x_i)} \cdot \prod_{m + 1\le j\le n - 1,j \not = i} (x - x_j) \\ &= f_0(x)\prod_{j = m + 1} ^ {n - 1} (x - x_j) + f_1(x)\prod_{j = 0} ^ m (x - x_j) \end{aligned}\]

其中 \(\prod_{j = l} ^ r (x - x_j)\) 是分治求 \(M(x)\) 的产物，直接计算即可。

点击查看代码

#include <bits/stdc++.h>
#define ll long long
#define ull unsigned ll
#define pir pair <ll, ll>
#define mkp make_pair
#define fi first
#define se second
#define pb push_back
using namespace std;
const ll maxn = 1e5 + 10, mod = 998244353;
ll power(ll a, ll b = mod - 2) {
	ll s = 1;
	while(b) {
		if(b & 1) s = s * a %mod;
		a = a * a %mod; b >>= 1;
	} return s;
}
ll pls(const ll x, const ll y) {return (x + y >= mod? x + y - mod : x + y);}
ll mus(const ll x, const ll y) {return (x < y? x + mod - y : x - y);}
void add(ll &x, const ll y) {x = (x + y >= mod? x + y - mod : x + y);}
void sub(ll &x, const ll y) {x = (x < y? x + mod - y : x - y);}
namespace Polygon {
	#define Poly vector <ll>
	void print(const Poly &A) {
		ll n = A.size();
		printf("size = %lld\n",n);
		for(ll i = 0; i < n; i++) printf("%lld ", A[i]);
		puts("");
	}
	const Poly operator + (const Poly &A, const Poly &B) {
		Poly C = A;
		if(B.size() > C.size()) C.resize(B.size());
		for(ll i = 0; i < B.size(); i++) add(C[i], B[i]);
		return C;
	}
	const Poly operator - (const Poly &A, const Poly &B) {
		Poly C = A;
		if(B.size() > C.size()) C.resize(B.size());
		for(ll i = 0; i < B.size(); i++) sub(C[i], B[i]);
		return C;
	}
	const Poly operator * (const ll k, const Poly &A) {
		Poly C = A;
		for(ll i = 0; i < C.size(); i++) C[i] = C[i] * k %mod;
		return C; 
	}
	ll rev[maxn << 2], tr;
	void Getrev(ll n) {
		if(tr == n) return;
		tr = n;
		for(ll i = 1; i < n; i++)
			rev[i] = (rev[i >> 1] >> 1) | (i & 1? n >> 1 : 0);
	}
	void ntt(ll *a, ll n) {
		Getrev(n);
		for(ll i = 1; i < n; i++)
			if(i < rev[i]) swap(a[i], a[rev[i]]);
		for(ll i = 1; i < n; i <<= 1) {
			ll g = power(3, (mod - 1) / (i << 1));
			for(ll j = 0; j < n; j += (i << 1)) {
				for(ll k = 0, t = 1; k < i; k++, t = t * g %mod) {
					ll x = a[j|k], y = a[i|j|k] * t %mod;
					a[j|k] = pls(x, y), a[i|j|k] = mus(x, y);
				}
			}
		}
	}
	ll a[maxn << 2], b[maxn << 2];
	void Mul(const Poly &A, const Poly &B, Poly &C) {
		ll n = A.size(), m = B.size(), k = n + m - 1;
		for(ll i = 0; i < n; i++) a[i] = A[i];
		for(ll i = 0; i < m; i++) b[i] = B[i];
		ll l = 1;
		while(l < k) l <<= 1;
		ntt(a, l), ntt(b, l);
		for(ll i = 0; i < l; i++) a[i] = a[i] * b[i] %mod;
		ntt(a, l);
		reverse(a + 1, a + l);
		ll Inv = power(l);
		C.resize(k);
		for(ll i = 0; i < l; i++) {
			if(i < k) C[i] = a[i] * Inv %mod;
			a[i] = b[i] = 0;
		}
	}
	Poly Mul(const Poly &A, const Poly &B) {
		Poly C;
		Mul(A, B, C);
		return C;
	}
	ll o[maxn], tt;
	Poly tmp, fw;
	void Getinv(const Poly &A, Poly &B) {
		ll n = A.size();
		while(n > 1) {
			o[++tt] = n;
			n = (n + 1) >> 1;
		}
		tmp.resize(1);
		tmp[0] = power(A[0]);
		o[tt + 1] = 1;
		for(ll i = tt; i; i--) {
			fw = tmp;
			Mul(fw, fw, fw);
			Poly _A; _A.resize(o[i]);
			for(ll j = 0; j < o[i]; j++) _A[j] = A[j];
			Mul(fw, _A, fw);
			fw.resize(o[i]);
			tmp.resize(o[i]);
			tmp = 2 * tmp - fw;
		}
		B = tmp;
	}
	Poly Getinv(const Poly &A) {
		Poly B;
		Getinv(A, B);
		return B;
	}
	void Mus_Mul(const Poly &A, Poly B, Poly &C) {
		reverse(B.begin(), B.end());
		Mul(A, B, C);
		ll n = A.size(), m = B.size();
		for(ll i = 0; i < n; i++)
			C[i] = C[i + m - 1];
		C.resize(A.size());
	}
	Poly Mus_Mul(const Poly &A, const Poly &B) {
		Poly C;
		Mus_Mul(A, B, C);
		return C;
	}
	void Dao(const Poly &A, Poly &B) {
		B.resize(A.size() - 1);
		for(ll i = 0; i < B.size(); i++)
			B[i] = A[i + 1] * (i + 1) %mod;
	}
	Poly Dao(Poly A) {
		for(ll i = 1; i < A.size(); i++)
			A[i - 1] = A[i] * i %mod;
		A.pop_back();
	}
	namespace MPE_FI {
		Poly Q[maxn << 4], A, res;
		ll n, m;
		void calc_Q(ll p, ll l, ll r) {
			if(l == r) {
				Q[p].resize(2);
				Q[p][0] = 1;
				Q[p][1] = mod - A[l];
				return;
			}
			ll mid = l + r >> 1;
			calc_Q(p << 1, l, mid);
			calc_Q(p << 1|1, mid + 1, r);
			Mul(Q[p << 1], Q[p << 1|1], Q[p]);
//			printf("Q[%lld, %lld]\n", l, r);
//			print(Q[p]);
		}
		void solve(ll p, ll l, ll r, Poly ret) {
			ret.resize(r - l + 1);
			if(l == r) {
				if(l < m) res[l] = ret[0];
				return;
			}
			ll mid = l + r >> 1;
			solve(p << 1, l, mid, Mus_Mul(ret, Q[p << 1|1]));
			solve(p << 1|1, mid + 1, r, Mus_Mul(ret, Q[p << 1]));
		}
		Poly solve(Poly F, Poly _A) {
			A = _A;
			n = F.size(), m = A.size();
			if(n < m) F.resize(m), n = m;
			else A.resize(n);
			calc_Q(1, 0, n - 1);
			res.resize(m);
			solve(1, 0, n - 1, Mus_Mul(F, Getinv(Q[1])));
			return res;
		}
		Poly M[maxn << 2], M_, G, X;
		Poly Y;
		void calc_M(ll p, ll l, ll r) {
			if(l == r) {
				M[p].resize(2);
				M[p][0] = mod - X[l];
				M[p][1] = 1;
				return;
			}
			ll mid = l + r >> 1;
			calc_M(p << 1, l, mid);
			calc_M(p << 1|1, mid + 1, r);
			Mul(M[p << 1], M[p << 1|1], M[p]);
		}
		Poly solve_FI(ll p, ll l, ll r) {
			if(l == r) {
				Poly r;
				r.pb(Y[l] * power(G[l]) %mod);
				return r;
			}
			ll mid = l + r >> 1;
			return Mul(solve_FI(p << 1, l, mid), M[p << 1|1]) +
			Mul(solve_FI(p << 1|1, mid + 1, r), M[p << 1]);
		}
		Poly solve_FI(const Poly &_X, const Poly &_Y) {
			X = _X, Y = _Y;
			ll n = X.size();
			calc_M(1, 0, n - 1);
//			print(M[1]);
			Dao(M[1], M_);
//			puts("M_");
//			print(M_);
			G = solve(M_, X);
//			puts("G");
//			print(G);
			return solve_FI(1, 0, n - 1);
		}
	};
	Poly Evaluation(const Poly &F, const Poly &A) {
		return MPE_FI::solve(F, A);
	}
	Poly Interpolation(const Poly &X, const Poly &Y) {
		return MPE_FI::solve_FI(X, Y);
	}
};
using Polygon::Mul;
using Polygon::Mus_Mul;
using Polygon::Getinv;
using Polygon::Dao;
using Polygon::Evaluation;
using Polygon::Interpolation;
ll n;
Poly X, Y;
int main() {
	scanf("%lld", &n);
	X.resize(n), Y.resize(n);
	for(ll i = 0; i < n; i++)
		scanf("%lld%lld", &X[i], &Y[i]);
	Poly ret = Interpolation(X, Y);
	for(ll i = 0; i < n; i++) printf("%lld ", ret[i]);
    return 0;
}