拉格朗日插值学习笔记

拉格朗日插值学习笔记

应用

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

思路

考虑给定的 \(n\) 个点的坐标表示为 \((x_i,y_i)\)，不难构造出如下函数：

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

那么此时只需要构造出符合要求的 \(g_i\) 即可。不难发现，为了使 \(f\) 符合条件，\(g_i\) 应该满足：

\[\forall n\in x,g_i(n)=\left\{\begin{matrix}1&n=x_i\\0&n\ne x_i\end{matrix}\right. \]

其中，\(x\) 表示所有 \(x_i\) 构成的集合。不难发现，当 \(n\ne x_i\) 时函数值为 \(0\)，说明该函数 \(g_i(x)\) 一定有因式 \((x-x_j)(i\ne j)\)，因此不难看出 \(g_i(x)\) 应当有因式 \(\prod_{j\ne i}(x-x_j)\)。又因为当 \(n=x_i\) 时，函数值为 \(0\)，不难构造出：

\[g_i(n)=\frac{\prod_{j\ne i}(n-x_j)}{\prod_{j\ne i}(x_i-x_j)} \]

如此就构造出了符合条件的函数，此时可以 \(O(n^2)\) 求解。在单次求解的时候可以采用。

不妨进行进一步转换，令 \(t_i=\frac{y_i}{\prod_{j\ne i}(x_i-x_j)}\)，不难看出，\(f(x)=t_i\prod_{j\ne i}(x-x_j)\)。不妨计算出 \(\prod_{i=1}^{n}(x-x_i)\) 后将对应项除掉，接着利用多项式乘法算出每一次项前的系数后，就可以按照多项式加法的方式 \(O(n)\) 求解。这样的方式适用于多组查询或需要求出函数对应项系数的情况。

代码

\(O(n^2)\) 写法

#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;
}

\(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;
}