Math Functions in AES and Sage

16/01/19 — capitol


AES, also known as Rijndael, performs its operations in a Galois field of 28. It’s a construct where multiplication, division, addition and subtraction have new meanings and implementations. We will look at how these work and how Sage have made our life easier by implementing them for us.

One way to view the numbers in the field is that they represent polynomials. 0x53 is 01010011 in binary, or the polynomial x6 + x4 + x + 1.

Addition and Subtraction

Since all operations in the field are reduced with modulo 2, addition and subtraction become XOR. This is easiest demonstrated with a few examples borrowed from Wikipedia.

p1 p2 p1 + p2 (normal algebra) p1 + p2 in GF(2n)
x3 + x + 1 x3 + x2 2x3 + x2 + x + 1 x2 + x + 1
x4 + x2 x6 + x2 x6 + x4 + 2x2 x6 + x4
x + 1 x2 + 1 x2 + x + 2 x2 + x
x3 + x x2 + 1 x3 + x2 + x + 1 x3 + x2 + x + 1
x2 + x x2 + x 2x2 + 2x 0

Or seen in Sage:

R.<x> = PolynomialRing(GF(2^8), 'x')
S.<a> = QuotientRing(R, R.ideal(x^8 + x^4 + x^3 + x + 1))

z3 = a^3 + a + 1
z4 = a^3 + a^2
print(z3 + z4)

which prints a^2 + a + 1.


This gets a lot more complicated. If we continue to look at the numbers as polynomials, then all multiplication in GF(28) happens modulo x8 + x4 + x3 + x + 1.

The reason that it’s not a straight modulo operation is that if it were then the multiplication of two numbers could result in 0, which wouldn’t work.

An example is 0x53 * 0xCA = 0x01, since:

(x6 + x4 + x + 1)(x7 + x6 + x3 + x)

= (x13 + x12 + x9 + x7) + (x11 + x10 + x7 + x5) + (x8 + x7 + x4 + x2) + (x7 + x6 + x3 + x)

= x13 + x12 + x11 + x10 + x9 + x8 + x6 + x5 + x4 + x3 + x2 + x

And x13 + x12 + x11 + x10 + x9 + x8 + x6 + x5 + x4 + x3 + x2 + x modulo x8 + x4 + x3 + x + 1 is 1.

We can also see this in Sage:

R.<x> = PolynomialRing(GF(2^8), 'x')
S.<a> = QuotientRing(R, R.ideal(x^8 + x^4 + x^3 + x + 1))

z1 = a^6 + a^4 + a + 1
z2 = a^7 + a^6 + a^3 + a
print(z1 * z2)

which prints 1.


Generally easy, since exponentiation is just repeated multiplication, but some of the numbers in AES have interesting properties. When you multiply those numbers with themselves they traverse all 255 non-zero numbers in the field. They are called generators and here is a list of them:

0x03 0x05 0x06 0x09 0x0b 0x0e 0x11 0x12 0x13 0x14 0x17 0x18 0x19 0x1a 0x1c 0x1e
0x1f 0x21 0x22 0x23 0x27 0x28 0x2a 0x2c 0x30 0x31 0x3c 0x3e 0x3f 0x41 0x45 0x46
0x47 0x48 0x49 0x4b 0x4c 0x4e 0x4f 0x52 0x54 0x56 0x57 0x58 0x59 0x5a 0x5b 0x5f
0x64 0x65 0x68 0x69 0x6d 0x6e 0x70 0x71 0x76 0x77 0x79 0x7a 0x7b 0x7e 0x81 0x84
0x86 0x87 0x88 0x8a 0x8e 0x8f 0x90 0x93 0x95 0x96 0x98 0x99 0x9b 0x9d 0xa0 0xa4
0xa5 0xa6 0xa7 0xa9 0xaa 0xac 0xad 0xb2 0xb4 0xb7 0xb8 0xb9 0xba 0xbe 0xbf 0xc0
0xc1 0xc4 0xc8 0xc9 0xce 0xcf 0xd0 0xd6 0xd7 0xda 0xdc 0xdd 0xde 0xe2 0xe3 0xe5
0xe6 0xe7 0xe9 0xea 0xeb 0xee 0xf0 0xf1 0xf4 0xf5 0xf6 0xf8 0xfb 0xfd 0xfe 0xff

We could write a small Sage program to print a complete loop for the generator 3:

R.<x> = PolynomialRing(GF(2^8), 'x')
S.<a> = QuotientRing(R, R.ideal(x^8 + x^4 + x^3 + x + 1))

z5 = a^0
z6 = a + 1

for i in range(16):
    s = ""
    for j in range(16):
        X1_bitvec = ''.join(map(str, z5.list()));
        X1 = Integer(str(X1_bitvec)[::-1], base=2);
        s += "0x%02X " % (X1)
        z5 *= z6

which prints this:

0x01 0x03 0x05 0x0F 0x11 0x33 0x55 0xFF 0x1A 0x2E 0x72 0x96 0xA1 0xF8 0x13 0x35 
0x5F 0xE1 0x38 0x48 0xD8 0x73 0x95 0xA4 0xF7 0x02 0x06 0x0A 0x1E 0x22 0x66 0xAA 
0xE5 0x34 0x5C 0xE4 0x37 0x59 0xEB 0x26 0x6A 0xBE 0xD9 0x70 0x90 0xAB 0xE6 0x31 
0x53 0xF5 0x04 0x0C 0x14 0x3C 0x44 0xCC 0x4F 0xD1 0x68 0xB8 0xD3 0x6E 0xB2 0xCD 
0x4C 0xD4 0x67 0xA9 0xE0 0x3B 0x4D 0xD7 0x62 0xA6 0xF1 0x08 0x18 0x28 0x78 0x88 
0x83 0x9E 0xB9 0xD0 0x6B 0xBD 0xDC 0x7F 0x81 0x98 0xB3 0xCE 0x49 0xDB 0x76 0x9A 
0xB5 0xC4 0x57 0xF9 0x10 0x30 0x50 0xF0 0x0B 0x1D 0x27 0x69 0xBB 0xD6 0x61 0xA3 
0xFE 0x19 0x2B 0x7D 0x87 0x92 0xAD 0xEC 0x2F 0x71 0x93 0xAE 0xE9 0x20 0x60 0xA0 
0xFB 0x16 0x3A 0x4E 0xD2 0x6D 0xB7 0xC2 0x5D 0xE7 0x32 0x56 0xFA 0x15 0x3F 0x41 
0xC3 0x5E 0xE2 0x3D 0x47 0xC9 0x40 0xC0 0x5B 0xED 0x2C 0x74 0x9C 0xBF 0xDA 0x75 
0x9F 0xBA 0xD5 0x64 0xAC 0xEF 0x2A 0x7E 0x82 0x9D 0xBC 0xDF 0x7A 0x8E 0x89 0x80 
0x9B 0xB6 0xC1 0x58 0xE8 0x23 0x65 0xAF 0xEA 0x25 0x6F 0xB1 0xC8 0x43 0xC5 0x54 
0xFC 0x1F 0x21 0x63 0xA5 0xF4 0x07 0x09 0x1B 0x2D 0x77 0x99 0xB0 0xCB 0x46 0xCA 
0x45 0xCF 0x4A 0xDE 0x79 0x8B 0x86 0x91 0xA8 0xE3 0x3E 0x42 0xC6 0x51 0xF3 0x0E 
0x12 0x36 0x5A 0xEE 0x29 0x7B 0x8D 0x8C 0x8F 0x8A 0x85 0x94 0xA7 0xF2 0x0D 0x17 
0x39 0x4B 0xDD 0x7C 0x84 0x97 0xA2 0xFD 0x1C 0x24 0x6C 0xB4 0xC7 0x52 0xF6 0x01 


When we have this table we can answer the question How many times must I multiply 3 with itself to get number x? by turning it inside out:

G = GF(2^8)
R.<x> = PolynomialRing(G, 'x')
S.<a> = QuotientRing(R, R.ideal(x^8 + x^4 + x^3 + x + 1))

z5 = a^0
z6 = a + 1

log_array = [None for v in range(256)]

for i in range(16):
    for j in range(16):
        X1_bitvec = ''.join(map(str, z5.list()));
        X1 = Integer(str(X1_bitvec)[::-1], base=2);
        z5 *= z6
        log_array[X1] = i * 16 + j;

for i in range(16):
    s = ""
    for j in range(16):
        if log_array[i * 16 + j] == None:
            s += "none "
            s += "0x%02X " % (log_array[i * 16 + j])

Which prints this table:

none 0xFF 0x19 0x01 0x32 0x02 0x1A 0xC6 0x4B 0xC7 0x1B 0x68 0x33 0xEE 0xDF 0x03 
0x64 0x04 0xE0 0x0E 0x34 0x8D 0x81 0xEF 0x4C 0x71 0x08 0xC8 0xF8 0x69 0x1C 0xC1 
0x7D 0xC2 0x1D 0xB5 0xF9 0xB9 0x27 0x6A 0x4D 0xE4 0xA6 0x72 0x9A 0xC9 0x09 0x78 
0x65 0x2F 0x8A 0x05 0x21 0x0F 0xE1 0x24 0x12 0xF0 0x82 0x45 0x35 0x93 0xDA 0x8E 
0x96 0x8F 0xDB 0xBD 0x36 0xD0 0xCE 0x94 0x13 0x5C 0xD2 0xF1 0x40 0x46 0x83 0x38 
0x66 0xDD 0xFD 0x30 0xBF 0x06 0x8B 0x62 0xB3 0x25 0xE2 0x98 0x22 0x88 0x91 0x10 
0x7E 0x6E 0x48 0xC3 0xA3 0xB6 0x1E 0x42 0x3A 0x6B 0x28 0x54 0xFA 0x85 0x3D 0xBA 
0x2B 0x79 0x0A 0x15 0x9B 0x9F 0x5E 0xCA 0x4E 0xD4 0xAC 0xE5 0xF3 0x73 0xA7 0x57 
0xAF 0x58 0xA8 0x50 0xF4 0xEA 0xD6 0x74 0x4F 0xAE 0xE9 0xD5 0xE7 0xE6 0xAD 0xE8 
0x2C 0xD7 0x75 0x7A 0xEB 0x16 0x0B 0xF5 0x59 0xCB 0x5F 0xB0 0x9C 0xA9 0x51 0xA0 
0x7F 0x0C 0xF6 0x6F 0x17 0xC4 0x49 0xEC 0xD8 0x43 0x1F 0x2D 0xA4 0x76 0x7B 0xB7 
0xCC 0xBB 0x3E 0x5A 0xFB 0x60 0xB1 0x86 0x3B 0x52 0xA1 0x6C 0xAA 0x55 0x29 0x9D 
0x97 0xB2 0x87 0x90 0x61 0xBE 0xDC 0xFC 0xBC 0x95 0xCF 0xCD 0x37 0x3F 0x5B 0xD1 
0x53 0x39 0x84 0x3C 0x41 0xA2 0x6D 0x47 0x14 0x2A 0x9E 0x5D 0x56 0xF2 0xD3 0xAB 
0x44 0x11 0x92 0xD9 0x23 0x20 0x2E 0x89 0xB4 0x7C 0xB8 0x26 0x77 0x99 0xE3 0xA5 
0x67 0x4A 0xED 0xDE 0xC5 0x31 0xFE 0x18 0x0D 0x63 0x8C 0x80 0xC0 0xF7 0x70 0x07 

One can notice that the first value is none. That is because multiplying two none-zero numbers never results in zero.


We can use the logarithm table to perform division, because if a = gx and b = gy, then a/b = gx/gy = gx-y = g(x-y) mod 255 where g is one of the generators we talked about earlier.

To take an example, 49 / 11 becomes 347 / 3104 by looking up the numbers 49 and 11 in the logarithm table above. 49 is 0x31 so row 4 column 2 which reads 0x2F or 47 in decimal. 11 is 0x0B so row 1 column 12 which reads 0x68 which is 104 in decimal.

347 / 3104 becomes 347 - 104 mod 255 or 3198. 198 is 0xC6 so by looking in row 13 column 7 of the exponent table we find the value 0x07, which is the answer.

0x07 could also be written as a polynomial: x2 + x + 1

Or we could use the built in division support in Sage:

R.<x> = PolynomialRing(GF(2^8), 'x')
S.<a> = QuotientRing(R, R.ideal(x^8 + x^4 + x^3 + x + 1))

z7 = a^3 + a + 1
z8 = x^5 + a^4 + 1
print(z8 / z7)

Which prints: a^2 + a + 1


By understanding how the basic math functions work we have now laid a foundation to try to understand more advanced analyses of AES, but that will be another blog post.