【算法竞赛学习笔记】拉格朗日插值-数学提升计划

本文详细介绍了拉格朗日插值法,包括基本原理、公式、优化方法以及在ACM竞赛和多项式计算中的应用。通过拉格朗日插值,可以在O(n^2)的时间复杂度内解决多项式插值问题,同时探讨了连续插值优化、重心拉格朗日插值以及快速插值等技术,提供了具体的例题和代码实现。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >


title : 拉格朗日插值
date : 2022-2-19
tags : ACM,数学,杂项
author : Linno

拉格朗日插值

简介

在平面直角坐标系中,n+1n+1n+1xxx坐标不同的点可以确认唯一的最高次为nnn的多项式。当我们要解决这些点求多项式的问题时,可以使用高斯消元来获得每一项的系数,但复杂度是O(n3)O(n^3)O(n3)的,而且往往会存在精度问题,而拉格朗日插值法可以在O(n2)O(n^2)O(n2)的复杂度内解决这一问题。

公式

假设多项式经过(x1,y1),(x2,y2),(x3,y3)(x_1,y_1),(x_2,y_2),(x_3,y_3)(x1,y1),(x2,y2),(x3,y3)三个点,那么可以假设

f(x)=y1(x−x2)(x−x3)(x1−x2)(x1−x3)+y2(x−x1)(x−x3)(x2−x1)(x2−x3)+y3(x−x1)(x−x2)(x3−x1)(x1−x2)f(x)=y_1\frac{(x-x_2)(x-x_3)}{(x_1-x_2)(x_1-x3)}+y_2\frac{(x-x_1)(x-x_3)}{(x_2-x_1)(x_2-x3)}+y_3\frac{(x-x_1)(x-x_2)}{(x_3-x_1)(x_1-x2)}f(x)=y1(x1x2)(x1x3)(xx2)(xx3)+y2(x2x1)(x2x3)(xx1)(xx3)+y3(x3x1)(x1x2)(xx1)(xx2)

这样我们直接把每个点代入,都只会剩下一项yyy,其正确性可以保证。

给出公式f(k)=∑i=0n∏i≠jk−x[j]x[i]−x[j]f(k)=\sum_{i=0}^n \prod_{i\ne j}\frac{k-x[j]}{x[i]-x[j]}f(k)=i=0ni=jx[i]x[j]kx[j],我们以此求出给定kkk时任一点的拟合yyy值。

优化

①关于x值连续时的优化

xix_ixi的取值是连续的时候,我们可以把式子优化成O(n)O(n)O(n)的复杂度。

xix_ixi换成i∈[0,n]i\in [0,n]i[0,n],新的式子就是f(k)=∑i=0nyi∏i≠jk−ji−jf(k)=\sum_{i=0}^ny_i\prod_{i\ne j}\frac{k-j}{i-j}f(k)=i=0nyii=jijkj

我们维护关于k的前缀积和后缀积prei=∏j=0i(k−j),sufi=∏j=in(k−j)pre_i=\prod_{j=0}^i(k-j),suf_i=\prod_{j=i}^n(k-j)prei=j=0i(kj),sufi=j=in(kj)

那么式子就变成了f(k)=∑i=0nyiprei−1∗sufi+1fac[i]∗fac[n−i]f(k)=\sum_{i=0}^n y_i\frac{pre_{i-1}*suf_{i+1}}{fac[i]*fac[n-i]}f(k)=i=0nyifac[i]fac[ni]prei1sufi+1

不过这一优化适用条件太特殊了,不常用。

②重心拉格朗日

gk=∏i=1nk−xi,ti=yi∏j≠ixi−xjg_k=\prod_{i=1}^nk-x_i,t_i=\frac{y_i}{\prod_{j\ne i}x_i-x_j}gk=i=1nkxi,ti=j=ixixjyi

这样式子就化成了f(k)=gk∑i=0ntik−xif(k)=g_k\sum_{i=0}^n\frac{t_i}{k-x_i}f(k)=gki=0nkxiti

如此一来我们每次加入一个点时只需要记录tit_iti即可。

③快速插值

快速插值法可以在任意情况下做到O(nlog2n)O(nlog^2n)O(nlog2n)的时间复杂度的插值。

对于原公式f(x)=∑i=0nyi∏j≠i(x−xj)∏j≠i(xi−xj)f(x)=\sum_{i=0}^ny_i\frac{\prod_{j\ne i}(x-x_j)}{\prod_{j\ne i}(x_i-x_j)}f(x)=i=0nyij=i(xixj)j=i(xxj),我们考虑如何快速求分子分母。

首先考虑分母:gi=∏j≠i(xi−xj)=lim⁡x→xi∏j=0n(x−xj)x−xi=(∏j=0n(x−xj)′)∣x=xig_i=\prod_{j\ne i}(x_i-x_j)=\lim_{x\to x_i}\frac{\prod_{j=0}^n(x-x_j)}{x-x_i}=(\prod_{j=0}^n(x-x_j)')|_{x=x_i}gi=j=i(xixj)=limxxixxij=0n(xxj)=(j=0n(xxj))x=xi

这样我们可以用分治求出∏j=0n(x−xj)\prod_{j=0}^n(x-x_j)j=0n(xxj),求导后对所有xix_ixi多点求值即可。

求分子是同样的步骤,不过是把某项xix_ixi换成特定的xxx值罢了。

例题

luogu P4781【模板】拉格朗日插值
//#pragma GCC optimize("Ofast", "inline", "-ffast-math")
//#pragma GCC target("avx,sse2,sse3,sse4,mmx")
#include<bits/stdc++.h>
#define inf 0x3f3f3f3f
#define int long long
using namespace std;
const int N=2006;
const int mod=998244353;

int fpow(int a,int b){
	int res=1;
	while(b){
		if(b&1) res=res*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return res%mod;
}

int n,k,x[N],y[N];

int fc(int kk){
	int ans=0;
	for(int i=1;i<=n;i++){
		int fz=y[i],fm=1;
		for(int j=1;j<=n;j++){
			if(j==i) continue;
			fz=fz*(kk-x[j])%mod;
			fm=fm*(x[i]-x[j])%mod;
		}
		ans=(ans+fz*fpow(fm,mod-2)%mod+mod)%mod;
	}
	return ans%mod;
}

signed main(){
	scanf("%lld%lld",&n,&k);
	for(int i=1;i<=n;i++) scanf("%lld%lld",x+i,y+i);
	printf("%lld",fc(k));
	return 0;
}
luoguP4593 [TJOI2018]教科书般的亵渎

题意是求∑i=1nik\sum_{i=1}^n i^ki=1nik,可以看出这是一个关于nnnk+1k+1k+1次多项式,可以插值解决。

n=0...k+1n=0...k+1n=0...k+1的函数值求出来,就可以做到O(k)O(k)O(k)或者O(klogk)O(klogk)O(klogk)n0n_0n0时的函数值。

其他求法:递推法、Bernoulli数、Stirling数、多项式差分

//#pragma GCC optimize("Ofast", "inline", "-ffast-math")
//#pragma GCC target("avx,sse2,sse3,sse4,mmx")
#include<bits/stdc++.h>
#define inf 0x3f3f3f3f
#define int long long
using namespace std;
const int N=66;
const int mod=1e9+7;

inline int read(){	int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9'){if(ch=='-') f=f*-1;ch=getchar();}while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}return x*f;}

int n,m,a[N],tot,inv[3601];


inline int fp(int a,int p){
	int res=1;
	while(p){
		if(p&1) res=res*a%mod;
		a=a*a%mod;
		p>>=1;
	}
	return res;
}

int x[N],y[N],fac[N],ifac[N],pre[N],suf[N];

inline int get(int nn,int mm){
	int lim=mm+1,ans=0;
	memset(y,0,sizeof(y));
	for(int i=1;i<=lim;i++) y[i]=(y[i]+y[i-1]+fp(i,mm))%mod;
	pre[0]=nn;suf[lim+1]=1;
	for(int i=1;i<=lim;i++) pre[i]=pre[i-1]*(nn-i)%mod;
	for(int i=lim;i>=1;i--) suf[i]=suf[i+1]*(nn-i)%mod;
	for(int i=0;i<=lim;i++){
		int up=pre[i-1]*suf[i+1]%mod*y[i]%mod,down=ifac[i]*ifac[lim-i]%mod;
		if((lim-i)&1) down=mod-down;
		ans=(ans+up*down%mod)%mod;
	}
	return ans;
}

inline void solve(){
	n=read();m=read();
	memset(a,0,sizeof(a));
	int ans=0;
	for(int i=1;i<=m;i++) a[i]=read();
	a[++m]=++n;
	sort(a+1,a+1+m);
	for(int i=1;i<=m;i++){
		for(int j=i;j<=m;j++) ans=(ans+get(a[j]-1,m)-get(a[j-1],m)+mod)%mod;
		for(int j=i+1;j<=m;j++) a[j]=(a[j]-a[i]+mod)%mod;
		a[i]=0;
	}
	printf("%lld\n",ans);
}

signed main(){
	inv[1]=1;for(int i=2;i<=3600;i++) inv[i]=(mod-mod/i)*inv[mod%i]%mod;
	fac[0]=1;for(int i=1;i<=60;i++) fac[i]=fac[i-1]*i%mod;
	ifac[60]=fp(fac[60],mod-2); 
	for(int i=60;i>=1;i--) ifac[i-1]=ifac[i]*i%mod;
	int t=read();
	while(t--) 	solve();
	return 0;
}

2019 ICPC 南昌邀请赛-Polynomial

已知f(x)f(x)f(x)是最高项为n次的多项式,给定f(0)...f(n)f(0)...f(n)f(0)...f(n) 多次询问∑i=LRf(i)mod  9999991\sum_{i=L}^Rf(i) \mod 9999991i=LRf(i)mod9999991

连续插值的多项式优化得f[k]=∑i=0nyi∏i≠jk−ji−jf[k]=\sum_{i=0}^ny_i\prod_{i\ne j}\frac{k-j}{i-j}f[k]=i=0nyii=jijkj,对于分母预处理阶乘和阶乘逆元,对于分子做前缀积和后缀积,化简得f[k]=∑i=0nyipre[i−1]∗suf[i+1]fac[i]∗fac[n−i][(n−i)&1?−1:1]f[k]=\sum_{i=0}^ny_i\frac{pre[i-1]*suf[i+1]}{fac[i]*fac[n-i]}[(n-i)\&1?-1:1]f[k]=i=0nyifac[i]fac[ni]pre[i1]suf[i+1][(ni)&1?1:1]

#include<bits/stdc++.h>
#define int long long
using namespace std;

const int N=1005;
const int MAXN=1e7+10;
const int mod=9999991;

int t,n,q,l,r;
int F[N],pre[N],suf[N],fac[N],ifac[N],sum[MAXN];

int fpow(int a, int b){
    int res=1;
    while(b){
        if(b&1) res=res*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return res;
}

void init(){
    fac[0]=1;
    for(int i=1;i<N;i++) fac[i]=fac[i-1]*i%mod;
    ifac[N-1]=fpow(fac[N-1],mod-2);
    for(int i=N-1;i>=1;i--) ifac[i-1]=ifac[i]*i%mod;
}

int cal(int *f, int k, int n){
    if(k<=n) return f[k];
    pre[0]=suf[n]=1;
    for(int i=1;i<=n;i++) pre[i]=pre[i-1]*(k-i+1)%mod;
    for(int i=n;i>=1;i--) suf[i-1]=suf[i]*(k-i)%mod;
    int ans=0;
    for(int i=0;i<=n;i++){
        int opt=(n-i)&1?-1:1;
        ans=(ans+opt*pre[i]%mod*suf[i]%mod*ifac[i]%mod*ifac[n-i]%mod*f[i]%mod+mod)%mod;
    }
    return f[k]=ans;
}


signed main(){
    init();
    scanf("%lld",&t);
    while(t--){
    	scanf("%lld%lld",&n,&q);
        for(int i=0;i<=n;i++) cin>>F[i];
        F[n+1]=cal(F,n+1,n);
        sum[0]=F[0];
        for(int i=1;i<=n+1;i++) sum[i]=(sum[i-1]+F[i])%mod;
        while(q--){
        	scanf("%lld%lld",&l,&r);
        	int ans=(cal(sum,r,n+1)-cal(sum,l-1,n+1)+mod)%mod;
        	printf("%lld\n",ans);
        }
    }
    return 0;
}

参考资料

https://www.cnblogs.com/zwfymqz/p/10063039.html

https://blog.youkuaiyun.com/Code92007/article/details/94412729

https://www.cnblogs.com/ywwyww/p/8511505.html

https://www.cnblogs.com/dcdcbigbig/p/9715638.html

https://blog.youkuaiyun.com/fztsilly/article/details/109529465

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

RWLinno

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值