超详细拉格朗日插值法 @ Wings            分类 ACM
发布于 星期五, 七月 30 日, 2021 年
更新于 星期五, 七月 30 日, 2021 年

本来是多项式里的内容, 学完后发现好多啊, 于是独立成文.

拉格朗日插值

拉格朗日公式如下:

$$L(x) = \sum_{i=0}^{n-1} y_i \frac{\prod \limits_{j \ne i} (x - x_j)}{\prod \limits_{j \ne i} (x_i - x_j)}$$

这个公式的意思是, 根据 $n$ 个点计算的 $n-1$ 次($n$ 次数界)多项式为 $L(x)$

有个大佬用通俗易懂的方法解释了这个公式, 我直接copy. (其实是我看不懂crt方法的证明)

下文引自 拉格朗日插值法入门 @crashed 大佬下标 $1 \sim n$, 我这里下标 $0 \sim n-1$, 注意一下.


它的思想非常简单, 大概可以理解为——硬性配凑!

怎么配凑呢?我们举个例子, 平面上有三个点, $(x_1,y_1),(x_2,y_2),(x_3,y_3)(x_1<x_2<x_3)$, 我们现在用它们仨插值.

根据小学基础知识, 我们知道, 这三个点肯定可以唯一确定一个二次函数.

那么我们就尝试找到它, 怎么找?拉格朗日想到了一个比较粗暴的方法——咱对于每个点都搞一个子函数 $f_i(x)$, 要求 $f_i(x)$ 在 $x = x_i$ 的时候得到 $1$, 在 $x = x_j(j \ne i)$ 的时候得到 $0$, 把 $n$ 个子函数凑起来, 得到的函数不就过了 $n$ 个点了!

也就是说, 我们要计算 $n$ 个子函数, 第 $i$ 个子函数为:

$$f_i(x)= \begin{cases} 1 & x=x_i\newline 0 & x=x_j(j\not=i)\newline I\ don’t\ care & otherwise \end{cases}$$

那么插值的结果就是:

$$f(x)=\sum_{i=1}^n y_if_i(x)$$

回到原问题上面来. 考虑构造 $f_1(x)$, 对于 $f_1(x_j) = 0 (j >1 )$ 的情况很好满足, 可以想到:

$$f_1(x) = (x - x2) (x - x3)$$

怎么让 $f_1(x_1) = 1$ 呢? 我们只需要把不用的丢掉就好:

$$f_1(x)=\frac{(x-x_2)(x-x_3)}{(x_1-x_2)(x_1-x_3)}$$

所以最后就有 $f_1(x)$ 长成下面这个亚子:

f1(x)

同理构造出 $f_2(x)$:

f2(x)

和 $f_3(x)$:

f3(x)

求和得到 $f(x)$:

$$f(x)=\frac{y_1(x-x_2)(x-x_3)}{(x_1-x_2)(x_1-x_3)}+\frac{y_2(x-x_1)(x-x_3)}{(x_2-x_1)(x_2-x_3)}+\frac{y_3(x-x_1)(x-x_2)}{(x_3-x_1)(x_3-x_2)}$$

f(x)

也可以很方便地推广到一般形式: 对于 $n$ 个点 $(x_1,y_1),(x_2,y_2),…,(x_n,y_n)(x_1<x_2<…<x_n)$, 设:

$$f_i(x)=\frac{\prod_{j\not=i}(x-x_j)}{\prod_{j\not=i}(x_i-x_j)}$$

那么插值结果就是:

$$f(x)=\sum_{i=1}^ny_if_i(x)=\sum_{i=1}^ny_i\times \frac{\prod_{j\not=i}(x-x_j)}{\prod_{j\not=i}(x_i-x_j)}$$


$\prod \limits_{j \ne i} (x - x_j)$ 决定了 $L(x)$ 的次数, 一共有 $n-1$ 个形如 $(x - k)$ (视 $x_j$ 为常数 $k$) 的东西相乘, 那么它就是 $n-1$ 次多项式($n$ 次数界多项式).

这个公式可以这样记录:

$$L(x) = \sum_{i=0}^{n-1} y_i \ell_i(x)$$ $$\ell_i(x) = \prod_{j=0,j\ne i}^{n-1} \frac{(x - x_j)}{(x_i - x_j)} = \frac{(x - x_0) \cdots (x - x_{i-1})(x - x_{i+1}) \cdots (x - x_{n-1})}{(x_i - x_0) \cdots (x_i - x_{i-1})(x_i - x_{i+1}) \cdots (x_i - x_{n-1})}$$

其中 $\ell_i(x)$ 称为拉格朗日基本多项式.

有了这个公式, 我们虽然不能直接求出多项式的系数, 但是能够直接求出 $L(x)$ 的值, 也就是 $n$ 个点确定的多项式 $L(x)$, 在任意点 $\xi$ 处的值 $L(\xi)$. 看这个公式, 两重循环, 所以复杂度是 $O(n^2)$ 的.

洛谷-4781【模板】拉格朗日插值

代码
int pls(LL a, LL b) {
	a += a < 0 ? P : 0;
	b += b < 0 ? P : 0;
	return a + b >= P ? a + b - P : a + b;
}

int mult(LL a, LL b) {
	a += a < 0 ? P : 0;
	b += b < 0 ? P : 0;
	return a * b >= P ? a * b % P : a * b;
}

int power(int a, int b = P-2) {
	int res = 1;
	while (b) {
		if (b&1)
			res = mult(res, a);
		a = mult(a, a);
		b >>= 1;
	}
	return res;
}

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

int main() {
	scanf("%d%d", &n, &k);
	for (int i = 1; i <= n; i++)
		scanf("%d%d", x + i, y + i);
	int L = 0;
	for (int i = 1; i <= n; i++) {
		int p = 1, q = 1;
		for (int j = 1; j <= n; j++)
			if (j != i) {
				p = mult(p, k - x[j]);
				q = mult(q, x[i] - x[j]);
			}
		L = pls(L, mult(y[i], mult(p, power(q))));
	}
	printf("%d\n", L);
	return 0;
}

计算多项式系数

虽然根据这个拉格朗日插值公式不能直接计算系数, 但是有间接计算的方法, 复杂度是 $O(n^2)$ 的. 算法导论上有习题. 先咕咕钴, 碰到再说.

当 x 连续时的优化

如果点值表示中 $n$ 个点的横座标 $x$ 是连续的, 那么可以把求 $L(\xi)$ 优化到 $O(n)$. 这也是拉格朗日插值法的常见考点(因为复杂度高于 $O(n\log n)$ 就可以直接 FFT 了, 要拉格朗日插值何用).

很明显, 这里主要优化的是 $\ell_i(x)$ 的计算.

$$\ell_i(x) = \frac{(x - x_0) \cdots (x - x_{i-1})(x - x_{i+1}) \cdots (x - x_{n-1})}{(x_i - x_0) \cdots (x_i - x_{i-1})(x_i - x_{i+1}) \cdots (x_i - x_{n-1})}$$

对不同的 $i$, 分子其实是有很多重复的部分的. 对于要求的值 $x = \xi$, $O(n)$ 预处理前缀积和后缀积, 即 $pre_i(\xi) = \prod_{j=0}^i (\xi - x_j)$, $suf_i(\xi) = \prod_{j=i}^{n-1} (\xi - x_j)$, 这样每次求 $\ell_i(\xi)$ 的分子, 就是 $O(1)$ 的了.

再考虑分母, 因为 $x_i$ 连续, 不妨设 $x_i < x_{i+1}$, 那么对于前一半, 即 $(x_i - x_0) \cdots (x_i - x_{i-1})$, 就是个阶乘 $(x_i - x_0)!$; 对于后一半来说, 先把每一项都取相反数, 变成 $(x_{n-1} - x_i) \cdots (x_{i+1} - x_i)$, 这也是个阶乘 $(x_{n-1} - x_i)!$, 然后再看对多少项取了相反数, 也就是说, 分母需要乘以 $(-1)^{n-1-i}$. 阶乘可以 $O(n)$ 预处理.

那么最后 $\ell_i(x)$ 可以 $O(1)$ 计算

$$\ell_i(x) = \frac{pre_{i-1}(x) \cdot suf_{i+1}(x)}{(-1)^{n-1-i}(x_i - x_0)! (x_{n-1} - x_i)!}$$

$L(x)$ 就可以 $O(n)$ 计算了.

(没有找到很板的板题)

多项式的前缀和

令 $S(x) = \sum_{i=0}^{x} A(i)$, 称 $S(x)$ 为 $A(x)$ 的(离散)前缀和. $S(x)$ 是一个多项式, 用点值表达的解释如下(其实是我看不懂数学推导):

假设有 $A(x)$ 是 $n$ 次数界多项式, 那么它可以由点 $(i, A(i)), i \in [0, n]$ 表示. 考虑 $S(x)$ 的点值表示为 $(i, S(i)) = (i, \sum_{j=0}^i A(j))$. 这样"对点做前缀和"得到的点, 就是点值表示的多项式 $S(x)$ 了.

需要注意的是, $n$ 次数界多项式的离散前缀和是一个 $n+1$ 次数界的多项式. 这里有一个简单理解(其实是我看不懂差分前缀和母函数者伯努力什么的):

在做前缀和的时候, 默认第一个点的前一个点($-1$?)的值是 $0$, 如果它不是 $0$, 那么答案显然不一样. 也就是说, 这里还有一个点. 这 $n+1$ 个点, 才能确定唯一一个 $n+1$ 次数界多项式.

计蒜客-40254 2019ICPC南昌邀请赛B.Polynomial

代码
int pre[MAXN], suf[MAXN], finv[MAXN], inv_neg_1, l[MAXN];

void init(int n) {
	int facn = 1;
	for (int i = 2; i <= n; i++)
		facn = mult(facn, i);
	finv[n] = power(facn);
	for (int i = n-1; ~i; i--)
		finv[i] = mult(i+1, finv[i+1]);
	inv_neg_1 = power(-1);
}

int Lagrange(int *y, int n, int x) {
	pre[0] = x;
	for (int i = 1; i < n; i++)
		pre[i] = mult(pre[i-1], x - i);
	suf[n] = 1;
	for (int i = n-1; ~i; i--)
		suf[i] = mult(suf[i+1], x - i);
	l[0] = mult(suf[1], finv[n-1]);
	if ((n-1)&1)
		l[0] = mult(l[0], inv_neg_1);
	for (int i = 1; i < n; i++) {
		l[i] = mult(mult(pre[i-1], suf[i+1]), mult(finv[i], finv[n-1-i]));
		if ((n-1-i)&1)
			l[i] = mult(l[i], inv_neg_1);
	}
	int L = 0;
	for (int i = 0; i < n; i++)
		L = pls(L, mult(y[i], l[i]));
	return L;
}

int T, n, m, y[MAXN], sy[MAXN];

int main() {
	scanf("%d", &T);
	init(1005);
	while (T--) {
		scanf("%d%d", &n, &m);
		for (int i = 0; i <= n; i++)
			scanf("%d", y + i);
		y[n+1] = Lagrange(y, n+1, n+1);
		sy[0] = y[0];
		for (int i = 1; i <= n+1; i++)
			sy[i] = pls(sy[i-1], y[i]);
		while (m--) {
			int l, r;
			scanf("%d%d", &l, &r);
			printf("%d\n", pls(Lagrange(sy, n+2, r), -Lagrange(sy, n+2, l-1)));
		}
	}
	return 0;
}

CF-622F The Sum of the k-th Powers

代码
int pre[MAXN], suf[MAXN], finv[MAXN], inv_neg_1, l[MAXN];
void init(int n) {
	int facn = 1;
	for (int i = 2; i <= n; i++) facn = mult(facn, i);
	finv[n] = power(facn);
	for (int i = n-1; ~i; i--) finv[i] = mult(i+1, finv[i+1]);
	inv_neg_1 = power(-1);
}
int Lagrange(int s, int *y, int n, int x) {
	pre[0] = x - s;
	for (int i = 1; i < n; i++) pre[i] = mult(pre[i-1], x - i - s);
	suf[n] = 1;
	for (int i = n-1; ~i; i--) suf[i] = mult(suf[i+1], x - i - s);
	l[0] = mult(suf[1], finv[n-1]);
	if ((n-1)&1) l[0] = mult(l[0], inv_neg_1);
	for (int i = 1; i < n; i++) {
		l[i] = mult(mult(pre[i-1], suf[i+1]), mult(finv[i], finv[n-1-i]));
		if ((n-1-i)&1) l[i] = mult(l[i], inv_neg_1);
	}
	int L = 0;
	for (int i = 0; i < n; i++) L = pls(L, mult(y[i], l[i]));
	return L;
}

int n, k, y[MAXN];

int main() {
	scanf("%d%d", &n, &k);
	init(k+2);
	y[0] = 0;
	for (int i = 1; i < k+2; i++)
		y[i] = pls(y[i-1], power(i, k));
	printf("%d\n", Lagrange(0, y, k+2, n));
	return 0;
}

拉格朗日基本多项式与多项式线性空间

一组 $\vec x$ 可以确定一组 $\vec l$, 那么对于不同的 $\vec y$, 可以看成是 $\ell_i$ 的不同系数. 所以可以把它看成多项式线性空间, 即 $L = \sum y_i \ell_i$. 好玩的东西在于, $\vec l$ 是线性无关的, 证明如下:

即证明, 仅有系数 $y_i$ 全为 $0$ 时, $L = 0$ (多项式线性空间下的 $\vec 0$ 等于数值 $0$). 设 $L = \sum y_i \ell_i = 0$, 有 $\forall \xi L(\xi) = 0$ , 因为 $L(x_i) = y_i$, 所以 $y_i = 0$. 证毕.

也就是说, 基本多项式 $\ell_i$ 一定是它张成的多项式线性空间的一组基底.

虽然目前不知道有什么用但是感觉好奇秒啊

并没有找到题, 但是我想到一个, 先不说, 看看之后能不能参与一波出题工作QwQ

重心拉格朗日插值法

改进版拉格朗日插值法

也称重心拉格朗日插值法第一型

对于已经计算过的点 $L(\varepsilon)$, 如果这个时候添加(或者删除)一个点值 $(x_k, y_k)$, 要重新求 $L(\varepsilon)$ 的话是 $O(n^2)$ 的. 对公式进行变形, 预处理一些东西, 能够做到 $O(n)$.

先对 $\ell_i(x)$ 变个形:

$$\begin{aligned} \ell_i(x) &= \frac{\prod_{j \ne i} (x - x_j)}{\prod_{j \ne i}(x_i - x_j)} \newline &= \frac{\prod_{j=0}^{n-1}(x - x_j)}{(x - x_i)} \cdot \frac{1}{\prod_{j \ne i}(x_i - x_j)} \end{aligned}$$

可以看到, 第一个分式的分子是和 $i$ 无关的东西, 令它为 $\ell(x) = \prod_{j=0}^{n-1}(x - x_j)$, 同时令 $w_i = \frac{1}{\prod_{j \ne i}(x_i - x_j)}$ (分母是 $\ell_i(x)$ 第二个分式的分母), 那么有

$$\ell_i(x) = \frac{\ell(x)w_i}{x - x_i}$$

其中, $w_i$ 称为重心权. 那么 $L(x)$ 就可以化为:

$$\begin{aligned} L(x) &= \sum_{i=0}^{n-1} y_i\ell_i(x) \newline &= \ell(x) \sum_{i=0}^{n-1}\frac{w_iy_i}{x - x_i} \end{aligned}$$

以增加点为例, 先 $O(n)$ 求 $\ell(x)$, 然后对 $O(n)$ 个 $w_i$ 在 $O(1)$ 内修改, 然后 $O(n)$ 求新增的一个 $w_k$, 最后 $O(n)$ 计算 $L(x)$.

重心拉格朗日插值法

也称重心拉格朗日插值法第二型

对一个 $n$ 次数界多项式 $g(x) \equiv 1$ 取 $n$ 个点, 进行改进版拉格朗日插值, 可以得到:

$$g(x) = \ell(x)\sum_{i=0}^{n-1} \frac{w_i}{x - x_i}$$

然后进行神奇的多项式除法, 得到

$$L(x) = \frac{L(x)}{g(x)} = \frac{\sum_{i=0}^{n-1}\frac{w_i}{x-x_i}y_i}{\sum_{i=0}^{n-1}\frac{w_i}{x-x_i}}$$

这个公式可以少维护一个 $\ell(x)$, 不过好像差不太多…

LOJ-165 拉格朗日插值

代码
void add(int x, int y, VI &X, VI &Y, VI &W) {
	int n = X.size();
	X.push_back(x);
	Y.push_back(y);
	W.push_back(1);
	for (int i = 0; i < n; i++)
		W[i] = mult(W[i], power(X[i] - x));
	for (int i = 0; i < n; i++)
		W[n] = mult(W[n], x - X[i]);
	W[n] = power(W[n]);
}

int Lagrange(VI &X, VI &Y, VI &W, int x) {
	int L = 0, g = 0, n = X.size();
	for (int i = 0; i < n; i++) {
		// 如果是已知点, 则直接求, 否则 $\frac{1}{x - x_i}$ 没有意义
		if (x == X[i]) return Y[i];
		int p = mult(W[i], power(x - X[i]));
		L = pls(L, mult(p, Y[i]));
		g = pls(g, p);
	}
	return mult(L, power(g));
}

int n;
VI X, Y, W;
int main() {
	scanf("%d", &n);
	while (n--) {
		int op, x, y;
		scanf("%d", &op);
		if (op == 1) {
			scanf("%d%d", &x, &y);
			add(x, y, X, Y, W);
		}
		else {
			scanf("%d", &x);
			printf("%d\n", Lagrange(X, Y, W, x));
		}
	}
	return 0;
}

留下昵称和邮箱, 可在第一时间获悉回复通知哦~

2021 FLAG

  • 找个妹子
  • 进计科
  • XCPC拿块金牌
  • 补全算法知识, 整全板子
  • 学会Web开发相关知识
  • 在服务器上搭建电子书库
  • 写个游戏并上线
  • 能弹一首曲子
  • 写首完整的曲子
  • 练习悠悠球

个人简介

我叫 Wings, 来自江西上饶, 目前人在西安, 是西电的一名学生.

常以 WingsWingsZengWingsWings的ID在各大小网站上游走, 一般来说, Wings不是我 😔, WingsZeng 一定是我 😊.

热爱算法, 喜欢钻研各种计算机技术.

业余爱好广泛, 只要不是文化课基本上都感兴趣😏.

开发/项目经历

  1. Android游戏 小墨滴的复仇 (弃坑)
  2. Android游戏 Circle Run (弃坑)
  3. Windows游戏 Snague (可能弃坑了吧)
  4. Python后端 Fathy' (可能弃坑了吧)

to be continued

教育经历

时间 学历 学校
2008-2014 小学 上饶市第十二小学
2014-2017 初中 上饶市第四中学
2017-2020 高中 上饶市第一中学
2020-2024 本科 西安电子科技大学
to be continued

比赛/竞赛经历

太久远太小的记不到了…

  1. 2017 国学竞赛初赛江西 没有分数或排名 二乙
  2. 2018 NOIP提高 258 省二
  3. 2019 CSP-S江西专场 145 省二
  4. 2019 数学竞赛初赛 70 没排名 (复赛打铁qaq)
  5. 2020 Gitee|Python贪吃蛇魔改大赛 可能是第四? 二等奖
  6. 2020 西电ACM训练基地熊猫杯 第四 银牌
  7. 2020 西安三校微软学生俱乐部Hackathon 和二等奖最后一名差0.5分 三等奖
  8. 2020 西电星火杯 三等奖
  9. 2020 西电ACM新生赛 第九 金牌
  10. 2020 ICPC 亚洲区域赛 济南站 132名 铜牌
  11. 2020-2021 第二届全国大学生算法设计与编程挑战赛(冬季赛) 924名 铜牌 (别骂了别骂了)
  12. 2020 ICPC 亚洲区域赛 昆明站 打星
  13. 2020 ICPC Asia-East Continent Final 签完到溜 打铁
  14. 西电"智能星"第一届自动驾驶小车比赛 第五 优胜奖|极速奖 本来可以冠军的别骂了别骂了
  15. 2021团体程序设计天体赛(CCCC) 个人二等奖
  16. 2021 西电 miniL CTF 优胜奖
  17. 2021 西电ACM校赛 第9名 金牌
  18. 2021 西电数模校赛 二等奖
  19. 2021 第15届IEEE 第48名
  20. 2021 CCPC 桂林站 打星

to be continued

爱好

技术

  • 算法
  • 独立游戏开发

游戏

  • Minecraft
  • Black Survival
  • I Wanna
  • Celeste
  • Life is Strange
  • Need for speed

运动

  • 篮球
  • 桌球
  • 乒乓球
  • 羽毛球
  • 慢跑

音乐

  • 吉他
  • 词曲
  • 流行

玩具

  • 魔方
    • 三阶速拧
    • 三阶盲拧
    • 高阶
  • yoyo球

追星

  • VAE
  • Benedict Cumberbatch