拉格朗日插值法
给定 \((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)\)，有
\[ \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} \]
求解完毕。
证明
存在性
显而易见，我们总能按公式构造出 \(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 用来表示有限域 Secp256k1.n 上的多项式
\[ y \equiv a_0x^0 + a_1x^1 + ... a_tx^t \pmod{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