在有限域上插值多项式 [英] Interpolate polynomial over a finite field

问题描述

我想在有限域的点上使用python插值多项式，并在该域中获得具有系数的多项式. 当前，我正在尝试使用SymPy并专门进行插值(来自sympy.polys.polyfuncs)，但是我不知道如何强制在特定gf中进行插值.如果没有，可以用另一个模块来完成吗?

I want to use python interpolate polynomial on points from a finite-field and get a polynomial with coefficients in that field. Currently I'm trying to use SymPy and specifically interpolate (from sympy.polys.polyfuncs), but I don't know how to force the interpolation to happen in a specific gf. If not, can this be done with another module?

我对Python实现/库感兴趣.

I'm interested in a Python implementation/library.

推荐答案

SymPy的拉格朗日多项式以一种残酷的直接方式.

SymPy's 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 in a brutally direct fashion.

通常，有限域GF(p ⁿ )的元素为

As usual, the elements of finite field GF(pⁿ) 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.

多项式由系数列表表示，首先是最高次数.例如，GF(3 ² )的元素是:

The polynomials are represented by lists of coefficients, highest degrees first. For example, the elements of GF(3²) are:

[], [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 ⁿ )评估给定的多项式.

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(pⁿ) at a point of GF(pⁿ).

请注意，构造函数用作G(3，2)，而不用作G(9)-素数及其幂分别提供.

Note that the constructor is used as G(3, 2), not as G(9), - the prime and its power are supplied separately.

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-它是基本多项式的线性组合，每个X的元素对应一个.每个基本多项式通过将(x-x_k)多项式相乘而获得，表示为[[1], K.sub([], x_k)].分母是一个标量，因此它甚至更容易计算.

Construction of interpolating polynomial.

The Lagrange polynomial is constructed for given x-values in list X and corresponding y-values in array Y. It is a linear combination of basis polynomials, one for each element of X. Each basis polynomial is obtained by multiplying (x-x_k) polynomials, represented as [[1], K.sub([], x_k)]. The denominator is a scalar, so it's even easier to compute.

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是多项式中的变量.

The pretty print is just to avoid some type-related decorations on the output. The polynomial is shown as [[1], [1, 1, 1], [1, 0]]. To help readability, I added a function to turn this in a more familiar form, with a symbol a being a generator of finite field, and x being the variable in the polynomial.

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]')

这个代数对象不是我们领域的多项式，这只是为了输出可读性.

This algebraic object is not a polynomial over our field, this is just for the sake of readable output.

作为一种替代选择，或只是另一种安全检查，可以使用

As an alternative, or just another safety check, one can use the lagrange_polynomial from Sage for the same data.

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

这篇关于在有限域上插值多项式的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！