拉格朗日插值法

发表于 阅读次数：

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

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

方法

已知点 (x0,y0)、(x1,y1)、…、(xt,yt)，拉格朗日插值法的思路是寻找多项式 lj(x)，使其在 x = xj 时的取值为 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 = xj 时，lj(xj) 的分子和分母相同，所以 lj(xj) = 1，当 x 取其他值时，其分子总是 0，所以结果为 0。

多项式 lj(x) 就像开关一样，可以让 yjlj(x) 满足只有在 x = xj 时的取值为 yj，而在其他点的取值都是 0。

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

L(x) = y0l0(x) + y1l1(x) + ... + ytlt(x)

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

例子

给定点 (x0,y0) = (1,350)、(x1,y1) = (2,770) 和 (x2,y2) = (3,1350)，求过这三个点的多项式。

先构造 lj(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)，有

$$ \begin{aligned} L(x) &= y_0l_0(x) + y_1l_1(x) + y_2l_2(x) \\[1em] &= 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) \\[1em] &= 80x^2 + 180x + 90 \end{aligned} $$

求解完毕。

证明

存在性

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

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

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

唯一性

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

因为 L1 和 L2 都过已知的 (t+1) 个点，所以两者的差 (L1−L2) 在这 (t+1) 个点的取值都是 0。

也就是说，(L1−L2) 一定是多项式 (x−x0)(x−x1)...(x−xt) 的倍数。

若 L1 − L2 ≠ 0，则其次数一定不小于 (t+1)。

而任意个最高次不大于 t 的多项式相减，结果多项式的最高次也不会大于 t，所以 (L1−L2) 的最高次也不大于 t。

矛盾。

所以 L1 − L2 = 0，即 L1 = L2。

代码

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

y ≡ a0x0 + a1x1 + ...atxt (mod  n)

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

对于方法 interpolate_evaluate，我们的目标是计算插值得到的多项式在某点的取值，所以并不需要知道 L(x) 的具体形式，可以直接将参数 x 的值带入公式计算最终结果。方法的内层循环用来计算 li(x) 的分子和分母的值，外层循环将每个 yili(x)（分子和分母）的值保存到数组 lagrange 中，最后对所有的 yili(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

参考