[浅谈] 拉格朗日插值

Introduce

给定 $n$ 个点，那么可以确定一个不超过 $n-1$ 项的多项式函数值。我们可以使用高斯消元，但是 $O(n^3)$ 的时间复杂度和精度误差难以接受。

Principle

我们考虑构造函数 $fi$ ，满足其在 $x=x_i$ 时函数值为 $1$ ，在 $x=x_j(j\ne i,j\in [1,n])$ 是 $0$，这很好构造： $fi(x)=\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}$ 。

然后再构造 $f(x)=\sum _{i=1}^n{y}_{i}fi(x)$ ，那么这个函数就经过这 $n$ 个点。

$f(x)=\sum _{i=1}^n{y}_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}$

$\text{P4781 【模板】拉格朗日插值}$

给定 $n$ 点确定多项式，求 $f(k)$ ，直接代入上面的公式。

点击查看代码

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=2e3+110,mod=998244353;
int read(){
	int x=0,f=1;char c=getchar();
	while(c>'9' || c<'0'){if(c=='-')f=-1;c=getchar();}
	while(c>='0' && c<='9'){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
	return x*f;
}
int n,k,x[N],y[N],ans;
int ksm(int x,int y){
	int sum=1;
	while(y){
		if(y&1)sum=(sum*x)%mod;
		y>>=1;
		x=(x*x)%mod;
	}
	return sum%mod;
}
signed main(){
	n=read(),k=read();
	for(int i=1;i<=n;i++)x[i]=read(),y[i]=read();
	for(int i=1;i<=n;i++){
		int tmp=y[i];
		for(int j=1;j<=n;j++)if(j!=i)tmp=tmp%mod*(k-x[j])%mod*ksm(x[i]-x[j],mod-2)%mod;
		ans=(ans+tmp)%mod;
	}
	printf("%lld\n",(ans+mod)%mod);
	return 0;
}

$x$ 连续时优化

如果 $x_i=i$ ：
$f(x)=\sum _{i=1}^n{y}_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}=\sum _{i=1}^n{y}_i\prod _{j\ne i}\frac{x-j}{i-j}=\sum _{i=1}^n{y}_i\frac{pre[i-1]suf[i+1]}{fac[i-1]fac[n-i](-1{)}^{n-i}}$

此时 $pre[i]=\prod _{j=1}^{i}x-j$
$suf[i]=\prod _{j=i}^{n}x-j$

时间复杂度 $O(n)$ 。

可能有锅

#include<bits/stdc++.h>
#define int long long

using namespace std;
const int N=2e6+110,mod=998244353;
int read(){
	int x=0,f=1;char c=getchar();
	while(c>'9' || c<'0'){if(c=='-')f=-1;c=getchar();}
	while(c>='0' && c<='9'){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
	return x*f;
}
int n,k,fac[N],pre[N],suf[N],y[N],ans;
int ksm(int x,int y){
	int sum=1;
	while(y){
		if(y&1)sum=(sum*x)%mod;
		y>>=1;
		x=(x*x)%mod;
	}
	return sum%mod;
}
void pre_opt(){
	fac[0]=pre[0]=suf[n+1]=1;
	for(int i=1;i<=n;i++)fac[i]=fac[i-1]*i%mod;
	for(int i=1;i<=n;i++)pre[i]=(pre[i-1]*(k-i))%mod;//,cout<<pre[i]<<" ";cout<<endl;
	for(int i=n;i>=1;i--)suf[i]=(suf[i+1]*(k-i))%mod;//,cout<<suf[i]<<" ";cout<<endl;
	return;
}
signed main(){
	n=read(),k=read()%mod;
	for(int i=1;i<=n;i++)y[i]=read()%mod;
	pre_opt();
	for(int i=1;i<=n;i++){
		int tmp=y[i];
		tmp=tmp*pre[i-1]%mod*suf[i+1]%mod;
		tmp=tmp*ksm(fac[i-1]*fac[n-i]%mod*ksm(-1,n-i)%mod,mod-2)%mod;
		ans=(ans+tmp)%mod;
	}
	printf("%lld\n",(ans+mod)%mod);
	return 0;
}

$\text{The Sum of the k-th Powers}$

求 $\sum _{i=1}^n{i}^k$ 。$n\le {10}^9,k\le {10}^6$
把这个式子展开其实是个 $k+1$ 项式。那么我们把它看成一个函数 $f(n)$ ，带入 $k+2$ 个点即可。选择连续的点可以做到 $O(n\log {10}^9)$ 。（逆元带来的 $\log$ ）。

最后提醒一句：分数取模不要漏模，多模，减法时注意 $+mod$

点击查看代码

#include<bits/stdc++.h>
#define int long long

using namespace std;
const int N=2e6+110,mod=1e9+7;
int read(){
	int x=0,f=1;char c=getchar();
	while(c>'9' || c<'0'){if(c=='-')f=-1;c=getchar();}
	while(c>='0' && c<='9'){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
	return x*f;
}
int n,k,fac[N],pre[N],suf[N],y,ans;
int ksm(int x,int y){
	int sum=1;
	while(y){
		if(y&1)sum=(sum*x)%mod;
		y>>=1;
		x=(x*x)%mod;
	}
	return sum%mod;
}
void pre_opt(){
	fac[0]=pre[0]=suf[k+3]=1;
	for(int i=1;i<=k+2;i++)fac[i]=fac[i-1]*i%mod;
	for(int i=1;i<=k+2;i++)pre[i]=(pre[i-1]*(n-i))%mod;//,cout<<pre[i]<<" ";cout<<endl;
	for(int i=k+2;i>=1;i--)suf[i]=(suf[i+1]*(n-i))%mod;//,cout<<suf[i]<<" ";cout<<endl;
	return;
}
signed main(){
	n=read(),k=read();
	pre_opt();
	for(int i=1;i<=k+2;i++){
		y=(y+ksm(i,k))%mod;int tmp=y;
		int a=y*pre[i-1]%mod*suf[i+1]%mod;
		int b=fac[i-1]*fac[m+3-i]%mod*ksm(-1,m+3-i)%mod;
		ans=(ans+a*ksm(b,mod-2)%mod+mod)%mod;
	}
	printf("%lld\n",(ans+mod)%mod);
	return 0;
}