拉格朗日插值学习笔记

\(\text{拉格朗日插值学习笔记}\)

插值是一种通过已知的、离散的数据点推算一定范围内的新数据点的方法。插值法常用于函数拟合中。我们知道一个 \(k\) 次多项式可以由 \(k+1\) 个点唯一确定，由此我们可以得到拉格朗日插值算法。

思路

我们现在给定 \(n\) 个点坐标为 \((x_i,y_i)\)，考虑构造函数：

\[f(x)=\sum_{i=1}^ny_ig_i(x) \]

那么只需要构造出符合条件的 \(g_i\)。考虑 \(g_i\) 的构造需要在 \(x\) 取值为 \(x_i\) 时为 \(1\)，\(x\) 取 \(x_j(i\neq j)\) 时为 \(0\)。那么容易得到的是 \(g_i(x)\) 有因式 \(\prod_{j\neq i}(x-x_j)\)。那么满足 \(x\) 取 \(x_i\) 时 \(g_i\) 为 \(1\) 是容易的，将上式除以 \(\prod_{j\neq i}(x_i-x_j)\) 即可。

于是构造出的 \(g_i(x)\) 如下：

\[g_i(x)=\prod_{j\neq i}\frac{x-x_j}{x_i-x_j} \]

那么求单点 \(x\) 的函数值 \(f(x)\) 时复杂度是 \(O(n^2)\) 的。

下面给出 拉格朗日插值 模板题的代码：

#include <bits/stdc++.h>
#define int long long
#define mod 998244353 
#define N 2005
using namespace std;
int n, k;
int x[N], y[N];
int ans;
int qpow(int x, int y) {
	int ans = 1;
	while (y) {
		if (y & 1) ans = ans * x % mod;
		x = x * x % mod;
		y >>= 1; 
	}
	return ans;
}

signed main() {
	ios::sync_with_stdio(0);
	cin.tie(0);
	cin >> n >> k;
	for (int i = 1; i <= n; i++) cin >> x[i] >> y[i];
	for (int i = 1; i <= n; i++) {
		int res = y[i];
		for (int j = 1; j <= n; j++)
			if (j != i) res = res * (k - x[j]) % mod * qpow(x[i] - x[j], mod - 2) % mod;
		ans = (ans + res + mod) % mod;
	}
	cout << ans << "\n";
	return 0;
}

求系数

当需要多次求解函数值时，或许 \(O(n^2)\) 的复杂度无法接受，于是我们考虑求出函数每一项的系数，从而 \(O(n)\) 求解函数值。我们进行转换，令 \(t_i=\dfrac{y_i}{\prod_{j\neq i}(x_i-x_j)}\)，那么 \(f(i)=\sum_{i=1}^nt_i\prod_{j\neq i}(x-x_j)\)。于是我们求出 \(\prod_{j=1}^n\) 展开后每一项的系数，然后除掉对应项即可。下面给出预处理以及查询的代码。（认为插值的是 \(n\) 次多项式，\(x,y\) 的范围为 \([0,n]\)）：

int x[N], y[N], t[N];
int f[N], g[N], fg[N];
#define __(i) memset(i, 0, sizeof i)
void init(int n) {
	__(t), __(f), __(g), __(fg);
	for (int i = 0; i <= n; i++) {
		t[i] = 1;
		for (int j = 0; j <= n; j++)
			if (i != j) t[i] = t[i] * (x[i] - x[j]) % mod;
		t[i] = y[i] * qpow(t[i], mod - 2) % mod;
	}
	fg[0] = 1;
	for (int i = 0; i <= n; i++) {
		for (int j = n; j > 0; j--)
			fg[j] = (fg[j - 1] + fg[j] * (mod - x[i]) % mod) % mod;
		fg[0] = fg[0] * (mod - x[i]) % mod;
	}
	for (int i = 0; i <= n; i++) {
		int inv = qpow(mod - x[i], mod - 2);
		g[0] = fg[0] * inv % mod;
		for (int j = 1; j <= n; j++) g[j] = (fg[j] - g[j - 1]) * inv % mod;
		for (int j = 0; j <= n; j++) f[j] = (f[j] + t[i] * g[j] % mod + mod) % mod;
	}
}
int lgl(int x, int n) {
	int ans = 0;
	for (int i = 0; i <= n; i++)
		ans = (ans + f[i] * qpow(x, i) % mod) % mod;
	return ans;
}