我想在有限域的点上使用 python 插值多项式，并得到一个具有该域系数的多项式。目前我正在尝试使用 SymPy 并专门插值（从sympy.polys.polyfuncs），但我不知道如何强制插值发生在特定的 gf 中。如果没有，这可以用另一个模块来完成吗？
编辑：我对 Python 实现/库感兴趣。
问问题
3294 次
1 回答
9
SymPy 的interpolating_poly不支持有限域上的多项式。但是在 SymPy 的底层有足够的细节来组合一个有限域类，并以一种直接的方式找到拉格朗日多项式的系数。
像往常一样，有限域 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。
GF 类，有限域
将算术实现为方法add, sub, mul, inv（乘法逆）。为方便测试，插值包括在 GF(p n ) 的点处使用 GF(p neval_poly ) 中的系数评估给定多项式。
请注意，构造函数用作 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
PolyRing 类，字段上的多项式
这一个更简单：它实现多项式的加法、减法和乘法，参考地面场进行系数运算。[::-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 中的给定 x 值和数组 Y 中的相应 y 值构建。它是基本多项式的线性组合，对于 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来自 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
于 2018-01-02T20:27:00.210 回答