在有限域上插值多项式
Interpolate polynomial over a finite field
我想在有限域的点上使用 python 插值多项式，并在该域中获得具有系数的多项式。
目前我正在尝试使用 SymPy 并专门进行插值（来自 sympy.polys.polyfuncs），但我不知道如何强制插值发生在特定的 gf 中。如果没有，可以用另一个模块来完成吗？
编辑：我对 Python implementation/library 感兴趣。
SymPy 的 interpolating_poly does not support polynomials over finite fields. But there are enough details under the hood of SymPy to put together a class for finite fields, and find the coefficients of Lagrange polynomial 以一种残酷直接的方式。
像往常一样，有限域GF(pn)的元素是小于n的represented by polynomials次，系数在GF(p)中。乘法以 n 次还原多项式为模进行，n 次多项式是在域构造时选择的。反演是用扩展欧几里德算法完成的。
多项式由系数列表表示，最高次数在前。比如GF(32)的元素是：
[], [1], [2], [1, 0], [1, 1], [1, 2], [2, 0], [2, 1], [2, 2]
空列表代表0。
Class GF，有限域
将算术实现为方法 add、sub、mul、inv（乘法逆）。为了方便测试，插值包括 eval_poly ，它在 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
Class 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
构造插值多项式。
Lagrange polynomial 是针对列表 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]')
这个代数对象不是我们领域的多项式，这只是为了可读输出。
贤者
作为一种替代方法，或者只是另一种安全检查，可以使用 Sage 中的 lagrange_polynomial 来获取相同的数据。
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
我想在有限域的点上使用 python 插值多项式，并在该域中获得具有系数的多项式。
目前我正在尝试使用 SymPy 并专门进行插值（来自 sympy.polys.polyfuncs），但我不知道如何强制插值发生在特定的 gf 中。如果没有，可以用另一个模块来完成吗？
编辑：我对 Python implementation/library 感兴趣。
SymPy 的 interpolating_poly does not support polynomials over finite fields. But there are enough details under the hood of SymPy to put together a class for finite fields, and find the coefficients of Lagrange polynomial 以一种残酷直接的方式。
像往常一样，有限域GF(pn)的元素是小于n的represented by polynomials次，系数在GF(p)中。乘法以 n 次还原多项式为模进行，n 次多项式是在域构造时选择的。反演是用扩展欧几里德算法完成的。
多项式由系数列表表示，最高次数在前。比如GF(32)的元素是：
[], [1], [2], [1, 0], [1, 1], [1, 2], [2, 0], [2, 1], [2, 2]
空列表代表0。
Class GF，有限域
将算术实现为方法 add、sub、mul、inv（乘法逆）。为了方便测试，插值包括 eval_poly ，它在 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
Class 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
构造插值多项式。
Lagrange polynomial 是针对列表 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]')
这个代数对象不是我们领域的多项式，这只是为了可读输出。
贤者
作为一种替代方法，或者只是另一种安全检查，可以使用 Sage 中的 lagrange_polynomial 来获取相同的数据。
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