Интерполировать многочлен по конечному полю

interpolating_poly SymPy не поддерживает полиномы над конечные поля. Но под капотом SymPy достаточно деталей, чтобы собрать класс для конечных полей и найти коэффициенты Лагранжа полиномиальный в грубой прямолинейной манере.

Как обычно, элементы конечного поля 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ⁿ) в точке 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, многочлены над полем

Этот проще: он реализует сложение, вычитание и умножение многочленов, обращаясь к основному полю для операций над коэффициентами. Существует много перестановок списков [::-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