在有限域上插值多项式

Question

在有限域上插值多项式

pythonsympypolynomial-mathgalois-fieldfinite-field

5

我想在有限域上使用Python插值多项式，并得到具有该域系数的多项式。目前我正在尝试使用SymPy，特别是插值（来自sympy.polys.polyfuncs），但我不知道如何强制插值发生在特定的gf中。如果不行，是否可以使用另一个模块完成此操作？注：我对Python实现/库很感兴趣。

- Avihu28

2个回答

2

我是 galois Python 库的作者。使用 lagrange_poly() 函数可以进行多项式插值。以下是一个简单的示例。

In [1]: import galois

In [2]: galois.__version__
Out[2]: '0.0.32'

In [3]: GF = galois.GF(3**5)

In [4]: x = GF.Random(10); x
Out[4]: GF([ 33,  58,  59,  21, 141, 133, 207, 182, 125, 162], order=3^5)

In [5]: y = GF.Random(10); y
Out[5]: GF([ 34, 239, 120, 170,  31, 165, 180,  79, 215, 215], order=3^5)

In [6]: f = galois.lagrange_poly(x, y); f
Out[6]: Poly(165x^9 + 96x^8 + 9x^7 + 111x^6 + 40x^5 + 208x^4 + 55x^3 + 17x^2 + 118x + 203, GF(3^5))

In [7]: f(x)
Out[7]: GF([ 34, 239, 120, 170,  31, 165, 180,  79, 215, 215], order=3^5)

有限域元素的显示可以更改为多项式或幂表示形式。

In [8]: GF.display("poly"); f(x)
Out[8]: 
GF([              α^3 + 2α + 1, 2α^4 + 2α^3 + 2α^2 + α + 2,
           α^4 + α^3 + α^2 + α,              2α^4 + 2α + 2,
                   α^3 + α + 1,                   2α^4 + α,
                   2α^4 + 2α^2,       2α^3 + 2α^2 + 2α + 1,
    2α^4 + α^3 + 2α^2 + 2α + 2, 2α^4 + α^3 + 2α^2 + 2α + 2], order=3^5)

In [9]: GF.display("power"); f(x)
Out[9]: 
GF([α^198, α^162, α^116, α^100, α^214, α^137, α^169,  α^95, α^175, α^175],
   order=3^5)

- Matt Hostetter

这个插值函数比普通的Python实现慢20倍，这可能是我做错了什么，但是测试运行结果确实如此。我本来希望它能快20倍呢。 - popeye

@popeye，你能否在 GitHub 上开一个问题（https://github.com/mhostetter/galois/issues），我会看一下？请添加原始的 Python 示例和 Galois 示例。 - Matt Hostetter

当然，我会找时间明天做的。 - popeye

在v0.1.2版本中增加了性能改进。 - Matt Hostetter

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user6655984 · Accepted Answer

SymPy的interpolating_poly不支持有限域上的多项式。但是SymPy底层有足够的细节可以组合成一个有限域类，并以一种非常直接的方式找到Lagrange polynomial的系数。

通常，有限域GF(pⁿ)的元素用次数小于n的多项式表示，其系数在GF(p)中。乘法是模n次降幂多项式的余数进行的，该多项式在构建域时选择。倒数是用扩展欧几里得算法完成的。

这些多项式由系数列表表示，最高次数优先。例如，GF(3²)的元素为：

[], [1], [2], [1, 0], [1, 1], [1, 2], [2, 0], [2, 1], [2, 2]

空列表表示0。

GF类，有限域

实现算术方法包括add、sub、mul、inv（乘法逆元）。为了方便插值测试，还包括eval_poly，它在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

类PolyRing，域上的多项式

这个类更简单：它实现了多项式的加减乘法，在系数上使用了底层域。由于SymPy约定以最高幂次开头来列出单项式，所以有很多列表翻转[::-1]。

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_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

作为另一种选择，或者只是另一个安全检查，可以使用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