拉格朗日插值学习笔记

拉格朗日插值学习笔记

用途

\(n+1\) 个横坐标互不相同的点可以唯一确定一个最高次项为 \(n+1\) 的多项式。

拉格朗日插值法可以在 \(n^2\) 的复杂度内根据这 \(n+1\) 个点求出多项式在某一点处的取值。

构造

\(f(k) = \sum_{i = 0}^{n} y_i \prod_{i \not = j} \frac{k - x_j}{x_i - x_j}\)。

如果我们把 \(k=x_p\) 带入的话，

当 \(i \neq p\) 时，后面连乘的部分会有一项的分母为 \(0\)，所以 \(y_i\) 不会作出贡献。

当 \(i=p\) 时，后面的部分分子和分母都相等，贡献的系数为 \(1\)。

所以得到的结果就是我们想要的 \(y_p\)。

如果题目中给出的点的横坐标是连续的，那么我们可以通过预处理一些前缀积后缀积来 \(O(1)\) 计算后面的部分，复杂度可以做到 \(O(n)\)。

例题

P4781 【模板】拉格朗日插值

模板题，\(n^2\) 插值即可。

代码

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<iostream>
#include<cmath>
#define rg register
template<typename T>void read(rg T& x){
	x=0;rg int fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	x*=fh;
}
const int maxn=2e3+5,mod=998244353;
inline int addmod(rg int now1,rg int now2){
	return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int delmod(rg int now1,rg int now2){
	return now1-=now2,now1<0?now1+mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
	return now1*=now2,now1>=mod?now1%mod:now1;
}
inline int ksm(rg int ds,rg int zs){
	rg int nans=1;
	while(zs){
		if(zs&1) nans=mulmod(nans,ds);
		ds=mulmod(ds,ds);
		zs>>=1;
	}
	return nans;
}
int n,k,x[maxn],y[maxn];
int main(){
	read(n),read(k);
	for(rg int i=1;i<=n;i++){
		read(x[i]),read(y[i]);
	}
	rg int nans=0,sum=0;
	for(rg int i=1;i<=n;i++){
		sum=1;
		for(rg int j=1;j<=n;j++){
			if(i==j) continue;
			sum=mulmod(sum,mulmod(delmod(k,x[j]),ksm(delmod(x[i],x[j]),mod-2)));
		}
		nans=addmod(nans,mulmod(y[i],sum));
	}
	printf("%d\n",nans);
	return 0;
}

CF622F The Sum of the k-th Powers

求 \(\sum_{i=1}^ni^k \bmod (10^9+7),1 \leq n \leq 10^{9},0 \leq k \leq 10^{6}\)。

\(\sum_{i=1}^ni^k\) 是一个 \(k+1\) 次多项式，预处理出前 \(k+2\) 项的值可以 \(O(n)\) 插值。

代码

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<iostream>
#include<cmath>
#define rg register
template<typename T>void read(rg T& x){
	x=0;rg int fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	x*=fh;
}
const int maxn=1e6+5,mod=1e9+7;
inline int addmod(rg int now1,rg int now2){
	return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int delmod(rg int now1,rg int now2){
	return now1-=now2,now1<0?now1+mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
	return now1*=now2,now1>=mod?now1%mod:now1;
}
inline int ksm(rg int ds,rg int zs){
	rg int nans=1;
	while(zs){
		if(zs&1) nans=mulmod(nans,ds);
		ds=mulmod(ds,ds);
		zs>>=1;
	}
	return nans;
}
int ny[maxn],jcc[maxn],rjcc[maxn],pre[maxn],suf[maxn],n,k,f[maxn],ans;
int main(){
	read(n),read(k);
	for(rg int i=1;i<=k+2;i++){
		f[i]=addmod(f[i-1],ksm(i,k));
	}
	if(n<=k+2){
		printf("%d\n",f[n]);
	} else {
		ny[1]=1;
		for(rg int i=2;i<maxn;i++) ny[i]=mulmod(mod-mod/i,ny[mod%i]);
		jcc[0]=1;
		for(rg int i=1;i<maxn;i++) jcc[i]=mulmod(jcc[i-1],ny[i]);
		rjcc[0]=1;
		for(rg int i=1;i<maxn;i++) rjcc[i]=mulmod(rjcc[i-1],mod-ny[i]);
		pre[0]=1;
		for(rg int i=1;i<=k+2;i++) pre[i]=mulmod(pre[i-1],n-i);
		suf[k+3]=1;
		for(rg int i=k+2;i>=1;i--) suf[i]=mulmod(suf[i+1],n-i);
		for(rg int i=1;i<=k+2;i++){
			rg int nans=mulmod(jcc[i-1],mulmod(rjcc[k+2-i],mulmod(pre[i-1],suf[i+1])));
			ans=addmod(ans,mulmod(nans,f[i]));
		}
		printf("%d\n",ans);
	}
	return 0;
}