luogu P5667 拉格朗日插值2 拉格朗日插值 多项式多点求值 NTT

LINK：P5667 拉格朗日插值2

给出了n个连续的取值的自变量的点值 求 f(m+1),f(m+2),...f(m+n).

如果我们直接把f这个函数给插值出来就变成了了多项式多点求值 这个难度好像有点大。

不妨 直接考虑拉格朗日插值。

设此时要求f(k) 那么则有 (sum_{i=0}^nf(i)frac{Pi_{i eq j}(k-j)}{Pi_{i eq j} (i-j)})

可以化简一下 (f(k)=sum_{i=0}^nf(i)frac{ Pi_{i eq j}(k-j) cdot (-1)^{k-i} }{fac_icdot fac_{n-i}})

其实还是可以化简的 (f(k)=sum_{i=0}^nf(i)frac{ Pi (k-j) cdot (-1)^{k-i} }{fac_icdot fac_{n-i}}cdot frac{1}{k-i})

(f(k)=sum_{i=0}^nf(i)frac{k!cdot (-1)^{k-i} }{fac_icdot fac_{n-i}}cdot frac{1}{(k-i)cdot (k-n-1)!})

再提出来一些项 (f(k)=k!cdot frac{1}{(k-n-1)!}sum_{i=0}^nf(i)frac{(-1)^{k-i} }{fac_icdot fac_{n-i}}cdot frac{1}{(k-i)})

容易发现这类似于卷积的形式。

不过 如果直接做卷积会出现问题 因为卷积的时候i<=j 而并非i<=n.

需要解决这个问题 比较容易想到 把j放到靠后的位置就能得到贡献了。

那么其实就是 所以维护一个长度为2n的序列进行卷积即可。

这个时候 最好的方法就是 把两个多项式写下来 看一下卷积的过程。

分析一下每一个位置都应该是什么数字即可。

注意 卡常的话中间有一个求2n长度逆元的东西 可以采用前缀积后缀积的方法来优化成O(n).

但是 由于多项式长度2^20 还是跑的很慢.. 勉强卡过。

码力还行 出错的地方是 数组开小了 开了1e6 少开了一点 wa了3,4发 我tcl.

const int MAXN=1100010,G=3;
int lim,n,m,N;
int rev[MAXN];
ll f[MAXN],b[MAXN],a[MAXN];
ll fac[MAXN],pre[MAXN],suf[MAXN],inv[MAXN];
inline int ksm(ll b,int p)
{
	ll cnt=1;
	while(p)
	{
		if(p&1)cnt=cnt*b%mod;
		p=p>>1;b=b*b%mod;
	}
	return cnt;
}
inline void NTT(ll *a,int op)
{
	rep(0,lim-1,i)if(i<rev[i])swap(a[i],a[rev[i]]);
	for(int len=2;len<=lim;len=len<<1)
	{
		int mid=len>>1;
		int wn=ksm(G,op==1?(mod-1)/len:mod-1-(mod-1)/len);
		for(int j=0;j<lim;j+=len)
		{
			ll d=1;
			for(int i=0;i<mid;++i)
			{
				ll x=a[i+j],y=a[i+j+mid]*d%mod;
				a[i+j]=(x+y)%mod;a[i+j+mid]=(x-y+mod)%mod;
				d=d*wn%mod;
			}
		}
	}
	if(op==-1)
	{
		int IN=ksm(lim,mod-2);
		rep(0,lim-1,i)a[i]=a[i]*IN%mod;
	}
}
int main()
{
	//freopen("1.in","r",stdin);
	get(n);get(m);fac[0]=1;
	rep(0,n,i)get(f[i]);
	rep(1,n,i)fac[i]=fac[i-1]*i%mod;
	inv[n]=ksm(fac[n],mod-2);
	fep(n-1,0,i)inv[i]=inv[i+1]*(i+1)%mod;
	rep(0,n,i)a[i]=(f[i]*(((n-i)&1)?-1:1)%mod*inv[i]%mod*inv[n-i]%mod+mod)%mod;
	pre[0]=1;N=n<<1|1;suf[N+1]=1;
	rep(1,N,i)pre[i]=pre[i-1]*(m-n+i-1)%mod;
	fep(N,1,i)suf[i]=suf[i+1]*(m+i-n-1)%mod;
	ll IN=ksm(pre[N],mod-2);
	rep(1,N,i)b[i]=pre[i-1]*suf[i+1]%mod*IN%mod;
	lim=1;while(lim<N+N)lim=lim<<1;
	rep(0,lim-1,i)rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
	NTT(a,1);NTT(b,1);
	rep(0,lim-1,i)a[i]=a[i]*b[i]%mod;
	ll w=1;rep(m-n,m,i)w=w*i%mod;
	b[n+1]=w;rep(n+2,N,i)b[i]=b[i-1]*(m+i-n-1)%mod*IN%mod*pre[i-n-2]%mod*suf[i-n]%mod;
	NTT(a,-1);rep(n+1,N,i)printf("%lld ",a[i]*b[i]%mod);return 0;
}

中间的地方看起来确实比较ex... 这也是没有办法的事情。