How do Käsper and Schwabe's Bitsliced AES Mixcolumns work?

The only way I see it possible to do the matrix-multiplication in the MixColumns operation of AES is by shifting the bits in the multiplied number, and then reduce with the polynomial if needed.

This can be done in constant time for a bitsliced implementation as shown below in C-like pseudo-code (the example is for multiplication by $3$ in $GF(2^8)$):

vector bit[8] = "bits from a bitsliced byte";
vector new_bit[8];

//First multiply with 2 in GF(2^8) with polynomial x^8+x^4+x^3+x+1:
new_bit[0] = {0} ^ bit[7];
new_bit[1] = bit[0] ^ bit[7];
new_bit[2] = bit[1];
new_bit[3] = bit[2] ^ bit[7];
new_bit[4] = bit[3] ^ bit[7];
new_bit[5] = bit[4];
new_bit[6] = bit[5];
new_bit[7] = bit[6];

//XOR M2 with original bits to get the result for multiplying with 3
for (int bit = 0; bit < 8; bit++){
new_bit[i] = new_bit[i] ^ bit[i];
}


However, this does not seem to be the way that Käsper and Schwabe do it in their bit-sliced AES MixColumns operation. Instead they are doing rotates. I have read their paper and looked at the layout of their bit-sliced state, but I'm unable to understand, how they arrive at those equations. If someone understands it, could they perhaps help a poor soul out there to understand it too? Just sending me in the direction for a paper that explains it more in depth would mean the world to me as well!

I'm currently trying to implement a bit-sliced mix columns implementation for a cipher other than AES (I'm trying to bit-slice the primates cipher: http://primates.ae/wp-content/uploads/primatesv1.02.pdf), and have currently made the "simple" MixColumns that I described. I'm curious in whether it can be done better, after having seen that they do it differently in their AES implementation. There's not the same symmetry in this MixColumns matrice, so I doubt a bit that it is possible.

EDIT: Here is a link to their AES implementation paper: http://www.chesworkshop.org/ches2009/presentations/01_Session_1/CHES2009_ekasper.pdf

Slide #8 in the presentation you linked to describes the way Käsper and Schwabe pack the bits of the AES data blocks into CPU registers.

According to the slide, what they're doing is processing eight 128-bit AES blocks in parallel, using eight 128-bit XMM registers to store them. They're not doing basic "naïve bitslicing", which would involve using 128 $n$-bit registers to process $n$ different 128-bit AES blocks in parallel, but rather exploiting the structure of AES to effectively do "byte-level bitslicing".

Conceptually, each AES block consists of a 4×4 matrix of 8-bit bytes, for a total of 4 × 4 × 8 = 128 bits. The standard packing order of these bits in AES input / output looks like this:

| Block 0:
|-------------------------------+-------------------------------+
|           Column 0            |           Column 1            |
|-------+-------+-------+-------|-------+-------+-------+-------|
| Row 0 | Row 1 | Row 2 | Row 3 | Row 0 | Row 1 | Row 2 | Row 3 | ...
|-------|-------|-------|-------|-------|-------|-------|-------|
| bits  | bits  | bits  | bits  | bits  | bits  | bits  | bits  |
| 0 - 7 | 0 - 7 | 0 - 7 | 0 - 7 | 0 - 7 | 0 - 7 | 0 - 7 | 0 - 7 |


That is, the bits of each input / output byte are stored consecutively, and the bytes within each block are sorted first by column and then by row.

Note that, even though standard non-bitsliced AES only processes one block of data at a time, I've included a block number at the top row of the diagram. This becomes relevant when comparing this standard packing order with the internal order used by Käsper and Schwabe, since their bitsliced AES implementation processes 8 blocks at the same time.

Specifically, whereas the canonical AES bit packing order illustrated above can be compactly summarized as "block > column > row > bit", the internal packing order used by Käsper and Schwabe is "bit > row > column > block":

| Register 0 (stores bit 0 of each byte):
|-------------------------------+-------------------------------+
|             Row 0             |             Row 1             |
|-------+-------+-------+-------|-------+-------+-------+-------|
| Col 0 | Col 1 | Col 2 | Col 3 | Col 0 | Col 1 | Col 2 | Col 3 | ...
|-------|-------|-------|-------|-------|-------|-------|-------|
| block | block | block | block | block | block | block | block |
| 0 - 7 | 0 - 7 | 0 - 7 | 0 - 7 | 0 - 7 | 0 - 7 | 0 - 7 | 0 - 7 |


Since they're processing 8 independent input blocks at a time, the bits from different input blocks never influence each other. In particular, this implies that, using the "reversed" representation above, the bits in each byte never influence each other, and so the entire AES encryption can be implemented using only byte shuffling and bitwise logic operations at byte (or larger) granularity.

Specifically, looking at slide #11 of the presentation, Käsper and Schwabe note that, in their AES implementation:

• ShiftRows is reduced to byte shuffling, for which the Intel SSSE3 instruction set they're using has a dedicated instruction (pshufb); and
• MixColumns can be implemented using a combination of byte shuffling and XOR.

I have not actually looked at Käsper and Schwabe's AES implementation code (and, in any case, I suck at reading Intel assembly code), but a basic naïve implementation of AES MixColumns using their bit packing would look something like this pseudo-C code:

// rotate each 32-bit block of state by 8 bits
void rotateColumns (uint128_t state[8]) {
for (int i = 0; i < 8; i++) {
uint128_t temp = state[i];
state[i] = ((temp & 0x00FFFFFF00FFFFFF00FFFFFF00FFFFFF) << 8)
| ((temp & 0xFF000000FF000000FF000000FF000000) >> 24);
}
}

// multiply each AES byte (spread across the 8 state registers)
// by 0x02 = x modulo the AES polynomial x^8 + x^4 + x^3 + x + 1
void doubleBytes (uint128_t state[8]) {
uint128_t temp = state[7];
for (int i = 0; i < 7; i++) state[i+1] = state[i];
state[0] = temp;
state[1] ^= temp;
state[3] ^= temp;
state[4] ^= temp;
}

void bitslicedMixColumns (uint128_t state[8]) {
uint128_t oldstate[8];
for (int i = 0; i < 8; i++) oldstate[i] = state[i];
rotateColumns(oldstate);
for (int i = 0; i < 8; i++) state[i] ^= oldstate[i];
doubleBytes(state);
for (int i = 0; i < 8; i++) state[i] ^= oldstate[i];
rotateColumns(oldstate);
for (int i = 0; i < 8; i++) state[i] ^= oldstate[i];
rotateColumns(oldstate);
for (int i = 0; i < 8; i++) state[i] ^= oldstate[i];
}


A less naïve implementation would of course inline the helper functions, unroll the loops, use registers instead of arrays, and replace the klugy byte rotation in rotateColumns e.g. with a pshufb instruction and the trivial shift loop in doubleBytes with register renaming.