拉格朗日插值法

给定 $(t+1)$ 个不同的点,过这些点且最高次不大于 $t$ 的多项式,有且只有一条。

本文将介绍如何使用拉格朗日插值法,求解这样的多项式。

方法

已知点 $(x_0, y_0)$、$(x_1, y_1)$、…、$(x_t, y_t)$,拉格朗日插值法的思路是寻找多项式 $l_j(x)$,使其在 $x = x_j$ 时的取值为 1,在其他点的取值都是 0。

这样的多项式很好构造。

$$
l_j(x) = \frac{(x-x_0)…(x-x_{j-1})(x-x_{j+1})…(x-x_t)}{(x_j-x_0)…(x_j-x_{j-1})(x_j-x_{j+1})…(x_j-x_t)}
$$

注意到,当 $x=x_j$ 时,$l_j(x_j)$ 的分子和分母相同,所以 $l_j(x_j)=1$,当 $x$ 取其他值时,其分子总是 0,所以结果为 0。

多项式 $l_j(x)$ 就像开关一样,可以让 $y_jl_j(x)$ 满足只有在 $x = x_j$ 时的取值为 $y_j$,而在其他点的取值都是 0。

基于此,可以继续构造多项式

$$
L(x) = y_0l_0(x) + y_1l_1(x) + … + y_tl_t(x)
$$

满足过已知的 $(t+1)$ 个点。

给定点 $(x_0, y_0) = (1, 350)$、$(x_1, y_1) = (2, 770)$ 和 $(x_2, y_2) = (3, 1350)$,求过这三个点的多项式。

先构造 $l_j(x)$,有

$$
l_0(x) = \frac{(x - x_1)(x - x_2)}{(x_0 - x_1)(x_0 - x_2)} = \frac{(x-2)(x-3)}{(1-2)(1-3)} = \frac{1}{2}x^2 - \frac{5}{2}x + 3
$$

$$
l_1(x) = \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} = \frac{(x-1)(x-3)}{(2-1)(2-3)} = -x^2 + 4x - 3
$$

$$
l_2(x) = \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} = \frac{(x-1)(x-2)}{(3-1)(3-2)} = \frac{1}{2}x^2 - \frac{3}{2}x + 1
$$

再构造 $L(x)$,有

$$
L(x) = \sum_{i=0}^{2}y_il_i(x) = 350(\frac{1}{2}x^2 - \frac{5}{2}x + 3) + 770(-x^2 + 4x - 3) + 1350(\frac{1}{2}x^2 - \frac{3}{2}x + 1) = 80x^2 + 180x + 90
$$

求解完毕。

证明

存在性

显而易见,我们总能按公式构造出 $l_j(x)$,从而构造出 $L(x)$。

注意到,$l_j(x)$ 的分子是 $t$ 个 1 次多项式相乘,分母是一个整数,所以 $l_j(x)$ 是一个最高次不大于 $t$ 的多项式。

任意个最高次不大于 $t$ 的多项式相加,结果多项式的最高次也不会大于 $t$,所以 $L(x)$ 是一个最高次不大于 $t$ 的多项式。

唯一性

假设这样的多项式不唯一,对任意两个最高次不大于 $t$ 的拉格朗日多项式 $L_1$ 和 $L_2$。

因为 $L_1$ 和 $L_2$ 都过已知的 $(t+1)$ 个点,所以两者的差 $(L_1 - L_2)$ 在这 $(t+1)$ 个点的取值都是 0。

也就是说,$(L_1 - L_2)$ 一定是多项式 $(x-x_0)(x-x_1)…(x-x_t)$ 的倍数。

若 $L_1 - L_2 \neq 0$,则其次数一定不小于 $(t+1)$。

而任意个最高次不大于 $t$ 的多项式相减,结果多项式的最高次也不会大于 $t$,所以 $(L_1 - L_2)$ 的最高次也不大于 $t$。

矛盾。

所以 $L_1 - L_2 = 0$,即 $L_1 = L_2$。

完整代码

polynomial.py

Polynomial用来表示有限域 Secp256k1.n 上的多项式

$$
y \equiv a_0x^0 + a_1x^1 + … a_tx^t \pmod{Secp256k1.n}
$$

并实现了多项式加法、乘法和求值等基本运算。

对于方法interpolate_evaluate,我们的目标是计算插值得到的多项式在某点的取值,所以并不需要知道 $L(x)$ 的具体形式,可以直接将参数 $x$ 的值带入公式计算最终结果。方法的内层循环用来计算 $l_i(x)$ 的分子和分母的值,外层循环将每个 $y_il_i(x)$(分子和分母)的值保存到数组lagrange中,最后对所有的 $y_il_i(x)$ 通分并求和,得到 $L(x)$ 的值。

def interpolate_evaluate(points: list, x: int) -> int:
"""Lagrange interpolate with the giving points, then evaluate y at x"""
if len(points) < 2:
raise ValueError('Lagrange interpolation requires at least 2 points')
# [(numerator, denominator) ...]
lagrange = [(0, 0)] * len(points)
# the product of all the denominator
denominator_product = 1
for i in range(len(points)):
numerator, denominator = 1, 1
for j in range(len(points)):
if j != i:
numerator *= (x - points[j][0])
denominator *= (points[i][0] - points[j][0])
lagrange[i] = (points[i][1] * numerator, denominator)
denominator_product *= denominator
numerator_sum = 0
for (numerator, denominator) in lagrange:
numerator_sum += numerator * denominator_product // denominator
return numerator_sum * modular_multiplicative_inverse(denominator_product, curve.n) % curve.n

参考