[学习笔记]拉格朗日&重心拉格朗日插值法
〇、前言 ¶
鸽了好久了，于 \(2021/4/5\) 看到重心 \(\tt Lagrange\) 插值法，惊讶于 \(\mathcal O(n)\) 的快速插值，想要把这篇完善一下。
壹、普通 Lagrange 插值法 ¶
对于一个 \(k\) 次多项式 \(f(x)\)，假如我们已经知道其中 \(k+1\) 个点的坐标 \((x_i,y_i)\)，那么我们并不需要知道这个多项式的系数（其实也是可解的）而直接求任意一个 \(x\) 对应的坐标。
假设我们已知 \(k+1\) 对 \((x_i,y_i)\)，现在我们要求 \(x\) 所对应的 \(y\)，有
\[y=\sum_{i=0}^k y_i\prod_{j=0,j\neq i}^k \frac{x-x_j}{x_i-x_j}
\]
直接求得 \(y\) 咯。
当然，如果你想要联立 \(k+1\) 个多项式高斯消元解出 \(f\)，再求 \(f(x)\)，这也是可以的 但是似乎拉格朗日比高斯消元好记
求原函数或单点差值均 \(\mathcal O(n^2)\).
贰、重心 Lagrange 插值法 ¶
考虑我们的插值法计算了很多重复的部分，不妨把式子变个型：
\[\begin{aligned}
y&=\sum_{i=0}^k y_i\prod_{j=0,j\neq i}^k \frac{x-x_j}{x_i-x_j} \\
&=\sum_{i=0}^ky_i\prod_{j=0,j\neq i}^k(x-x_j)\prod_{p=0,p\neq i}^k{1\over x_i-x_p} \\
&=\sum_{i=0}^k\prod_{j=0,j\neq i}^k(x-x_j)\prod_{p=0,p\neq i}^k{y_i\over x_i-x_p} \\
&=\prod_{j=0}^k(x-x_j)\left (\sum_{i=0}^k{y_i\over \prod_{p=0,p\neq i}^k(x_i-x_p)}\times {1\over x-x_i}\right)
\end{aligned}
\]
最后一步的变形，为了让 \(\prod\) 脱离 \(i\) 的限制，补上了一个 \(\large {1\over x-x_i}\)，再在后面乘上分数。
稍微整理一下，若
\[l(x)\overset\Delta=\prod_{i=0}^k(x-x_i) \\
w_i\overset\Delta={y_i\over \prod_{p=0,p\neq j}^k(x_i-x_j)}
\]
那么最后，得到差值过程：
\[y=l(x)\sum_{i=0}^k{w_i\over x-x_i}
\]
预处理仍然是 \(\mathcal O(n^2)\)，但是单次差值变成 \(\mathcal O(n)\)，求原函数仍然是 \(\mathcal O(n^2)\).
叁、参考代码 ¶
还没有时间打重心 \(\tt lagrange\) 的插值法，目前只打了 \(\mathcal O(n^2)\) 的。
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
template<class T>inline T readin(T x){
x=0; int f=0; char c;
while((c=getchar())<'0' || '9'<c) if(c=='-') f=1;
for(x=(c^48); '0'<=(c=getchar()) && c<='9'; x=(x<<1)+(x<<3)+(c^48));
return f? -x: x;
}
const int maxn=2000;
const int mod=998244353;
inline int qkpow(int a, int n){
int ret=1;
for(; n>0; n>>=1, a=1ll*a*a%mod)
if(n&1) ret=1ll*ret*a%mod;
return ret;
}
int x[maxn+5], y[maxn+5], n, k;
inline void input(){
n=readin(1), k=readin(1);
for(int i=1; i<=n; ++i)
x[i]=readin(1), y[i]=readin(1);
}
inline void getval(int k){
int ans=0;
for(int i=1; i<=n; ++i){
int up=1, down=1;
for(int j=1; j<=n; ++j) if(j!=i){
up=1ll*up*(k-x[j])%mod;
down=1ll*down*(x[i]-x[j])%mod;
}
ans=(ans+1ll*y[i]*up%mod*qkpow(down, mod-2)%mod)%mod;
}
printf("%d\n", (ans+mod)%mod);
}
signed main(){
input();
getval(k);
return 0;
}
