拉格朗日插值法

参考资料

---

参考资料

拉格朗日插值法入门
拉格朗日插值学习小结

简介

插值

在离散数据的基础上补插连续的函数，使得这条连续函数经过所有离散数据点，这个过程就叫插值。

其意义在于：
插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。

拉格朗日插值法

其核心思想在于直接硬性配凑，先分成不同的子函数，在合并成一个新的函数，很多时候只需要知道 \(O\)。

引入

模板题

给定 \(n\) 个点 \((x_i, y_i)\)，满足所有 \(x_i\) 互不相同，则存在唯一次数不超过 \(n-1\) 的多项式 \(y=f(x)\) 经过这些点。要求确定该多项式并计算 \(f(k) \bmod 998244353\)。
\(n \le 2\times 10^3,\ 1 \le x_i, y_i, k < 998244353\)，\(x_i\) 两两不同。

解法

基础

我们对每个点构造一个子函数，令其只在该点取值为 1，否则为 0。
也就是说，我们要计算 $ n $ 个子函数，第 $ i $ 个子函数为：

\[f_i(x) = \begin{cases} 1 & x = x_i \\ 0 & x = x_j (j \neq i) \\ Who \ \ cares & \text{otherwise} \end{cases} \]

那么插值的结果就是：

\[f(x) = \sum_{i=1}^n y_i f_i(x) \]

对于一般形式：对于 $ n $ 个点 \((x_1, y_1), (x_2, y_2), \dots, (x_n, y_n)\) (\(x_1 < x_2 < \cdots < x_n\))，设：

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

那么插值结果就是：

\[f(k) = \sum_{i=1}^n y_i f_i(k) = \sum_{i=1}^n y_i \times \frac{\prod_{j \neq i}(k - x_j)}{\prod_{j \neq i}(x_i - x_j)} \]

直接算上面这个求和式子就可以得到 \(O(n^2)\) 的做法，我们可以通过一些小小的变形变成 \(O(n^2)\) 预处理 \(O(n)\) 查询，具体看重心拉格朗日差值。

code

ll ans=0;
for(int i=1;i<=n;i++){
	ll p=1,q=1;
	for(int j=1;j<=n;j++) if(j!=i) {
		p=p*(k-x[j]+mod)%mod;
		q=q*(x[i]-x[j]+mod)%mod;
	}
	ans=(ans + y[i]%mod*p%mod*inv(q)%mod)%mod;
}

拉插其实是将 \(k+1\) 个点值表示成 \(k\) 次多项式，所以如果你想求一个 \(k\) 次多项式的值，我们需要知道 \(k+1\) 个点值。

拓展

当 \(x_i\) 连续

我们将原来的式子替换为

\[f(k)=\sum_{i=1}^n y_i\times {pre_{i-1}\times suf_{i+1} \over \prod_{j\neq i}(x_i-x_j)} \]

其中，\(pre_i,suf_i\) 分别表示 \(k-x_j\) 的前缀积和后缀积。

由于 \(x_i\) 是连续的，所以 \(x_i-x_{i-1}\) 差值固定，因此可以进行预处理。比如说，常见的对于 \(x_i-x_{i-1}=1\)，我们就可以预处理阶乘及其逆元。把式子改写成

\[f(k)=\sum_{i=1}^n y_i\times (-1)^{n-i} {pre_{i-1}\times suf_{i+1} \over fac_{i-1}\times fac_{n-i}} \]

然后我们就得到了 \(O(n)\) 的做法。

重心拉格朗日差值法

动态加点，可以做到 \(O(n)\) 修改，\(O(n)\) 查询。
重写原来的式子

\[f(k) = \sum_{i=1}^n y_i \times \frac{\prod_{j \neq i}(k - x_j)}{\prod_{j \neq i}(x_i - x_j)} \]

令 \(g=\prod_j (k-x_j)\)，于是可以将原来的式子改写成

\[f(k)=g\sum_{i=1}y_i\times {1\over (k-x_i)w_i} \]

其中 \(w_i={\prod_{j\neq i} (x_i-x_j)}\)，每次插入一个新点是，这个是可以 \(O(n)\) 动态维护的。

code

for(int i=1;i<=n;i++){ //O(n^2)预处理 w 
	w[i]=1;
	for(int j=1;j<i;j++) w[j]=w[j]*(x[j]-x[i]+mod)%mod;
	for(int j=1;j<i;j++) w[i]=w[i]*(x[i]-x[j]+mod)%mod;
}
ll ans=0,g=1;
for(int i=1;i<=n;i++) g=g*(k-x[i]+mod)%mod;
for(int i=1;i<=n;i++) ans=(ans + y[i]*inv((k-x[i]+mod)%mod*w[i]%mod)%mod)%mod;
cout<<g*ans%mod<<"\n"; 

应用

一、多项式差值

暂时不会，可以参考一下资料 1。

二、求自然数幂和

除了拉插还有其他很多其他方法。由于这里是拉插文章，这里就讲拉插做法。
我们可以通过一些神秘方法（比如说差分/扰动法）知道一个重要结论：\(\sum_{i=1}^n i^k\) 可以表示成一个关于 \(n\) 的 \(k+1\) 次多项式。
现在问题变成我们要算出这个 \(k+1\) 次多项式的在 \(n\) 这个位置的单点值，我们要用拉插的前提条件是要先得出 \(k+2\) 个点值，于是我们算出 \(1\sim k+2\) 的点值，带入拉插式子中，又由于我们选择的横坐标是连续的，所以我们可以预处理从而做到 \(O(k)\) 的复杂度，当然由于算点值带个 \(\log\)，所以总复杂度是 \(O(k\log k)\)。

补充一下：其实可以做到 \(O(k)\)。
因为 \(i^k\) 这个是积性函数，所以只用算出质数位置上的 \(p^k\)，然后用线性筛算出来，前缀和一下即可，复杂度 \(O({k\over \ln k}\log_2 k)\)。

扰动法证明

我们从二项式展开入手：

\[(i+1)^{k+1} - i^{k+1} = \sum_{j=0}^{k} \binom{k+1}{j} i^j \]

现在，对等式两边同时求和（\(i\) 从 1 到 \(n\)）：

\[\sum_{i=1}^n \left[ (i+1)^{k+1} - i^{k+1} \right] = \sum_{i=1}^n \sum_{j=0}^k \binom{k+1}{j} i^j \]

左边（裂项相消）：

这是一个经典的 telescoping sum（裂项相消）：

\[(2^{k+1} - 1^{k+1}) + (3^{k+1} - 2^{k+1}) + \cdots + ((n+1)^{k+1} - n^{k+1}) = (n+1)^{k+1} - 1 \]

右边（交换求和顺序）：

我们把 \(\sum_{i=1}^n\) 和 \(\sum_{j=0}^k\) 交换顺序：

\[\sum_{j=0}^k \binom{k+1}{j} \sum_{i=1}^n i^j = \sum_{j=0}^k \binom{k+1}{j} S_j(n) \]

把右边展开，并把 \(j=k\) 的项单独提出来：

\[\binom{k+1}{k} S_k(n) + \sum_{j=0}^{k-1} \binom{k+1}{j} S_j(n) = (n+1)^{k+1} - 1 \]

因为 \(\binom{k+1}{k} = k+1\)，所以：

\[(k+1) S_k(n) + \sum_{j=0}^{k-1} \binom{k+1}{j} S_j(n) = (n+1)^{k+1} - 1 \]

移项，解出 \(S_k(n)\)：

\[S_k(n) = \frac{1}{k+1} \left[ (n+1)^{k+1} - 1 - \sum_{j=0}^{k-1} \binom{k+1}{j} S_j(n) \right] \]

结论观察：

现在我们用数学归纳法来看这个公式：

1. 假设对于所有 \(j < k\)，\(S_j(n)\) 都是 \(j+1\) 次多项式。
2. 那么 \(\sum_{j=0}^{k-1} \cdots S_j(n)\) 这个求和里，最高次项来自 \(S_{k-1}(n)\)，它是 \(k\) 次多项式。
3. 等式右边最前面的 \((n+1)^{k+1}\) 展开后，最高次项是 \(n^{k+1}\)（\(k+1\) 次）。
4. 后面减去的求和部分最高只有 \(k\) 次。
5. 整体除以 \(k+1\) 后，最高次项依然是 \(n^{k+1}\)。

因此，\(S_k(n)\) 必然是一个 \(k+1\) 次多项式。

三、优化DP

注意事项

1. 重心拉格朗日差值法在 \(x_i\) 连续时表现不是最佳。

2. 拉插其实是将 \(k+1\) 个点值表示成 \(k\) 次多项式，所以如果你想求一个 \(k\) 次多项式的值，我们需要知道 \(k+2\) 个点值。

题目

P4781 【模板】拉格朗日插值 \(\checkmark\)

其实就是上面的引入。

P5667 拉格朗日插值2

题目大意

给定不超过 \(n\) 次多项式的 \(n+1\) 个点值 \(f(0),f(1),\cdots,f(n)\) 和正整数 \(m\)，求 \(f(m),f(m+1),\cdots,f(m+n)\) 对 998244353 取模。
\(n\le 1.6\times 10^5,m\le 10^8\)。

解题思路

这个需要掌握 NTT，暂时不会。

CF622F The Sum of the k-th Powers \(\checkmark\)

就是第二个应用。

P4593 [TJOI2018] 教科书般的亵渎 \(\checkmark\)

题目大意

有血量分别为 \(1\sim n\) 中某些值的怪物，其中 \(m\) 个血量没有怪物。每次操作：所有存活怪物血量减 \(1\)，若有怪物死亡则立即再次释放（连锁，但连锁不额外计次也不额外得分）。杀死所有怪物所需的手动操作次数记为 \(k\)。每次手动操作时，对每个存活怪物，若其操作前血量为 \(x\)，则获得 \(x^k\) 分。求总得分 \(\bmod 10^9+7\)。

\(T\le 10,\ n\le 10^{13},\ m\le 50,\ 1\le a_i<n\)。

解题思路

容易发现 \(k=m+1-len\)，\(len\) 是从 \(n\) 开始有连续的空缺。

然后就是计算 \(\sum_{i=1}^n i^k\) 再减去空缺位置的贡献，注意：答案是每次操作贡献的答案累加，所以每次都要重新计算一下，实现自己推一推就行了，具体也可以看看代码。计算就是应用二。复杂度就是 \(O(Tm^2)\)。这篇文章会 hack 掉不勤奋取模。