拉格朗日插值学习笔记

数学上来说，拉格朗日插值法可以给出一个恰好穿过二维平面上若干个已知点的多项式函数.

概述

根据小学知识我们可以知道,平面上任意$n+1$可以确定一个$n$次多项式$F(x)$.

现在我们已知这$n+1$个点和一个数$k$,求经过这$n+1$个点的多项式$F(x)$在$x=k$时的函数值$F(k)$.

一个显然的小学做法是利用待定系数法列出方程,然后利用高斯消元法$\mathcal O(n^3)$的解决.

但是拉格朗日插值可以做到$\mathcal O(n^2)$,在特殊情况下甚至可以$\mathcal O(n)$

拉格朗日插值

拉格朗日插值的精髓在于被构造出来的拉格朗日插值多项式

设已经给出的$n+1$个点的坐标分别为$(x_0,y_0), (x_1,y_1) \dots , (x_n,y_n)$

那么我们构造多项式:

$$\ell_i(k) = \prod_{j=0, i \neq j}^n \frac {k-x_j}{x_i - x_j}$$

容易发现:

$$\ell_i(x_j) = 1 \quad i = j$$

$$\ell_i(x_j) = 0 \quad i \neq j$$

为了满足$\ell_i(x_i) = y_i$,我们考虑构造$F(x)$:

$$F(x) = \sum_{i=0}^n \ell_i(x)y_i$$

然后可以得出$F(x)$的系数表达式,直接带入$k$计算即可.

Code

inline int lagrange(int *x, int *y, int k){
	int ans = 0;
	for(int i = 0; i < n; ++i) {
		int ell = 1;
		for(int j = 0; j < n; ++j) {
			if(i == j) continue;
			ell = 1ll * ell * (k - x[j]) % mod * quick_power(x[i] - x[j], mod - 2) % mod;
		}
		ans = (ans + 1ll * y[i] * ell) % mod;
	}
	return ans;
}

时间复杂度$\mathcal O(n^2)$.

特殊情况

当我们已知的这$n+1$个点取值连续的时候, 即已知$(0,y_1),(1,y_2) \dots (n,y_n)$时,

那么我们把他们带入上面构造的$F(x)$:

$$F(x) = \sum_{i=0}^n y_i \prod_{j=0,j \neq i}^{n} \frac {x-j}{i-j}$$

考虑构造:

$$P_i = \prod_{j=0}^i (x-x_j), \quad S_i = \prod_{j=i}^n (x-x_j)$$

带回原式可得:

$$F(x) = \sum_{i=0}^n \frac {P_{i-1} S_{i+1}}{(-1)^{n-i} \Gamma(i) \Gamma(n-i)}$$

时间复杂度$\mathcal O(n)$

例题

这里我们以LOJ2578「TJOI2018」为例.

可以发现我们使用的亵渎张数就是$m+1$

那么最后的贡献一定是:

$$\sum \sum_{i=1}^n i^{m+1} - a_i^{m+1}$$.

后面那一部分是空缺部分的贡献, 我们可以直接快速幂暴力算一下,

于是主要考虑前面一半.

首先这个式子显然是一个自变量为$n$的$m+2$次多项式,那么我们只需要暴力算出前$m+3$项,然后直接拉格朗日插值求解即可.

注意到这里我们求出的点的取值连续,所以复杂度为$\mathcal O(Tm^2)$.

参考资料

• 「算法笔记」拉格朗日插值 - By Siyuan

THE END