拉格朗日插值法

https://aaron67.cc/2020/10/29/lagrange-interpolate/

给定 $(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_j{l}_j(x)$ 满足只有在 $x=x_j$ 时的取值为 $y_j$，而在其他点的取值都是 0。

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

$$L(x)=y_0{l}_0(x)+y_1{l}_1(x)+...+y_t{l}_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}_i{l}_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\ne 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_0{x}^0+a_1{x}^1+...a_t{x}^t(\mathrm{m}\mathrm{o}\mathrm{d}Secp256k1.n)$$

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

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

参考

• 维基百科
• 拉格朗日插值法