拉格朗日插值法

大橘猫喵～ · 2026-09-21 19:11:36 · 算法·理论

-2. 引言

相信我们小时候都做过找规律题目，比如 1,3,5,7\dots。

相信我们肯定脱口而出下一个肯定是 9！

但并不一定，下一个可以是任何实数，比如 33。

-1. 我会瞪眼！

比如一个数列：1,3,5,7,33。

发现可以构造为 P_x=(x-1)(x-2)(x-3)(x-4)+2x-1。

此时 P_1=1,P_2=3,P_3=5,P_4=7,P_5=33。

0. 寻找规律

我们发现 (x-1) 在 x=1 时为 0，这一项 (x-1)(x-2)(x-3)(x-4) 就全是 0 了。

所以 1\le x\le 4 时，满足这个条件。

我们发现可以专门构造一些 (x-x_i) 来使某些项 =0。

1. 正题：拉格朗日插值法

给定 n 个点：(x_1,y_1),(x_2,y_2),(x_3,y_3),\dots,(x_n,y_n)。

其中 x_i 互不相同，我们尝试构造：

首先构造基函数 L_i(x)，L_i(x)=\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}。

可能有些同学看不懂这个式子，展开来就是（这里写 n=4 的情况）：

L_1(x)=\frac{(x-x_2)(x-x_3)(x-x_4)}{(x_1-x_2)(x_1-x_3)(x_1-x_4)}

L_2(x)=\frac{(x-x_1)(x-x_3)(x-x_4)}{(x_2-x_1)(x_2-x_3)(x_2-x_4)}

L_3(x)=\frac{(x-x_1)(x-x_2)(x-x_4)}{(x_3-x_1)(x_3-x_2)(x_3-x_4)}

L_4(x)=\frac{(x-x_1)(x-x_2)(x-x_3)}{(x_4-x_1)(x_4-x_2)(x_4-x_3)}

那这个函数有什么用呢，他有如下性质：

当 i=j 时，L_i(x_j)=1。

当 i\ne j 时，L_i(x_j)=0。

然后答案就是：

P(x)=\sum^{n}_{i=1}y_iL_i(x)

因为当 x=x_j 时，只有 L_j(x_j)=1，其他都是 0，所以 P(x_j)=y_j。

2. 如何用 C++ 实现呢

给出一道例题

先看代码：

::::success[代码]

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int mod=998244353;
int n,k,ans,x[2005],y[2005];
int pw(int a,int b){ // 快速幂
    int res=1;
    while(b){
        if(b&1) res=res*a%mod;
        a=a*a%mod,b>>=1;
    }return res;
}int inv(int a){return pw(a,mod-2);} // 逆元
signed main(){
    cin>>n>>k;
    for(int i=1;i<=n;i++) cin>>x[i]>>y[i];
    for(int i=1;i<=n;i++){
        int z=y[i]%mod,m=1;
        for(int j=1;j<=n;j++){
            if(i==j) continue;
            z=z*(k-x[j])%mod; // 分子
            m=m*(x[i]-x[j])%mod; // 分母
        }ans=((ans+z*inv(m)%mod)+mod)%mod; // 分子/分母的和（即为答案）
    }cout<<(ans+mod)%mod;
}

::::

显然，这个代码就是模拟一遍刚刚的流程，时间复杂度是 \mathcal{O}(n^2)。

3. 拓展

连续点值优化

当插值节点的横坐标是连续整数时（例如 x_i=1,2,\dots,n），拉格朗日插值可以大幅加速。

看看这个式子，也就是 P(x)：

P(x)=\sum^{n}_{i=1}y_i\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}

我们发现了什么？分母可以预处理，分子也可以用前缀积和后缀积在 \mathcal{O}(1) 时间内得到！

我们先令 x_i=i（方便之后推式子，反正如果不是就把 i 移个几位就好了），则：

f(k)=\sum_{i=1}^{n} y_{i} \frac{\text { pre}_{i-1} \text{ suf}_{i+1}}{(i-1)!(-1)^{n-i}(n-i)!}

我们来仔细看一下这个式子：

首先他本质上是把 \prod_{j\ne i} 优化掉，另外 \operatorname{pre}_{i}=\prod_{j=1}^{i}(k-j),\operatorname{suf}_{i}=\prod_{j=i}^{n}(k-j)
先看分子，是 k-x_j，因为 x 一直不变（在这个 f(k) 里），所以只要前后缀乘一下就可以表示出 (k-1)\times(k-2)\times\dots\times(k-i-1)\times(k-i+1)\times(k-i+2)\times\dots\times(k-n)。
再看分母，观察发现这其实就是阶乘的形式，就是 (i-1)!(n-i)!，但是要调一下符号，就是 (-1)^{n-i}。

重心拉格朗日插值

重心拉格朗日插值将插值多项式写为一种更稳定的形式，适合动态加点的场景。

也就是换一种写法：定义重心权 \omega_i 为 \frac{1}{\prod_{j\ne i}(x_i-x_j)}。

即 \frac{1}{\omega_i} 就是 \prod_{j\ne i}(x_i-x_j)，它出现在 P(x) 的分母中：

P(x)=\sum^{n}_{i=1}y_i \textcolor{red}{\prod_{j\ne i}} \frac{x-x_j}{\textcolor{red}{(x_i-x_j)}}

令 \ell(x)=\prod_{i=1}^n(x-x_i)，这是所有节点构成的多项式。

然后我们试试重写基函数：

L_i(x)=\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}=\frac{\prod_{j\ne i}(x-x_j)}{\prod_{j\ne i}(x_i-x_j)}=\omega_i\prod_{j\ne i}(x-x_j)=\omega_i \frac{\ell(x)}{x-x_i}

把基函数带入最后的多项式：

P(x)=\sum^{n}_{i=1}y_i\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}\Rightarrow P(x)=\ell(x)\sum_{i=1}^n \frac{\omega_i y_i}{x-x_i}

就大功告成了！

第二种形式（再优化！）

考虑用同样的节点插值常数函数 1。因为插值多项式唯一，所以：

1=\sum^{n}_{i=1}L_i(x)=\ell(x)\sum_{i=1}^n \frac{\omega_i}{x-x_i}

于是 \ell(x)=\frac{1}{\sum^n_{i=1}\frac{\omega_i}{x-x_i}}。

然后再带回原来的式子，得到：

f(x)=\frac{\sum^n_{i=1} \frac{\omega_i y_i}{x-x_i}}{\sum^n_{i=1} \frac{\omega_i}{x-x_i}}

我们发现他不用算 \ell(x) 了！

求值时只需要计算两个和，每个和都是 O(n) 的时间复杂度。

使用完重心拉插后发现单次求值变成 O(n) 了，只有预处理是 O(n^2)，还可以动态添加节点。

动态添加节点

假设我们现在已经有了 n 个节点 (x_1,y_1),(x_2,y_2),(x_3,y_3),\dots,(x_n,y_n)，往里面添加新点 (x_{n+1},y_{n+1})。

然后更新重心权，旧点：

\omega_i^′=\frac{1}{\prod_{j\ne i,j\le n+1}(x_i-x_j)}=\frac{\omega_i}{x_i-x_{n+1}}

显然，更新一个旧点的时间复杂度是 O(1)。

新点（n+1）：

\omega_{n+1}=\frac{1}{\prod^n_{j=1}(x_{n+1}-x_j)}

因为只有一个点，所以哪怕处理是 O(n) 的也没事。

之后求值仍然是 O(n)。这比重新计算所有 \omega_i 的 O(n^2) 快得多。

总结一下：

重心拉插的公式（第二种形式）：

f(x)=\frac{\sum^n_{i=1} \frac{\omega_i y_i}{x-x_i}}{\sum^n_{i=1} \frac{\omega_i}{x-x_i}}

其中 \omega_i=\frac{1}{\prod_{j\ne i}(x_i-x_j)}。

End.