插值

留个脚本
就是用来求解

$\left\{\begin{array}{c} s+a_1x_1+a_2x_1^2+\cdots+a_{k-1}x_1^{k-1}=f(x_1) \\ s+a_1x_2+a_2x_2^2+\cdots+a_{k-1}x_2^{k-1}=f(x_2) \\ \vdots \\ s+a_1x_n+a_2x_n^2+\cdots+a_{k-1}x_n^{k-1}=f(x_n) \end{array}\right.$

已知$(x_1,f(x_1)),(x_2,f(x_2)),\cdots,(x_n,f(x_n))$,这样的n组数据

拉格朗日插值法

代码：

# 定义有限域或实数域
R.<x> = PolynomialRing(QQ)

# 已知点
points = [(0, 1), (1, 3), (2, 2)]

# 使用内置插值函数
f = R.lagrange_polynomial(points)

print("插值多项式 f(x) =", f)
print("验证 f(1) =", f(1))

快速拉格朗日插值

代码：

import multiprocessing as mp

p = 1041231053

X = []
Y = []
for line in open(r'D:\CTF_challs\py\crypto\2024\XYCTF 2024\encoded.txt', 'r').read().strip().split('\n'):
    x, y = line.split(' ')
    X.append(int(x))
    Y.append(int(y))

K = GF(p)
R = PolynomialRing(K, 'x')

def compZ(X):
    x = R.gen()
    Z = K(1)
    for xk in X:
        Z *= (x-xk)
    return Z

def comp(X, Y, Xother):
    Z = compZ(Xother)
    Y = [y/Z(x) for x, y in zip(X, Y)]
    return Y, Z

def solve(X, Y):
    n = len(Y)
    print("Solving for", n, "points...")

    # just use lagrange interpolation if the degree is small enough
    if n <= 10:
        return R.lagrange_polynomial(list(zip(X, Y)))

    nhalf = n // 2

    X1 = X[:nhalf]
    Y1 = Y[:nhalf]
    X2 = X[nhalf:]
    Y2 = Y[nhalf:]

    # parallelize the computation of the two halves
    if nhalf > 10000:
        with mp.Pool(2) as pool:
            result1 = pool.apply_async(comp, (X1, Y1, X2))
            result2 = pool.apply_async(comp, (X2, Y2, X1))

            Y1, Z2 = result1.get()
            Y2, Z1 = result2.get()
    else:
        Y1, Z2 = comp(X1, Y1, X2)
        Y2, Z1 = comp(X2, Y2, X1)

    # solve recursively
    f1 = solve(X1, Y1)
    f2 = solve(X2, Y2)

    # put it back together
    return f1*Z2 + f2*Z1

def test():
    Xt = X[:1000]
    Yt = Y[:1000]
    f = solve(Xt, Yt)
    for x, y in zip(Xt, Yt):
        assert f(x) == y

test()

f = solve(X, Y)

with open("111.bmp","wb") as tt:
    tt.write(bytearray(f.coefficients(sparse=False)[:-1]))

ProductTree(最优版本)：

from sage.rings.generic import ProductTree

F.<x> = Zmod(N, is_field=True)[]
# 1. 从 arr 中分离出插值节点 (xs) 和对应的函数值 (ys)
# arr 的结构是 [(i^e, result), ...]，所以 item[0] 是 x 坐标，item[1] 是 y 坐标
xs = [item[0] for item in arr]
ys = [item[1] for item in arr]

# 2. 使用 ProductTree 构建树并进行插值
# x 是你之前定义的 F.<x> 中的变量
poly = ProductTree([x - xi for xi in xs]).interpolation(ys)

Ben-Or/Tiwari算法

$\left(\begin {array}{c} 1 & x_1 & x_1^2 & \cdots & x_1^{k-1} \\ 1 & x_2 & x_2^2 & \cdots & x_2^{k-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^{k-1} \end{array} \right) \times \left(\begin {array}{c} s \\ a_1 \\ \vdots \\ a_{k-1} \end{array} \right) = \left(\begin {array}{c} f(x_1) \\ f(x_2) \\ \vdots \\ f{x_n} \end{array} \right)$

V = matrix(x)
assert rank(V) == t

s = -vector(y)
L = V.solve_right(s)

PR.<z> = PolynomialRing(F)
Lambda = z^t + sum(F(L[i]) * z^i for i in range(t))
R = Lambda.roots(multiplicities=False)
B = [m.log(g) for m in R]

V = matrix([[m^i for m in R] for i in range(t)])
b = vector(v[:t])

a = V.solve_right(b)
a0 = a[B.index(min(B))]

参考：
NSS 4th|N1gh7ma12e的小站
其他加密算法|Lazzaro