多项式进阶操作

多点求值#

问题：给定一个 $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=\left[\begin{array}{cccc}{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{array}\right]$。

那么我们要求的答案为 $Y_i=\sum _{j=0}^{n-1}{a}_i^j{f}_j$，即 $Y=AF$。

考虑转置原理，我们考虑求 $Y^{\prime }=A^{T}F$，有 $Y_i^{\prime }=\sum _{j=0}^{n-1}{a}_j^i{f}_j$。

考虑 $Y_{0..n-1}^{\prime }$ 的生成函数：

$$\begin{array}{rl}\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}^i{a}_j^i\\ & =\sum _{j=0}^{n-1}{f}_j\sum _{i=0}(a_{j}x{)}^i\\ & =\sum _{j=0}^{n-1}{f}_j\cdot {\displaystyle \frac{1}{1-a_{j}x}}\end{array}$$

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

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

$${\displaystyle \frac{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_{i}x$。

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

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

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

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

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

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

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

至于倒着计算，我们先对 $F$ 先与 ${\displaystyle \frac{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 }^{2}n)$ 的复杂度完成了多项式多点求值操作。

点击查看代码

#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{array}{rl}f(x)& =\sum _{i=0}^{n-1}{y}_i\prod _{j\ne i}{\displaystyle \frac{x-x_j}{x_i-x_j}}\\ & =\sum _{i=0}^{n-1}{y}_i\cdot \frac{1}{\prod _{j\ne i}(x_i-x_j)}\cdot \prod _{j\ne i}(x-x_j)\end{array}$$

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

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

那么 $\prod _{j\ne i}(x_i-x_j)={\displaystyle \frac{\underset{x\to x_i}{lim}M(x)}{\underset{x\to x_i}{lim}x-x_i}}$。

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

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

对于最后一部分 $\prod _{j\ne i}(x-x_j)$，我们考虑分治解决：

设 $m=\lfloor {\displaystyle \frac{n-1}{2}}\rfloor$。

$$\begin{array}{rl}f(x)& =\sum _{i=0}^m\frac{y_i}{M^{\prime }(x_i)}\cdot \prod _{j\ne i}(x-x_j)+\sum _{i=m+1}^{n-1}\frac{y_i}{M^{\prime }(x_i)}\cdot \prod _{j\ne i}(x-x_j)\\ & =(\prod _{j=m+1}^{n-1}(x-x_j))\sum _{i=0}^m\frac{y_i}{M^{\prime }(x_i)}\cdot \prod _{0\le j\le m,j\ne i}(x-x_j)+(\prod _{j=0}^m(x-x_j))\sum _{i=m+1}^{n-1}\frac{y_i}{M^{\prime }(x_i)}\cdot \prod _{m+1\le j\le n-1,j\ne 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{array}$$

其中 $\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;
}