LG5050 多项式多点求值 和 LG5158 多项式快速插值

LG5050 多项式多点求值

给定一个 \(n\) 次多项式 \(f(x)\) ,现在请你对于 \(i \in [1,m]\) ,求出 \(f(a_i)\)

答案对 \(998244353\) 取模。

\(n,m∈[1,64000]\)

题解

我们将要求值的点均分成两份,构造多项式\(P_0(x)=\prod\limits_{i=1}^{\lfloor\frac n 2\rfloor}(x-x_i)\)\(P_1(x)=\prod\limits_{i=\lfloor\frac n 2\rfloor+1}^{n}(x-x_i)\)

显然,对于\(\forall i\in[1,\lfloor\frac n 2 \rfloor]\),有\(P_0(x_i)=0\)\(P_1\)同理。

我们假设多项式\(D(x),R(x)\)满足:\(D(x)\)是一个\(n-\lfloor\frac n 2\rfloor\)次多项式,\(R(x)\)是一个次数不超过\(\lfloor\frac n 2\rfloor-1\)的多项式,且\(A(x)=P_0(x)D(x)+R(x)\)

那么对于\(\forall i\in[1,\lfloor\frac n 2 \rfloor]\),有\(A(x_i)=R(x_i)\)\(P_1\)同理可得。

由于\(R(x)\)的次数是\(A(x)\)的次数的一半,所以我们把一个规模为\(n\)的问题,用\(O(n\log n)\)的复杂度(多项式取模,详见多项式除法),转化为两个规模为\(\frac n 2\)的问题。

分治计算即可,由主定理得时间复杂度\(O(n\log^2 n)\)

求每次的\(P_0(x),P_1(x)\),可以开始用一次分治NTT预处理。此处时间复杂度也为\(O(n\log^2 n)\)

typedef vector<int> polynomial;
void num_trans(polynomial&a,int dir){
	int lim=a.size();
	static vector<int> rev,w[2];
	if(rev.size()!=lim){
		rev.resize(lim);
		int len=log2(lim);
		for(int i=0;i<lim;++i) rev[i]=rev[i>>1]>>1|(i&1)<<(len-1);
		for(int dir=0;dir<2;++dir){
			w[dir].resize(lim);
			w[dir][0]=1,w[dir][1]=fpow(g[dir],(mod-1)/lim);
			for(int i=2;i<lim;++i) w[dir][i]=mul(w[dir][i-1],w[dir][1]);
		}
	}
	for(int i=0;i<lim;++i)if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int step=1;step<lim;step<<=1){
		int quot=lim/(step<<1);
		for(int i=0;i<lim;i+=step<<1){
			int j=i+step;
			for(int k=0;k<step;++k){
				int t=mul(w[dir][quot*k],a[j+k]);
				a[j+k]=add(a[i+k],mod-t),a[i+k]=add(a[i+k],t);
			}
		}
	}
	if(dir){
		int ilim=fpow(lim,mod-2);
		for(int i=0;i<lim;++i) a[i]=mul(a[i],ilim);
	}
}
polynomial poly_inv(polynomial a,int n){
	polynomial b(1,fpow(a[0],mod-2));
	if(n==1) return b;
	int lim=2;
	for(;lim<n;lim<<=1){
		polynomial a1(a.begin(),a.begin()+lim);
		a1.resize(lim<<1),num_trans(a1,0);
		b.resize(lim<<1),num_trans(b,0);
		for(int i=0;i<lim<<1;++i) b[i]=mul(add(2,mod-mul(a1[i],b[i])),b[i]);
		num_trans(b,1),b.resize(lim);
	}
	a.resize(lim<<1),num_trans(a,0);
	b.resize(lim<<1),num_trans(b,0);
	for(int i=0;i<lim<<1;++i) b[i]=mul(add(2,mod-mul(a[i],b[i])),b[i]);
	num_trans(b,1),b.resize(n);
	return b;
}
polynomial poly_div(polynomial f,polynomial g){
	int n=f.size()-1,m=g.size()-1;
	reverse(g.begin(),g.end()),g.resize(n-m+1),g=poly_inv(g,n-m+1);
	reverse(f.begin(),f.end()),f.resize(n-m+1);
	int lim=1<<int(ceil(log2(((n-m)<<1)+1)));
	f.resize(lim),num_trans(f,0);
	g.resize(lim),num_trans(g,0);
	for(int i=0;i<lim;++i) f[i]=mul(f[i],g[i]);
	num_trans(f,1),f.resize(n-m+1);
	return reverse(f.begin(),f.end()),f;
}
polynomial poly_mod(polynomial f,polynomial g){
	int n=f.size()-1,m=g.size()-1;
	polynomial q=poly_div(f,g);
	int lim=1<<int(ceil(log2(n+1)));
	q.resize(lim),num_trans(q,0);
	g.resize(lim),num_trans(g,0);
	for(int i=0;i<lim;++i) q[i]=mul(q[i],g[i]);
	num_trans(q,1);
	for(int i=0;i<m;++i) f[i]=add(f[i],mod-q[i]);
	return f.resize(m),f;
}

co int N=64000+1;
int a[N];
polynomial s[N<<2];

#define lc (x<<1)
#define rc (x<<1|1)
polynomial poly_mul(polynomial a,polynomial b){
	int n=a.size()-1,m=b.size()-1;
	int lim=1<<int(ceil(log2(n+m+1)));
	a.resize(lim),num_trans(a,0);
	b.resize(lim),num_trans(b,0);
	for(int i=0;i<lim;++i) a[i]=mul(a[i],b[i]);
	num_trans(a,1),a.resize(n+m+1);
	return a;
}
void build(int x,int l,int r){
	if(l==r){
		s[x].resize(2);
		s[x][0]=mod-a[l],s[x][1]=1;
		return;
	}
	int mid=(l+r)>>1;
	build(lc,l,mid),build(rc,mid+1,r);
	s[x]=poly_mul(s[lc],s[rc]);
}
void solve(int x,int l,int r,polynomial f){
	if(l==r){
		printf("%d\n",f[0]);
		return;
	}
	int mid=(l+r)>>1;
	solve(lc,l,mid,poly_mod(f,s[lc])),solve(rc,mid+1,r,poly_mod(f,s[rc]));
}
int main(){
	int n=read<int>(),m=read<int>();
	if(!m) return 0;
	polynomial f(n+1);
	for(int i=0;i<=n;++i) read(f[i]);
	for(int i=1;i<=m;++i) read(a[i]);
	build(1,1,m);
	if(n>=m) f=poly_mod(f,s[1]);
	solve(1,1,m,f);
	return 0;
}

卡常数做法

定义减法卷积\(\text{mul}(F(x),G(x))=\sum_i\sum_{j\leq i}f_ig_jx^{i-j}\)

注意到\(F(a_i)=[x^0]\text{mul}(F(x),\frac{1}{1-a_ix})\)

减法卷积有性质\(\text{mul}(F(x),G(x)H(x))=\text{mul}(\text{mul}(F(x),G(x)),H(x))\),证明考虑\(i-(j+k)=i-j-k\)

于是类似分治,根节点传入\(F_{1,m}(x)=\text{mul}(F(x),\frac{1}{\prod_{i=1}^m1-a_ix})\)

\([l,r]\)递归的时候,\([l,\text{mid}]\)传入\(F_{l,\text{mid}}(x)=\text{mul}(F_{l,r}(x),\prod_{i=\text{mid}+1}^r(1-a_ix))\)\([\text{mid}+1,r]\)传入\(F_{\text{mid+1,r}(x)}=\text{mul}(F_{l,r}(x),\prod_{i=l}^{\text{mid}}(1-a_ix))\)

区间\([l,r]\)的多项式只用保留到\(x^{r-l+1}\),因为之后的项不可能通过减法卷积贡献到\([x^0]\)

时间复杂度\(O(n\log^2 n)\),因为没有多项式取模所以常数小。

CO int N=1<<17;
int omg[2][N],rev[N];

void NTT(poly&a,int dir){
	int lim=a.size(),len=log2(lim);
	for(int i=0;i<lim;++i){
		int r=rev[i]>>(17-len);
		if(i<r) swap(a[i],a[r]);
	}
	for(int i=1;i<lim;i<<=1)
		for(int j=0;j<lim;j+=i<<1)for(int k=0;k<i;++k){
			int t=mul(omg[dir][N/(i<<1)*k],a[j+i+k]);
			a[j+i+k]=add(a[j+k],mod-t),a[j+k]=add(a[j+k],t);
		}
	if(dir){
		int ilim=fpow(lim,mod-2);
		for(int i=0;i<lim;++i) a[i]=mul(a[i],ilim);
	}
}
poly operator~(poly a){
	int n=a.size();
	poly b={fpow(a[0],mod-2)};
	if(n==1) return b;
	a.resize(1<<(int)ceil(log2(n)));
	for(int lim=2;lim<2*n;lim<<=1){
		poly c(a.begin(),a.begin()+lim);
		c.resize(lim<<1),NTT(c,0);
		b.resize(lim<<1),NTT(b,0);
		for(int i=0;i<lim<<1;++i) b[i]=mul(2+mod-mul(c[i],b[i]),b[i]);
		NTT(b,1),b.resize(lim);
	}
	return b.resize(n),b;
}
poly operator*(poly a,poly b){
	if(a.size()<=20 or b.size()<=20){
		int n=a.size()-1,m=b.size()-1;
		a.resize(n+m+1);
		for(int i=n;i>=0;--i){
			for(int j=m;j>=1;--j) a[i+j]=add(a[i+j],mul(a[i],b[j]));
			a[i]=mul(a[i],b[0]);
		}
		return a;
	}
	int n=a.size()+b.size()-1,lim=1<<(int)ceil(log2(n));
	a.resize(lim),NTT(a,0);
	b.resize(lim),NTT(b,0);
	for(int i=0;i<lim;++i) a[i]=mul(a[i],b[i]);
	NTT(a,1),a.resize(n);
	return a;
}
poly operator/(poly a,poly b){
	if(a.size()<=20 or b.size()<=20){
		int n=a.size()-1,m=b.size()-1;
		for(int i=0;i<=n;++i){
			for(int j=min(i,m);j>=1;--j) a[i-j]=add(a[i-j],mul(a[i],b[j]));
			a[i]=mul(a[i],b[0]);
		}
		return a;
	}
	int n=a.size()-1,m=b.size()-1;
	if(n<m) b.resize(n+1),m=n;
	reverse(b.begin(),b.end());
	int lim=1<<(int)ceil(log2(n+m+1));
	a.resize(lim),NTT(a,0);
	b.resize(lim),NTT(b,0);
	for(int i=0;i<lim;++i) a[i]=mul(a[i],b[i]);
	NTT(a,1);
	for(int i=0;i<=n;++i) a[i]=a[i+m];
	return a.resize(n+1),a;
}

int a[N];
poly g[N];

#define lc (x<<1)
#define rc (x<<1|1)
#define mid ((l+r)>>1)
void build(int x,int l,int r){
	if(l==r) {g[x]={1,mod-read<int>()}; return;}
	build(lc,l,mid),build(rc,mid+1,r);
	g[x]=g[lc]*g[rc];
}
void solve(int x,int l,int r,CO poly&f){
	if(l==r) {printf("%d\n",f[0]); return;}
	poly lf=f/g[rc];lf.resize(mid-l+2);
	solve(lc,l,mid,lf);
	poly rf=f/g[lc];rf.resize(r-mid+1);
	solve(rc,mid+1,r,rf);
}
#undef lc
#undef rc
#undef mid

int main(){
	omg[0][0]=1,omg[0][1]=fpow(3,(mod-1)/N);
	omg[1][0]=1,omg[1][1]=fpow(omg[0][1],mod-2);
	rev[0]=0,rev[1]=1<<16;
	for(int i=2;i<N;++i){
		omg[0][i]=mul(omg[0][i-1],omg[0][1]);
		omg[1][i]=mul(omg[1][i-1],omg[1][1]);
		rev[i]=rev[i>>1]>>1|(i&1)<<16;
	}
	int n=read<int>(),m=read<int>();
	poly f(n+1);
	for(int i=0;i<=n;++i) read(f[i]);
	build(1,1,m);
	g[1].resize(n+1),f=f/~g[1],f.resize(m+1);
	solve(1,1,m,f);
	return 0;
}

LG5158 多项式快速插值

给出 \(n\) 个点 \((X_i,Y_i)\)。求一个 \(n-1\) 次的多项式 \(f(x)\),使得 \(f(X_i)= Y_i\mod{998244353}\)

\(1⩽n⩽100000\)

题解

思路是优化拉格朗日插值。先看看原式

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

分离常数

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

来考虑一下\({ y_i\over \prod_{i \not = j}{x_i - x_j}}\)这个东西,上面是个常数,那么只考虑下面。如果我们设多项式\(g(x)=\prod_{i=1}^n (x-x_i)\),那么下面那个东西就是\({g(x_i)\over (x-x_i)}\)

这分子分母全为\(0\)怎么求?根据洛必达法则,如果

\[\lim_{x\to a}f(x)=0,\lim_{x\to a}g(x)=0 \]

则有

\[\lim_{x\to a}{f(x)\over g(x)}=\lim_{x\to a}{f'(x)\over g'(x)} \]

那么我们代入之后可以发现\({g(x_i)\over (x-x_i)}=g'(x_i)\)

先分治NTT算出\(g\),然后对\(g'\)多点求值把每个点处的值算出来就好了

那么接下来我们就考虑分治,设\(f_{l,r}\)表示\((x_l,y_l)\sim(x_r,y_r)\)这些点算出来的答案,则有

\[\begin{aligned} f_{l,r}&=\sum_{i = l}^{r}{ y_i\over g'(x_i)} \prod_{j=l,i \not = j}^r(x - x_j)\\&=\prod_{j=mid+1}^r(x - x_j)\sum_{i = l}^{mid}{ y_i\over g'(x_i)} \prod_{j=l,i \not = j}^{mid}(x - x_j)+\prod_{j=l}^{mid}(x - x_j)\sum_{i = mid+1}^{r}{ y_i\over g'(x_i)} \prod_{j=mid+1,i \not = j}^{r}(x - x_j)\\&=\prod_{j=mid+1}^r(x - x_j)f_{l,mid}+\prod_{j=l}^{mid}(x - x_j)f_{mid+1,r}\\ \end{aligned} \]

时间复杂度为\(O(n\log^2n)\)。这题卡常,需要小范围多项式乘法暴力才能过。

typedef vector<int> polynomial;
void num_trans(polynomial&a,int dir){
	int lim=a.size();
	static vector<int> rev,w[2];
	if(rev.size()!=lim){
		rev.resize(lim);
		int len=log2(lim);
		for(int i=0;i<lim;++i) rev[i]=rev[i>>1]>>1|(i&1)<<(len-1);
		for(int dir=0;dir<2;++dir){
			static co int g[2]={3,332748118};
			w[dir].resize(lim);
			w[dir][0]=1,w[dir][1]=fpow(g[dir],(mod-1)/lim);
			for(int i=2;i<lim;++i) w[dir][i]=mul(w[dir][i-1],w[dir][1]);
		}
	}
	for(int i=0;i<lim;++i)if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int step=1;step<lim;step<<=1){
		int quot=lim/(step<<1);
		for(int i=0;i<lim;i+=step<<1){
			int j=i+step;
			for(int k=0;k<step;++k){
				int t=mul(w[dir][quot*k],a[j+k]);
				a[j+k]=add(a[i+k],mod-t),a[i+k]=add(a[i+k],t);
			}
		}
	}
	if(dir){
		int ilim=fpow(lim,mod-2);
		for(int i=0;i<lim;++i) a[i]=mul(a[i],ilim);
	}
}
polynomial poly_inv(polynomial a,int n){
	polynomial b(1,fpow(a[0],mod-2));
	if(n==1) return b;
	int lim=2;
	for(;lim<n;lim<<=1){
		polynomial a1(a.begin(),a.begin()+lim);
		a1.resize(lim<<1),num_trans(a1,0);
		b.resize(lim<<1),num_trans(b,0);
		for(int i=0;i<lim<<1;++i) b[i]=mul(add(2,mod-mul(a1[i],b[i])),b[i]);
		num_trans(b,1),b.resize(lim);
	}
	a.resize(lim<<1),num_trans(a,0);
	b.resize(lim<<1),num_trans(b,0);
	for(int i=0;i<lim<<1;++i) b[i]=mul(add(2,mod-mul(a[i],b[i])),b[i]);
	num_trans(b,1),b.resize(n);
	return b;
}
polynomial poly_div(polynomial f,polynomial g){
	int n=f.size()-1,m=g.size()-1;
	reverse(g.begin(),g.end()),g.resize(n-m+1),g=poly_inv(g,n-m+1);
	reverse(f.begin(),f.end()),f.resize(n-m+1);
	int lim=1<<int(ceil(log2(((n-m)<<1)+1)));
	f.resize(lim),num_trans(f,0);
	g.resize(lim),num_trans(g,0);
	for(int i=0;i<lim;++i) f[i]=mul(f[i],g[i]);
	num_trans(f,1),f.resize(n-m+1);
	return reverse(f.begin(),f.end()),f;
}
polynomial poly_mod(polynomial f,polynomial g){
	int n=f.size()-1,m=g.size()-1;
	polynomial q=poly_div(f,g);
	int lim=1<<int(ceil(log2(n+1)));
	q.resize(lim),num_trans(q,0);
	g.resize(lim),num_trans(g,0);
	for(int i=0;i<lim;++i) q[i]=mul(q[i],g[i]);
	num_trans(q,1);
	for(int i=0;i<m;++i) f[i]=add(f[i],mod-q[i]);
	return f.resize(m),f;
}

polynomial poly_mul(polynomial a,polynomial b){
	int n=a.size()-1,m=b.size()-1;
	if((LL)n*m<=5000){ // edit 1
		polynomial c(n+m+1);
		for(int i=0;i<=n;++i)
			for(int j=0;j<=m;++j) c[i+j]=add(c[i+j],mul(a[i],b[j]));
		return c;
	}
	int lim=1<<int(ceil(log2(n+m+1)));
	a.resize(lim),num_trans(a,0);
	b.resize(lim),num_trans(b,0);
	for(int i=0;i<lim;++i) a[i]=mul(a[i],b[i]);
	num_trans(a,1),a.resize(n+m+1);
	return a;
}
polynomial poly_add(polynomial a,polynomial b){
	if(a.size()<b.size()) swap(a,b);
	for(int i=0;i<b.size();++i) a[i]=add(a[i],b[i]);
	return a;
}

co int N=100000+1;
#define lc (x<<1)
#define rc (x<<1|1)
namespace eval{
	int n,m,a[N],ans[N];
	polynomial f,s[N<<2];
	
	void build(int x,int l,int r){
		if(l==r){
			s[x].resize(2);
			s[x][0]=mod-a[l],s[x][1]=1;
			return;
		}
		int mid=(l+r)>>1;
		build(lc,l,mid),build(rc,mid+1,r);
		s[x]=poly_mul(s[lc],s[rc]);
	}
	void solve(int x,int l,int r,polynomial f){
		if(l==r){
			ans[l]=f[0];
			return;
		}
		int mid=(l+r)>>1;
		solve(lc,l,mid,poly_mod(f,s[lc])),solve(rc,mid+1,r,poly_mod(f,s[rc]));
	}
	void main(){
		build(1,1,m);
		if(n>=m) f=poly_mod(f,s[1]);
		solve(1,1,m,f);
	}
}

int X[N],Y[N];
polynomial g[N<<2],f[N<<2];
void build(int x,int l,int r){
	if(l==r){
		g[x].resize(2);
		g[x][0]=mod-X[l],g[x][1]=1;
		return;
	}
	int mid=(l+r)>>1;
	build(lc,l,mid),build(rc,mid+1,r);
	g[x]=poly_mul(g[lc],g[rc]);
}
void solve(int x,int l,int r){
	if(l==r){
		f[x].assign(1,mul(Y[l],fpow(eval::ans[l],mod-2)));
		return;
	}
	int mid=(l+r)>>1;
	solve(lc,l,mid),solve(rc,mid+1,r);
	f[x]=poly_add(poly_mul(f[lc],g[rc]),poly_mul(f[rc],g[lc]));
}
int main(){
//	freopen("LG5158.in","r",stdin),freopen("LG5158.out","w",stdout);
	int n=read<int>();
	for(int i=1;i<=n;++i) read(X[i]),read(Y[i]);
	build(1,1,n);
	eval::n=n-1,eval::f.resize(n);
	for(int i=0;i<=eval::n;++i) eval::f[i]=mul(g[1][i+1],i+1);
	eval::m=n;
	for(int i=1;i<=eval::m;++i) eval::a[i]=X[i];
	eval::main();
	solve(1,1,n);
	for(int i=0;i<n;++i) printf("%d ",f[1][i]);
	return 0;
}

posted on 2020-07-23 20:10  autoint  阅读(110)  评论(0编辑  收藏

导航