多项式的多点求值和快速插值

前言

常数好大

---

多点求值

描述

给一个 \(n-1\) 次的多项式 \(f(x)\) 和 \(m\) 个位置 \(x_0,x_1,..,x_{m-1}\)，求\(f(x_0),f(x_1),..,f(x_{m-1})\)

暴力

做法

不难发现

\[f(x_0)\equiv f(x)\pmod{(x-x_0)}\]

因为设 \(f(x)=d(x)(x-x_0)+r\)，显然有\(f(x_0)=r\)

于是我们可以分治，设

\[ \begin{align} L(x)&=\prod_{i=0}^{\lfloor\frac{m}{2}\rfloor}(x-x_i) \\ R(x)&=\prod_{i=\lfloor\frac{m}{2}\rfloor+1}^{m-1}(x-x_i) \end{align} \]

对于 \(0\le i\le \lfloor\frac{m}{2}\rfloor\)，有 \(f(x_i)=(f\ mod\ L)(x_i)\)

对于 \(\lfloor\frac{m}{2}\rfloor<i<m\)，有 \(f(x_i)=(f\ mod\ R)(x_i)\)

可以把次数减少一半

求出每个 \(L(x)\) 和 \(R(x)\) 可以用分治FFT，复杂度和后面的取模相同

\(n,m\) 同阶，总复杂度 \(\mathcal O(m\log^2m)\)

---

插值

描述

给 \(n\) 个点 \((x_0,y_0),(x_1,y_1),..,(x_{n-1},y_{n-1})\)，求一个 \(n-1\) 次多项式 \(f(x)\) 满足 \(f(x_i)=y_i\)

暴力

有 \(\mathcal O(n^2)\) 拉格朗日插值，这里不展开

快速做法

还是考虑拉格朗日插值

\[ f(x)= \sum_{i=0}^{n-1} \frac {\prod_{j\ne i}(x-x_j)} {\prod_{j\ne i}(x_i-x_j)} y_i \]

先计算每个 \(\prod_{j\ne i}(x_i-x_j)\)

设

\[ \begin{align} M(x)&=\prod_{i=0}^{n-1} (x-x_i) \\ g_i(x)&=\frac{M(x)}{x-x_i} \end{align} \]

那么

\[\prod_{j\ne i}(x_i-x_j)=g_i(x_i)\]

分子分母都是 \(0\)，根据洛必达法则

\[\prod_{j\ne i}(x_i-x_j)=g_i(x_i)=M'(x_i)\]

于是可以使用多点求值在 \(\mathcal O(n\log^2n)\) 的时间内求得

设 \(z_i=\frac{y_i}{\prod_{j\ne i}(x_i-x_j)}\)，问题就是求 \(f(x)=\sum_{i=0}^{n-1}z_i\prod_{j\ne i}(x-x_j)\)

可以分治FFT，同样地设

\[ \begin{align} L(x)&=\prod_{i=0}^{\lfloor\frac{n}{2}\rfloor}(x-x_i) \\ R(x)&=\prod_{i=\lfloor\frac{n}{2}\rfloor+1}^{n-1}(x-x_i) \end{align} \]

于是

\[ f(x)= R(x) \sum_{i=0}^{\lfloor\frac{n}{2}\rfloor}z_i \prod_{0\le j\le \lfloor\frac{n}{2}\rfloor,j\ne i} (x-x_j) + L(x) \sum_{i=\lfloor\frac{n}{2}\rfloor+1}^{n-1} \prod_{\lfloor\frac{n}{2}\rfloor<j<n,j\ne i} (x-x_j) \]

总复杂度还是 \(\mathcal O(n\log^2n)\)

---

例题

NFLSOJ #53. 多项式多点求值和多点插值

模板

这是一个不完整的多项式板子.

之后会放一个比较全的

在这里挑战多项式

#include<cstdio>
#include<algorithm>
#include<ctype.h>
#include<string.h>
#include<math.h>
#include<vector>

using namespace std;
#define ll long long

inline char read() {
	static const int IN_LEN = 1000000;
	static char buf[IN_LEN], *s, *t;
	return (s==t?t=(s=buf)+fread(buf,1,IN_LEN,stdin),(s==t?-1:*s++):*s++);
}
template<class T>
inline void read(T &x) {
	static bool iosig;
	static char c;
	for (iosig=false, c=read(); !isdigit(c); c=read()) {
		if (c == '-') iosig=true;
		if (c == -1) return;
	}
	for (x=0; isdigit(c); c=read()) x=((x+(x<<2))<<1)+(c^'0');
	if (iosig) x=-x;
}
const int OUT_LEN = 10000000;
char obuf[OUT_LEN], *ooh=obuf;
inline void print(char c) {
	if (ooh==obuf+OUT_LEN) fwrite(obuf, 1, OUT_LEN, stdout), ooh=obuf;
	*ooh++=c;
}
template<class T>
inline void print(T x) {
	static int buf[30], cnt;
	if (x==0) print('0');
	else {
		if (x<0) print('-'), x=-x;
		for (cnt=0; x; x/=10) buf[++cnt]=x%10+48;
		while(cnt) print((char)buf[cnt--]);
	}
}
inline void flush() { fwrite(obuf, 1, ooh - obuf, stdout); }

int n, m;
vector<int> x, y, z;

namespace Poly{
	const int P = 998244353;

	vector<int> ans;// for Evaluate()
	vector<vector<int>> p;// for Evaluate() & Interpolate()

	inline int Pow(ll x, int y=P-2){ // x^y
		int ans=1;
		for(; y; y>>=1, x=x*x%P) if(y&1) ans=ans*x%P;
		return ans;
	}
	inline int Ge(int x){ int n=1; while(n<=x) n<<=1; return n;}
	inline int Mod(int x){ return x<P?x:x-P;}
	inline void NTT(vector<int> &f, int g, int n){
		f.resize(n);
		for(int i=0, j=0; i<n; ++i){
			if(i>j) swap(f[i], f[j]);
			for(int k=n>>1; (j^=k)<k; k>>=1);
		}
		vector<int> w(n>>1);
		for(int i=1; i<n; i<<=1){
			for(int j=w[0]=1, w0=(g==1?Pow(3, (P-1)/i/2):Pow(Pow(3, (P-1)/i/2))); j<i; ++j) w[j]=(ll)w[j-1]*w0%P;
			for(int j=0; j<n; j+=i<<1){
				for(int k=j; k<j+i; ++k){
					int t=(ll)f[k+i]*w[k-j]%P;
					f[k+i]=Mod(f[k]-t+P);
					f[k]=Mod(f[k]+t);
				}
			}
		}
		if(g==-1) for(int i=0, I=Pow(n); i<n; ++i) f[i]=(ll)f[i]*I%P;
	}
	inline vector<int> Add(const vector<int> &f, const vector<int> &g){
		vector<int> ans=f;
		for(unsigned i=0; i<f.size(); ++i) (ans[i]+=g[i])%=P;
		return ans;
	}
	inline vector<int> Mul(const vector<int> &f, const vector<int> &g){// f*g
		vector<int> F=f, G=g;
		int p=Ge(f.size()+g.size()-2);
		NTT(F, 1, p), NTT(G, 1, p);
		for(int i=0; i<p; ++i) F[i]=(ll)F[i]*G[i]%P;
		NTT(F, -1, p);
		return F.resize(f.size()+g.size()-1), F;
	}
	inline vector<int> PolyInv(const vector<int> &f, int n=-1){// 1/f
		if(n==-1) n=f.size();
		vector<int> ans;
		if(n==1) return ans.push_back(Pow(f[0])), ans;
		ans=PolyInv(f, (n+1)/2);
		vector<int> tmp(&f[0], &f[0]+n);
		int p=Ge(n*2-2);
		NTT(tmp, 1, p), NTT(ans, 1, p);
		for(int i=0; i<p; ++i) ans[i]=(2-(ll)ans[i]*tmp[i]%P+P)*ans[i]%P;
		NTT(ans, -1, p);
		return ans.resize(n), ans;
	}
	inline void PolyDiv(const vector<int> &a, const vector<int> &b, vector<int> &d, vector<int> &r){//a=d*b+r
		if(b.size()>a.size()) return (void)(d.clear(), r=a);

		vector<int> A=a, B=b, iB;
		int n=a.size(), m=b.size();
		reverse(A.begin(), A.end()), reverse(B.begin(), B.end());
		B.resize(n-m+1), iB=PolyInv(B, n-m+1);
		d=Mul(A, iB);
		d.resize(n-m+1), reverse(d.begin(), d.end());

		r=Mul(b, d);
		for(int i=0; i<m-1; ++i) r[i]=(P+a[i]-r[i])%P;
		r.resize(m-1);
	}
	inline vector<int> Derivative(const vector<int> &a){ // a'
		vector<int> ans;
		ans.resize(a.size()-1);
		for(unsigned i=1; i<a.size(); ++i) ans[i-1]=(ll)a[i]*i%P;
		return ans;
	}
	void Evaluate_Interpolate_Init(int l, int r, int t, const vector<int> &a){
		if(l==r) return p[t].clear(), p[t].push_back(P-a[l]), p[t].push_back(1);
		int mid=(l+r)/2, k=t<<1;
		Evaluate_Interpolate_Init(l, mid, k, a), Evaluate_Interpolate_Init(mid+1, r, k|1, a);
		p[t]=Mul(p[k], p[k|1]);
	}
	inline void Evaluate(int l, int r, int t, const vector<int> &f, const vector<int> &a){
		if(r-l+1<=512){
			for(int i=l; i<=r; ++i){
				int x=0, j=f.size(), a1=a[i], a2=(ll)a[i]*a[i]%P, a3=(ll)a[i]*a2%P, a4=(ll)a[i]*a3%P, a5=(ll)a[i]*a4%P, a6=(ll)a[i]*a5%P, a7=(ll)a[i]*a6%P, a8=(ll)a[i]*a7%P;
				while(j>=8)
				x=((ll)x*a8+(ll)f[j-1]*a7+(ll)f[j-2]*a6+(ll)f[j-3]*a5+(ll)f[j-4]*a4+(ll)f[j-5]*a3+(ll)f[j-6]*a2+(ll)f[j-7]*a1+f[j-8])%P, j-=8;
				while(j--) x=((ll)x*a[i]+f[j])%P;
				ans.push_back(x);
			}
			return;
		}
		vector<int> tmp;
		PolyDiv(f, p[t], tmp, tmp);
		Evaluate(l, (l+r)/2, t<<1, tmp, a), Evaluate((l+r)/2+1, r, t<<1|1, tmp, a);
	}
	inline vector<int> Evaluate(const vector<int> &f, const vector<int> &a, int flag=-1){// f(a_i)
		if(flag==-1) p.resize(a.size()<<2), Evaluate_Interpolate_Init(0, a.size()-1, 1, a);
		ans.clear(), Evaluate(0, a.size()-1, 1, f, a);
		return ans;
	}
	vector<int> Interpolate(int l, int r, int t, const vector<int> &x, const vector<int> &f){
		if(l==r){
			vector<int> ans;
			return ans.push_back(f[l]), ans;
		}
		int mid=(l+r)/2, k=t<<1;
		return Add(Mul(Interpolate(l, mid, k, x, f), p[k|1]), Mul(Interpolate(mid+1, r, k|1, x, f), p[k]));
	}
	inline vector<int> Interpolate(const vector<int> &x, const vector<int> &y){// (x_i,y_i)
		int n=x.size();
		p.resize(n<<2), Evaluate_Interpolate_Init(0, n-1, 1, x);
		vector<int> f=Evaluate(Derivative(p[1]), x, 0);
		for(int i=0; i<n; ++i) f[i]=(ll)y[i]*Pow(f[i])%P;
		return Interpolate(0, n-1, 1, x, f);
	}
}
using namespace Poly;

int main() {
	read(n), x.resize(n), y.resize(n);
	for(int i=0; i<n; ++i) read(x[i]), read(y[i]);
	read(m), z.resize(m);
	for(int i=0; i<m; ++i) read(z[i]);

	x=Evaluate(Interpolate(x, y), z);
	for(int i:x) print(i), print(' ');
	return flush(), 0;
}