拉格朗日插值

合集 - 学习笔记(52)

1.随机化算法2024-02-272.组合数学2024-02-273.概率和期望2024-02-274.树链剖分2024-02-275.矩阵乘法与矩阵快速幂2024-02-276.最短路2024-02-277.卡特兰数、Prüfer 序列、BSGS2024-02-278.分块和莫队2024-02-279.AC 自动机2024-02-2710.平衡树2024-02-2711.基础字符串算法2024-02-2712.基础数论2024-02-2713.Miller-Rabin 与 Pollard-Rho2024-08-2314.广义后缀自动机2024-06-2515.后缀自动机 SAM2024-06-2216.回文自动机 PAM2024-06-1817.Manacher2024-06-1618.后缀数组2024-06-0819.01 分数规划2024-06-0220.网络流2024-05-0521.整体二分2024-05-0522.cdq 分治2024-05-0423.点分治2024-04-2124.虚树2024-04-1425.扫描线2024-04-1126.笛卡尔树2024-04-0627.基环树2024-04-0228.树哈希2024-03-0929.莫比乌斯反演2024-11-2830.二分图2024-11-0731.朱刘算法2024-10-2932.最大团2024-10-2933.杂项2024-08-2434.杜教筛2024-12-01

35.拉格朗日插值2024-12-10

36.线段树综合2024-12-1737.可持久化数据结构2024-12-2438.K-D Tree2024-12-2539.Burnside 引理与 Polya 定理2024-12-2740.线性基2024-12-3041.替罪羊树01-0942.LCT01-1843.插头 dp01-2044.原根01-2045.多项式乘法01-2246.斯特林数01-2347.二项式反演与斯特林反演01-2348.Min-Max 容斥01-2449.辛普森积分法01-2550.Min_25 筛02-0351.凸包02-1652.2-SAT03-03

收起

1 引入

众所周知，给定二维平面上的任意 $n$ 个点，一定会有一个不超过 $n-1$ 次的多项式与之对应。那么假如给出 $n$ 个点，需要找到一个多项式 $f(x)$ 使其对应，这个过程就是多项式插值。

而拉格朗日插值就是求解多项式插值的一种插值法。

2 拉格朗日插值

2.1 朴素解法

假设我们现在给出的 $n$ 个点坐标为 $(x_i,y_i)$，考虑构造出如下函数：

$$f(x)=\sum _{i=1}^n{y}_i{g}_i(x)$$

其中 $g_i(x)$ 需要满足当 $x$ 取值为 $x_i$ 时值为 $1$，取值为其他 $x_j(i\ne j)$ 时值为 $0$。考虑因式定理，当 $x=x_j(i\ne j)$ 时该函数值为 $0$，说明 $g_i(x)$ 有因式 $(x-x_j)(i\ne j)$。也就是说，$g_i(x)$ 有因式 $\prod _{j\ne i}(x-x_j)$。

这个时候我们已经可以满足取值为其他 $x_j(i\ne j)$ 时值为 $0$ 这个条件了，考虑满足 $x$ 取值为 $x_i$ 时值为 $1$ 这个条件。实际上这个很简单，我们给上面那个因式除以 $\prod _{j\ne i}(x_i-x_j)$ 即可，所以最后构造出的 $g_i(x)$ 如下：

$$g_i(x)={\displaystyle \frac{\prod _{j\ne i}(x-x_j)}{\prod _{j\ne i}(x_i-x_j)}}$$

此时假如我们要查询一个单点 $x$ 的值 $f(x)$，直接套用上面公式求解显然是 $O(n^2)$ 的。

模板题：【模板】拉格朗日插值，代码如下：

#include <bits/stdc++.h>
#define int long long

using namespace std;

const int Maxn = 2e5 + 5;
const int Inf = 2e9;
const int Mod = 998244353;

int n, k, x[Maxn], y[Maxn];

int qpow(int a, int b) {
	int res = 1;
	while(b) {
		if(b & 1) res = res * a % Mod;
		a = a * a % Mod;
		b >>= 1;
	}
	return res;
}

int lagrange(int x[], int y[], int X) {
	int ans = 0;
	for(int i = 1; i <= n; i++) {
		int res = 1;
		for(int j = 1; j <= n; j++) {
			if(i == j) continue;
			res = res * (X - x[j]) % Mod * qpow(x[i] - x[j], Mod - 2) % Mod;
		}
		ans = (ans + res * y[i] % Mod + Mod) % Mod;
	}
	return ans;
}

signed main() {
	ios::sync_with_stdio(0);
	cin.tie(0), cout.tie(0);
	cin >> n >> k;
	for(int i = 1; i <= n; i++) {
		cin >> x[i] >> y[i];
	}
	cout << lagrange(x, y, k) << '\n';
	return 0;
}

2.2 线性求值解法

上面的朴素解法单次查询的复杂度是 $O(n^2)$ 的，显然不优，因为我们一定可以将构造出来的式子写成一个多项式，单次求解的复杂度应该是 $O(n)$ 的。所以问题在于怎么去求出多项式中每一项的系数。

考虑到上面构造出来的多项式本质上是：

$$f(x)=\sum _{i=1}^n{y}_i\times {\displaystyle \frac{\prod _{j\ne i}(x-x_j)}{\prod _{j\ne i}(x_i-x_j)}}$$

考虑令 $t_i={\displaystyle \frac{y_i}{\prod _{j\ne i}(x_i-x_j)}}$，此时多项式转化为 $f(x)=\sum _{i=1}^n{t}_i\times \prod _{j\ne i}(x-x_j)$。考虑先求出 $\prod _{j=1}^n(x-x_j)$ 这个多项式展开后每一项的系数，然后对于每一个 $x_i$ 将 $(x-x_i)$ 除掉即可。这样求出的多项式系数乘上 $t_i$ 之后就是其对应的贡献了。

这里说一下怎样直接除以多项式 $(x-x_i)$。

我们先看常数项 $a_0$，显然它的值就是 $\prod x_i$，想要实现除掉多项式 $(x-x_i)$ 的效果直接给 $a_0$ 除以 $-x_i$ 即可。

然后看一次项 $a_1$，我们同样直接给它除以一个 $-x_i$，但是此时会发现一个问题，即 $a_1$ 中本身有一项为 $\prod _{j\ne i}{x}_j$，我们需要先将这个值减去才能除以 $-x_i$。不过好巧不巧，我们发现此时的 $a_0$ 就是这个值，所以新的 $a_1$ 就是 $a_1$ 减去 $a_0$ 后再除以 $-x_i$。

对于后面的每一项都是同理，因此可以得出 ${a_i}^{\prime }={\displaystyle \frac{a_i-{a_{i-1}}^{\prime }}{-x_i}}$。

代码如下：

#include <bits/stdc++.h>
#define int long long

using namespace std;

const int Maxn = 2e5 + 5;
const int Inf = 2e9;
const int Mod = 998244353;

int n, k, x[Maxn], y[Maxn];

int qpow(int a, int b) {
	int res = 1;
	while(b) {
		if(b & 1) res = res * a % Mod;
		a = a * a % Mod;
		b >>= 1;
	}
	return res;
}

int t[Maxn], f[Maxn], g[Maxn], fg[Maxn];

void init(int x[], int y[]) {
	for(int i = 1; i <= n; i++) {
		t[i] = 1;
		for(int j = 1; j <= n; j++) {
			if(i == j) continue;
			t[i] = t[i] * (x[i] - x[j]) % Mod;
		}
		t[i] = qpow(t[i], Mod - 2) * y[i] % Mod;//求出 t[i]
	}
	fg[1] = 1;
	for(int i = 1; i <= n; i++) {
		for(int j = n; j >= 1; j--) fg[j] = (fg[j - 1] + fg[j] * (Mod - x[i]) % Mod) % Mod;//先求出所有 (x - x[i]) 的积
	}
	for(int i = 1; i <= n; i++) {
		int inv = qpow(Mod - x[i], Mod - 2);
		for(int j = 1; j <= n; j++) {
			g[j] = (fg[j] - g[j - 1]) * inv % Mod;//利用上面的式子
			f[j] = (f[j] + g[j] * t[i] % Mod + Mod) % Mod;//乘上 t[i] 之后贡献给对应系数
		}
	}
}

int lagrange(int X) {
	int ans = 0, res = 1;
	for(int i = 1; i <= n; i++) {//多项式求值
		ans = (ans + res * f[i] % Mod) % Mod;
		res = res * X % Mod;
	}
	return ans;
}

signed main() {
	ios::sync_with_stdio(0);
	cin.tie(0), cout.tie(0);
	cin >> n >> k;
	for(int i = 1; i <= n; i++) {
		cin >> x[i] >> y[i];
	}
	init(x, y);
	cout << lagrange(k) << '\n';
	return 0;
}

2.3 横坐标连续的插值

在某些题目中，我们获得的 $x_i$ 坐标是连续的整数，这个时候我们有一种 $O(n)$ 预处理，$O(n)$ 求解的做法。

考虑回到构造出的函数 $f(x)=\sum _{i=1}^n{y}_i\times {\displaystyle \frac{\prod _{j\ne i}(x-x_j)}{\prod _{j\ne i}(x_i-x_j)}}$，此时它应该变成 $f(x)=\sum _{i=1}^n{y}_i\times {\displaystyle \frac{\prod _{j\ne i}(x-j)}{\prod _{j\ne i}(i-j)}}$。考虑求出 $(x-i)$ 的前缀积和后缀积，分子就应该变成 $pr{e}_{i-1}\times su{f}_{i+1}$。然后再看分母，发现它其实就是两个阶乘而已，所以分母可以写作 $(i-1)!\times (n-i)!\times (-1{)}^{n-i}$。因此插值公式就可以写作：

$$f(x)=\sum _{i=1}^n{y}_i\times {\displaystyle \frac{pr{e}_{i-1}\times su{f}_{i+1}}{(i-1)!\times (n-i)!\times (-1{)}^{n-i}}}$$

此时即可 $O(n)$ 求解。

代码如下：

int y[Maxn], f[Maxn], g[Maxn], pre[Maxn], suf[Maxn];

int qpow(int a, int b) {
	int res = 1;
	while(b) {
		if(b & 1) res = res * a % Mod;
		a = a * a % Mod;
		b >>= 1;
	}
	return res;
}

void init(int n) {
	f[0] = 1;
	for(int i = 1; i <= n; i++) f[i] = f[i - 1] * i % Mod;
	g[n] = qpow(f[n], Mod - 2);
	for(int i = n - 1; i >= 0; i--) g[i] = g[i + 1] * (i + 1) % Mod;
}

int lagrange(int n, int X) {
	pre[0] = suf[n + 1] = 1;
	for(int i = 1; i <= n; i++) pre[i] = pre[i - 1] * (X - i) % Mod;
	for(int i = n; i >= 1; i--) suf[i] = suf[i + 1] * (X - i) % Mod;
	int res = 0;
	for(int i = 1; i <= n; i++) {
		int num = y[i] * pre[i - 1] % Mod * suf[i + 1] % Mod * g[i - 1] % Mod * g[n - i] % Mod;
		if((n - i) & 1) res = (res - num + Mod) % Mod;
		else res = (res + num) % Mod;
	}
	return res;
}

下面来看几道经典例题。

3 例题

例 1： [湖北省队互测2014] 一个人的数论

题意： 令 $f_d(n)$ 表示所有小于 $n$ 且与 $n$ 互质的正整数的 $d$ 次方之和，求 $f_d(n)$。

发现 $n$ 过于巨大，因此我们想要构造出一个积性函数，这样就可以拆成每一个质因子去求解了。首先进行如下变换：

$$\begin{array}{rl}{f}_k(n)=& \sum _{i=1}^n{i}^k\cdot [gcd(i,n)=1]\\ =& \sum _{i=1}^n{i}^k\sum _{d\mid gcd(i,n)}\mu (d)\\ =& \sum _{d\mid n}\mu (d)\sum _{d\mid i}{i}^k\\ =& \sum _{d\mid n}\mu (d)\cdot d^k\sum _{i=1}^{\frac{n}{d}}{i}^k\end{array}$$

发现整个式子是一个狄利克雷卷积的形式，但是后面自然数的幂次方和并不是一个积性函数。这个时候运用拉格朗日插值，熟知自然数的 $k$ 次方和可以写作一个 $k+1$ 次多项式，所以只需要给定 $k+2$ 个点值即可求出该多项式。

于是设该多项式为 $g(x)=\sum _{i=0}^{k+1}{a}_i{x}^i$，则可得：

$$\begin{array}{rl}{f}_k(n)=& \sum _{d\mid n}\mu (d)\cdot d^k\sum _{i=1}^{\frac{n}{d}}{i}^k\\ =& \sum _{d\mid n}\mu (d)\cdot d^k\sum _{i=0}^{k+1}{a}_i({\displaystyle \frac{n}{d}}{)}^i\\ =& \sum _{i=0}^{k+1}{a}_i\sum _{d\mid n}\mu (d)\cdot d^k\cdot ({\displaystyle \frac{n}{d}}{)}^i\\ =& \sum _{i=0}^{k+1}{a}_i\cdot n^i\sum _{d\mid n}\mu (d)\cdot d^{k-i}\end{array}$$

这个时候发现后面的和式就是一个积性函数了。由于 $\mu (d)$ 在 $d$ 有平方因子的时候值为 $0$，所以 $d$ 的每种质因子都应只出现一次。于是令 $x=\prod p_i$，则后面和式的条件应该化为 $d\mid x$。设 $F(n)=\sum _{d\mid n}\mu (d)\cdot d^{k-i}$，则有：

$$\begin{array}{rl}{f}_k(n)=& \sum _{i=0}^{k+1}{a}_i\cdot n^i\cdot F(x)\\ =& \sum _{i=0}^{k+1}{a}_i\cdot n^i\prod _{j=1}^{w}F(p_j)\\ =& \sum _{i=0}^{k+1}{a}_i\cdot n^i\prod _{j=1}^w(1-p_j^{k-i})\end{array}$$

直接求解该式子即可，复杂度 $O(k^2+kw)$。

通过上例发现，拉格朗日插值的一个经典用途就是求出自然数的幂次方和。

例 2： [JLOI2016] 成绩比较

我们考虑 dp。由于最后一定会有 $k$ 个人被碾压，所以剩下的 $n-1-k$ 个人都至少有一科考的比 B 高。所以设 $dp(i,j)$ 表示枚举到第 $i$ 科，还有 $j$ 个人没有满足 “至少有一科考的比 B 高” 的条件。枚举当前有多少个人考的比 B 高，则可得如下转移方程：

$$dp(i,j)=\sum _{x=0}^{n}dp(i-1,j+x)(\genfrac{}{}{0ex}{}{j+x}{x})(\genfrac{}{}{0ex}{}{n-1-k-j-x}{r_i-1-x})\times \sum _{p=1}^{u_i}(u_i-p{)}^{r_i-1}{p}^{n-r_i}$$

首先从 $j+x$ 个人中任选 $x$ 个，让他们考的比 B 高，则这些人就满足了要求。此时我们会有 $n-1-k-j-x$ 已经满足了要求的人，从中任选几个填满 $r_i-1$ 即可。

发现这个式子的后半部分无法直接求解，因为 $u_i\le {10}^9$。但此时不难发现后面的和式中，每一项都是一个关于 $u_i$ 的 $n-1$ 次多项式，求总和后就是一个关于 $u_i$ 的 $n$ 次多项式。所以预处理出前 $n+1$ 个点值，然后拉格朗日插值求解即可。复杂度 $O(n^3)$。

例 3： [集训队互测 2012] calc

首先将有序的 $a$ 转化成选出一个无序的集合，算出总权值后乘上 $n!$ 即为最后答案。

考虑对于无序集合怎么做，不难想到设 $dp(i,j)$ 表示在 $[1,j]$ 中选了 $i$ 个数的权值总和，则有转移方程：

$$dp(i,j)=dp(i,j-1)+dp(i-1,j-1)\times i$$

发现此时直接暴力做是 $O(kn)$ 的，肯定无法通过。此时我们将 $i$ 看作常数，$j$ 看作变量，$dp(i,j)$ 实际上可以看作是一个关于 $j$ 的多项式。那么现在的问题就在于这个多项式的次数，于是设 $dp(i,j)$ 是一个 $g(i)$ 次多项式。根据上面的转移式，可以得到 $dp(i,j)-dp(i,j-1)=dp(i-1,j-1)\times i$。

左边是多项式差分的形式，所以差分过后应该会得到一个 $g(i)-1$ 次多项式，而右边本来是一个 $g(i-1)$ 次多项式，再乘上 $i$ 后就变为了 $g(i-1)+1$ 次多项式，所以会有如下关系式：

$$g(i)-1=g(i-1)+1$$

所以可以得到 $g(i)=g(i-1)+2$ 这个递推式。又由于 $dp(0,0)=1$ 导出 $g(0)=0$，所以 $g(i)=2i$。

也就是说，$dp(i,j)$ 是一个关于 $j$ 的 $2i$ 次多项式，只需要知道 $2i+1$ 个点值即可拉格朗日插值求出 $dp(n,k)$。