拉格朗日插值法

简介

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

推导

假设现在我们得到了n+1n+1n+1 个点：(x0,y0),(x1,y1),…,(xn,yn)\left(x_{0}, y_{0}\right),\left(x_{1}, y_{1}\right), \ldots,\left(x_{n}, y_{n}\right)(x0,y0),(x1,y1),…,(xn,yn)，需要求过这些点的nnn次多项式。

给出拉格朗日基本多项式：ℓi(x)=∏j=0,j!=inx−xjxi−xj\ell_{i}(x)=\prod_{j=0, j != i}^{n} \frac{x-x_{j}}{x_{i}-x_{j}} ℓi(x)=j=0,j!=i∏nxi−xjx−xj
可以看出这个基本多项式构造十分巧妙，它具有ℓi(xi)=1\ell_{i}\left(x_{i}\right)=1ℓi(xi)=1，并且ℓi(xj)=0,∀i!=j\ell_{i}\left(x_{j}\right)=0, \forall i !=jℓi(xj)=0,∀i!=j的性质。那么，接着构造出这个nnn次多项式：P(x)=∑i=0nyiℓi(x)P(x)=\sum_{i=0}^{n} y_{i} \ell_{i}(x) P(x)=i=0∑nyiℓi(x)
根据基本多项式的性质，我们可以知道P(xi)=yiP\left(x_{i}\right)=y_{i}P(xi)=yi，也就是经过了这n+1n+1n+1个点。通过简单的多项式乘法和多项式除法就可以在O(n2)O(n^2)O(n2)的时间求出这个多项式的系数表达。

举例：假设给出的三个点为(1,3)(2,7)(3,13)(1,3)(2,7)(3,13)(1,3)(2,7)(3,13)，那么：
P(x)=3(x−2)(x−3)(1−2)(1−3)+7(x−1)(x−2)(2−1)(2−3)+13(x−1)(x−2)(3−1)(3−2)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)}P(x)=3(1−2)(1−3)(x−2)(x−3)+7(2−1)(2−3)(x−1)(x−2)+13(3−1)(3−2)(x−1)(x−2)

当x连续时的优化

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

先给出上文推导出的拉格朗日差值公式：f(k)=∑i=0nyi∏i!=jk−x[j]x[i]−x[j]f(k)=\sum_{i=0}^{n} y_{i} \prod_{i != j} \frac{k-x[j]}{x[i]-x[j]} f(k)=i=0∑nyii!=j∏x[i]−x[j]k−x[j]
对于分子来说，我们维护出关于kkk的前缀积和后缀积，也就是：prei=∏j=0ik−j,sufi=∏j=ink−jp r e_{i}=\prod_{j=0}^{i} k-j , s u f_{i}=\prod_{j=i}^{n} k-j prei=j=0∏ik−j,sufi=j=i∏nk−j
对于分母来说，观察发现这其实就是阶乘的形式，我们用 fac[i]fac[i]fac[i] 表示 i!i!i!,
式子变为：f(k)=∑i=0nyiprei−1∗sufi+1fac[i]∗fac[N−i]f(k)=\sum_{i=0}^{n} y_{i} \frac{p r e_{i-1} * s u f_{i+1}}{f a c[i] * f a c[N-i]} f(k)=i=0∑nyifac[i]∗fac[N−i]prei−1∗sufi+1
注意:分母可能会出现符号问题，也就是说，当N−iN−iN−i为奇数时，分母应该取负号

重心拉格朗日插值法

再来看一下前面的式子f(k)=∑i=0nyi∏i!=jk−x[j]x[i]−x[j]f(k)=\sum_{i=0}^{n} y_{i} \prod_{i != j} \frac{k-x[j]}{x[i]-x[j]} f(k)=i=0∑nyii!=j∏x[i]−x[j]k−x[j]设g=∏i=1nk−x[i]g=\prod_{i=1}^{n} k-x[i]g=∏i=1nk−x[i]，则：
f(k)=g∑i=0n∏i≠jyi(k−x[i])(x[i]−x[j])f(k)=g \sum_{i=0}^{n} \prod_{i \neq j} \frac{y_{i}}{(k-x[i])(x[i]-x[j])} f(k)=gi=0∑ni̸=j∏(k−x[i])(x[i]−x[j])yi设ti=yi∏j!=ixi−xjt_{i}=\frac{y_{i}}{\prod_{j != i} x_{i}-x_{j}}ti=∏j!=ixi−xjyi
f(k)=g∑i=0nti(k−x[i])f(k)=g \sum_{i=0}^{n} \frac{t_{i}}{(k-x[i])} f(k)=gi=0∑n(k−x[i])ti
这样每次新加入一个点的时候只需要计算它的 tititi 即可

应用

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

模板

#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