Avi*_*u28 4 python sympy polynomial-math galois-field finite-field
我想在有限域的点上使用python插值多项式，并在该域中获取具有系数的多项式。当前，我正在尝试使用SymPy并专门进行插值（来自sympy.polys.polyfuncs），但是我不知道如何强制在特定gf中进行插值。如果没有，可以使用其他模块来完成吗？
编辑：我对Python实现/库感兴趣。
小智 5
SymPy的interpolating_poly不支持有限域上的多项式。但是，在SymPy的内部有足够的细节可以组合一个有限域类，并以一种残酷直接的方式找到Lagrange多项式的系数。
通常，有限域GF（p n）的元素由小于n 的多项式表示，系数在GF（p）中。模是在场构造时选择的次数为n的递减多项式的模上完成的。用扩展的欧几里得算法进行反演。
多项式由系数列表表示，首先是最高次数。例如，GF（3 2）的元素是：
[], [1], [2], [1, 0], [1, 1], [1, 2], [2, 0], [2, 1], [2, 2]
空列表表示0。
工具算术为方法add，sub，mul，inv（乘法逆）。用于测试的便利插包括eval_poly其评估在GF（P具有系数给定的多项式Ñ在GF（p的点）ñ）。
请注意，该构造函数用作G（3，2），而不用作G（9）-素数及其幂分别提供。
import itertools
from functools import reduce
from sympy import symbols, Dummy
from sympy.polys.domains import ZZ
from sympy.polys.galoistools import (gf_irreducible_p, gf_add, \
gf_sub, gf_mul, gf_rem, gf_gcdex)
from sympy.ntheory.primetest import isprime
class GF():
def __init__(self, p, n=1):
p, n = int(p), int(n)
if not isprime(p):
raise ValueError("p must be a prime number, not %s" % p)
if n <= 0:
raise ValueError("n must be a positive integer, not %s" % n)
self.p = p
self.n = n
if n == 1:
self.reducing = [1, 0]
else:
for c in itertools.product(range(p), repeat=n):
poly = (1, *c)
if gf_irreducible_p(poly, p, ZZ):
self.reducing = poly
break
def add(self, x, y):
return gf_add(x, y, self.p, ZZ)
def sub(self, x, y):
return gf_sub(x, y, self.p, ZZ)
def mul(self, x, y):
return gf_rem(gf_mul(x, y, self.p, ZZ), self.reducing, self.p, ZZ)
def inv(self, x):
s, t, h = gf_gcdex(x, self.reducing, self.p, ZZ)
return s
def eval_poly(self, poly, point):
val = []
for c in poly:
val = self.mul(val, point)
val = self.add(val, c)
return val
这比较简单：它实现多项式的加，减和乘法，是指对系数进行运算的基础字段。[::-1]由于SymPy约定从最高功率开始列出单项式，因此存在很多列表反转。
class PolyRing():
def __init__(self, field):
self.K = field
def add(self, p, q):
s = [self.K.add(x, y) for x, y in \
itertools.zip_longest(p[::-1], q[::-1], fillvalue=[])]
return s[::-1]
def sub(self, p, q):
s = [self.K.sub(x, y) for x, y in \
itertools.zip_longest(p[::-1], q[::-1], fillvalue=[])]
return s[::-1]
def mul(self, p, q):
if len(p) < len(q):
p, q = q, p
s = [[]]
for j, c in enumerate(q):
s = self.add(s, [self.K.mul(b, c) for b in p] + \
[[]] * (len(q) - j - 1))
return s
所述拉格朗日多项式被构造成用于在列表X和对应于阵列Y. y值是基础多项式，一个用于每个基础多项式由相乘得到X的每个元素的线性组合给定的x值(x-x_k)的多项式，表示为[[1], K.sub([], x_k)]。分母是一个标量，因此它甚至更容易计算。
def interp_poly(X, Y, K):
R = PolyRing(K)
poly = [[]]
for j, y in enumerate(Y):
Xe = X[:j] + X[j+1:]
numer = reduce(lambda p, q: R.mul(p, q), ([[1], K.sub([], x)] for x in Xe))
denom = reduce(lambda x, y: K.mul(x, y), (K.sub(X[j], x) for x in Xe))
poly = R.add(poly, R.mul(numer, [K.mul(y, K.inv(denom))]))
return poly
K = GF(2, 4)
X = [[], [1], [1, 0, 1]] # 0, 1, a^2 + 1
Y = [[1, 0], [1, 0, 0], [1, 0, 0, 0]] # a, a^2, a^3
intpoly = interp_poly(X, Y, K)
pprint(intpoly)
pprint([K.eval_poly(intpoly, x) for x in X]) # same as Y
漂亮的印刷品只是为了避免输出中与类型相关的装饰。多项式显示为[[1], [1, 1, 1], [1, 0]]。为了提高可读性，我添加了一个函数来将其转换为更熟悉的形式，其中的符号a是有限域的生成器，并且x是多项式中的变量。
def readable(poly, a, x):
return Poly(sum((sum((c*a**j for j, c in enumerate(coef[::-1])), S.Zero) * x**k \
for k, coef in enumerate(poly[::-1])), S.Zero), x)
所以我们可以做
a, x = symbols('a x')
print(readable(intpoly, a, x))
并得到
Poly(x**2 + (a**2 + a + 1)*x + a, x, domain='ZZ[a]')
这个代数对象不是我们领域内的多项式，这仅仅是出于可读输出的目的。
作为一种选择，或仅是另一种安全检查，您可以lagrange_polynomial对相同的数据使用from Sage。
field = GF(16, 'a')
a = field.gen()
R = PolynomialRing(field, "x")
points = [(0, a), (1, a^2), (a^2+1, a^3)]
R.lagrange_polynomial(points)
输出： x^2 + (a^2 + a + 1)*x + a