How to perform addition and multiplication in F_{2^8}

Question

I want to perform addition and multiplication in F_{2^8}

I currently have this code which seems to work for add but doesn't work for multiply; the issue seems to be that when I modulo by 100011011 (which represents x^8 + x^4 + x^3 + x + 1), it doesn't seem to do it. Another idea would be to use numpy.polynomial but it isn't as intuitive.

def toBinary(self, n):
    return ''.join(str(1 & int(n) >> i) for i in range(8)[::-1])

def add(self, x, y):
    """
     "10111001" + "10010100" = "00101101"
    """

    if len(x)<8:
        self.add('0'+x,y)
    elif len(y)<8:
        self.add(x,'0'+y)

    try:
        a = int(x,2); b = int(y,2)
        z = int(x)+int(y)
        s = ''
    
        for i in str(z):
            if int(i)%2 == 0:
                s+='0'
            else:
                s+='1'
    except:
        return '00000000'

    return s
    
def multiply(self, x, y):
    """
    "10111001" * "10010100" = "10110010"
    """
    if len(x)<8:
        self.multiply('0'+x,y)
    elif len(y)<8:
        self.multiply(x,'0'+y)

    result = '00000000'

    result = '00000000'
    while y!= '00000000' :
        print(f'x:{x},y:{y},result:{result}')
        if int(y[-1]) == 1 :
            result = self.add(result ,x)
            y = self.add(y, '00000001')
        x = self.add(self.toBinary(int(x,2)<<1),'100011011')
        y = self.toBinary(int(y,2)>>1) #b = self.multiply(b,inverse('00000010'))
    return result

Answer 1

Python example for add (same as subtract), multiply, divide, and inverse. Assumes the input parameters are 8 bit values, and there is no check for divide by 0.

def add(x, y):                  # add is xor
    return x^y

def sub(x, y):                  # sub is xor
    return x^y

def mpy(x, y):                  # mpy two 8 bit values
    p = 0b100011011             # mpy modulo x^8+x^4+x^3+x+1
    m = 0                       # m will be product
    for i in range(8):
        m = m << 1
        if m & 0b100000000:
            m = m ^ p
        if y & 0b010000000:
            m = m ^ x
        y = y << 1
    return m

def div(x, y):                   # divide using inverse
    return mpy(x, inv(y))        #  (no check for y = 0)

def inv(x):                      # x^254 = 1/x
    p=mpy(x,x)                   # p = x^2
    x=mpy(p,p)                   # x = x^4
    p=mpy(p,x)                   # p = x^(2+4)
    x=mpy(x,x)                   # x = x^8
    p=mpy(p,x)                   # p = x^(2+4+8)
    x=mpy(x,x)                   # x = x^16
    p=mpy(p,x)                   # p = x^(2+4+8+16)
    x=mpy(x,x)                   # x = x^32
    p=mpy(p,x)                   # p = x^(2+4+8+16+32)
    x=mpy(x,x)                   # x = x^64
    p=mpy(p,x)                   # p = x^(2+4+8+16+32+64)
    x=mpy(x,x)                   # x = x^128
    p=mpy(p,x)                   # p = x^(2+4+8+16+32+64+128)
    return p

print hex(add(0b01010101, 0b10101010))      # returns 0xff
print hex(mpy(0b01010101, 0b10101010))      # returns 0x59
print hex(div(0b01011001, 0b10101010))      # returns 0x55

---

For GF(2^n), both add and subtract are XOR. This means multiplies are carryless and divides are borrowless. The X86 has a carryless multiply for XMM registers, PCLMULQDQ. Divide by a constant can be done with carryless multiply by 2^64 / constant and using the upper 64 bits of the product. The inverse constant is generated using a loop for borrowless divide.

The reason for this is GF(2^n) elements are polynomials with 1 bit coefficients, (the coefficients are elements of GF(2)).

For GF(2^8), it would be simpler to generate exponentiate and log tables. Example C code:

#define POLY (0x11b)
/* all non-zero elements are powers of 3 for POLY == 0x11b */
typedef unsigned char BYTE;
/* ... */
static BYTE exp2[512];
static BYTE log2[256];
/* ... */
static void Tbli()
{
int i;
int b;
    b = 0x01;                   /* init exp2 table */
    for(i = 0; i < 512; i++){
        exp2[i] = (BYTE)b;
        b = (b << 1) ^ b;       /*  powers of 3 */
        if(b & 0x100)
            b ^= POLY;
    }

    log2[0] = 0xff;             /* init log2 table */
    for(i = 0; i < 255; i++)
        log2[exp2[i]] = (BYTE)i;
}
/* ... */
static BYTE GFMpy(BYTE m0, BYTE m1) /* multiply */
{
    if(0 == m0 || 0 == m1)
        return(0);
    return(exp2[log2[m0] + log2[m1]]);
}
/* ... */
static BYTE GFDiv(BYTE m0, BYTE m1) /* divide */
{
    if(0 == m0)
        return(0);
    return(exp2[log2[m0] + 255 - log2[m1]]);
}

Answer 2

I created a Python package galois that extends NumPy arrays over finite fields. Working with GF(2^8) is quite easy, see my below example.

In [1]: import galois

In [2]: GF = galois.GF(2**8, irreducible_poly="x^8 + x^4 + x^3 + x + 1")

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

# Your original values from your example
In [4]: a = GF(0b10111001); a
Out[4]: GF(185, order=2^8)

In [5]: b = GF(0b10010100); b
Out[5]: GF(148, order=2^8)

In [6]: c = a * b; c
Out[6]: GF(178, order=2^8)

# You can display the result as a polynomial over GF(2)
In [7]: GF.display("poly");

# This matches 0b10110010
In [8]: c
Out[8]: GF(x^7 + x^5 + x^4 + x, order=2^8)

You can work with arrays too.

In [12]: a = GF([1, 2, 3, 4]); a
Out[12]: GF([1, 2, 3, 4], order=2^8)

In [13]: b = GF([100, 110, 120, 130]); b
Out[13]: GF([100, 110, 120, 130], order=2^8)

In [14]: a * b
Out[14]: GF([100, 220, 136,  62], order=2^8)

It's open source, so you can review all the code. Here's a snippet of multiplication in GF(2^m). All of the inputs are integers. Here's how to perform the "polynomial multiplication" using integers with characteristic 2.

    def _multiply_calculate(a, b, CHARACTERISTIC, DEGREE, IRREDUCIBLE_POLY):
        """
        a in GF(2^m), can be represented as a degree m-1 polynomial a(x) in GF(2)[x]
        b in GF(2^m), can be represented as a degree m-1 polynomial b(x) in GF(2)[x]
        p(x) in GF(2)[x] with degree m is the irreducible polynomial of GF(2^m)

        a * b = c
            = (a(x) * b(x)) % p(x) in GF(2)
            = c(x)
            = c
        """
        ORDER = CHARACTERISTIC**DEGREE

        # Re-order operands such that a > b so the while loop has less loops
        if b > a:
            a, b = b, a

        c = 0
        while b > 0:
            if b & 0b1:
                c ^= a  # Add a(x) to c(x)

            b >>= 1  # Divide b(x) by x
            a <<= 1  # Multiply a(x) by x
            if a >= ORDER:
                a ^= IRREDUCIBLE_POLY  # Compute a(x) % p(x)

        return c

The same example runs as follows.

In [72]: _multiply_calculate(0b10111001, 0b10010100, 2, 8, 0b100011011)
Out[72]: 178

In [73]: bin(_multiply_calculate(0b10111001, 0b10010100, 2, 8, 0b100011011))
Out[73]: '0b10110010'