拉格朗日插值法总结模板(1~n)

拉格朗日插值法

---

简介

在数值分析中，拉格朗日插值法是以法国十八世纪数学家约瑟夫·路易斯·拉格朗日命名的一种多项式插值方法。许多实际问题中都用函数来表示某种内在联系或规律，而不少函数都只能通过实验和观测来了解。如对实践中的某个物理量进行观测，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日（插值）多项式。

推导

假设现在我们得到了$n+1$ 个点：$\left(x_0,y_0\right),\left(x_1,y_1\right),\dots ,\left(x_n,y_n\right)$，需要求过这些点的$n$次多项式。

给出拉格朗日基本多项式：$${\ell }_i(x)=\prod _{j=0,j!=i}^n\frac{x-x_j}{x_i-x_j}$$
可以看出这个基本多项式构造十分巧妙，它具有${\ell }_i\left(x_i\right)=1$，并且${\ell }_i\left(x_j\right)=0,\forall i!=j$的性质。那么，接着构造出这个$n$次多项式：$$P(x)=\sum _{i=0}^n{y}_i{\ell }_i(x)$$
根据基本多项式的性质，我们可以知道$P\left(x_i\right)=y_i$，也就是经过了这$n+1$个点。通过简单的多项式乘法和多项式除法就可以在$O(n^2)$的时间求出这个多项式的系数表达。

举例：假设给出的三个点为$(1,3)(2,7)(3,13)$，那么：
$$P(x)=3\frac{(x-2)(x-3)}{(1-2)(1-3)}+7\frac{(x-1)(x-2)}{(2-1)(2-3)}+13\frac{(x-1)(x-2)}{(3-1)(3-2)}$$

当x连续时的优化

在绝大多数题目中我们需要用到的 $x_i$ 的取值都是连续的，这样的话我们可以把上面的算法优化到$O(n)$复杂度。

先给出上文推导出的拉格朗日差值公式：$$f(k)=\sum _{i=0}^n{y}_i\prod _{i!=j}\frac{k-x[j]}{x[i]-x[j]}$$
对于分子来说，我们维护出关于$k$的前缀积和后缀积，也就是：$$pr{e}_i=\prod _{j=0}^{i}k-j,su{f}_i=\prod _{j=i}^{n}k-j$$
对于分母来说，观察发现这其实就是阶乘的形式，我们用 $fac[i]$ 表示 $i!$,
式子变为：$$f(k)=\sum _{i=0}^n{y}_i\frac{pr{e}_{i-1}\ast su{f}_{i+1}}{fac[i]\ast fac[N-i]}$$
注意:分母可能会出现符号问题，也就是说，当$N-i$为奇数时，分母应该取负号

重心拉格朗日插值法

再来看一下前面的式子$$f(k)=\sum _{i=0}^n{y}_i\prod _{i!=j}\frac{k-x[j]}{x[i]-x[j]}$$设$g={\prod }_{i=1}^{n}k-x[i]$，则：
$$f(k)=g\sum _{i=0}^n\prod _{i̸=j}\frac{y_i}{(k-x[i])(x[i]-x[j])}$$设$t_i=\frac{y_i}{{\prod }_{j!=i}{x}_i-x_j}$
$$f(k)=g\sum _{i=0}^n\frac{t_i}{(k-x[i])}$$
这样每次新加入一个点的时候只需要计算它的 $ti$ 即可

应用

拉格朗日插值法求自然数幂的前缀和

模板

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod = 1e9 + 7;
const int N = 5e4 + 10;

int t, e;
ll n, k;

//求过 (xi, yi) 的 e-1 次的多项式
ll inv[N];          //逆元
ll a[N];            //1 到 e，e个点，  也就是yi，视题目而定
ll pre[N], suf[N];  //前缀积 后缀积
ll fac[N];          //阶乘的逆元

ll q_pow(ll x, ll y){
    ll ans = 1;
    while(y){
        if(y&1)
            ans = ans * x % mod;
        x = x * x % mod, y >>= 1;
    }
    return ans;
}

ll solve(ll x, int e){         // 1 到 e，e个点插值
    if(x <= e)        //直接返回
        return a[x];
    x %= mod;
    ll ans = 0;
    pre[0] = suf[e+1] = 1;    //边界
    for(ll i = 1; i <= e; i++){      //(x-i) 前缀积
        pre[i] = pre[i - 1] * (x - i) % mod;    
    }
    for (ll i = e; i >=1 ; i--){   //(x-i) 后缀积
        suf[i] = suf[i + 1] * (x - i) % mod;    
    }
    for(ll i = 1; i <= e; i++){
        ll f = fac[i - 1] * fac[e - i] % mod; //从 0 到 e时，fac[i-1] 变为 fac[i];
        f = (e - i) & 1 ? -f : f;           //判断正负
        (ans += a[i] * f % mod * pre[i - 1] % mod * suf[i + 1]) %= mod;
    }
    ans += ans < 0 ? mod : 0;
    return ans;
}

int main()
{
    for(ll i = 1; i < N; i++){     //预处理逆元
        inv[i] = q_pow(i, mod - 2);
    }
    fac[0] = 1;
    for(int i = 1; i < N; i++){     //预处理阶乘逆元
        fac[i] = fac[i - 1] * inv[i] % mod;
    }
    cout << solve(n, e) << endl;
    return 0;
}

参考

https://www.cnblogs.com/zwfymqz/p/10063039.html#_label1_0
https://riteme.site/blog/2017-3-18/lagrange-interpolation.html
https://www.cnblogs.com/ECJTUACM-873284962/p/6833391.html#_labelTop
插值详解：https://www.cnblogs.com/duye/p/8671820.html