拉格朗日插值学习笔记

拉格朗日插值学习笔记

应用

众所周知,在平面直角坐标系中,对于任意的 n 个点,都一定有一个不超过 n1 次的函数与之相对应。拉格朗日插值适用于求解这 n 个点对应的函数。

思路

考虑给定的 n 个点的坐标表示为 (xi,yi),不难构造出如下函数:

f(x)=i=1nyigi(x)

那么此时只需要构造出符合要求的 gi 即可。不难发现,为了使 f 符合条件,gi 应该满足:

nx,gi(n)={1n=xi0nxi

其中,x 表示所有 xi 构成的集合。不难发现,当 nxi 时函数值为 0,说明该函数 gi(x) 一定有因式 (xxj)(ij),因此不难看出 gi(x) 应当有因式 ji(xxj)。又因为当 n=xi 时,函数值为 1,不难构造出:

gi(n)=ji(nxj)ji(xixj)

如此就构造出了符合条件的函数,此时可以 O(n2) 求解。在单次求解的时候可以采用。### O(n2) 写法

代码

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

const int N=2e3+5,M=998244353;

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

ll QuickPow(ll a,int b){
    ll res=1;
    while(b>0){
        if(b&1)res=res*a%M;
        a=a*a%M;
        b>>=1;
    }
    return res;
}

ll Lagrange(int x[],int y[],int X){
    ll res=0;
    for(int i=1;i<=n;i++){
        ll l_i=1;
        for(int j=1;j<=n;j++){
            if(i!=j)l_i=l_i*(X-x[j])%M*QuickPow(x[i]-x[j],M-2)%M;
        }
        res=(res+l_i*y[i]%M)%M;
    }
    return res;
}

int main(){
    scanf("%d%d",&n,&k);
    for(int i=1;i<=n;i++)scanf("%d%d",&x[i],&y[i]);
    printf("%lld",(Lagrange(x,y,k)+M)%M);
    return 0;
}

那么可能有人就要说了,我要是需要多次求解,那我不就似了吗,没错,你似了。所以不妨进行进一步转换,令 ti=yiji(xixj),不难看出,f(x)=i=1ntiji(xxj)。不妨计算出 i=1n(xxi) 后将对应项除掉,接着利用多项式乘法算出每一次项前的系数后,就可以按照多项式加法的方式 O(n2) 预处理,O(n) 求解。这样的方式适用于多组查询或需要求出函数对应项系数的情况。

代码

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

const int N=2e3+5,M=998244353;

int n,k;
ll t[N],x[N],y[N],g[N],fs[N],f[N];

ll QuickPow(ll a,int b){
    ll res=1;
    while(b>0){
        if(b&1)res=res*a%M;
        a=a*a%M;
        b>>=1;
    }
    return res;
}

void Init(){
    for(int i=0;i<n;i++){
        t[i]=1;
        for(int j=0;j<n;j++){
            if(i==j)continue;
            t[i]=t[i]*(x[i]-x[j])%M;
        }
        t[i]=QuickPow(t[i],M-2)*y[i]%M;
    }
    fs[0]=1;
    for(int i=0;i<n;i++){
        for(int j=n-1;j>0;j--)fs[j]=(fs[j-1]+fs[j]*(M-x[i])%M)%M;
        fs[0]=fs[0]*(M-x[i])%M;
    }
    for(int i=0;i<n;i++){
        ll inv=QuickPow(M-x[i],M-2);
        g[0]=fs[0]*inv%M;
        for(int j=1;j<n;j++)g[j]=(fs[j]-g[j-1])*inv%M;
        for(int j=0;j<n;j++)f[j]=(f[j]+t[i]*g[j]%M)%M;
    }
    return ;
}

ll Lagrange(ll X){
    ll res=0,Pow=1;
    for(int i=0;i<n;i++){
        res=(res+Pow*f[i]%M)%M;
        Pow=Pow*X%M;
    }
    return res;
}

int main(){
    scanf("%d%d",&n,&k);
    for(int i=0;i<n;i++)scanf("%lld%lld",&x[i],&y[i]);
    Init();
    printf("%lld",(Lagrange(k)+M)%M);
    return 0;
}

但也许题目中给出的 n 不接受 n2 的预处理,那么你可以在题目不强制给定点坐标的题里大展身手。更具体的,我们可以利用一些特殊值,即尝试代入连续的点值,那么此时 xi 是连续的,分母就变成了类似 i=1j1ii=1n+1j(i) 的形式,这种情况只需要预处理阶乘就可以 O(1) 求解,分子也可以通过 O(n) 预处理前缀和后缀积达到 O(1) 求解的目的,那么单次求解的复杂度是 O(n)。但是这种求解方式常数较大,因为前缀后缀积也需要在每次求解时重新处理,并且要求能够快速求出点值对应的函数值,所以推荐在复杂度允许的情况下选择第二种做法。下方代码已求解函数 f(x)=i=1xik 为例。

代码

void Init(){
	f[0]=1;
	pw[0]=0;
	for(int i=1;i<=k+2;i++){
		f[i]=f[i-1]*i%P;
		pw[i]=(pw[i-1]+QuickPow(i,k,P))%P;
	}
    return ;
}

ll Lagrange(ll x,const int M){
	pre[0]=nxt[k+3]=1;
	for(int i=1;i<=k+2;i++)pre[i]=pre[i-1]*(x-i)%M;
	for(int i=k+2;i>=1;i--)nxt[i]=nxt[i+1]*(x-i)%M;
	ll res=0,y=0;
	for(int i=1;i<=k+2;i++)res=(res+(((k-i)&1)?-1:1)*pw[i]*pre[i-1]%M*nxt[i+1]%M*QuickPow(f[i-1]*f[k+2-i],M-2,M))%M;
	return res;
}
posted @   DycIsMyName  阅读(32)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· [翻译] 为什么 Tracebit 用 C# 开发
· 腾讯ima接入deepseek-r1,借用别人脑子用用成真了~
· Deepseek官网太卡,教你白嫖阿里云的Deepseek-R1满血版
· DeepSeek崛起:程序员“饭碗”被抢,还是职业进化新起点?
· RFID实践——.NET IoT程序读取高频RFID卡/标签
点击右上角即可分享
微信分享提示