[学习笔记]拉格朗日&重心拉格朗日插值法

目录

• 〇、前言 ¶
• 壹、普通 Lagrange 插值法 ¶
• 贰、重心 Lagrange 插值法 ¶
• 叁、参考代码 ¶

〇、前言 ¶

鸽了好久了，于 $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;
}