拉格朗日插值法

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

給定 ( t + 1 ) (t+1) (t+1) 個不同的點，過這些點且最高次不大于 t t t 的多項式，有且隻有一條。

本文将介紹如何使用拉格朗日插值法，求解這樣的多項式。

方法

已知點 ( x 0 , y 0 ) (x_0, y_0) (x0​,y0​)、 ( x 1 , y 1 ) (x_1, y_1) (x1​,y1​)、…、 ( x t , y t ) (x_t, y_t) (xt​,yt​)，拉格朗日插值法的思路是尋找多項式 l j ( x ) l_j(x) lj​(x)，使其在 x = x j x = x_j x=xj​ 時的取值為 1，在其他點的取值都是 0。

這樣的多項式很好構造。

l j ( x ) = ( 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 ) 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)} lj​(x)=(xj​−x0​)...(xj​−xj−1​)(xj​−xj+1​)...(xj​−xt​)(x−x0​)...(x−xj−1​)(x−xj+1​)...(x−xt​)​

注意到，當 x = x j x=x_j x=xj​ 時， l j ( x j ) l_j(x_j) lj​(xj​) 的分子和分母相同，是以 l j ( x j ) = 1 l_j(x_j)=1 lj​(xj​)=1，當 x x x 取其他值時，其分子總是 0，是以結果為 0。

多項式 l j ( x ) l_j(x) lj​(x) 就像開關一樣，可以讓 y j l j ( x ) y_jl_j(x) yj​lj​(x) 滿足隻有在 x = x j x = x_j x=xj​ 時的取值為 y j y_j yj​，而在其他點的取值都是 0。

基于此，可以繼續構造多項式

L ( x ) = y 0 l 0 ( x ) + y 1 l 1 ( x ) + . . . + y t l t ( x ) L(x) = y_0l_0(x) + y_1l_1(x) + ... + y_tl_t(x) L(x)=y0​l0​(x)+y1​l1​(x)+...+yt​lt​(x)

滿足過已知的 ( t + 1 ) (t+1) (t+1) 個點。

給定點 ( x 0 , y 0 ) = ( 1 , 350 ) (x_0, y_0) = (1, 350) (x0​,y0​)=(1,350)、 ( x 1 , y 1 ) = ( 2 , 770 ) (x_1, y_1) = (2, 770) (x1​,y1​)=(2,770) 和 ( x 2 , y 2 ) = ( 3 , 1350 ) (x_2, y_2) = (3, 1350) (x2​,y2​)=(3,1350)，求過這三個點的多項式。

先構造 l j ( x ) l_j(x) lj​(x)，有

l 0 ( x ) = ( x − x 1 ) ( x − x 2 ) ( x 0 − x 1 ) ( x 0 − x 2 ) = ( x − 2 ) ( x − 3 ) ( 1 − 2 ) ( 1 − 3 ) = 1 2 x 2 − 5 2 x + 3 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 l0​(x)=(x0​−x1​)(x0​−x2​)(x−x1​)(x−x2​)​=(1−2)(1−3)(x−2)(x−3)​=21​x2−25​x+3

l 1 ( x ) = ( x − x 0 ) ( x − x 2 ) ( x 1 − x 0 ) ( x 1 − x 2 ) = ( x − 1 ) ( x − 3 ) ( 2 − 1 ) ( 2 − 3 ) = − x 2 + 4 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 l1​(x)=(x1​−x0​)(x1​−x2​)(x−x0​)(x−x2​)​=(2−1)(2−3)(x−1)(x−3)​=−x2+4x−3

l 2 ( x ) = ( x − x 0 ) ( x − x 1 ) ( x 2 − x 0 ) ( x 2 − x 1 ) = ( x − 1 ) ( x − 2 ) ( 3 − 1 ) ( 3 − 2 ) = 1 2 x 2 − 3 2 x + 1 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 l2​(x)=(x2​−x0​)(x2​−x1​)(x−x0​)(x−x1​)​=(3−1)(3−2)(x−1)(x−2)​=21​x2−23​x+1

再構造 L ( x ) L(x) L(x)，有

L ( x ) = ∑ i = 0 2 y i l i ( x ) = 350 ( 1 2 x 2 − 5 2 x + 3 ) + 770 ( − x 2 + 4 x − 3 ) + 1350 ( 1 2 x 2 − 3 2 x + 1 ) = 80 x 2 + 180 x + 90 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(x)=i=0∑2​yi​li​(x)=350(21​x2−25​x+3)+770(−x2+4x−3)+1350(21​x2−23​x+1)=80x2+180x+90

求解完畢。

證明

存在性

顯而易見，我們總能按公式構造出 l j ( x ) l_j(x) lj​(x)，進而構造出 L ( x ) L(x) L(x)。

注意到， l j ( x ) l_j(x) lj​(x) 的分子是 t t t 個 1 次多項式相乘，分母是一個整數，是以 l j ( x ) l_j(x) lj​(x) 是一個最高次不大于 t t t 的多項式。

任意個最高次不大于 t t t 的多項式相加，結果多項式的最高次也不會大于 t t t，是以 L ( x ) L(x) L(x) 是一個最高次不大于 t t t 的多項式。

唯一性

假設這樣的多項式不唯一，對任意兩個最高次不大于 t t t 的拉格朗日多項式 L 1 L_1 L1​ 和 L 2 L_2 L2​。

因為 L 1 L_1 L1​ 和 L 2 L_2 L2​ 都過已知的 ( t + 1 ) (t+1) (t+1) 個點，是以兩者的差 ( L 1 − L 2 ) (L_1 - L_2) (L1​−L2​) 在這 ( t + 1 ) (t+1) (t+1) 個點的取值都是 0。

也就是說， ( L 1 − L 2 ) (L_1 - L_2) (L1​−L2​) 一定是多項式 ( x − x 0 ) ( x − x 1 ) . . . ( x − x t ) (x-x_0)(x-x_1)...(x-x_t) (x−x0​)(x−x1​)...(x−xt​) 的倍數。

若 L 1 − L 2 ≠ 0 L_1 - L_2 \neq 0 L1​−L2​​=0，則其次數一定不小于 ( t + 1 ) (t+1) (t+1)。

而任意個最高次不大于 t t t 的多項式相減，結果多項式的最高次也不會大于 t t t，是以 ( L 1 − L 2 ) (L_1 - L_2) (L1​−L2​) 的最高次也不大于 t t t。

沖突。

是以 L 1 − L 2 = 0 L_1 - L_2 = 0 L1​−L2​=0，即 L 1 = L 2 L_1 = L_2 L1​=L2​。

完整代碼

polynomial.py

類

Polynomial

用來表示有限域 Secp256k1.n 上的多項式

y ≡ a 0 x 0 + a 1 x 1 + . . . a t x t ( m o d S e c p 256 k 1. n ) y \equiv a_0x^0 + a_1x^1 + ... a_tx^t \pmod{Secp256k1.n} y≡a0​x0+a1​x1+...at​xt(modSecp256k1.n)

并實作了多項式加法、乘法和求值等基本運算。

對于方法

interpolate_evaluate

，我們的目标是計算插值得到的多項式在某點的取值，是以并不需要知道 L ( x ) L(x) L(x) 的具體形式，可以直接将參數 x x x 的值帶入公式計算最終結果。方法的内層循環用來計算 l i ( x ) l_i(x) li​(x) 的分子和分母的值，外層循環将每個 y i l i ( x ) y_il_i(x) yi​li​(x)（分子和分母）的值儲存到數組

lagrange

中，最後對所有的 y i l i ( x ) y_il_i(x) yi​li​(x) 通分并求和，得到 L ( x ) L(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
           

參考

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