拉格朗日插值法学习笔记

拉格朗日插值法学习笔记

​ 在数值分析中，拉格朗日插值法是以法国十八世纪数学家约瑟夫·拉格朗日命名的一种多项式插值方法。许多实际问题中都用函数来表示某种内在联系或规律，而不少函数都只能通过实验和观测来了解。如对实践中的某个物理量进行观测，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日（插值）多项式。数学上来说，拉格朗日插值法可以给出一个恰好穿过二维平面上若干个已知点的多项式函数。拉格朗日插值法最早被英国数学家爱德华·华林于1779年发现，不久后（1783年）由莱昂哈德·欧拉再次发现。1795年，拉格朗日在其著作《师范学校数学基础教程》中发表了这个插值方法，从此他的名字就和这个方法联系在一起。

​ 给定 \(n + 1​\)个横坐标不相同的点，可以唯一确定一个 \(n​\) 次的多项式，我们可以直接列方程利用高斯消元求解，时间复杂度\(O(n^3)​\)。

​ 而另一种方法拉格朗日插值可以在\(O(n^2)\)的时间复杂度内完美的求解出答案，可以说是十分完美了。

拉格朗日插值法

为了接下来的讲解方便，我们引进一些函数和变量来辅助讲解。

拉格朗日插值法基函数（后文内简称为基函数）\(g(x)​\)：

下面给出\(g(x)​\)的定义式：\(g(x)=\Pi_{i=1,i!=x}^{n}\frac{k-X[j]}{X[x]-X[j]}​\)

根据拉格朗日插值法，\(f(k)=\sum_{i=1}^{n}y_i*g(i)​\)

于是，我们就可以在\(O(n^2)​\)的复杂度内求出\(f(k)​\)。

我们来看一个例子：

对于\(f(x)\),我们有\(f(1)=1,f(2)=2,f(3)=3\)，求\(f(100)\)。

我们先求出\(g(1,2,3)​\)：

\(g(1)=\frac{(100-2)*(100-3)}{(1-2)*(1-3)}=4753​\)，

\(g(2)=\frac{(100-1)*(100-3)}{(2-1)*(2-3)}=-9603​\)，

\(g(3)=\frac{(100-1)*(100-2)}{(3-1)*(3-2)}=4851​\)。

我们有\(f(100)=g(1)*1+g(2)*2+g(3)*3=4753-9603*2+4851*3=100\)。

同时，我们将\(f(k)​\)暴力展开有：

\(f(k)=1*\frac{(k-2)*(k-3)}{(1-2)*(1-3)}+2*\frac{(k-1)*(k-3)}{(2-1)*(2-3)}+3*\frac{(k-1)*(k-2)}{(3-1)*(3-2)}​\),

我们发现若我们将\(x_i\)带入时，其他项中都有\(x_i-x_i\)的项，必定为\(0\)，所以这个正确性是很显然的。

这是一道模板题：洛谷P4781

My code

#include <bits/stdc++.h>

using namespace std;

#define int long long
#define reg register
#define clr(a,b) memset(a,b,sizeof a)
#define Mod(x) (x>=mod)&&(x-=mod)
#define abs(a) ((a)<0?-(a):(a))
#define debug(x) cerr<<#x<<"="<<x<<endl;
#define debug2(x,y) cerr<<#x<<"="<<x<<" "<<#y<<"="<<y<<endl;
#define debug3(x,y,z) cerr<<#x<<"="<<x<<" "<<#y<<"="<<y<<" "<<#z<<"="<<z<<endl;
#define rep(a,b,c) for(reg int a=(b),a##_end_=(c); a<=a##_end_; ++a)
#define ret(a,b,c) for(reg int a=(b),a##_end_=(c); a<a##_end_; ++a)
#define drep(a,b,c) for(reg int a=(b),a##_end_=(c); a>=a##_end_; --a)
#define erep(i,G,x) for(int i=(G).Head[x]; i; i=(G).Nxt[i])
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#pragma GCC optimize(3,"Ofast","inline")

inline int Read(void) {
	int res = 0, f = 1;
	char c;
	while (c = getchar(), c < 48 || c > 57)if (c == '-')f = 0;
	do res = (res << 3) + (res << 1) + (c ^ 48);
	while (c = getchar(), c >= 48 && c <= 57);
	return f ? res : -res;
}

template<class T>inline bool Min(T &a, T const&b) {
	return a > b ? a = b, 1 : 0;
}
template<class T>inline bool Max(T &a, T const&b) {
	return a < b ? a = b, 1 : 0;
}

const int N=2e3+5,M=1e5+5,K=10,mod=998244353,INF=2e9;

bool MOP1;

namespace LP {
	int n,k,A[N],B[N];
	inline int Pow(int x,int y) {
		int res=1;
		while(y) {
			if(y&1)res=res*x%mod;
			x=x*x%mod,y>>=1;
		}
		return res;
	}
	inline int Inv(int x) {
		return Pow(x,mod-2);
	}
	inline void solve(void) {
		n=Read(),k=Read();
		rep(i,1,n)A[i]=Read(),B[i]=Read();
		int Ans=0;
		rep(i,1,n) {
			int L1=1,L2=1;
			rep(j,1,n)if(i^j) {
				L1=L1*(k-A[j])%mod;
				L2=L2*(A[i]%mod-A[j]%mod)%mod;
			}
			Ans+=(L1*Inv(L2)%mod)*B[i]%mod,Mod(Ans);
			Ans+=mod,Mod(Ans);
		}
		printf("%lld\n",Ans);
	}
}

bool MOP2;

void _main(void) {
	LP::solve();
}

signed main() {
	_main();
	return 0;
}

在\(x_i\)连续时的做法

根据拉格朗日插值法，我们有\(g(x)=\Pi_{i=1,i!=x}^{n}\frac{k-X[j]}{X[x]-X[j]}​\),\(f(k)=\sum_{i=1}^{n}y_i*g(i)​\)

关于\(f(k)\)的\(O(n)\)的复杂度已经无法在进行优化了，复杂度主要堆积在\(g(x)\)的求值上，我们需要优化\(g(x)\)的求值过程。

我们先用\(i​\)来代替\(x_i​\)可得：\(g(x)=\Pi_{i=1,i!=x}^{n}\frac{k-j}{x-j}​\)，

仔细看一下，我们发现了什么呢？

先看一下分子，我们可以利用前缀和的思想，维护一个前缀积和后缀积。

即\(pre_i=\Pi_{j=0}^{i}k-i,nxt_i=\Pi_{j=i+1}^{n}k-j\)。

我们发现分子就变成了\(pre_{i-1}*nxt_{i+1}​\)，可以在\(O(1)​\)的复杂度内求出。

在来看一下分母，这不就是一个阶乘的形式吗？

我们有分母\(=Frac[i]*Frac[n-i]*(-1)^{n-i}​\)。

所以我们对于\(g(i)​\)有\(g(i)=\frac {pre_{i-1}*nxt_{i+1}}{Frac[i]*Frac[n-i]*(-1)^{n-i}}​\)

于是我们就有了\(f(x)=\sum_{i=1}^{n}{\frac {pre_{i-1}*nxt_{i+1}}{Frac[i]*Frac[n-i]*(-1)^{n-i}}}\)。

于是我们就可以在\(O(n)\)的复杂度内解出了。。。

重心拉格朗日插值法

先咕咕咕了