拉格朗日插值法

考试中 插值往往很重要 如一个多项式很难求的时候 我们经常带入数值生成点值 最终进行插值插出来。

不过GAUSS消元的复杂度太高 一般采用拉格朗日插值法。

n+1个x坐标不同的点可以确定唯一的最高为n次的多项式。

设该多项式为f(x) 第i个点的坐标为(xi,yi)

如果我们要在k点处取值 那么 f(k)=\(\sum_{i=0}^{n}y_i\Pi_{i\neq j} \frac{k-x_j}{x_i-x_j}\)

证明的话 可以带入一下x_i看一下 其余的 可以自己感性理解。

LINK：模板题

值得一提的是 x切记要两两不同 不然插出来的式子就不一定是一个n次多项式了。

const ll MAXN=200005;
ll n,k;
ll x[MAXN],y[MAXN];
inline ll ksm(ll b,ll p){ll cnt=1;while(p){if(p&1)cnt=cnt*b%mod;b=b*b%mod;p=p>>1;}return cnt;}
signed main()
{
	freopen("1.in","r",stdin);
	get(n);get(k);
	rep(1,n,i)get(x[i]),get(y[i]);
	ll ans=0;
	rep(1,n,i)
	{
		ll w1=1;
		ll w2=1;
		rep(1,n,j)if(i!=j)w1=w1*(k-x[j])%mod,w2=w2*(x[i]-x[j])%mod;
		w2=ksm(w2,mod-2);
		ans=(ans+w2*w1%mod*y[i])%mod;
	}
	putl((ans+mod)%mod);
	return 0;
}

但是一般我们自己插值出来一个多项式的时候 可以取一段连续的点值 然后 预处理一些东西 就可以在nlogn的时间内求出某个点的点值。

算法瓶颈在于求逆元。

LINK：教科书般的亵渎

仔细观察 可以发现需要求 m+1次值。但是k是固定的。k=m+1.

可以发现每次求得都是形似\(\sum_{i=1}^ni^k\)的和。

人类智慧告诉我们i^k这个形式 是以i为自变量 k+1次多项式。

直接求和 插值出来即可。

注意要减掉不合法方案 当然这里是我们自己插值 所以可以采取nlogn的做法。

值得注意的是 \(sum_{i=1}^n i^k\)这个东西是从1开始的 所以不能从0开始插值。

我们的k为m+1 次数为m+2 需要求出来的项的个数为m+3 不要弄混了。

const ll MAXN=60;
ll T,n,m;
ll fac[MAXN],pre[MAXN],suf[MAXN];
ll ans,a[MAXN],f[MAXN],inv[MAXN],inv2;
inline ll ksm(ll b,ll p){ll cnt=1;while(p){if(p&1)cnt=cnt*b%mod;b=b*b%mod;p=p>>1;}return cnt;}
inline ll calc(ll n)//求\sum{i=0}^n i^{m+1}
{
	ll k=n%mod;
	pre[0]=1;suf[m+4]=1;
	rep(1,m+3,i)pre[i]=pre[i-1]*(k-i)%mod;
	fep(m+3,1,i)suf[i]=suf[i+1]*(k-i)%mod;
	ll cnt=0;
	rep(1,m+3,i)
	{
		ll w1=pre[i-1];
		ll w2=suf[i+1];
		ll op=((m-i+3)&1)?-1:1;
		cnt=(cnt+op*f[i]*w1%mod*w2%mod*inv[i-1]%mod*inv[m+3-i]%mod)%mod;
	}
	return (cnt+mod)%mod;
}
signed main()
{
	freopen("1.in","r",stdin);
	get(T);fac[0]=1;inv2=ksm(2,mod-2);
	rep(1,55,i)fac[i]=fac[i-1]*i%mod;
	inv[55]=ksm(fac[55],mod-2);
	fep(54,0,i)inv[i]=(inv[i+1]*(i+1))%mod;
	while(T--)
	{
		get(n);get(m);
		rep(1,m,i)get(a[i]);
		sort(a+1,a+1+m);
		rep(1,m+3,i)f[i]=(f[i-1]+ksm(i,m+1))%mod;
		ans=calc(n);
		rep(1,m,i)ans=(ans-ksm(a[i],m+1))%mod;
		rep(1,m,i)ans=(ans+calc(n-a[i]))%mod;
		rep(1,m,i)fep(i-1,1,j)ans=(ans-ksm(a[i]-a[j],m+1))%mod;
		putl((ans+mod)%mod);
	}
	return 0;
}