AES in Mathematica

Disclaimer I: The following code is for educational purposes only. Its sole function is to show how the AES algorithm works. It is not by any means intended for real life applications. Do not use it to encrypt your personal data! If you fiddle around and break something important, this is your own fault.

Disclaimer II: If you are interested in following the code examples presented here, you will need a working version of Wolfram Mathematica. In case you don’t have one, there is a 30-day trial available at http://www.wolfram.com/mathematica/trial. If you own a Raspberry Pi, there are operating systems that ship with Mathematica, e.g. Raspbian, which will work perfectly fine.

What we are going to do

Using Mathematica, I will implement the AES (Rijndael) encryption working on a single block of 128 bits using a 256 bit Key. This is a prototypic, laboratory-like implementation of the algorithm. It is designed for conceptual purposes and addresses beginners and non-professionals.

What this is not

This implementation, though mathematically correct, is not recommended for real life applications. It lacks performance (I will embrace improvements submitted to me) and is probably vulnerable to various side channel attacks. Again, use it only to encrypt non-sensitive data! It is not suited for production process.

Overview

I will present the algorithm “chronologically”, i.e. in the sense of the actual workflow. Concepts are introduced “on the fly” as we will need them. If you are new to the subject, it is highly recommended to read this post from the beginning to the end. Impatiently jumping to the middle of the document and working your way from there will give you a hard time. If you read through all the warnings above and are still interested, I will now continue and stop bothering you with further disclaimers, assuming you know what you’re doing.

The first fact that you need to absorb is that AES (with the parameters described above) will operate in 15 rounds. Sometimes people refer to 14 rounds. I never understood how the counting works. For each round we will need a round key of the same size as our plain text, i.e, 128 bits or 16 bytes. The round keys are calculated from the 256 bit (or 32 byte) key, the “password”, if you will. Well, “password” is the wrong name here, but for the sake of simplicity you can imagine that it is something like that. If you are interested you need to search for “key derivation functions” to get some insight on how a password is related to a key. Whatever, let’s stick to notion of a “key” and show how the round keys are obtained in the next section.

Key expansion

Let us start with an example of a 32 byte key. If it were a password, I have some advice how to choose a secure one. For definiteness, we type into Mathematica

key = "66OlSO8L7KoW44awcg2xHJ9X1FbOoF4z"

It’s more convenient to work with numbers rather than letters, so we need the ASCII character codes, which simply get by

keyCC = ToCharacterCode[key]

The resulting output (and also the value of keyCC) is a sequence of 32 numbers

{54, 54, 79, 108, 83, 79, 56, 76, 55, 75, 111, 87, 52, 52, 97, 119, 99, 103, 50, 120, 72, 74, 57, 88, 49, 70, 98, 79, 111, 70, 52, 122}

If you want to view it in hexadecimal representation type BaseForm[keyCC, 16], but we wont need that. As I have told you before, the key needs to be expand in 15 chunks of 16 bytes. For convenience we create a table with 4 lines and 60 columns:

keytable = ConstantArray[0, {4, 60}]

The first round key is simply the key itself, and it is written “top-down then left-right” into
keytable.

For[j = 1, j <= 8, j++,
  For[i = 1, i <= 4, i++,
    keytable[[i]][[j]] = keyCC[[4*(j-1)+i]]
  ]
]

The For-loops work exactly as in your usual programming languages. Well, arrays start with index 1 instead of 0, but that’s nothing to be afraid of.
If you want to view the result,

keytable //MatrixForm
 54, 83,  55,  52,  99, 72, 49, 111, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
 54, 79,  75,  52, 103, 74, 70,  70, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
 79, 56, 111,  97,  50, 57, 98,  52, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
108, 76,  87, 119, 120, 88, 79, 122, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0

That was easy, right? Now we need to fill the remaining 52 rows of this table. This is “usually” (with a few exceptions described below) done for a given row by taking the row directly left of it and the row 8 places left of it. The latter two are XORed (entry by entry) and the result gives the values for the row in question. For example tenth row = ninth row XOR second row. Well, we need to have the ninth row first, for which a special treatment is needed. So now let’s talk about the exceptions:

  1. For the 9th, 17th, 25th, …, 57th row, i.e. those with an index of the form 9+multiple of 8: let’s say the row in question has position i. Then
    ith row = (i-8)th row XOR sbox[RotateLeft[(i-1)th row]] XOR Rcon[(i-1)/8]

    Shit just got real! You should know what XOR does. I tell you in a second what sbox, RotateLeft and Rcon do.
  2. For the 13th, 21st, 29th, …, 53rd row, i.e. those with an index of the form 9+multiple of 4 and not in the above category: let’s call the row index i again. In contrast to the usual key schedule
    ith row = (i-8)th row XOR sbox((i-1)th row)
    I.e. we beat the left neighbor row with an additional sbox.

Now I need to supply definitions. The sbox is not a new gaming console. It is called a substitution box that interchanges bytes according to a specified table. For example, it replaces 0 by 99 (decimal), 1 by 124. The complete implementation is

sboxsub = {99, 124, 119, 123, 242, 107, 111, 197, 48, 1, 103, 43, 254, 215, 171, 118, 202, 130, 201, 125, 250, 89, 71, 240, 173, 212, 162, 175, 156, 164, 114, 192, 183, 253, 147, 38, 54, 63, 247, 204, 52, 165, 229, 241, 113, 216, 49, 21, 4, 199, 35, 195, 24, 150, 5, 154, 7, 18, 128, 226, 235, 39, 178, 117, 9, 131, 44, 26, 27, 110, 90, 160, 82, 59, 214, 179, 41, 227, 47, 132, 83, 209, 0, 237, 32, 252, 177, 91, 106, 203, 190, 57, 74, 76, 88, 207, 208, 239, 170, 251, 67, 77, 51, 133, 69, 249, 2, 127, 80, 60, 159, 168, 81, 163, 64, 143, 146, 157, 56, 245, 188, 182, 218, 33, 16, 255, 243, 210, 205, 12, 19, 236, 95, 151, 68, 23, 196, 167, 126, 61, 100, 93, 25, 115, 96, 129, 79, 220, 34, 42, 144, 136, 70, 238, 184, 20, 222, 94, 11, 219, 224, 50, 58, 10, 73, 6, 36, 92, 194, 211, 172, 98, 145, 149, 228, 121, 231, 200, 55, 109, 141, 213, 78, 169, 108, 86, 244, 234, 101, 122, 174, 8, 186, 120, 37, 46, 28, 166, 180, 198, 232, 221, 116, 31, 75, 189, 139, 138, 112, 62, 181, 102, 72, 3, 246, 14, 97, 53, 87, 185, 134, 193, 29, 158, 225, 248, 152, 17, 105, 217, 142, 148, 155, 30, 135, 233, 206, 85, 40, 223, 140, 161, 137, 13, 191, 230, 66, 104, 65, 153, 45, 15, 176, 84, 187, 22}
sbox[i_] := sboxsub[[i+1]

Next, RotateLeft simply cuts the first element out of the row and appends it to the end. To be clear, after the rotation the former second element stands now on the first place, the former third is now second etc. Since we are talking about rows, it should rather be called RotateUp, or something like that. Well, RotateLeft is a built-in symbol of Mathematica, so we stick to it.

Finally, Rcon returns a row whose only nonzero component is the first one and contains powers of 2:

rcontab = {{1, 0, 0, 0}, {2, 0, 0, 0}, {4, 0, 0, 0}, {8, 0, 0, 0}, {16, 0, 0, 0}, {32, 0, 0, 0}, {64, 0, 0, 0}}
Rcon[i_] := rcontab[[i]]

For AES-256 it is sufficient to have 7 “Rcon-vectors”. There is some more mathematical background about polynomials over a finite field, but I leave that for you to read and it’s not terribly complicated. Alright. Let’s walk through an example.

Key expansion example and round keys

Remember, the first 8 rows were for free. So we need to start with row number 9, which needs exceptional treatment. According to the above formula, we need Rcon[1]={1, 0, 0, 0} and

RotateLeft[8th row] = RotateLeft[{111, 70, 52, 122}] = {70, 52, 122, 111}

from which we obtain

sbox({70, 52, 122, 111}) = {90, 24, 218, 168}

And then XORing with the first row and Rcon[1]

{54, 54, 79, 108} XOR {90, 24, 218, 168} XOR {1, 0, 0, 0} = {109, 46, 148, 196}

Ok, the 9th row is {109, 46, 148, 196}.

Now the 10th row. It is non-exceptional and therefore

{109, 46, 148, 196} XOR {83, 79, 56, 76} = {62, 97, 173, 136}

The 11th and 12th ones are also regular and given by

{62, 97, 173, 136} XOR {55, 75, 111,  87} = { 9, 42, 194, 223}
{ 9, 42, 194, 223} XOR {52, 52,  97, 119} = {61, 30, 163, 168}

respectively.

For the 13th we encounter an exception of type “2”

sbox({61, 30, 163, 168}) XOR {99, 103, 50, 120} = {39, 114, 10, 194} XOR {99, 103, 50, 120} = {68, 21, 56, 186}

The 14-16th are again regular, the 17th is an exception of type “1” etc. Got the idea? good. So lets make Mathematica do it for us,

Clear[i, t];
For[i = 9, i <= 60, i++,
  t = keytable[[1;;4, i-1]];
  If[Mod[i-1, 8] == 0,
    t = BitXor[sbox/@(RotateLeft[t]), Rcon[(i-1)/8]];
    ,
    If[(Mod[i-1, 8] == 4), t = sbox/@t];
  ];
  t = BitXor[keytable[[1;;4, i-8]], t];
  keytable[[1, i]] = t[[1]];
  keytable[[2, i]] = t[[2]];
  keytable[[3, i]] = t[[3]];
  keytable[[4, i]] = t[[4]];
]

The whole result reads

keytable //MatrixForm
 54, 83,  55,  52,  99, 72, 49, 111, 109,  62,   9,  61,  68,  12,  61,  82, 160, 158, 151, 170, 232, 228, 217, 139,  78, 208,  71, 237, 189,  89, 128,  11, 123, 171, 236,   1, 193, 152,  24,  19, 196, 111, 131, 130, 210,  74,  82,  65, 235, 132,   7, 133,  69,  15,  93,  28, 203,  79, 72, 205
 54, 79,  75,  52, 103, 74, 70,  70,  46,  97,  42,  30,  21,  95,  25,  95, 117,  20,  62,  32, 162, 253, 228, 187, 196, 208, 238, 206,  41, 212,  48, 139, 125, 173,  67, 141, 116, 160, 144,  27, 126, 211, 144,  29, 208, 112, 224, 251,  85, 134,  22,  11, 251, 139, 107, 144, 197,  67, 85,  94
 79, 56, 111,  97,  50, 57, 98,  52, 149, 173, 194, 163,  56,   1,  99,  87, 155,  54, 244,  87,  99,  98,   1,  86,  33,  23, 227, 180, 238, 140, 141, 219, 160, 183,  84, 224,  15, 131,  14, 213, 164,  19,  71, 167,  83, 208, 222,  11, 236, 255, 184,  31, 147,  67, 157, 150, 125, 130, 58,  37
108, 76,  87, 119, 120, 88, 79, 122, 196, 136, 223, 168, 186, 226, 173, 215, 196,  76, 147,  59,  88, 186,  23, 192, 249, 181,  38,  29, 252,  70,  81, 145, 210, 103,  65,  92, 182, 240, 161,  48, 175, 200, 137, 213, 181,  69, 228, 212,  44, 228, 109, 184, 217, 156, 120, 172, 176,  84, 57, 129

OK, now this table contains our 16-byte round keys successively from top to bottom and left to right (in that order). We can conveniently address them by

keylistflat = Flatten@Transpose@keytable
RoundKey[i_] := keylistflat[[1+(i-1)*16;;i*16]]

Then for example

RoundKey[1]
RoundKey[2]
RoundKey[3]

gives

{ 54,  54,  79, 108, 83, 79,  56,  76, 55, 75, 111,  87,  52, 52,  97, 119}
{ 99, 103,  50, 120, 72, 74,  57,  88, 49, 70,  98,  79, 111, 70,  52, 122}
{109,  46, 149, 196, 62, 97, 173, 136,  9, 42, 194, 223,  61, 30, 163, 168}

There are 15 of them as mentioned.

Encryption

Now that we have the round keys, we can start the AES algorithm to encrypt a message block of 128 bits or 16 bytes. Let’s pick one and convert it to numbers

mes = "gooby has cansur";
mesCC = ToCharacterCode[mes]
{103, 111, 111, 98, 121, 32, 104, 97, 115, 32, 99, 97, 110, 115, 117, 114}

Now the AES algorithm does some whacky stuff to this message. We will call the result after each round a “state”. The first state or the first round is easy,

State[1] = BitXor[mesCC, RoundKey[1]]
{81, 89, 32, 14, 42, 111, 80, 45, 68, 107, 12, 54, 90, 71, 20, 5}

Rounds 2-14 are more involved. We need some additional definitions

ShiftRows[state_] :=  Function[x, state[[x]]]/@{1, 6, 11, 16, 5, 10, 15, 4, 9, 14, 3, 8, 13, 2, 7, 12}

does some funny permutation of the state.

Next we need to define an odd looking multiplication, called GaloisMultiply. We will define it in terms of simple prototypes. If I have number x, then the Galois multiplication by 1 is simply x. That sounds reasonable, right? If I have number x which is less than 128, then the Galois multiplication by 2 is simply 2x. If x happens to be greater than or equal to 128, then the Galois multiplication by 2 is defined as 2x XOR 283. In doing so, the highest bit of 2x which would exceed the one-byte range is removed by this and the rest is scrambled with the other bits of 283. Now GaloisMultiply by 3, 5, 7 etc. is defined as 3x = x XOR GaloisMultiply[2,x], 5x = x XOR GaloisMultiply[4,x], 7x = x XOR GaloisMultiply[6,x] and so on, where GaloisMultiply[4,x] = GaloisMultiply[2, GaloisMultiply[2, x]], GaloisMultiply[6,x] = GaloisMultiply[2, GaloisMultiply[3, x]]. We will need that up to multiplications by 14. One can build a lookup table for this. In Mathematica we simply use the concept of memorizing functions:

Clear[GaloisMultiply]
GreaterEqual128[x_] = x >= 128;
Less128[x_] = x < 128;
GaloisMultiply[1, x_] := x;
GaloisMultiply[2, x_?GreaterEqual128] := GaloisMultiply[2, x] = BitXor[2 x, 283];
GaloisMultiply[2, x_?Less128] := GaloisMultiply[2, x] = 2 x;
GaloisMultiply[3, x_] := GaloisMultiply[3, x] = BitXor[x, GaloisMultiply[2, x]];
GaloisMultiply[4, x_] := GaloisMultiply[2, GaloisMultiply[2, x]];
GaloisMultiply[5, x_] := BitXor[x, GaloisMultiply[4, x]];
GaloisMultiply[6, x_] := GaloisMultiply[2, GaloisMultiply[3, x]];
GaloisMultiply[7, x_] := BitXor[x, GaloisMultiply[6, x]];
GaloisMultiply[8, x_] := GaloisMultiply[2, GaloisMultiply[2, GaloisMultiply[2, x]]];
GaloisMultiply[9, x_] := GaloisMultiply[9, x] = BitXor[x, GaloisMultiply[8, x]];
GaloisMultiply[10, x_] := GaloisMultiply[2, GaloisMultiply[5, x]];
GaloisMultiply[11, x_] := GaloisMultiply[11, x] = BitXor[x, GaloisMultiply[10, x]];
GaloisMultiply[12, x_] := GaloisMultiply[2, GaloisMultiply[6, x]];
GaloisMultiply[13, x_] := GaloisMultiply[13, x] = BitXor[x, GaloisMultiply[12, x]];
GaloisMultiply[14, x_] := GaloisMultiply[14, x] = GaloisMultiply[2, GaloisMultiply[7, x]];

I invite you to play around with this definition and find out the properties of this product.

Next, we define an operation called MixColumns which implements a matrix multiplication in the Galois sense. The state will be divided in 4 chunks of 4 bytes each, and on each of these chunks we act with a matrix multiplication in the following form

MixColumn[col_] := { BitXor[GaloisMultiply[2, col[[1]]], GaloisMultiply[3, col[[2]]], col[[3]], col[[4]]],
                     BitXor[col[[1]], GaloisMultiply[2, col[[2]]], GaloisMultiply[3, col[[3]]], col[[4]]],
                     BitXor[col[[1]], col[[2]], GaloisMultiply[2, col[[3]]], GaloisMultiply[3, col[[4]]]],
                     BitXor[GaloisMultiply[3, col[[1]]], col[[2]], col[[3]], GaloisMultiply[2, col[[4]]]] }

If one interprets XOR as “+”, then the first line would read something like 2*col[[1]]+3*col[[2]]+col[[3]]+*col[[4]], which is just the first component of your usual matrix-vector multiplication in four dimensions. As I said this is done in blocks to the whole state:

MixColumns[state_] := Flatten@{ MixColumn[state[[1;;4]]], MixColumn[state[[5;;8]]], MixColumn[state[[9;;12]]], MixColumn[state[[13;;16]]] }

Now finally, for i= 2,…,14

State[i_] := BitXor[ MixColumns@(ShiftRows@(sbox/@State[i-1])), RoundKey[i] ]

and the last round is without column mixing

State[15] := BitXor[(ShiftRows@(sbox/@State[14])), RoundKey[15] ]

…aaand we’re done. State[15] is the cipher key! FYI

State[15]
{205, 95, 203, 120, 35, 143, 230, 63, 225, 53, 189, 238, 178, 44, 132, 203}

Encryption example

The previous section contained a lot of material, which you have to absorb. I will walk you through our accompanying example. So where were we? We had State[1]={81, 89, 32, 14, 42, 111, 80, 45, 68, 107, 12, 54, 90, 71, 20, 5}. To obtain the 2nd state we need to  sbox the first one, which yields sbox(State[1])={209, 203, 183, 171, 229, 168, 83, 216, 27, 127, 254, 5, 190, 160, 250, 107}. This is followed by the ShiftRows operation which gives {209, 168, 254, 107, 229, 127, 250, 171, 27, 160, 183, 216, 190, 203, 83, 5}. After that comes the MixColumns, resulting in {207, 232, 35, 232, 1, 165, 147, 252, 162, 90, 189, 145, 119, 195, 220, 75}. Last comes the XOR with the second round key from above. The result is

State[2]
{172, 143, 17, 144, 73, 239, 170, 164, 147, 28, 223, 222, 24, 133, 232, 49}

Perfect! Only thirteen more rounds to go. One arrives at

State[14]
{165, 30, 234, 5, 184, 55, 241, 13, 183, 39, 121, 92, 107, 144, 140, 177}

The result of the final round is quoted above!

Encryption code

Let us collect all relevant code in this section and generalize it to arbitrary keys and message blocks.

sboxsub = {99, 124, 119, 123, 242, 107, 111, 197, 48, 1, 103, 43, 254, 215, 171, 118, 202, 130, 201, 125, 250, 89, 71, 240, 173, 212, 162, 175, 156, 164, 114, 192, 183, 253, 147, 38, 54, 63, 247, 204, 52, 165, 229, 241, 113, 216, 49, 21, 4, 199, 35, 195, 24, 150, 5, 154, 7, 18, 128, 226, 235, 39, 178, 117, 9, 131, 44, 26, 27, 110, 90, 160, 82, 59, 214, 179, 41, 227, 47, 132, 83, 209, 0, 237, 32, 252, 177, 91, 106, 203, 190, 57, 74, 76, 88, 207, 208, 239, 170, 251, 67, 77, 51, 133, 69, 249, 2, 127, 80, 60, 159, 168, 81, 163, 64, 143, 146, 157, 56, 245, 188, 182, 218, 33, 16, 255, 243, 210, 205, 12, 19, 236, 95, 151, 68, 23, 196, 167, 126, 61, 100, 93, 25, 115, 96, 129, 79, 220, 34, 42, 144, 136, 70, 238, 184, 20, 222, 94, 11, 219, 224, 50, 58, 10, 73, 6, 36, 92, 194, 211, 172, 98, 145, 149, 228, 121, 231, 200, 55, 109, 141, 213, 78, 169, 108, 86, 244, 234, 101, 122, 174, 8, 186, 120, 37, 46, 28, 166, 180, 198, 232, 221, 116, 31, 75, 189, 139, 138, 112, 62, 181, 102, 72, 3, 246, 14, 97, 53, 87, 185, 134, 193, 29, 158, 225, 248, 152, 17, 105, 217, 142, 148, 155, 30, 135, 233, 206, 85, 40, 223, 140, 161, 137, 13, 191, 230, 66, 104, 65, 153, 45, 15, 176, 84, 187, 22}
sbox[i_] := sbox[i] = sboxsub[[i+1]];

rcontab = {{1, 0, 0, 0}, {2, 0, 0, 0}, {4, 0, 0, 0}, {8, 0, 0, 0}, {16, 0, 0, 0}, {32, 0, 0, 0}, {64, 0, 0, 0}}
Rcon[i_] := Rcon[i] = rcontab[[i]];

KeyExpand256[k_String] := KeyExpand256[k] = Module[{kcc, kt, i, j, t, t2, t3},
  kcc = ToCharacterCode[k];
  kt = ConstantArray[0, {4, 60}];
  For[j = 1, j <= 8, j++,
    For[i = 1, i <= 4, i++,
      kt[[i]][[j]] = kcc[[4*(j-1)+i]]
    ]
  ]; 
  For[i = 9, i <= 60, i++,
    t = kt[[1;;4, i-1]];
    t2 = t;
    If[Mod[i-1, 8] == 0,
      t2 = BitXor[sbox/@(RotateLeft[t]), Rcon[(i-1)/8]];
    ,
      If[(Mod[i-1, 8] == 4), t3 = sbox/@t2; t2 = t3;];
    ];
    t = BitXor[kt[[1;;4, i-8]], t2];
    kt[[1, i]] = t[[1]];
    kt[[2, i]] = t[[2]];
    kt[[3, i]] = t[[3]];
    kt[[4, i]] = t[[4]];
  ];
  Flatten[Transpose[kt]]
]

RoundKey[i_, list_] := list[[1+(i-1)*16;;i*16]]

ShiftRows[state_] :=  Function[x, state[[x]]]/@{1, 6, 11, 16, 5, 10, 15, 4, 9, 14, 3, 8, 13, 2, 7, 12}

GreaterEqual128[x_] = x >= 128;
Less128[x_] = x < 128;
GaloisMultiply[1, x_] := x;
GaloisMultiply[2, x_?GreaterEqual128] := GaloisMultiply[2, x] = BitXor[2 x, 283];
GaloisMultiply[2, x_?Less128] := GaloisMultiply[2, x] = 2 x;
GaloisMultiply[3, x_] := GaloisMultiply[3, x] = BitXor[x, GaloisMultiply[2, x]];
GaloisMultiply[4, x_] := GaloisMultiply[2, GaloisMultiply[2, x]]; GaloisMultiply[5, x_] := BitXor[x, GaloisMultiply[4, x]];
GaloisMultiply[6, x_] := GaloisMultiply[2, GaloisMultiply[3, x]]; GaloisMultiply[7, x_] := BitXor[x, GaloisMultiply[6, x]];
GaloisMultiply[8, x_] := GaloisMultiply[2, GaloisMultiply[2, GaloisMultiply[2, x]]];
GaloisMultiply[9, x_] := GaloisMultiply[9, x] = BitXor[x, GaloisMultiply[8, x]]; GaloisMultiply[10, x_] := GaloisMultiply[2, GaloisMultiply[5, x]];
GaloisMultiply[11, x_] := GaloisMultiply[11, x] = BitXor[x, GaloisMultiply[10, x]];
GaloisMultiply[12, x_] := GaloisMultiply[2, GaloisMultiply[6, x]];
GaloisMultiply[13, x_] := GaloisMultiply[13, x] = BitXor[x, GaloisMultiply[12, x]];
GaloisMultiply[14, x_] := GaloisMultiply[14, x] = GaloisMultiply[2, GaloisMultiply[7, x]];

MixColumn[col_] := { BitXor[GaloisMultiply[2, col[[1]]], GaloisMultiply[3, col[[2]]], col[[3]], col[[4]]],
                     BitXor[col[[1]], GaloisMultiply[2, col[[2]]], GaloisMultiply[3, col[[3]]], col[[4]]],
                     BitXor[col[[1]], col[[2]], GaloisMultiply[2, col[[3]]], GaloisMultiply[3, col[[4]]]],
                     BitXor[GaloisMultiply[3, col[[1]]], col[[2]], col[[3]], GaloisMultiply[2, col[[4]]]] }
MixColumns[state_] := Flatten@{ MixColumn[state[[1;;4]]], MixColumn[state[[5;;8]]], MixColumn[state[[9;;12]]], MixColumn[state[[13;;16]]] }

State[i_, prev_, list_] := BitXor[ MixColumns@(ShiftRows@(sbox/@prev)), RoundKey[i, list] ]

Is256Key[k_] := (Head[k] == String) && (StringLength[k] == 32);
Is128Block[b_] := (Head[b] == String) && (StringLength[b] == 16);

aes[m_?Is128Block, k_?Is256Key] := aes[m, k] = Module[{mcc, kcc, kt, kt2, i, j, t, t2, t3, state, state2},
  mcc = ToCharacterCode[m];
  kt = KeyExpand256[k];
  state = BitXor[mcc, RoundKey[1, kt]];
  state2 = State[2, state, kt];
  state = State[3, state2, kt];
  state2 = State[4, state, kt];
  state = State[5, state2, kt];
  state2 = State[6, state, kt];
  state = State[7, state2, kt];
  state2 = State[8, state, kt];
  state = State[9, state2, kt];
  state2 = State[10, state, kt];
  state = State[11, state2, kt];
  state2 = State[12, state, kt];
  state = State[13, state2, kt];
  state2 = State[14, state, kt];
  state = BitXor[(ShiftRows@(sbox/@state2)), RoundKey[15, kt] ];
  FromCharacterCode[state]
]

We took the opportunity to define a function aes[m, k] which spits out the AES cipher for a message m encrypted under the key k, both formatted as strings. In our example

aes["gooby has cansur", "66OlSO8L7KoW44awcg2xHJ9X1FbOoF4z"]
Í_Ëx#æ?á5½î\.b2,„Ë

Cool! We have turned the plain text into gibberish. That is what encryption is roughly about! Now we need to address the reverse scenario: turning the gibberish into plain text, occasionally known as ‘decryption’.

Decryption

Decryption is easy if you understood encryption. One basically needs to undo the steps that led to the cipher text.

  • The round keys are obtained just as before, they will just be applied in reverse order, i.e. starting with the 15th.
  • Certainly we need the inverse sbox. You can copy it from the code I provide below, or use Mathematica to find it: (Transpose@SortBy[Transpose@({sboxsub, Sort@sboxsub}), First]) [[2]]
  • The inverse of ShiftRows can be easily found.
  • For MixColumn(s) we need an inverse matrix. It is implemented below.

Let’s do it:

dsboxsub = {82, 9, 106, 213, 48, 54, 165, 56, 191, 64, 163, 158, 129, 243, 215, 251, 124, 227, 57, 130, 155, 47, 255, 135, 52, 142, 67, 68, 196, 222, 233, 203, 84, 123, 148, 50, 166, 194, 35, 61, 238, 76, 149, 11, 66, 250, 195, 78, 8, 46, 161, 102, 40, 217, 36, 178, 118, 91, 162, 73, 109, 139, 209, 37, 114, 248, 246, 100, 134, 104, 152, 22, 212, 164, 92, 204, 93, 101, 182, 146, 108, 112, 72, 80, 253, 237, 185, 218, 94, 21, 70, 87, 167, 141, 157, 132, 144, 216, 171, 0, 140, 188, 211, 10, 247, 228, 88, 5, 184, 179, 69, 6, 208, 44, 30, 143, 202, 63, 15, 2, 193, 175, 189, 3, 1, 19, 138, 107, 58, 145, 17, 65, 79, 103, 220, 234, 151, 242, 207, 206, 240, 180, 230, 115, 150, 172, 116, 34, 231, 173, 53, 133, 226, 249, 55, 232, 28, 117, 223, 110, 71, 241, 26, 113, 29, 41, 197, 137, 111, 183, 98, 14, 170, 24, 190, 27, 252, 86, 62, 75, 198, 210, 121, 32, 154, 219, 192, 254, 120, 205, 90, 244, 31, 221, 168, 51, 136, 7, 199, 49, 177, 18, 16, 89, 39, 128, 236, 95, 96, 81, 127, 169, 25, 181, 74, 13, 45, 229, 122, 159, 147, 201, 156, 239, 160, 224, 59, 77, 174, 42, 245, 176, 200, 235, 187, 60, 131, 83, 153, 97, 23, 43, 4, 126, 186, 119, 214, 38, 225, 105, 20, 99, 85, 33, 12, 125};
dsbox[i_] := dsbox[i] = dsboxsub[[i+1]];

dShiftRows[state_] := Function[x, state[[x]]] /@ {1, 14, 11, 8, 5, 2, 15, 12, 9, 6, 3, 16, 13, 10, 7, 4}

dMixColumn[col_] := { BitXor[GaloisMultiply[14, col[[1]]],  GaloisMultiply[11, col[[2]]], GaloisMultiply[13, col[[3]]], GaloisMultiply[9, col[[4]]]],
                      BitXor[GaloisMultiply[9, col[[1]]], GaloisMultiply[14, col[[2]]], GaloisMultiply[11, col[[3]]], GaloisMultiply[13, col[[4]]]],
                      BitXor[GaloisMultiply[13, col[[1]]], GaloisMultiply[9, col[[2]]], GaloisMultiply[14, col[[3]]], GaloisMultiply[11, col[[4]]]],
                      BitXor[GaloisMultiply[11, col[[1]]], GaloisMultiply[13, col[[2]]], GaloisMultiply[9, col[[3]]], GaloisMultiply[14, col[[4]]]]
                     }
dMixColumns[state_] := Flatten@{ dMixColumn[state[[1;;4]]], dMixColumn[state[[5;;8]]], dMixColumn[state[[9;;12]]], dMixColumn[state[[13;;16]]] }

dState[i_, prev_, list_] := dsbox/@(dShiftRows@(dMixColumns@(BitXor[prev, RoundKey[i, list]])))

daes[m_?Is128Block, k_?Is256Key] := daes[m, k] = Module[{mcc, kt, kt2, i, j, t, t2, t3, state, state2},
  mcc = ToCharacterCode[m];
  kt = KeyExpand256[k]; 
  state = dsbox/@(dShiftRows@(BitXor[mcc, RoundKey[15, kt] ])) ;
  state2 = dState[14, state, kt];
  state = dState[13, state2, kt];
  state2 = dState[12, state, kt];
  state = dState[11, state2, kt];
  state2 = dState[10, state, kt];
  state = dState[9, state2, kt];
  state2 = dState[8, state, kt];
  state = dState[7, state2, kt];
  state2 = dState[6, state, kt];
  state = dState[5, state2, kt];
  state2 = dState[4, state, kt];
  state = dState[3, state2, kt];
  state2 = dState[2, state, kt];
  state = BitXor[state2, RoundKey[1, kt]];
  FromCharacterCode[state]
 ]

That was it! Let’s try if it works:

daes["Í_Ëx#æ?á5½î\.b2,„Ë", "66OlSO8L7KoW44awcg2xHJ9X1FbOoF4z"]
gooby has cansur

Now I will think of a funny application for this code! Cheers. (Update: the code is put into action in this post.)

Final remarks

  • There are other variants of AES/Rijndael varying in key- and blocksizes. The en- and decryption algorithms are slightly modified. I leave it to you as a homework to look at the details.
  • The question “How do I encrypt arbitrary data, not just 16-byte blocks?” will be addressed elsewhere. (Update: in this post we start to scratch on the surface of the topic).
  • If you have questions or suggestions about the Mathematica code here, feel free to post them in the comments section.
Advertisements

About goobypl5

pizza baker, autodidact, particle physicist
This entry was posted in Encryption, Mathematica and tagged , , , . Bookmark the permalink.

2 Responses to AES in Mathematica

  1. Pingback: Block mode of operation: why ECB is bad | randomgooby

  2. Pingback: How to get a job with 250 lines of code | randomgooby

Share your thoughts

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s