Lagrange 插值

有一个 \(n\) 次多项式函数，现在你知道它 \(n + 1\) 处的点值，你该怎么还原这个多项式函数？

普通的高斯消元是 \(\mathcal{O}(n^{3})\) 的，但是这实在是太慢了，有没有什么更快的方法？

假设这 \(n + 1\) 个点分别是 \((x_{1}, y_{1}), (x_{2}, y_{2}), \cdots, (x_{n + 1}, y_{n + 1})\)，那么其实构造一个满足 \(x = x_{i}\) 时 \(f(x) = y_{i}\) 的函数即可。(默认所有的 \(x_{i}\) 互不相同。)

考虑对于每一个点设计一个函数 \(l_{i}(x)\) 使得 \(x = x_{i}\) 时 \(l_{i}(x) = 1\)，当 \(x = x_{j}(i \not = j)\) 时 \(l_{i}(x) = 0\)（这里不是 \(x \not = x_{i}\) 时 \(l_{i}(x) = 0\) 哦，要注意区分）。

此时就可以得到：

\[f(x) = \sum\limits_{i = 1}^{n + 1}l_{i}(x)y_{i} \]

接下来就是构造 \(l_{i}(x)\) 的事情了，首先 \(x = x_{j}(i \not = j)\) 时 \(l_{i}(x) = 0\)，那么很大可能 \(l_{i}(x)\) 里面会有 \((x - x_{j})(j \not = i)\) 这个因式，我们先把它丢进去，此时 \(l_{i}(x) = \prod\limits_{j \not = i}(x - x_{j})\)。

接下来考虑 \(x = x_{i}\) 时 \(l_{i}(x) = 1\) 这个条件，最直接的想法是给它丢很多的分母将原来的分子一一约掉，所以把 \(\dfrac{1}{\prod\limits_{j \not = i}(x_{i} - x_{j})}\) 丢进去，这样就得到了最终的 \(l_{i}(x)\)：

\[l_{i}(x) = \prod\limits_{j \not = i}\dfrac{x - x_{j}}{x_{i} - x_{j}} \]

这时候我们随便代入一个值进去，发现计算 \(f(x)\) 的时间复杂度只有 \(\mathcal{O}(n^{2})\)！不仅如此，当 \(x_{i}\) 连续时，计算方法还有进一步优化的空间！

具体地，不妨设 \(x_{1 \sim n + 1}\) 分别为 \(1 \sim n + 1\)，那么有：

\[l_{i}(x) = \prod\limits_{j \not = i}\dfrac{x - j}{i - j} \]

代入到 \(f(x)\) 中：

\[f(x) = \sum\limits_{i = 1}^{n + 1}y_{i}\dfrac{\left[\prod\limits_{j < i}(x - j)\right]\left[\prod\limits_{j > i}(x - j)\right]}{(i - 1)!(n + 1 - i)!(-1)^{n + 1 - i}} \]

令 \(pre_{i} = \prod\limits_{j \leqslant i}(x - j), suf_{i} = \prod\limits_{j \geqslant i}(x - j)\)，这两个东西可以递推预处理，而阶乘逆元同样可以预处理，此时计算 \(f(x)\) 的复杂度就降到了 \(\mathcal{O}(n)\)！

大部分时候我们都会使用到后面这种形式。

放一个我自己的板子，里面用到了我自己写的 modint：

modint lagrange(int n, modint* x, modint* y, modint k) {
	modint ret = 0, num, den;
	for(int i = 1; i <= n; ++i) {
		num = y[i], den = 1;
		for(int j = 1; j <= n; ++j) {
			if(i != j) {
				num *= (modint)k - x[j];
				den *= x[i] - x[j];
			}
		}
		ret += num / den;
	}
	return ret;
}
modint lagrange(int n, modint* y, modint k) {
	static modint fct[1000005], inv[1000005], pre[1000005], suf[1000005];
	if(k.x <= n) return y[k.x];
	modint ret = 0;
	fct[0] = 1, pre[0] = 1, suf[n + 1] = 1;
	for(int i = 1; i <= n; ++i) {
		pre[i] = pre[i - 1] * (k - i);
		fct[i] = fct[i - 1] * i;
	}
	inv[n] = fct[n].ksm(mod - 2);
	for(int i = n; i >= 1; --i) {
		suf[i] = suf[i + 1] * (k - i);
		inv[i - 1] = inv[i] * i;
	}
	for(int i = 1; i <= n; ++i) {
		ret += y[i] * pre[i - 1] * suf[i + 1] * inv[i - 1] * inv[n - i] * ((n - i) & 1 ? -1 : 1);
	}
	return ret;
}