Interpolate polynomial over a finite field

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.

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.

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]

The empty list represents 0.

Class GF, finite fields

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ⁿ).

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

Class PolyRing, polynomials over a field

This one is simpler: it implements addition, subtraction, and multiplication of polynomials, referring to the ground field for operations on coefficients. There is a lot of list reversals [::-1] because of SymPy's convention to list monomials starting with highest powers.

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

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

Example of usage:

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

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)

So we can do

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

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

Sage

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)

Output: x^2 + (a^2 + a + 1)*x + a

I'm the author of the galois Python library. Polynomial interpolation can be performed with the lagrange_poly() function. Here's a simple example.

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)

The finite field element display may be changed to either the polynomial or power representation.

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)

Class GF, finite fields

Class PolyRing, polynomials over a field

Construction of interpolating polynomial.

Example of usage:

Sage

Recommended topics

Hot tags