拉格朗日插值

拉格朗日插值

任给定 \(F\) 中 \(2n+2\) 个数 \(x_1，x_2，…，x_{n+1}，y_1，y_2，…，y_{n+1}，\) 其中 \(x_1，x_2，…x_{n+1}\) 互不相同，则存在唯一的次数不超过n的多项式 \(f(x)\) ，满足 \(f(x_i)=y_i\) \((i=1，2，…，n+1)\) ，这里

\[f(x)=\sum_{i=1}^{n} y_i\cdot\prod_{j\ne i}^{1\le j\le n}\frac{x-x_j}{x_i-x_j} \]

叫做拉格朗日插值公式。

公式推导

对于一个 \(n-1\) 次多项式，知道 \(n\) 个点
设函数 \(\iota_i(x_i)\) \(i\in [1,n]\)
满足：

1. \(\iota_i\) 必须是 \(n-1\)次多项式
2. \( \iota_i(x_j)=\begin{cases} 1 & i=j \\ 0 & i \ne j \end{cases} \)

于是 \(y_i \iota_i (x)\) 可以保证，在 \(x_i\)点处取值为 \(y_i\)，其余 \(n-1\) 个点处取值为 \(0\)

那么

\[f(x) = \sum_{i=1}^{n} y_i \iota_i (x) \]

容易构造

\[\iota_i(x)=\prod_{j \ne i}^{1 \le j \le n}\frac{x - x_j}{x_i - x_j} \]

代入 \(\iota_i\) ，可得拉格朗日插值公式

\[f(x) = \sum_{i=1}^{n} y_i\cdot\prod_{j\ne i}^{1\le j\le n}\frac{x - x_j}{x_i - x_j} \]

其中 \(\iota_i\) 即为拉格朗日基本多项式，又称插值基函数

il int solve(){
    int as=0;
    for(ri int i=1;i<=n;++i){
        int s=1;
        for(ri int j=1;j<=n;++j){
            if(i==j) continue;
            s=s*(X-x[j])%mod*inv(x[i]-x[j])%mod;
        }
        as=(as+s*y[i]%mod)%mod;
    }
    return (as+mod)%mod;
}

code

#include<bits/stdc++.h>
#define il inline
#define cs const
#define ri register
#define int long long 
using namespace std;
namespace Q{
    il int rd(){
        ri int x=0,f=1;ri char c=getchar();
        while(c<'0'||c>'9') f=(c=='-')?-1:1,c=getchar();
        while(c>='0'&&c<='9') x=x*10+(c^48),c=getchar();
        return x*f; 
    }
    il void wt(int x){
        if(x<0) x=-x,putchar('-');
        if(x>=10) wt(x/10);
        return putchar(x%10+48),void();
    }
}
cs int mod=998244353,N=2e3+5;
int n,k,x[N],y[N];
il int qpow(int a,int b){
    int as=1;
    while(b){
        if(b&1) as=as*a%mod;
        a=a*a%mod,b>>=1;
    }
    return as;
}
signed main(){
    n=Q::rd(),k=Q::rd();
    for(ri int i=1;i<=n;++i){
        x[i]=Q::rd(),y[i]=Q::rd();
    }
    int as=0;
    for(ri int i=1;i<=n;++i){
        int s=1;
        for(ri int j=1;j<=n;++j){
            if(i==j) continue;
            s=s*(k-x[j])%mod*qpow(x[i]-x[j],mod-2)%mod;
        }
        as=(as+s*y[i]%mod+mod)%mod;
    }
    Q::wt(as);
    return 0;
}

横坐标是连续整数

在 \(x\) 连续时，可以做到 \(O(n)\) 插值

\[\begin{align*} f(x) &= \sum_{i=1}^{n} y_i\cdot\prod_{j\ne i}^{1\le j\le n}\frac{x - x_j}{x_i - x_j} \\ &= \sum_{i=1}^{n} y_i\cdot\prod_{j\ne i}^{1\le j\le n} \frac{x - j}{i - j} \end{align*} \]

考虑快速计算 \(\prod_{j\ne i}^{1\le j\le n} \frac{x - j}{i - j}\)

对于分子，维护关于 \(x\) 的前后缀积

\[pre_i = \prod_{j=1}^{i} x-j \space , \space nxt_i = \prod_{j=i}^{n} x-j \]

对于分母，预处理阶乘
则对于每个 \(i\) ，分母为

\[(i-1)! \cdot (-1)^{k+2-i} \cdot (k+2-i)! \]

即

\[\iota_i = (-1)^{k+2-i} \cdot \frac{pre_i \cdot nxt_i}{(i-1)! \cdot (k+2-i)!} \]

\(O(n)\) 预处理，对于每个 \(i\) ，可以 \(O(1)\) 得到 \(\iota_i\)

il void init(cs int k,cs int n){
    per[0]=nxt[k+3]=fct[0]=s[0]=1;
    for(ri int i=1;i<=k+2;++i){
        per[i]=1ll*per[i-1]*(n-i)%mod;
        fct[i]=1ll*fct[i-1]*i%mod;
    }
    for(ri int i=k+2;i;--i){
        nxt[i]=1ll*nxt[i+1]*(n-i)%mod;
    }
    for(ri int i=1;i<=k+2;++i){
        ss[i]=((k-i)&1)?(mod-fct[k+2-i]):fct[k+2-i];
        ss[i]=1ll*ss[i]*fct[i-1]%mod;
        s[i]=1ll*s[i-1]*ss[i]%mod;
    }
    inv[k+2]=qpow(s[k+2]);
    for(ri int i=k+1;~i;--i){
        inv[i]=1ll*inv[i+1]*ss[i+1]%mod;
        ss[i+1]=1ll*s[i]*inv[i+1]%mod;
        s[i+1]=1ll*per[i]*nxt[i+2]%mod;
    }
    return;
}

点击查看代码

#include<bits/stdc++.h>
#define il inline
#define cs const
#define ri register
using namespace std;
cs int mod=1e9+7,N=1e6+5;
int fct[N],per[N],nxt[N],inv[N],s[N],ss[N];
il int qpow(int a,int b=mod-2,int p=mod){
    ri int as=1;
    while(b>0){
        if(b&1) as=1ll*as*a%p;
        a=1ll*a*a%p,b>>=1;
    }
    return as;
}
il void init(cs int k,cs int n){
    per[0]=nxt[k+3]=fct[0]=s[0]=1;
    for(ri int i=1;i<=k+2;++i){
        per[i]=1ll*per[i-1]*(n-i)%mod;
        fct[i]=1ll*fct[i-1]*i%mod;
    }
    for(ri int i=k+2;i;--i){
        nxt[i]=1ll*nxt[i+1]*(n-i)%mod;
    }
    for(ri int i=1;i<=k+2;++i){
        ss[i]=((k-i)&1)?(mod-fct[k+2-i]):fct[k+2-i];
        ss[i]=1ll*ss[i]*fct[i-1]%mod;
        s[i]=1ll*s[i-1]*ss[i]%mod;
    }
    inv[k+2]=qpow(s[k+2]);
    for(ri int i=k+1;~i;--i){
        inv[i]=1ll*inv[i+1]*ss[i+1]%mod;
        ss[i+1]=1ll*s[i]*inv[i+1]%mod;
        s[i+1]=1ll*per[i]*nxt[i+2]%mod;
    }
    return;
}
signed main(){
    int n,k,as=0;cin>>n>>k;init(k,n);
    for(ri int i=1,y=0,a,b;i<=k+2;++i){
        y=(y+qpow(i,k))%mod;
        as=(as+1ll*y*s[i]%mod*ss[i]%mod)%mod;//1ll !!!!!
        // a=1ll*per[i-1]*nxt[i+1]%mod;
        // b=1ll*fct[i-1]*((k-i)&1?-1ll:1ll)*fct[k+2-i]%mod;
        // as=(as+1ll*y*a%mod*qpow(b)%mod)%mod;
    }
    cout<<(as%mod+mod)%mod;
    return 0;
}

重心拉格朗日插值法

对于式子

\[f(x) = \sum_{i=1}^{n} y_i\cdot\prod_{j\ne i}^{1\le j\le n}\frac{x - x_j}{x_i - x_j} \]

设 \( g=\prod_{i=1}^{n} x-x_i \)

\[f(x) = g \cdot \sum_{i=1}^{n} \prod_{j\ne i}^{1\le j\le n} \frac{y_i}{(x-x_i)(x_i-x_j)} \]

设 \( t_i=\frac{y_i}{\prod_{j\ne i}^{1\le j\le n} x_i-x_j} \)

\[f(x) =g\cdot \sum_{i=1}^n\frac{t_i}{x-x_i} \]

这样每次新加入一个点的时候只需要 \(O(n)\) 计算它的 \(t_i\) 即可

这样就可以 \(O(n^2)\) 计算 \(n\) 个点的插值

edit