关于拉格朗日插值

基于拉格朗日插值的多项式求和

拉格朗日插值. 假设 \(P\) 是一个度数为 \(D\) 的多项式。如果我们知道 \(P\) 在 \(D+1\) 个不同位置的取值 \((x_1,P(x_1)),\dots, (x_D,P(x_D))\)，那么我们可以唯一地把 \(P\) 写为

\[P(x) = \sum_{i=1}^{D+1} P(x_i) \frac{\prod_{j\ne i} (x-x_j)}{\prod_{j\ne i}(x_i-x_j)}. \]

link

点击查看代码

#include<bits/stdc++.h>
using namespace std;
int n,k;
long long x[2005],y[2005];
const long long md=998244353;
inline long long pwr(long long x,long long y){
	long long res=1;
	while(y){
		if(y&1)res=res*x%md;
		x=x*x%md;y>>=1;
	}return res;
}
inline long long Lagrange(int t){
	long long res=0;
	for(int i=1;i<=n;i++){
		long long tmp=y[i];
		for(int j=1;j<=n;j++){
			if(i!=j)tmp=tmp*(x[j]-t)%md*pwr(x[j]-x[i],md-2)%md;
		}res=(res+tmp)%md;
	}return (res+md)%md;
}
int main(){
	scanf("%d%d",&n,&k);
	for(int i=1;i<=n;i++)scanf("%lld%lld",&x[i],&y[i]);
	printf("%lld",Lagrange(k));

	return 0;
}

多项式的部分和. 假设 \(P:\mathbb{N}\to \mathbb{R}\) 可以表示为一个度数为 \(D\) 的多项式，那么我们考虑函数 \(Q:\mathbb{N}\to \mathbb{R}\)

\[Q(n) = \sum_{i\le n} P(i). \]

可以证明 \(Q\) 可以表示为一个度数为 \(D+1\) 的多项式，并且这个表示方法是唯一的。

如果已知 \(Q\) 在 \(0,1,\dots, D\) 这些输入上的输出，那么对于任意给定的数 \(n\)，可以在 \(O(D)\) 时间内计算出 \(Q(n)\)。算法：按照拉格朗日插值公式，只需要对每个 \(i\in \{0,\dots, D\}\) 计算 \(\prod_{j\ne i} (n - j)\) 和 \(\prod_{j\ne i} (i - j)\)。对于前者，我们只需要计算一个数列 \((n-0),(n-1),\dots, (n-D)\)，然后取一个前缀积和一个后缀积即可。对于后者，可以发现它就是 \((i-1)!\cdot (D-i)!\cdot (-1)^{D-i}\)。因此，整个计算可以在 \(O(D)\) 时间内完成。

题1. 给定 \(N,K\)，计算 \(\sum_{i=0}^{N} i^K\). 保证 \(K\le 10^5\).

题1 加强版. 给定 \(N,K,q\)，计算 \(\sum_{i=0}^N i^K q^i\). 保证 \(K\le 10^5\).

link

点击查看代码

#include<bits/stdc++.h>
using namespace std;
int n,k;
int fac[1000005],inv[1000005];
const int md=1e9+7;
inline void init(){
	fac[0]=fac[1]=inv[0]=inv[1]=1;
	for(int i=2;i<=k+1;i++)fac[i]=1ll*fac[i-1]*i%md;
	for(int i=2;i<=k+1;i++)inv[i]=1ll*(md-md/i)*inv[md%i]%md;
	for(int i=2;i<=k+1;i++)inv[i]=1ll*inv[i]*inv[i-1]%md;
}
int f[1000005];
inline int pwr(int x,int y){
	int res=1;
	while(y){
		if(y&1)res=1ll*res*x%md;
		x=1ll*x*x%md;y>>=1;
	}return res;
}
int pre[1000005],suf[1000005];
inline int Lagrange(int x){
	int res=0;
	for(int i=1;i<=k+2;i++)pre[i]=1ll*pre[i-1]*(x-i)%md;
	for(int i=k+2;i>0;i--)suf[i]=1ll*suf[i+1]*(x-i)%md;
	for(int i=1;i<=k+2;i++){
		int tmp=((k-i)&1?-1ll:1ll)*inv[i-1]*inv[k+2-i]%md;
		res=(res+1ll*f[i]*pre[i-1]%md*suf[i+1]%md*tmp)%md;
	}return (res+md)%md;
}
int main(){
	scanf("%d%d",&n,&k);init();pre[0]=suf[k+3]=1;
	for(int i=1;i<=k+2;i++)f[i]=(f[i-1]+pwr(i,k))%md;
	printf("%d",Lagrange(n));

	return 0;
}