Euclidean algorithm for polynomials in GF(2^8)

Question

I am trying to create an Euclidean algorithm (to solve Bezout's Relation) for 2 polynomials in the GF(2^8).

I currently have this code for my different operations

class ReedSolomon:   
    gfSize = 256
    genPoly = 285 
    log = [0]*gfSize
    antilog = [0]*gfSize


    def _genLogAntilogArrays(self):
        self.antilog[0] = 1
        self.log[0] = 0
        self.antilog[255] = 1
        for i in range(1,255):
            self.antilog[i] = self.antilog[i-1] << 1 
            if self.antilog[i] >= self.gfSize:
                self.antilog[i] = self.antilog[i] ^ self.genPoly
            self.log[self.antilog[i]] = i

def __init__(self):
        self._genLogAntilogArrays()

def _galPolynomialDivision(self,dividend, divisor):
        result = dividend.copy()
        for i in range(len(dividend) - (len(divisor)-1)):
            coef = result[i]
            if coef != 0:
                for j in range(1, len(divisor)):
                    if divisor[j] != 0: 
                        result[i + j] ^= self._galMult(divisor[j], coef) # équivalent result[i + j] += -divisor[j] * coef car dans un champ GF(2) addition <=> substraction <=> XOR

        remainderIndex = -(len(divisor)-1)
        return result[:remainderIndex], result[remainderIndex:]

def _galMultiplicationPolynomiale(self, x,y):
        result = [0]*(len(x)+len(y)-1)
        for i in range(len(x)):
            for j in range(len(y)):
                result[i+j] ^= self._galMult(x[i],y[j])
        return result

def _galMult(self,x,y):
        if ((x == 0) or (y == 0)):
            val = 0
        else:
            val = self.antilog[(self.log[x] + self.log[y])%255]
        return val

def _galPolynomialAddition(self, a, b):
        polSum = [0] * max(len(a), len(b))
        for index in range(0, len(a)):
            polSum[index + len(polSum) - len(a)] = a[index]
        for index in range(0, len(b)):
            polSum[index + len(polSum) - len(b)] ^= b[index]
        return (polSum)

And here is my euclidean algorithm :

def _galEuclideanAlgorithm(self,a,b):
        r0 = a.copy()
        r1 = b.copy()

        u0 = [1]
        u1 = [0]

        v0 = [0]
        v1 = [1]
        while max(r1) != 0:
            print(r1)
            q,r = self._galPolynomialDivision(r0,r1)
            r0 = self._galPolynomialAddition(self._galMultiplicationPolynomiale(q,r1),r)
            r1,r0 = self._galPolynomialAddition(r0,self._galMultiplicationPolynomiale(q,r1)),r1.copy()
            u1,u0 = self._galPolynomialAddition(u0,self._galMultiplicationPolynomiale(q,u1)),u1.copy()
            v1,v0 = self._galPolynomialAddition(v0,self._galMultiplicationPolynomiale(q,v1)),v1.copy()
        return r1,u1,v1

I don't understand my issue where my algorithm is looping, here is my remainder output with my tests:

rs = ReedSolomon()
a = [1,15,7,8,0,11]

b = [1,0,0,0,0,0,0]

print(rs._galEuclideanAlgorithm(b,a))

#Console output
'''
[1, 15, 7, 8, 0, 11]
[0, 0, 82, 37, 120, 11, 105]
[1, 15, 7, 8, 0, 11]
[0, 0, 82, 37, 120, 11, 105]
[1, 15, 7, 8, 0, 11]
[0, 0, 82, 37, 120, 11, 105]
[1, 15, 7, 8, 0, 11]
'''

I know it might seem like I'm throwing some code just expecting an answer, but I'm genuinely searching for the error.

Thanks in advance !

Answer 1

I created a Python package called galois that does this. galois extends NumPy arrays to operate over Galois fields. The code is written in Python but JIT compiled with Numba for speed. In addition to array arithmetic, it also supports polynomials over Galois fields. ...And Reed-Solomon codes are implemented too :)

The Extended Euclidean Algorithm to solve the Bezout identity for two polynomials in GF(2^8) would be solved this way. Below is an abbreviated chunk of source code. You can see my full source code here.

def poly_egcd(a, b):
    field = a.field
    zero = Poly.Zero(field)
    one = Poly.One(field)

    r2, r1 = a, b
    s2, s1 = one, zero
    t2, t1 = zero, one

    while r1 != zero:
        q = r2 / r1
        r2, r1 = r1, r2 - q*r1
        s2, s1 = s1, s2 - q*s1
        t2, t1 = t1, t2 - q*t1

    # Make the GCD polynomial monic
    c = r2.coeffs[0]  # The leading coefficient
    if c > 1:
        r2 /= c
        s2 /= c
        t2 /= c

    return r2, s2, t2

And here is a complete example using the galois library and the polynomials from your example. (I'm assuming the highest-degree coefficient is first?)

In [1]: import galois                                                                                                                                                                          

In [2]: GF = galois.GF(2**8)                                                                                                                                                                   

In [3]: print(GF.properties)                                                                                                                                                                   
GF(2^8):
  characteristic: 2
  degree: 8
  order: 256
  irreducible_poly: x^8 + x^4 + x^3 + x^2 + 1
  is_primitive_poly: True
  primitive_element: x

In [4]: a = galois.Poly([1,15,7,8,0,11], field=GF); a                                                                                                                                          
Out[4]: Poly(x^5 + 15x^4 + 7x^3 + 8x^2 + 11, GF(2^8))

In [5]: b = galois.Poly([1,0,0,0,0,0,0], field=GF); b                                                                                                                                          
Out[5]: Poly(x^6, GF(2^8))

In [6]: d, s, t = galois.poly_egcd(a, b); d, s, t                                                                                                                                              
Out[6]: 
(Poly(1, GF(2^8)),
 Poly(78x^5 + 7x^4 + 247x^3 + 74x^2 + 152, GF(2^8)),
 Poly(78x^4 + 186x^3 + 45x^2 + x + 70, GF(2^8)))

In [7]: a*s + b*t == d                                                                                                                                                                         
Out[7]: True