Mix Column

The "Mix Column" operates on the rows of the 4x4 state matrix.

It is a non-linear operation that involves multiplying each column by a fixed matrix and reducing the result modulo a fixed polynomial.

The fixed matrice is the following :

+-------------+
| 02 03 01 01 |
| 01 02 03 01 |
| 01 01 02 03 |
| 03 01 01 02 |
+-------------+

Then the operation is the following :

Fixed matrice        User data          Result
+-------------+     +-------------+     +-------------+
| 02 03 01 01 |     | 00 01 02 03 |     | r1 r2 r3 r4 |
| 01 02 03 01 |     | 11 12 13 10 |     | ...         |
| 01 01 02 03 |  *  | 22 23 20 21 |  =  |      ?      |
| 03 01 01 02 |     | 33 30 31 32 |     |             |
+-------------+     +-------------+     +-------------+

Each row of the result are made by the result of the operation of the fixed matrice row times each column : 

r1 = ( 02 * 00 ) ⊕ ( 03 * 11 ) ⊕ ( 01 * 22 ) ⊕ ( 01 * 33 )
r2 = ( 02 * 01 ) ⊕ ( 03 * 12 ) ⊕ ( 01 * 23 ) ⊕ ( 01 * 30 ) 
r3 = ( 02 * 02 ) ⊕ ( 03 * 13 ) ⊕ ( 01 * 20 ) ⊕ ( 01 * 31 )
etc. 

Because of all operation are a Galois Field multiplication, it is not only a times 2 or 3

Only 3 operators are used into this operation :

• 1

• 2

• 3

It make easier to explain a process to calculate the result of the mix column operation, but keep in mind that it look easy but it is not.

Here is only a limited explanation that do not cover the entire mathematics behavior of this maths tool

All operations are made with the binary forms of the number.

• Multiply by 1 is just like standard decimal operation, a * 1 = a

• Multiply by 2

a * 2 = a << 1 (<< is left shift, 1 is the number of shift done, pad on with 0's)

Because we work into GF(2^8), the number can't have more than 8 bits, so if there is more after the bit shift, then the number must be xored with 0b0011011

This value is calculated using the irreducible Polynomial theorem that it will not be explain here.

Here is an example with 0xd4 * 2

d4 = 11010100
d4 * 2 = 11010100 * 10 
       = 110101000
110101000 and 100000000 = 1
d4 * 2 = (110101000 % 100000000) ⊕ 0011011 
       = 10110011

The xor operation is done only because d4 * 02 >= 100000000

• Multiply by 3

3 = 11 (in binary form) 
11 = 10 ⊕ 01

It mean that : 

a * 3 = (bin(a) * 10) ⊕ (bin(a) * 01)

bin(a) * 01 = a * 1
            = a
                        
bin(a) * 10 = a * 2

So : 
a * 3 = a * 2 ⊕ a

Here is an example with 0xd4 * 3

d4 = 11010100
d4 * 2 =  10110011 ( such as explained before )
d4 * 3 = d4 * 2 ⊕ d4 
       = 10110011 ⊕ 11010100
       = 1100111

Python implementation

import copy

def galois_mult(a, b):
    """
    Multiplication in the Galois field GF(2^8).
    """
    p = 0
    hi_bit_set = 0
    for i in range(2):
        if b & 1 == 1: p ^= a
        hi_bit_set = a & 0x80
        a <<= 1
        if hi_bit_set == 0x80: a ^= 0x1b
        b >>= 1
    return p % 256

def mixcolumn(column):
    """
    Mix one column by by considering it as a polynomial and performing
    operations in the Galois field (2^8).
    """
    # XOR is addition in this field
    temp = copy.copy(column) # Store temporary column for operations
    column[0] = galois_mult(temp[0], 2) ^ galois_mult(temp[1], 3) ^ \
                galois_mult(temp[2], 1) ^ galois_mult(temp[3], 1)
    column[1] = galois_mult(temp[0], 1) ^ galois_mult(temp[1], 2) ^ \
                galois_mult(temp[2], 3) ^ galois_mult(temp[3], 1)
    column[2] = galois_mult(temp[0], 1) ^ galois_mult(temp[1], 1) ^ \
                galois_mult(temp[2], 2) ^ galois_mult(temp[3], 3)
    column[3] = galois_mult(temp[0], 3) ^ galois_mult(temp[1], 1) ^ \
                galois_mult(temp[2], 1) ^ galois_mult(temp[3], 2)

def mixcolumns(state):
    """
    Perform a mixing operation which operates on the columns of the states,
    combining the four bytes in each column.
    """
    for i in range(len(state)):
        # Create column from the corresponding array positions
        column = copy.copy(state[i])

        # Mix the extracted column
        mixcolumn(column)

        # Set the new column in the state
        state[i] = column
    return state

Optimized

# learned from http://cs.ucsb.edu/~koc/cs178/projects/JT/aes.c
xtime = lambda a: (((a << 1) ^ 0x1B) & 0xFF) if (a & 0x80) else (a << 1)


def mix_single_column(a):
    # see Sec 4.1.2 in The Design of Rijndael
    t = a[0] ^ a[1] ^ a[2] ^ a[3]
    u = a[0]
    a[0] ^= t ^ xtime(a[0] ^ a[1])
    a[1] ^= t ^ xtime(a[1] ^ a[2])
    a[2] ^= t ^ xtime(a[2] ^ a[3])
    a[3] ^= t ^ xtime(a[3] ^ u)
    return a


def mix_columns(s):
    for i in range(4):
        s[i] = mix_single_column(s[i])
    return s


def inv_mix_columns(s):
    # see Sec 4.1.3 in The Design of Rijndael
    for i in range(4):
        u = xtime(xtime(s[i][0] ^ s[i][2]))
        v = xtime(xtime(s[i][1] ^ s[i][3]))
        s[i][0] ^= u
        s[i][1] ^= v
        s[i][2] ^= u
        s[i][3] ^= v
    print(s)
    return mix_columns(s)