拉格朗日插值法

发表于 更新于 分类于 阅读次数：

这里只介绍关于多项式的拉格朗日插值法，对于一般函数的拟合当做一个多项式就好了。

插点值

已知多项式 \(f(x)\) 的 \(n\) 个点值 \((x_i, y_i = f(x_i))\) ，求 \(f(k)\) 。

拉格朗日插值法的思路在于： 对于每个 \((x_i, y_i)\) 找到 \(L_i(x)\) 使得 \(L_i(x_i) = y_i, L_i(x_j) = 0\) ， 其中 \(x_j\) 是已知的 \(x\) 中任意一个不等于 \(x_i\) 的 \(x\) 。

而由 \(L\) 的定义可知，\(f(k) = \sum_{i=1}^n L_i(k)\) 。
代入上式即可求解。

现在问题在于构造 \(L\) 。 下面的 \(L\) 可以满足定义：

\[L_i(x) = y_i \cdot \prod\limits_{j \ne i}\frac{x_j-x}{x_j-x_i}\]

代入可得这对于任意 \(x_i\) 可以使得 \(L_i(x_i) = y_i, L_i(x_j) = 0\)

复杂度 \(O(n^2)\) 。

特别地，如果已知的点值横坐标是一段连续的整数，该做法可以优化到线性，略。

插系数

已知多项式 \(f(x)\) 的 \(n\) 个点值，求 \(f(x)\) 的第 \(k\) 次项系数 。

先考虑插出一项 \(x^k\) 的系数。

首先依次考虑每个 \(L_i(x)\) 的 \(x^k\) 系数，最后累加即可。

注意到分母和 \(y_i\) 是常数可以直接算，考虑提出来，

\[ L_i(x) = (y_i \cdot \prod\limits_{j \ne i}\frac{1}{x_j-x_i}) (\prod\limits_{j \ne i} (x_j - x)) \]

左边的常数是一堆逆元相乘再乘上 \(y_i\) ，逆元往往不是 \(O(1)\) 计算的，
事实上可以先把分母中所有 \(x-x_j\) 乘起来再一起求逆元，这样就只需要计算 1 次逆元，此时一般可以忽略逆元对复杂度的影响。

那么只需要算分子部分，即上式右边的 \(\prod\) 的 \(x^k\) 系数即可。

假设把这个连乘暴力拆开，其实 \(x^k\) 的系数就是在之中选 \(n - k - 1\) 个常数。
DP 预处理 \(pre(i, k)\) 表示在 \(\prod\limits_{j \leq i} (x - x_j)\) 中选 k 个常数的系数和，
同理 \(suf(i, k)\) 表示在 \(\prod\limits_{j \geq i} (x - x_j)\) 中选 k 个常数的系数和。
那么枚举 \(p + q = n - k - 1\) 求 \(pre(i - 1, p)\) 乘 \(suf(i + 1, q)\) 的和就是上式的 \(x^k\) 系数了。

时间复杂度 \(O(n^2)\) 。

参考实现（未经过测试）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
/* 插出 L_i(x) 的 x^k 系数 */
int la(int i, int k) {
	int A = 1;
	for(int j = 1; j <= n; j ++)
		if(i != j)
			A *= X[i] - X[j];
	int B = 0;
	for(int l = 0; l <= n - k - 1; l ++)
		B += pre[i - 1][l] * suf[i + 1][n - k - 1 - l];
	return B / A * X[i];
}

/* 初始化 suf 和 pre */
void init() {
	for(int i = 0; i <= n; i ++) {
		pre[i][0] = 1;
		for(int j = 1; j <= i; j ++)
			pre[i][j] = pre[i - 1][j] - pre[i - 1][j - 1] * X[i];
	}

	for(int i = n + 1; i; i --) {
		suf[i][0] = 1;
		for(int j = 1; j <= n - i + 1; j ++)
			suf[i][j] = suf[i + 1][j] - suf[i + 1][j - 1] * X[i];
	}
}

插多项式

已知 \(n-1\) 次多项式 \(f(x)\) 的 \(n\) 个点值，求 \(f(x)\) 的所有系数 。

直接将所有暴力多项式相乘计算所有 \(L_i(x)\) ，或者使用 \(n\) 次上述的插系数，复杂度 \(O(n^3)\) ，难以接受。

点值比较特殊的情况下可以使用 FFT 或者 NTT ，但有失一般性。

沿用上述分离常数的方法，考虑 \(O(n)\) 求出每个 \(L_i(x)\) 的所有系数：

\[L_i(x) = (y_i \cdot \prod\limits_{j \ne i}\frac{1}{x_j-x_i}) (\prod\limits_{j \ne i} (x_j - x))\]

左边的常数还是同样地处理，不同的是对于右边的多项式现在要求的不是某一项的系数而是所有系数。

对于右边的多项式部分，可以看做：

\[\frac{\prod_{j=1}^n (x_j - x)}{x_i - x}\]

其分子是个与 \(i\) 无关的多项式，可以 \(O(n^2)\) 暴力预处理，而其分母是个简单的一次二项式。
那么可以用短除法来进行多项式除二项式，复杂度 \(O(n)\) 。

总时间复杂度为 \(O(n^2)\) 。

参考实现（未经过测试）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
void la(int n) {
	tmp[0] = 1; // tmp 表示 n 个二项式相乘的多项式
	for(int i = 1; i <= n; i ++) {
		tmp[i] = 0;
		for(int j = i; j; j --)
			tmp[j] = tmp[j - 1] - X[i] * tmp[j];
		tmp[0] *= - X[i];
	}

	for(int i = 1; i <= n; i ++) {
		ll A = 1; // A 表示 L_i 的常数部分
		for(int j = 1; j <= n; j ++)
			if(i != j)
				A *= X[i] - X[j];
		A = Y[i] / A;

		tmp2[n] = tmp[n]; // tmp2 表示将 tmp 除以一个二项式的结果
		for(int j = n; j; j --) {
			ll t = tmp2[j];
			get[j - 1] += t * A; // get 表示 f(x) 的系数
			tmp2[j - 1] = tmp[j - 1] + t * X[i];
		}

		// 此时应有 tmp2[0] = 0 ，否则说明除法有余数
		assert(tmp2[0] == 0);
	}
}