pythonsympypolynomial-mathgalois-fieldfinite-field# Interpolate polynomial over a finite field

### Class GF, finite fields

### Class PolyRing, polynomials over a field

### Construction of interpolating polynomial.

### Example of usage:

### Sage

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?

Edit: I'm interested in a Python implementation/library.

Solution

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^{n}) 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^{2}) are:

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

The empty list represents 0.

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^{n}) at a point of GF(p^{n}).

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

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

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

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.

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`

- Python Jinja2 LaTeX Table
- Getting attributes of a class
- How can I print many significant figures in Python?
- How to allow list append() method to return the new list
- Calculate Last Friday of Month in Pandas
- Python type hint for Iterable[str] that isn't str
- How to iterate over a list in chunks
- How to exit the entire application from a Python thread?
- Running shell command and capturing the output
- How do I pass a variable by reference?
- Convert range(r) to list of strings of length 2 in python
- How can I get the start and end dates for each week?
- how to use send_message() in python-telegram-bot
- Python conditional replacement based on element type
- How can I count the number of items in an arbitrary iterable (such as a generator)?
- Find longest consecutive range of numbers in list
- Insert text in braces with asyncpg
- How does one put a link / url to the web-site's home page in Django?
- How to determine if a path is a subdirectory of another?
- Custom Keybindings for Ipython terminal
- FastAPI asynchronous background tasks blocks other requests?
- How to make sure that information from one file is duplicated into several text documents, without specific lines
- Installing a Python environment with Anaconda
- sklearn pipeline model predicting same results for all input
- Brew command not found after installing Anaconda Python
- How to get an XPath from selenium webelement or from lxml?
- Pipe PuTTY console to Python script
- How to align the axes of a figure in matplotlib?
- Persist ParentDocumentRetriever of langchain
- How to reset index in a pandas dataframe?