拉格朗日插值2

X.拉格朗日插值2

从这题开始,拉格朗日插值就逐渐同多项式同流合污了。

我们列出式子:

\[f(m+i)=\sum\limits_{j=0}^nf(j)\prod\limits_{k\neq j}\dfrac{m+i-k}{j-k} \]

借鉴前面的思想,我们将它转成了

\[f(m+i)=\sum\limits_{j=0}^nf(j)\dfrac{\Big(\dfrac{(m+n+i)!}{(m+i-1)!\times(m+i-j)}\Big)}{j!\times(n-j)!\times(-1)^{n-j}} \]

我们将其整理一下,将只与\(i\)有关、只与\(j\)有关和与\(i,j\)都有关的项分开,就得到了

\[f(m+i)=\dfrac{(m+n+i)!}{(m+i-1)!}\times\sum\limits_{j=0}^n\dfrac{f(j)}{j!\times(n-j)!\times(-1)^{n-j}}\times\dfrac{1}{m+i-j} \]

显然\(\sum\)内部的一大坨是一个卷积的形式,可以NTT优化。

于是我们令\(F(x)=\begin{cases}\dfrac{f(x)}{x!\times(n-x)!\times(-1)^{n-x}}\ (x\leq n)\\0\ (x>n)\end{cases}\)

再令\(G(x)=\dfrac{1}{m-n+x}\)

计算\(FG\),我们会发现\(f(m+i)=\dfrac{(m+n+i)!}{(m+i-1)!}\times(FG)(n+i)\),可以通过暴力展开证明。

则复杂度为\(O(n\log n)\)

代码:

#include<bits/stdc++.h>
using namespace std;
const int N=1<<20;
const int mod=998244353;
const int G=3;
int n,m,fac[N],inv[N],rev[N],f[N],g[N],all,FAC[N];
int ksm(int x,int y){
	int rt=1;
	for(;y;x=(1ll*x*x)%mod,y>>=1)if(y&1)rt=(1ll*rt*x)%mod;
	return rt;
}
void NTT(int *a,int tp,int LG){
	int lim=(1<<LG),invlim=ksm(lim,mod-2);
	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(LG-1));
	for(int i=0;i<lim;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
	for(int md=1;md<lim;md<<=1){
		int rt=ksm(G,(mod-1)/(md<<1));
		if(tp==-1)rt=ksm(rt,mod-2);
		for(int stp=md<<1,pos=0;pos<lim;pos+=stp){
			int w=1;
			for(int i=0;i<md;i++,w=(1ll*w*rt)%mod){
				int x=a[pos+i],y=(1ll*w*a[pos+md+i])%mod;
				a[pos+i]=(x+y)%mod;
				a[pos+md+i]=(x-y+mod)%mod;
			}
		}
	}
	if(tp==-1)for(int i=0;i<lim;i++)a[i]=(1ll*a[i]*invlim)%mod;
}
int A[N],B[N];
void mul(int *a,int *b,int *c,int LG){//using: Array A and B
	int lim=(1<<LG);
	for(int i=0;i<lim;i++)A[i]=B[i]=0;
	for(int i=0;i<(lim>>1);i++)A[i]=a[i],B[i]=b[i];
	NTT(A,1,LG),NTT(B,1,LG);
	for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
	NTT(A,-1,LG);
	for(int i=0;i<lim;i++)c[i]=A[i];
}
int main(){
	scanf("%d%d",&n,&m),fac[0]=FAC[0]=1;
	for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%mod;
	inv[n]=ksm(fac[n],mod-2);
	for(int i=n-1;i>=0;i--)inv[i]=1ll*inv[i+1]*(i+1)%mod;
	for(int i=m;i>=m-n;i--)FAC[0]=1ll*FAC[0]*i%mod;
	for(int i=m+1;i<=m+n;i++)FAC[i-m]=1ll*FAC[i-m-1]*i%mod*ksm(i-n-1,mod-2)%mod;
//	for(int i=0;i<=n;i++)printf("%d ",FAC[i]);puts("");
	for(int i=0;i<=n;i++)scanf("%d",&f[i]),f[i]=1ll*f[i]*inv[i]%mod*inv[n-i]%mod,f[i]=((n-i)&1?(mod-f[i])%mod:f[i]);
	for(int i=0;i<=(n<<1);i++)g[i]=ksm(m-n+i,mod-2);
	while((1<<all)<=n)all++;
	mul(f,g,f,all+2);
	for(int i=0;i<=n;i++)printf("%d ",1ll*f[n+i]*FAC[i]%mod);
	return 0;
} 

原文地址:https://www.cnblogs.com/Troverld/p/14607925.html