Find inverse polynomial in python in GF2

I'm fairly new to Python and I have a question related to the polynomials. Lets have a high-degree polynomial in GF(2), for example : x^n + x^m + ... + 1, where n, m could be up to 10000. I need to find inverse polynomial to this one. What will be the fastest way to do that in Python (probably using numpy) ?

Thanks

Solution

You have basically two options:

Use sympy.polys.galoistools, Sympy is an excellent python library for symbolic mathematics. However, at least in my experience, if the polynomial has a very high degree (which is common in error-correction coding theory) sympy is very slow.
```
from sympy.polys.galoistools import gf_gcdex, gf_strip

def gf_inv(a, irr_poly):  # irriducible polynomial 

    return gf_gcdex(gf_strip(a), irr_poly, 2 , ZZ)[0]
```

This solution is generical for all GF(2^p) fields e.g. GF(2)/(x^p+1).

If your polynomial has a very high degree sympy is too slow and you have to develop your own Extended Euclidean Algorithm routine and all the other functions required to perform this algorithm.

Personally since i had to invert polynomial of in GF(2) of very high degree (in the ten-thousands !) too i have developed myself the EEA using numpy but only for GF2[x] polynomials:

def gf2_xgcd(b, a):
"""Perform Extended Euclidean Algorithm over GF2

Given polynomials ``b`` and ``a`` in ``GF(p)[x]``, computes polynomials
``s``, ``t`` and ``h``, such that ``h = gcd(f, g)`` and ``s*b + t*a = h``.
The typical application of EEA is solving polynomial diophantine equations and findining multiplicative inverse.


Parameters
----------
b : ndarray (uint8 or bool) or list
    b polynomial's coefficients.
a : ndarray (uint8 or bool) or list
    a polynomial's coefficients.
Returns
-------
y2 : ndarray of uint8
     s polynomial's coefficients.
x2 : ndarray of uint8
     t polynomial's coefficients.
b : ndarray of uint8
    h polynomial's coefficients.

Notes
-----
Rightmost element in the arrays is the leading coefficient of the polynomial.
In other words, the ordering for the coefficients of the polynomials is like the one used in MATLAB while
in Sympy, for example, the leftmost element is the leading coefficient.

Examples
========

>>> x = np.array([1, 1, 1, 1, 1, 0, 1, 0, 1], dtype="uint8")
>>> y = np.array([1, 0, 1], dtype="uint8")
>>> gf2_xgcd(x,y)
(array([0, 1, 1, 1], dtype=uint8),
 array([1, 1], dtype=uint8),
 array([1], dtype=uint8))

"""

x1 = np.array([1], dtype="uint8")
y0 = np.array([1], dtype="uint8")

x0 = np.array([], dtype="uint8")
y1 = np.array([], dtype="uint8")

while True:

    q, r = gf2_div(b, a)

    b = a

    if not r.any():
        break

    a = r

    if not (q.any() and x1.any()):  # if q is zero or x1 is zero
        x2 = x0
    elif not x0.any():  # if x0 is zero
        x2 = mul(x1, q)
    else:
        mulres = mul(x1, q)

        x2 = gf2_add(x0, mulres)

    if not (q.any() and y1.any()):
        y2 = y0
    elif not y0.any():
        y2 = mul(y1, q)
    else:
        mulres = mul(y1, q)

        y2 = gf2_add(y0, mulres)

    # update
    y0 = y1
    x0 = x1
    y1 = y2
    x1 = x2

return y2, x2, b

Of course to perform the EEA one needs firstly to define functions that allow to perform addition, multiplication and division on the GF2[x] field.

In fact i have developed a python module which is available in GitHub: pyGF2 which can efficiently compute polynomial arithmetic over GF2[x] including polynomial inversion orders of magnitude faster than Sympy. The EEA code here is taken directly from that repo.