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

LG5050 多项式多点求值

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

答案对 (998244353) 取模。

(n,m∈[1,64000])

题解

我们将要求值的点均分成两份，构造多项式(P_0(x)=prodlimits_{i=1}^{lfloorfrac n 2 floor}(x-x_i))，(P_1(x)=prodlimits_{i=lfloorfrac n 2 floor+1}^{n}(x-x_i))。

显然，对于(forall iin[1,lfloorfrac n 2 floor])，有(P_0(x_i)=0)。(P_1)同理。

我们假设多项式(D(x),R(x))满足：(D(x))是一个(n-lfloorfrac n 2 floor)次多项式，(R(x))是一个次数不超过(lfloorfrac n 2 floor-1)的多项式，且(A(x)=P_0(x)D(x)+R(x))。

那么对于(forall iin[1,lfloorfrac n 2 floor])，有(A(x_i)=R(x_i))。(P_1)同理可得。

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

分治计算即可，由主定理得时间复杂度(O(nlog^2 n))。

求每次的(P_0(x),P_1(x))，可以开始用一次分治NTT预处理。此处时间复杂度也为(O(nlog^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
",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;
}

卡常数做法

定义减法卷积( ext{mul}(F(x),G(x))=sum_isum_{jleq i}f_ig_jx^{i-j})

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

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

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

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

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

时间复杂度(O(nlog^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
",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_imod{998244353})。

(1⩽n⩽100000)

题解

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

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

分离常数

[f(x)=sum_{i = 1}^{n}{ y_iover prod_{i ot = j}{x_i - x_j}} prod_{i ot = j}(x - x_j) ]

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

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

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

则有

[lim_{x o a}{f(x)over g(x)}=lim_{x o 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))这些点算出来的答案，则有

[egin{aligned} f_{l,r}&=sum_{i = l}^{r}{ y_iover g'(x_i)} prod_{j=l,i ot = j}^r(x - x_j)\&=prod_{j=mid+1}^r(x - x_j)sum_{i = l}^{mid}{ y_iover g'(x_i)} prod_{j=l,i ot = j}^{mid}(x - x_j)+prod_{j=l}^{mid}(x - x_j)sum_{i = mid+1}^{r}{ y_iover g'(x_i)} prod_{j=mid+1,i ot = 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(nlog^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;
}