SymPy's 插值多边形不支持有限域上的多项式。但是 SymPy 背后有足够的细节来组合有限域的类，并找到拉格朗日多项式以一种残酷直接的方式。
As usual, the elements of finite field GF(pn) are represented by polynomials of degree less than n, with coefficients in GF(p). Multiplication is done modulo a reducing polynomial of degree n, which is selected at the time of field construction. Inversion is done with extended Euclidean algorithm.
The polynomials are represented by lists of coefficients, highest degrees first. For example, the elements of GF(32) are:
[], [1], [2], [1, 0], [1, 1], [1, 2], [2, 0], [2, 1], [2, 2]
空列表代表0。
GF 类，有限域
Implements arithmetics as methods add, sub, mul, inv (multiplicative inverse). For convenience of testing interpolation includes eval_poly which evaluates a given polynomial with coefficients in GF(pn) at a point of GF(pn).
请注意，构造函数用作 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
插值多项式的构造。
The 拉格朗日多项式是针对列表 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))
and get
Poly(x**2 + (a**2 + a + 1)*x + a, x, domain='ZZ[a]')
这个代数对象不是我们领域的多项式，这只是为了可读的输出。
Sage
作为另一种选择，或者只是另一种安全检查，可以使用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)
Output: x^2 + (a^2 + a + 1)*x + a
