Writeup by poustouflan for BadAES

Table of contents

Disclaimer: This is a rough markdown retranscription of my PDF writeup I originally sent, which is better formatted, and contains full computation results for better understanding. You may want to look for the PDF version of this writeup on FCSC Discord server.

1. Challenge Description

Quick overview The challenge was about an altered AES whitebox where all linear operations were replaced with lookup tables. In particular:

XORs in the key schedule are replaced with a $256\times256$ lookup table named XK;
XOR in AddRoundKey is replaced with a $256\times256$ lookup table named X0;
XORs in MixColumns are replaced with two $256\times256$ lookup tables, X1 for the inner XORs and X2 for the outer XORs;

Additionally, standard lookup tables were also altered:

The SBox used in SubBytes is replaced with an other 256-bytes SBox named S;
The SBox used in the key schedule is replaced with an other 256-bytes SBox named SK;
The lookup tables used for the Gallois multiplication by 2 and 3 are replaced with 256-bytes lookup tables M2 and M3 respectively, and an other lookup table M1 is added for the other elements MixColumns.

All of these lookup tables are directly loaded from a binary file included in the challenge:

with open("badaes.bin", "rb") as fp:
    S = fp.read(2**8)
    SK = fp.read(2**8)
    XK = []
    for _ in range(2**8):
        XK.append(fp.read(2**8))
    X0 = []
    for _ in range(2**8):
        X0.append(fp.read(2**8))
    X1 = []
    for _ in range(2**8):
        X1.append(fp.read(2**8))
    X2 = []
    for _ in range(2**8):
        X2.append(fp.read(2**8))
    M1 = fp.read(2**8)
    M2 = fp.read(2**8)
    M3 = fp.read(2**8)

Attack model In the challenge, we get access to a random plaintext-ciphertext pair, as well as half of the key used:

with open("key.bin", "rb") as fp:
    key = b"FCSC{{" + fp.read(8) + b"}}"
E = BADAES(key)
p = os.urandom(16)
c = E.encrypt(p)

The same key is then directly used to encrypt the flag with standard AES:

enc_key = SHA256.new(key).digest()
E = AES.new(enc_key, AES.MODE_ECB)
flag = open("flag.txt", "rb").read()
enc = E.encrypt(flag.ljust(80, b"\x00"))

The goal is thus to recover the missing 8 bytes of the key using our plaintext-ciphertext pair in order to decrypt the flag.

Let’s now see the cipher in full, starting with the key schedule:

Key Schedule The key schedule is entirely done in the __init__ method of the class:

RCON = [1, 2, 4, 8, 16, 32, 64, 128, 27, 54]

def __init__(self, k):
    self.rk = [k[i:i+4] for i in range(0, 16, 4)]
    for i in range(4, 44):
        temp = self.rk[i-1][:]
        if i % 4 == 0:
            temp = temp[1:] + temp[:1]
            temp = [self.SK[x] for x in temp]
            temp[0] = self.XK[temp[0]][self.RCON[i//4-1]]
        rk = [self.XK[self.rk[i-4][j]][temp[j]] for j in range(4)]
        self.rk.append(rk)
    self.rk = [self.rk[i:i+4] for i in range(0, 44, 4)]

Figure 1 depicts the key schedule for the two first keys, where $XK_x[y]=XK[y][x]$. This process is iterated 10 times, with the $x$ in $XK$ taking all values in RCON iteratively.

AddRoundKey AddRoundKey is similar to standard AES, except XOR has been replaced with a $256\times256$ lookup table, named X0.

def ARK(self, s, r):
    for c in range(4):
        for ro in range(4):
            s[ro][c] = self.X0[self.rk[r][c][ro]][s[ro][c]]
    return s

SubBytes SubBytes is again very similar to standard AES, except the SBox used is not the same one.

def SB(self, s):
    for i in range(4):
        for j in range(4):
            s[i][j] = self.S[s[i][j]]
    return s

ShiftRows This time absolutely no difference with standard AES.

def SR(self, s):
    for i in range(1, 4):
        s[i] = s[i][i:] + s[i][:i]
    return s

MixColumns MixColumns is probably the step that has been the most affected. The first difference compared to standard Rijndael MixColumns is that the multiplication of elements in $GF(2^8)$ is replaced with 256-bytes lookup tables: M1, M2 and M3. The second difference is that elements in the columns are then combined using $256\times256$ lookup tables instead of simple XORs: first $(M_i[a_0], M_j[a_1])$ and $(M_i[a_2], M_j[a_3])$ are combined using X1, then the two resulting bytes are combined using X2.

Figure 2 depicts this algorithm for a single byte of output.

def MC(self, s):
    for c in range(4):
        a0, a1, a2, a3 = [s[r][c] for r in range(4)]
        col = [
            self.X2[self.X1[self.M2[a0]][self.M3[a1]]][self.X1[self.M1[a2]][self.M1[a3]]],
            self.X2[self.X1[self.M1[a0]][self.M2[a1]]][self.X1[self.M3[a2]][self.M1[a3]]],
            self.X2[self.X1[self.M1[a0]][self.M1[a1]]][self.X1[self.M2[a2]][self.M3[a3]]],
            self.X2[self.X1[self.M3[a0]][self.M1[a1]]][self.X1[self.M1[a2]][self.M2[a3]]]
        ]
        for r in range(4):
            s[r][c] = col[r]
    return s

Encrypt Finally, the encryption function combines these exactly as you would expect from AES:

def encrypt(self, p):
    s = [
        [p[r+4*c] for c in range(4)]
        for r in range(4)
    ]
    s = self.ARK(s, 0)
    for r in range(1, 10):
        s = self.SB(s)
        s = self.SR(s)
        s = self.MC(s)
        s = self.ARK(s, r)
    s = self.SB(s)
    s = self.SR(s)
    s = self.ARK(s, 10)
    return bytes([s[r][c] for c in range(4) for r in range(4)])

Since we have 64 out of the 128 bits of the key, an exhaustive key search for the key would be a $2^{64}$ attack, which might actually be affordable with a big enough cluster of GPUs. However I am poor so that’s why we are going to see how to reduce greatly the time needed to solve this challenge.

2. Introduction on Statistical Analysis

When we are faced against a random-looking lookup table, it is not easy to determine whether the lookup table is actually random or has any kind of interesting pattern, structures or properties that could make the challenge easier. In order to detect the most common structures, we often use a set of metrics to see if anything looks out of the ordinary, without really knowing in advance which one will give any result.

Three years ago, I was faced with a similar challenge, and I ended up making a tool that performs usual cryptographical analysis and summarize any statistical anomalies it could find: https://github.com/PoustouFlan/SUnbox. This tool is far from being complete and still a work in progress that I left over years ago, but it was useful multiple times during this challenge.

This section describes a set of tools in linear and differential cryptanalysis that I will use during my solve. Feel free to skip to the next section if you are already familiar with differential and linear cryptanalysis. For the purpose of this writeup, I will explain in details the kind of analysis that is done, so that you can follow from scratch without any other tool than what is already bundled in SageMath. However, most visuals here are generated using SUnbox.

Differential Analysis One common thing to check in a SBox is its differential properties. This tells us how likely it is that a specific difference in input to an SBox will yield a specific difference in output. For an $m\times n$ SBox (meaning, the SBox takes m bits in input, and outputs n bits), we encode this number in a $2^m\times 2^n$ table named DDT (Difference Distribution Table):

DDT_S[\Delta_{in}][\Delta_{out}] = 2^m P(S[x\oplus\Delta_{in}] = S[x]\oplus\Delta_{out})

(the $2^m$ factor is simply to keep everything integers).

This table is quite easy to compute in $O(2^{m+n} + 2^{2m})$. For an $8\times8$ SBox:

def difference_distribution_table(S):
    """
    Returns the Difference Distribution Table (DDT) for a m x n SBox
    (meaning, m bits in input, n bits in output).
    DDT[din][dout] corresponds to the probability P[S(x ^ din) == S(x) ^ dout].
    Values are multiplied by 2^m to remain integers.
    """
    m = (len(S) - 1).bit_length()
    n = max(S).bit_length()
    nrows = 1 << m
    ncols = 1 << n
    DDT = [[0] * ncols for _ in range(nrows)]
    for x in range(nrows):
        sx = S[x]
        for din in range(nrows):
            dout = sx ^ S[x ^ din]
            DDT[din][dout] += 1
    return DDT

The same can be done for our $256\times256$ lookup tables, which can be seen as a $16\times8$ SBox (taking $2\times8$ bits in input, and 8 bits in output). This gives a $2^{16}\times2^8$ table, which can be computed in about $2^{32}$ operations (may take some minutes in Python).

Of course all of this is built-in in SageMath and you could just call Sbox.difference_distribution_table from the module sage.crypto.sbox, but I like to explain everything in details.

A simple metric you usually want to check for is the differential uniformity $\Delta_S$ of an SBox, which is simply the maximal entry in its DDT (for $\Delta_{in}\ne0$). I will not spill all the math here, but one could reasonably expect a differential uniformity of 10 or 12 for an $8\times8$ random SBox, and a differential uniformity of about 376 for a $16\times8$ SBox. (This can either be computed empirically, or by modeling the entries of the DDT as a binomial distribution $2 \mathcal{B} (2^m-1, 2^{-n})$ ) ¹.

[^1]: O’Connor. On the distribution of characteristics in bijective mappings.

This value can be directly computed in SageMath using the method differential_uniformity, or here you can simply iterate through the DDT and check for the maximal entry (skipping the first element). Checking the differential uniformity of all lookup tables involved in this challenge gives the following results:

$\Delta_S = 12$, $\Delta_{SK} = 12$, $\Delta_{M_1} = 10$, $\Delta_{M_2} = 12$, $\Delta_{M_3} = 12$
$\Delta_{XK} = 904$, $\Delta_{X_0} = 420$, $\Delta_{X_1} = 936$, $\Delta_{X_2} = 904$

For the $8\times8$ SBoxes, this is very reasonable. Regarding the $16\times8$ lookup tables, a differential uniformity above 900 is clearly not what you would get from a random table. Since we only have one plaintext-ciphertext pair, we will not be able to perform any kind of differential statistical attacks, but at least it tells us that there is something to look for in the $16\times8$ lookup tables.

Linear Analysis An other kind of statistical analysis that is often done is the linear analysis. Here, we check how likely is it that the linear combination of some output bits can be correlated to the linear combination of some input bits. For a $m\times n$ SBox, this can be summarized in a $2^m\times 2^n$ table that we usually call LAT (Linear Approximation Table):

LAT_S[\nabla_{in}][\nabla_{out}] = 2^m\left(P[\nabla_{in}\cdot x = \nabla_{out}\cdot S(x)] - \frac{1}{2}\right)

(the subtraction and the $2^m$ factor are again simply to keep everything as integers).

This is a little more tedious to compute faster than the naive way, but there is a trick involving a Fast Walsh-Hadamard Transform that allows to compute a column of the LAT (so the list of probabilities for a specific output mask) quickly. We compute for all inputs in the SBox its output and the value of the linear combination with our output mask. If we encode these results as 1 or -1, the Walsh-Hadamard Transform of this array yield our desired probability for each input mask. Using a Fast Walsh-Hadamard Transform, we can get each column in $O(m2^m)$, which allows to compute the entire table in $O(m2^{n+m})$:

def walsh_spectrum(array):
    """
    In-place Fast Walsh-Hadamard Transform of array
    """
    h = 1
    n = len(array)
    while h < n:
        for i in range(0, n, h*2):
            for j in range(i, i+h):
                x = array[j]
                y = array[j+h]
                array[j] = x + y
                array[j+h] = x - y
        h *= 2

def linear_approximation_table(S):
    """
    Returns the Linear Approximation Table (LAT) for a m x n SBox.
    LAT[mask_in][mask_out] corresponds to the probability P[mask_in . x == mask_out . S(x)], where . denotes a vector dot product.
    Absolute bias scale is used, therefore the actual value corresponds to the probability 1/2, multiplied by 2^m.
    """
    m = (len(S) - 1).bit_length()
    n = max(S).bit_length()
    nrows = 1 << m
    ncols = 1 << n
    LAT = [[None] * ncols for _ in range(nrows)]
    for mask_out in range(ncols):
        row = []
        for x in range(nrows):
            row.append(1 - (((mask_out & S[x]).bit_count() & 1)*2))
        walsh_spectrum(row)
        for mask_in in range(nrows):
            LAT[mask_in][mask_out] = row[mask_in] // 2
    return LAT

Again, most people don’t need that and can resort to use the method linear_approximation_table, built-in in SageMath. But it’s a cool trick, I like it, and now you know why it’s somewhat fast to compute. This table, however, represents absolute scaled bias, and contains both positive and negative entries. The commonly used metric here is the linearity ${L}(S)$ of an SBox S, which corresponds to twice the maximal entry in the LAT (in absolute value). You can get it yourself, or by simply calling the method SBox.linearity in SageMath.

Anyway, know that the linearity for a random $8\times8$ SBox is expected to be around 68 or 72, and around 1384 for a random $16\times8$ SBox. Let’s see what we get for our SBoxes:

$\mathcal{L}(S) = 80$, $\mathcal{L}(SK) = 76$, $\mathcal{L}(M_1) = 68$, $\mathcal{L}(M_2) = 76$, $\mathcal{L}(M_3) = 68$
$\mathcal{L}(XK) = 2704$, $\mathcal{L}(X_0) = 1392$, $\mathcal{L}(X_1) = 1840$, $\mathcal{L}(X_2) = 1760$

Again, for the $8\times8$ SBoxes, this is very reasonable. Regarding the $16\times8$ lookup tables, we get the same impression that $X_0$ is the only table which actually looks random. Since we only have one plaintext-ciphertext pair, we will not be able to perform any kind of linear statistical attacks, but at least it confirms that there is clearly something to look for in most of the $16\times8$ lookup tables.

Autocorrelation Table Finally, a third property that we may be interested in is a kind of mix of both previous ones: given a specific difference $\Delta_{in}$ in input, can we correlate some specific bits $\nabla_{out}$ of their output? For an $m\times n$ SBox, we encode this number in a $2^m\times 2^n$ table named ACT (Autocorrelation Table):

ACT_S[\Delta_{in}][\nabla_{out}] = 2^m\left(P[\nabla_{out}\cdot(S[x]\oplus S[x\oplus\Delta_{in}])=0] - \frac{1}{2}\right)

(again, the subtraction ans the $2^m$ factor are simply to keep everything as integers). This can be computed using basically the same trick as for the Linear Approximation Table, but this time taking the Walsh-Hadamard Transform of the Difference Distribution Table:

def hadamard_matrix(n):
    """
    Generates a hadamard matrix of size 2^n
    """
    if n == 0:
        return np.array([[1]])
    else:
        h_n_minus_1 = hadamard_matrix(n-1)
        h_n = np.vstack((
            np.hstack((h_n_minus_1, h_n_minus_1)),
            np.hstack((h_n_minus_1, -h_n_minus_1))
        ))
        return h_n

def autocorrelation_table(S):
    """
    Returns the Autocorrelation Table (ACT) for this SBox.
    ACT[din][mask_out] corresponds to the probability P[mask_out . (S(x) ^ S(x ^ din)) == 0], where . denotes a vector dot product.
    Absolute bias scale is used, therefore the actual value corresponds to the probability 1/2, multiplied by 2^m.
    """
    ddt = np.matrix(difference_distribution_table(S))
    n = max(S).bit_length()
    had = hadamard_matrix(n)
    ACT = ddt * had
    return ACT.tolist()

This is of course built-in in sage under the autocorrelation_table method of the class SBox in sage.crypto.sbox. We may again check for the maximal entry (in absolute value) in this table, denoted as the absolute indicator $\mathcal{AC}_S$ of the SBox. Here are the values we get:

$\mathcal{AC} _ S = 104$, $\mathcal{AC} _ {SK} = 96$ , $\mathcal{AC} _ {M_1} = 104$ , $\mathcal{AC} _ {M_2} = 96$ , $\mathcal{AC} _ {M_3} = 112$
$\mathcal{AC} _ {XK} = 3648$ , $\mathcal{AC} _ {X_0} = 2768$ , $\mathcal{AC} _ {X_1} = 2768$ , $\mathcal{AC} _ {X_2} = 2664$

Generally, we expect a random $8\times8$ SBox to have an absolute indicator around 104 or 112, and we expect a random $16\times8$ lookup table to have an absolute indicator around 2128. For our $8\times8$ SBoxes, we get again very reasonable values. However, yet again the $16\times8$ lookup tables appears to have some intricate hidden inner structures, and this time even $X_0$ fails the test.

3. Resolution

So, given a quick look at the $16\times8$ lookup tables, we know there is probably something hidden in there. From a lot of papers on AES whitebox cryptography, it is a common practice to encode the linear maps f in the cipher (such as the XORs) by using two random bijections A and B and replacing $f(x)$ with $A\circ f \circ B(x)$ ². This is even more likely considering the fact that those tables are replacing exactly the linear operations in our alternative AES.

[^2]: James A. Muir. A tutorial on white-box AES.

Testing XOR-isotopy If there exists three bijections A, B, C such that a function $f(x,y)=A(B(x)\oplus C(y))$, we will call $f$ a XOR-isotope. There is a way to check if our lookup tables are constructed this way. The method I am going to present is heavily inspired by the technique formalized by Billet, Gilbert, and Ech-Chatbi in their seminal 2004 attack on a white-box AES ³.

[^3]: Olivier Billet, Henri Gilbert, and Charaf Ech-Chatbi. Cryptanalysis of a white box AES implementation.

Suppose that our table $f$ is defined as $f(x,y)=A(B(x) \ast C(y))$ for some bijection A, B, C and some inner function $\ast$ . Then we can define $f_y(x)=A(B(x) \ast C_y)$ and compute its inverse $f_y^{-1}(x)=B^{-1}(A^{-1}(x)\ast C_y^{-1})$ . If we go forward with $f_{y_1}$ for some $y_1$ and backward with $f_{y_2}$ for some other $y_2$, we get this:

\begin{aligned}
f_{y_1}(f_{y_2}^{-1}(x)) &= f_{y_1}(B^{-1}(A^{-1}(x)*C_{y_2}^{-1})) \\
                         &= A(B(B^{-1}(A^{-1}(x)*C_{y_2}^{-1}))*C_{y_1}) \\
                         &= A(A^{-1}(x)*C_{y_2}^{-1}*C_{y_1})
\end{aligned}

Let us note $P_{y_1,y_2} = f_{y_1}\circ f_{y_2}^{-1}$ the function we just described, and let $\Delta_{y_1,y_2} = C_{y_2}^{-1} \ast C_{y_1}$ . Therefore we have $P_{y_1,y_2}(x) = A(A^{-1}(x) \ast \Delta_{y_1,y_2})$ . We can now use this function to check if $\ast$ exists and is an involution by looking at $P_{y_1,y_2}^2(x)$:

\begin{aligned}
P_{y_1,y_2}(P_{y_1,y_2}(x)) &= P_{y_1,y_2}(A(A^{-1}(x)\ast\Delta_{y_1,y_2})) \\
                            &= A(A^{-1}(A(A^{-1}(x)\ast\Delta_{y_1,y_2}))\ast\Delta_{y_1,y_2}) \\
                            &= A(A^{-1}(x)\ast\Delta_{y_1,y_2}\ast\Delta_{y_1,y_2})
\end{aligned}

If $\ast$ is involutive (so, in $\mathbb{F}_2^8$ this would mean it is a XOR), we will always get:

\begin{aligned}
P_{y_1,y_2}(P_{y_1,y_2}(x)) &= A(A^{-1}(x)\oplus\Delta_{y_1,y_2}\oplus\Delta_{y_1,y_2}) \\
&= A(A^{-1}(x)) = x
\end{aligned}

So, we have a simple test to check this. Let us build $P_{y_1,y_2}$ for all $(y_1,y_2)$ for our tables, and if $P_{y_1,y_2}^2$ is always the identity, it means our $16\times8$ lookup tables are simply XORs in disguise.

def compose(p1, p2):
    return [p1[p2[x]] for x in range(256)]

def invert(f):
    inverse_map = [-1] * 256
    for x, y in enumerate(f):
        inverse_map[y] = x
    return inverse_map

identity = list(range(256))

def test_xor_isotope(S):
    for y2 in range(256):
        fy2_inv = invert(S[y2])
        for y1 in range(256):
            P = compose(S[y1], fy2_inv)
            P2 = compose(P, P)
            if P2 != identity:
                return False
    return True

print(test_xor_isotope(XK)) # True
print(test_xor_isotope(X0)) # True
print(test_xor_isotope(X1)) # True
print(test_xor_isotope(X2)) # True

Decoding the isotopy More than testing for isotopy, we can also recover from our functions the bijections used for the isotopy, up to an affine transformation, again using the methodology formalized in the BGE attack ³. Recall our composition from earlier, but let us fix a baseline by arbitrarily choosing $y_2=0$.

We define the set of permutations $Q_y(x)$ as:

Q_y(x) = f_y \circ f_0^{-1}(x) = A(A^{-1}(x) \oplus \Delta_y)

where $\Delta_y = C(y) \oplus C(0)$.

Because C is a bijection, the set of values $\Delta_y$ covers all 256 elements of the vector space $\mathbb{F}_2^8$. Therefore, the set of permutations $\mathcal{Q} = \{Q_y | y \in \mathbb{F}_2^8\}$ forms an abelian group under function composition, and this group is strictly isomorphic to $(\mathbb{F}_2^8, \oplus)$.

To recover the outer encoding $A$ , we need to reconstruct this isomorphism explicitly. Since $\mathbb{F} _ 2^8$ is an 8-dimensional vector space over $\mathbb{F} _ 2$, we can find a basis of 8 independent permutations $\{Q_{b_0}, …, Q_{b_7}\}$ from our set $\mathcal{Q}$. Any permutation $Q_y$ in the group can be generated by a unique composition of these basis elements.

Let us define a mapping $\phi: \mathbb{F} _ 2^8 \rightarrow \mathcal{Q}$ by taking the binary representation of an integer $k = \sum _ {i=0}^7 k _ i 2^i$ and composing the corresponding basis elements:

\phi(k) = \bigcirc_{i=0}^7 (Q_{b_i})^{k_i}

Because of the structure of $\mathcal{Q}$, we know that for any integer $k$ , the resulting permutation takes the form $\phi(k)(x) = A(A^{-1}(x) \oplus L(k))$, for some linear map $L$ . If we evaluate this mapping at $x=0$, we get $\phi(k)(0) = A(A^{-1}(0) \oplus L(k))$. Since internal encodings are unique up to affine equivalence, let’s simply set $A’(k) = \phi(k)(0)$.

Once the outer encoding $A’$ is recovered, isolating $B$ and $C$ becomes trivial. Since we want $f_0(x) = A’(B’(x) \oplus C’(0))$, we can simply set $C’(0) = 0$ giving us:

B'(x) = A'^{-1}(f_0(x))

For $C’$, we simply find the integer index $k$ such that $\phi(k)$ matches the permutation $Q_y$ we generated earlier:

C'(y) = \phi^{-1}(Q_y)

This process allows us to recover three bijections $A’$ , $B’$ , $C’$ such that $f(x, y) = A’(B’(x) \oplus C’(y))$ , which will greaty simplify our analysis of the cipher.

def recover_isotopy(S):
    f0_inv = invert(S[0])
    Q = [compose(S[y], f0_inv) for y in range(256)]
    phi = [identity]
    for Qb in Q:
        if Qb in phi:
            continue
        phi.extend([compose(Qa, Qb) for Qa in phi])
        if len(phi) == 256:
            break
    A = [phi[k][0] for k in range(256)]
    B = compose(invert(A), S[0])
    C = [phi.index(Q[y]) for y in range(256)]
    # Sanity Check
    for y in range(256):
        for x in range(256):
            assert S[y][x] == A[B[x] ^ C[y]]
    return A, B, C

Symmetric Isotopies There is one further simplification we may do for some of our lookup tables. Indeed, it is quite easy to check which of them are actually symmetric:

def is_symmetric(S):
    for y in range(256):
        for x in range(256):
            if S[y][x] != S[x][y]:
                return False
    return True

print(is_symmetric(XK)) # True
print(is_symmetric(X0)) # False
print(is_symmetric(X1)) # True
print(is_symmetric(X2)) # True

This means most of our $256\times256$ lookup tables are symmetric. Since they are symmetric and XOR-isotopes, it means we can actually encode them using only two bijections. Indeed, we can observe that $B(x) \oplus C(x) = B(y) \oplus C(y)$. This literally means $B \oplus C$ is a constant, that we may call $\Delta$. Let’s define $A’(x) = A(x \oplus \Delta)$. We can now define f using only two bijections:

f(x, y) = A'(B(x) \oplus B(y))

Let’s edit our code to recover these clean isotopies:

def recover_isotopy(S):
    #... SNIP
    if is_symmetric(S):
        D = B[0] ^ C[0]
        A = [A[x ^ D] for x in range(256)]
        C = B
    # Sanity Check
    for y in range(256):
        for x in range(256):
            assert S[y][x] == A[B[x] ^ C[y]]
    return A, B, C

XKA, XKB, _ = recover_isotopy(XK)
X0A, X0B, X0C = recover_isotopy(X0)
X1A, X1B, _ = recover_isotopy(X1)
X2A, X2B, _ = recover_isotopy(X2)

Decoding the Key Schedule Now that we have an alternative and simpler representation of XK, let’s start with simplifying the key schedule. If we replace $XK[x, y]$ with our new representation $XK_A[XK_B[x] \oplus XK_B[y]]$, we get the diagram showing apparent XORs.

One thing we notice is that $XK_A$ is almost always followed by $XK_B$ of some kind. Therefore, we may attempt to get rid of most of the $XK_B$ in the key schedule by doing the following substitutions:

Replace all $XK_A$ with $XK_A’ = XK_B \circ XK_A$
Replace all SK with $SK’ = XK_B \circ SK \circ XK_B^{-1}$
Add an initial layer of $XK_B$ on the initial key;
Also apply $XK_B$ directly on the round constants RCON;
Add a final layer of $XK_B^{-1}$ before the output round keys, to compensate for the extra $XK_B$ added.

As a result, we get an equivalent but much simpler key schedule. Let us compute the new SBoxes we get from this simplification:

XKAp = compose(XKB, XKA)
XKBinv = invert(XKB)
SKp = compose(compose(XKB, SK), XKBinv)

Oh, did you notice? $XK_A’$ is the identity permutation! We can simply get rid of it in our key schedule, and what we get is now exactly the AES key schedule, up to a bijection of the input and output bytes, and a different SBox used.

Decoding AddRoundKey Now, what about all the other steps of the cipher that involved some $256\times256$ lookup tables? Let’s start with AddRoundKey. If we look closely at how AddRoundKey interacts with the other steps, we can perform the following substitutions to remove $X0_R$ completely and turn AddRoundKey back into a simple XOR:

Apply a substitution with $X0_B$ directly on the plaintext;
Replace S with $S’ = S \circ X0_A$;
Replace MC with $MC’ = X0_B \circ MC$
Apply $X0_C$ on each round key;
Apply a final substitution with $X0_A$ to get the ciphertext.

What’s more, if we recall from the key expansion, all of the round keys are substituted by $XK_B^{-1}$ once before entering $X0_C$. We can also directly combine them into a single SBox $XR = X0_C \circ XK_B^{-1}$.

Sp = compose(S, X0A)
XR = compose(X0C, XKBinv)

Notice how $XR$ is now simply a XOR with 0xEB? We can get rid of it entirely by incorporating this XOR directly into the preliminary $X0_B$ for the first round, and into the output of $MC’$ for all other rounds!

Decoding MixColumns Last but not least, let’s see how much we can simplify MixColumns this time. Let’s naturally expand $X_1[y][x]$ and $X_2[y][x]$ into our new representation $X1_A[X1_B[y] \oplus X1_B[x]]$ and $X2_A[X2_B[y] \oplus X2_B[x]]$ respectively. If we consider this, the additionnal $X0_B$, and the additional 0xEB we got from the decoding of AddRoundKey, we can perform the following substitution to have an equivalent but simpler MC’':

Replace all $M_i$ with $M_i’ = X1_B \circ M_i$;
Replace all $X_1$ with an XOR, and add a new SBox $X_1’ = X2_B \circ X1_A$ after it;
Also replace $X_2$ with an XOR, and add a new SBox $X_2’ = (X0_B \circ X2_A) \oplus 0xEB$ after it.

def xor(S, b):
    return [y ^ b for y in S]

M1p = compose(X1B, M1)
M2p = compose(X1B, M2)
M3p = compose(X1B, M3)
X1p = compose(X2B, X1A)
X2p = xor(compose(X0B, X2A), 0xEB)

Getting rid of $X_i’$ This is already getting nicer, but it’s only the beginning. If we wanted it to look exactly like the original MixColumns, we would need to get rid of the SBoxes $X_1’$ and $X_2’$ as well. This is something we could do if those were affine SBoxes, meaning the output bits of the SBoxes were some specific linear combinations of the input bits (which would therefore distribute well over the XORs before them). Fortunately, we now have a simple way to check this:

print(differential_uniformity(X1p)) # 256
print(linearity(X1p)) # 256
print(differential_uniformity(X2p)) # 256
print(linearity(X2p)) # 256

Good news! They are all perfectly affine. We thus have $X_1’[x] = A_1[x] \oplus B_1$ and $X_2’[x] = A_2[x] \oplus B_2$, where $A_1$ and $A_2$ are linear transformation. We are therefore allowed to distribute the linear transformations over the XORs back to the matrices $M_i’$, for instance:

\begin{aligned}
X_2'[X_1[x] \oplus X_1'[y]] &= A_2[X'_1[x] \oplus X'_1[y]] \oplus B_2 \\
&= A_2 \circ X_1'[x] \oplus A_2 \circ X_1'[y] \oplus B_2
\end{aligned}

And we can do the same for $X’_1$ :

\begin{aligned}
X_1[M_i'[x] \oplus M'_j[y]] &= A_1[M'_i[x] \oplus M'_j[y]] \oplus B_1 \\
&= A_1 \circ M_i'[x] \oplus A_1 \circ M_j'[y] \oplus B_1
\end{aligned}

Distributing everything, we get:

$$MC’(a_0,a_1,a_2,a_3) = A_2A_1M’{i_0}[a_0] \oplus A_2A_1M’{i_1}[a_1] \oplus A_2A_1M’{i_2}[a_2] \oplus A_2A_1M’{i_3}[a_3] \oplus B_2$$

We can extract $A_i$ and $B_i$ from $X’ _ i[x] = A _ i[x] \oplus B_i$ where $A _ i$ is linear very simply since we then have $B _ i = X’ _ i[0]$ , and $A _ i = X’ _ i \oplus B _ i$ . We can then perform the substitutions $M’’ _ i = A _ 2 A _ 1 M’ _ i$ and get back to a MixColumns that is starting to look like the standard Rijndael MixColumns.

def decompose_affine(S):
    B = S[0]
    A = xor(S, B)
    return A, B

A1, B1 = decompose_affine(X1p)
A2, B2 = decompose_affine(X2p)

M1s = compose(A2, compose(A1, M1p))
M2s = compose(A2, compose(A1, M2p))
M3s = compose(A2, compose(A1, M3p))

Getting rid of $M_1$ We can go one last step further in the reduction of MixColumns by getting rid of both $M_1’’$ and $B_2$ entirely. If we look back at the entire cipher, we first notice we can easily get rid of $B_2$ by integrating it into the SBox $S’$ that comes right after in the cipher.

For $M_1’’$ what we can do is simply compose all of the matrices with $M_1’’^{-1}$ in input of MixColumns, and compose $M_1’’$ back sooner in the cipher. The neat thing is that this goes through ShiftRows without any issue, so we can compose $S’$ with $M_1’’$ and it will do the job.

Using the following substitutions in the whole cipher, we get again something simplified yet equivalent:

Replace $M’’ _ 2$ in $MC’$ with $M^{(3)} _ 2 = M’’ _ 2 \circ M’’^{-1} _ 1 \quad$
Replace $M’’ _ 3$ in $MC’$ with $M^{(3)} _ 3 = M’’ _ 3 \circ M’’^{-1} _ 1 \quad$
Get rid of $M’’ _ 1$ entirely in $MC'$
Get rid of the $\oplus B _ 2$ at the end of $MC'$
Compensate for the added $M’’^{-1} _ 1$ in the input of $MC’$ and the removed $\oplus B _ 2$ at its end by replacing the SBox $S’$ with $S’’[x] = M’’ _ 1 \circ S[x \oplus B _ 2]$
Compensate for the missing $\oplus B _ 2$ at the end of the last $MC’’$ by replacing the final $X0 _ A$ with $X0’ _ A[x] = X0 _ A[x \oplus B _ 2]$
Compensate for the added $\oplus B _ 2$ on the input from the first SBox $S’’$ by replacing $X0 _ B$ with $X0’ _ B = X0 _ B \oplus B _ 2$

M1_inv = invert(M1s)
M2t = compose(M2s, M1_inv)
M3t = compose(M3s, M1_inv)
Ss = compose(M1s, [Sp[x ^ B2] for x in range(256)])
X0Ap = [X0A[x ^ B2] for x in range(256)]
X0Bp = xor(X0B, B2 ^ 0xEB)

One can notice that $M^{(3)} _ 2$ and $M^{(3)} _ 3$ are now affine.

The result The whole cipher is therefore strictly equivalent to the following one, which resembles very closely to the standard AES, except for a substitution of the inputs and outputs, and different SBoxes used for SubBytes and the key schedule:

class BADAES:
    def __init__(self, k):
        # Added permutation of the key
        RCON = [XKB[i] for i in [1, 2, 4, 8, 16, 32, 64, 128, 27, 54]]
        k = [XKB[x] for x in k]

        # AES Key Schedule
        self.rk = [[k[j] for j in range(i, i+4)] for i in range(0, 16, 4)]
        for i in range(4, 44):
            temp = self.rk[i-1][:]
            if i % 4 == 0:
                temp = temp[1:] + temp[:1]
                temp = [SK[x] for x in temp]
                temp[0] = temp[0] ^ RCON[i//4-1]
            rk = [self.rk[i-4][j] ^ temp[j] for j in range(4)]
            self.rk.append(rk)
        self.rk = [self.rk[i:i+4] for i in range(0, 44, 4)]

    def ARK(self, s, r):
        for c in range(4):
            for ro in range(4):
                s[ro][c] ^= self.rk[r][c][ro]
        return s

    def SB(self, s):
        for i in range(4):
            for j in range(4):
                s[i][j] = S[s[i][j]]
        return s

    def SR(self, s):
        for i in range(1, 4):
            s[i] = s[i][i:] + s[i][:i]
        return s

    def MC(self, s):
        for c in range(4):
            a0, a1, a2, a3 = [s[r][c] for r in range(4)]
            col = [
                M2[a0] ^ M3[a1] ^ a2 ^ a3,
                a0 ^ M2[a1] ^ M3[a2] ^ a3,
                a0 ^ a1 ^ M2[a2] ^ M3[a3],
                M3[a0] ^ a1 ^ a2 ^ M2[a3]
            ]
            for r in range(4):
                s[r][c] = col[r]
        return s

    def encrypt(self, p):
        # Added permutation of the plaintext
        p = [X0B[x] for x in p]
        
        # Standard AES (Except for the SBoxes)
        s = [[p[r+4*c] for c in range(4)] for r in range(4)]
        s = self.ARK(s, 0)
        for r in range(1, 10):
            s = self.SB(s)
            s = self.SR(s)
            s = self.MC(s)
            s = self.ARK(s, r)
        s = self.SB(s)
        s = self.SR(s)
        s = self.ARK(s, 10)
        
        # Added permutation of the ciphertext
        s = [X0A[x] for x in s]
        return bytes([s[r][c] for c in range(4) for r in range(4)])

4. The weakness

The whole cipher is now equivalent to a standard AES. The bijections $X0_A$ on the ciphertext, $X0_B$ on the plaintext and $XK_B$ on the key and round constants can be inverted easily for a given plaintext-ciphertext pair, so they do not impact the security of the cipher at all. The whole cipher relies solely on the non-linearity of the two new SBoxes $S$ and $SK$.

Doing a quick statistical test using the tools we saw in the first section on S and SK reveals something very promising:

$\Delta _ S = 32$, $\mathcal{L} _ S = 128$, $\mathcal{CA} _ S = 256$
$\Delta _ {SK} = 32$, $\mathcal{L} _ {SK} = 128$, $\mathcal{CA} _ {SK} = 256$

Their differential uniformity is already quite high, their linearity is even higher, but most importantly, their absolute indicator is maximal! Indeed, if we take a closer look at their Linear Approximation Tables and their Autocorrelation Tables, we find some very interesting patterns:

Those green horizontal and vertical lines with bright red spot at their crossing reveals what we denote as linear structures. For an $8\times8$ SBox $S$, an entry $ACT _ S[\Delta _ {in}][\nabla _ {out}] = 256$ (the maximal value) means we always have $\nabla _ {out}\cdot(S[x]\oplus S[x\oplus\Delta _ {in}])=0$ for all x.

In our case, this happens exactly for all $(\Delta _ {in}, \nabla _ {out})$ in $U \times V$, where:

$U = \{0x00, 0x03, 0x08, 0x0b, 0x24, 0x27, 0x2c, 0x2f, 0x81, 0x82, 0x89, 0x8a, 0xa5, 0xa6, 0xad, 0xae\}$
$V = \{0x00, 0x10, 0x24, 0x34, 0x40, 0x50, 0x64, 0x74, 0x83, 0x93, 0xa7, 0xb7, 0xc3, 0xd3, 0xe7, 0xf7\}$

These both form vector subspaces of dimension 4 of $\mathbb{F} _ 2^8$:

$U = Span((0x03, 0x08, 0x24, 0x81))$
$V = Span((0x10, 0x24, 0x40, 0x83))$

What does it mean for an SBox to have such linear structures? By definition, since the maximal entries occur for $(\Delta _ {in}, \nabla _ {out}) \in U \times V$, if we take any $u \in U$ and any $v \in V$, we will have $v \cdot (S[x] \oplus S[x \oplus u]) = 0$ for all $x$. This means $v$ is always orthogonal with $(S[x] \oplus S[x \oplus u])$.

Since this holds for all $v \in V$, the difference $S[x] \oplus S[x \oplus u]$ must be a value that is orthogonal to the entirety of V. We call the set of elements orthogonal to all elements of V the orthogonal complement of V, denoted as $V^{\perp}$.

Vc = []
for x in range(256):
    for v in V:
        if (v & x).bit_count() % 2 == 1:
            break
    else:
        Vc.append(x)
print(Vc) # {0x00, 0x03, 0x08, 0x0b, ..., 0xa5, 0xad, 0xae}

Computing and printing $V^{\perp}$ reveals an interesting surprise, we have $V^{\perp} = U$! So, for any $x \in \mathbb{F}_2^8$ and any $u \in U$, we know that the difference $S[x] \oplus S[x \oplus u]$ will always be orthogonal to all elements in V. By definition of $V^{\perp}$, this translates to a very powerful property:

S[x] \oplus S[x \oplus u] \in U

This is what we call a linear structure. If you change the input of the SBox by an element of U, the output will also change by an element of U. The same property also applies perfectly to the key schedule SBox SK, which has the same linear structures.

It is this specific invariant subspace U that perfectly aligns with the linear layers of the cipher and will allow us to break it. To understand how this vulnerability completely shatters the security of the cipher, we need to look at how this property interacts with the concept of a quotient space.

Let us consider the vector space $\mathbb{F}_2^8$. Because U is a subspace of $\mathbb{F}_2^8$, we can construct the quotient space $\mathbb{F}_2^8/U$. In simple terms, the quotient space groups the 256 elements of $\mathbb{F}_2^8$ into sets called cosets. Two elements x and y belong to the same coset if and only if their difference belongs to U (i.e., $x \oplus y \in U$). When this is true, we say that x and y are equivalent modulo U, denoted as:

x \equiv y \pmod U

Since U contains 16 elements, the quotient space $\mathbb{F}_2^8/U$ contains exactly $256/16=16$ distinct cosets. To compute with these cosets practically, we choose a single representative element for each coset. The set of these 16 representatives is called a transversal, which we will denote as T. We can then define an explicit projection map $\pi_U: \mathbb{F}_2^8 \rightarrow T$, which takes any byte and returns the unique representative of its coset.

Here is how we can build this in Python:

# Projection map to T
PU = [-1] * 256

for t in range(256):
    if PU[t] != -1:
        continue
    for u in U:
        PU[u ^ t] = t

Now, let us revisit the SBox property $S[x] \oplus S[x \oplus u] \in U$. In terms of quotient spaces, this states that if you take an element x, and you add any $u \in U$ to it, the output $S[x \oplus u]$ will only differ from $S[x]$ by some element in U. Mathematically, if $x \equiv y \pmod U$, then $S[x] \equiv S[y] \pmod U$.

This is devastating for the SBox. It means that the projection of the output of the SBox only depends on the projection of the input. The SBox induces a perfectly well-defined function on the 16 elements of the quotient space!

This property holds for both S and SK. But what about the rest of the cipher?

AddRoundKey is just a XOR, it preserves the quotient space.
ShiftRows only moves bytes around, so it preserves the quotient space byte-wise.
MixColumns multiplies some bytes by $M_2$ and $M_3$. For the quotient space to be preserved, the linear matrices $M_2$ and $M_3$ must map elements of U to elements of U. If you test this, you will find it is strictly true: U is an invariant subspace of both $M_2$ and $M_3$!

Because every single operation in the cipher preserves the equivalence classes modulo U, the entire AES encryption function $E_k(p)$ induces a well-defined “quotient cipher” $\tilde{E}_{\tilde{k}}(\tilde{p})$ operating exclusively on $(\mathbb{F}_2^8/U)^{16}$ instead of $(\mathbb{F}_2^8)^{16}$.

If two keys $k_1$ and $k_2$ belong to the same cosets byte-by-byte (meaning $k_1 \equiv k_2 \pmod U$), then for any plaintext p, their ciphertexts will also belong to the same cosets byte-by-byte $(c_1 \equiv c_2 \pmod U)$.

5. The Attack

We are missing 8 bytes of the key. A standard brute-force attack would require testing $256^8 = 2^{64}$ keys, which is computationally out of reach. However, thanks to the quotient space vulnerability, we can split this brute-force attack into two sequential phases of $16^8 = 2^{32}$ operations each, completely breaking the cipher in a matter of minutes.

Stage 1: Quotient Space Evaluation Instead of guessing the exact values of the 8 missing key bytes, we only guess their projections in $\mathbb{F}_2^8/U$. Since there are only 16 cosets, there are $16^8 = 2^{32}$ possible projected keys. For each guess k, we construct a full key using the representatives from T and run the encryption on our known plaintext. We then project the resulting ciphertext and compare it to the projection of our known target ciphertext.

E_{\tilde{k}}(p) \equiv c_{target} \pmod U

When this equivalence holds for all 16 bytes of the state, we have found the correct quotient key. This first phase requires $2^{32}$ encryptions.

Stage 2: Exact Fiber Inversion Once Stage 1 is complete, we know the exact coset for each of the 8 missing key bytes. Let $t_i \in T$ be the correct representative found for the i-th missing key byte. We now know that the true key byte $k_i$ must be of the form:

k_i = t_i \oplus u_i

for some unknown $u_i \in U$. Since U contains exactly 16 elements, there are again $16^8 = 2^{32}$ possible exact keys within this specific “fiber” (the set of keys that project to our known quotient key). We simply iterate over all possible combinations of $u_i \in U$, compute the exact key k, and encrypt the plaintext. This time, we check for a strict equality:

E_k(p) = c_{target}

When the exact match is found, we have recovered the missing half of the key! This second phase requires an additional $2^{32}$ AES encryptions.

Implementation Here is my C++ implementation of the attack, which finds the corresponding key for a given plaintext-ciphertext pair in about a minute.

#include <atomic>
#include <chrono>
#include <cstdint>
#include <omp.h>
#include <string>
#include <vector>

std::vector<uint8_t> hex_to_bytes(const std::string &hex) {
  std::vector<uint8_t> bytes;
  for (size_t i = 0; i < hex.length(); i += 2) {
    bytes.push_back((uint8_t)strtol(hex.substr(i, 2).c_str(), NULL, 16));
  }
  return bytes;
}
const std::string P_hex = "5aea1a4940f81959...56a1b99424afe97a";
const std::string F_hex = "6eb9011a067eda81...20e3f20a4f9330ce";
const std::string M2_hex = "c08883cb5d151e56...fab2b9f1672f246c";
const std::string M3_hex = "40090148d99098d1...82cbc38a1b525a13";
const std::string S_hex = "c7ba17465263e1db...81ab87825d10bedf";
const std::string RP_hex = "ebeae9efe3fbcbab...d68335b9dc5fd773";
const std::string SK_hex = "473a97c6d2e3615b...012b0702dd903e5f";
const std::string P_known_hex = "24d2e9d44b2b01a50f3d5ca80c4f57ab";
const std::string C_known_hex = "3fe8937a109576303437c9f1dd8a2a0d";
inline void key_schedule(const uint8_t *rk0, const uint8_t *SK,
                         const std::vector<uint8_t> &RCON, uint8_t *rk) {
  for (int i = 0; i < 16; i++)
    rk[i] = rk0[i];
  for (int i = 4; i < 44; i++) {
    uint8_t temp[4];
    temp[0] = rk[(i - 1) * 4 + 0];
    temp[1] = rk[(i - 1) * 4 + 1];
    temp[2] = rk[(i - 1) * 4 + 2];
    temp[3] = rk[(i - 1) * 4 + 3];
    if (i % 4 == 0) {
      uint8_t t0 = SK[temp[1]] ^ RCON[i / 4 - 1];
      uint8_t t1 = SK[temp[2]];
      uint8_t t2 = SK[temp[3]];
      uint8_t t3 = SK[temp[0]];
      temp[0] = t0;
      temp[1] = t1;
      temp[2] = t2;
      temp[3] = t3;
    }
    for (int j = 0; j < 4; j++) {
      rk[i * 4 + j] = rk[(i - 4) * 4 + j] ^ temp[j];
    }
  }
}
inline void encrypt_to_ark10(const uint8_t *rk, const uint8_t *P_mapped,
                             const uint8_t *S, const uint8_t *M2,
                             const uint8_t *M3, uint8_t *state) {
  for (int i = 0; i < 16; i++)
    state[i] = P_mapped[i];
  uint8_t new_state[16];
  const uint8_t SR_MAP[16] = {0, 5,  10, 15, 4,  9, 14, 3,
                              8, 13, 2,  7,  12, 1, 6,  11};
  for (int r_idx = 0; r_idx < 10; r_idx++) {
    for (int i = 0; i < 16; i++)
      state[i] ^= rk[r_idx * 16 + i];
    for (int i = 0; i < 16; i++)
      state[i] = S[state[i]];
    for (int i = 0; i < 16; i++)
      new_state[i] = state[SR_MAP[i]];
    for (int c = 0; c < 4; c++) {
      uint8_t a0 = new_state[c * 4 + 0];
      uint8_t a1 = new_state[c * 4 + 1];
      uint8_t a2 = new_state[c * 4 + 2];
      uint8_t a3 = new_state[c * 4 + 3];
      state[c * 4 + 0] = M2[a0] ^ M3[a1] ^ a2 ^ a3;
      state[c * 4 + 1] = a0 ^ M2[a1] ^ M3[a2] ^ a3;
      state[c * 4 + 2] = a0 ^ a1 ^ M2[a2] ^ M3[a3];
      state[c * 4 + 3] = M3[a0] ^ a1 ^ a2 ^ M2[a3];
    }
  }
  for (int i = 0; i < 16; i++)
    state[i] ^= rk[160 + i];
}
int main() {
  omp_set_num_threads(16);
  std::vector<uint8_t> P = hex_to_bytes(P_hex);
  std::vector<uint8_t> F = hex_to_bytes(F_hex);
  std::vector<uint8_t> M2 = hex_to_bytes(M2_hex);
  std::vector<uint8_t> M3 = hex_to_bytes(M3_hex);
  std::vector<uint8_t> S = hex_to_bytes(S_hex);
  std::vector<uint8_t> RP = hex_to_bytes(RP_hex);
  std::vector<uint8_t> SK = hex_to_bytes(SK_hex);
  std::vector<uint8_t> P_known = hex_to_bytes(P_known_hex);
  std::vector<uint8_t> C_known = hex_to_bytes(C_known_hex);
  std::vector<uint8_t> F_inv(256), RP_inv(256);
  for (int i = 0; i < 256; i++) {
    F_inv[F[i]] = i;
    RP_inv[RP[i]] = i;
  }
  uint8_t rcon_base[10] = {1, 2, 4, 8, 16, 32, 64, 128, 27, 54};
  std::vector<uint8_t> RCON(10);
  for (int i = 0; i < 10; i++)
    RCON[i] = RP[rcon_base[i]];
  uint8_t Va[16] = {0,   3,   8,   11,  36,  39,  44,  47,
                    129, 130, 137, 138, 165, 166, 173, 174};
  bool is_Va[256] = {false};
  for (int i = 0; i < 16; i++)
    is_Va[Va[i]] = true;
  std::vector<uint8_t> T;
  bool seen[256] = {false};
  for (int i = 0; i < 256; i++) {
    if (!seen[i]) {
      T.push_back(i);
      for (int j = 0; j < 16; j++)
        seen[i ^ Va[j]] = true;
    }
  }
  std::vector<uint8_t> P_mapped(16), C_target(16);
  for (int c = 0; c < 4; c++) {
    for (int r = 0; r < 4; r++) {
      P_mapped[c * 4 + r] = P[P_known[c * 4 + r]];
      C_target[c * 4 + r] = F_inv[C_known[c * 4 + r]];
    }
  }
  uint8_t fixed_rk0[16] = {0};
  uint8_t prefix[] = {'F', 'C', 'S', 'C', '{', '{'};
  uint8_t suffix[] = {'}', '}'};
  for (int i = 0; i < 6; i++)
    fixed_rk0[i] = RP[prefix[i]];
  for (int i = 14; i < 16; i++)
    fixed_rk0[i] = RP[suffix[i - 14]];
  // Stage 1
  printf("[*] Stage 1...\n");
  std::atomic<bool> found1(false);
  uint8_t quotient_key[8] = {0};
  auto start_time = std::chrono::high_resolution_clock::now();
#pragma omp parallel for schedule(dynamic, 65536)
  for (uint64_t guess = 0; guess < 0x100000000ULL; guess++) {
    if (found1.load(std::memory_order_relaxed))
      continue;
    if ((guess & 0x3FFFFFF) == 0 && guess > 0) {
      printf(" [.] Progress: %.2f%%\n", (double)guess / 42949672.96);
    }
    uint8_t rk0[16];
    for (int i = 0; i < 16; i++)
      rk0[i] = fixed_rk0[i];
    rk0[6] = T[(guess >> 28) & 0xF];
    rk0[7] = T[(guess >> 24) & 0xF];
    rk0[8] = T[(guess >> 20) & 0xF];
    rk0[9] = T[(guess >> 16) & 0xF];
    rk0[10] = T[(guess >> 12) & 0xF];
    rk0[11] = T[(guess >> 8) & 0xF];
    rk0[12] = T[(guess >> 4) & 0xF];
    rk0[13] = T[guess & 0xF];
    uint8_t rk[176];
    key_schedule(rk0, SK.data(), RCON, rk);
    uint8_t state[16];
    encrypt_to_ark10(rk, P_mapped.data(), S.data(), M2.data(), M3.data(),
                     state);
    bool matched = true;
    for (int i = 0; i < 16; i++) {
      if (!is_Va[state[i] ^ C_target[i]]) {
        matched = false;
        break;
      }
    }
    if (matched) {
      found1.store(true, std::memory_order_relaxed);
      for (int i = 0; i < 8; i++)
        quotient_key[i] = rk0[6 + i];
      printf("[+] Stage 1 Done!\n");
    }
  }
  // Stage 2
  printf("\n[*] Stage 2...\n");
  std::atomic<bool> found2(false);
  uint8_t exact_key[8] = {0};
#pragma omp parallel for schedule(dynamic, 65536)
  for (uint64_t delta_guess = 0; delta_guess < 0x100000000ULL; delta_guess++) {
    if (found2.load(std::memory_order_relaxed))
      continue;
    if ((delta_guess & 0x3FFFFFF) == 0 && delta_guess > 0) {
      printf(" [.] Progress: %.2f%%\n", (double)delta_guess / 42949672.96);
    }
    uint8_t rk0[16];
    for (int i = 0; i < 16; i++)
      rk0[i] = fixed_rk0[i];
    rk0[6] = quotient_key[0] ^ Va[(delta_guess >> 28) & 0xF];
    rk0[7] = quotient_key[1] ^ Va[(delta_guess >> 24) & 0xF];
    rk0[8] = quotient_key[2] ^ Va[(delta_guess >> 20) & 0xF];
    rk0[9] = quotient_key[3] ^ Va[(delta_guess >> 16) & 0xF];
    rk0[10] = quotient_key[4] ^ Va[(delta_guess >> 12) & 0xF];
    rk0[11] = quotient_key[5] ^ Va[(delta_guess >> 8) & 0xF];
    rk0[12] = quotient_key[6] ^ Va[(delta_guess >> 4) & 0xF];
    rk0[13] = quotient_key[7] ^ Va[delta_guess & 0xF];
    uint8_t rk[176];
    key_schedule(rk0, SK.data(), RCON, rk);
    uint8_t state[16];
    encrypt_to_ark10(rk, P_mapped.data(), S.data(), M2.data(), M3.data(),
                     state);
    bool matched = true;
    for (int i = 0; i < 16; i++) {
      if (state[i] != C_target[i]) {
        matched = false;
        break;
      }
    }
    if (matched) {
      found2.store(true, std::memory_order_relaxed);
      for (int i = 0; i < 8; i++)
        exact_key[i] = rk0[6 + i];
    }
  }
  auto end_time = std::chrono::high_resolution_clock::now();
  std::chrono::duration<double> diff = end_time - start_time;
  printf("\n[+] Stage 2 Done. Key found in %.2f seconds.\n", diff.count());
  printf("[!] Master key bytes: ");
  for (int i = 0; i < 8; i++) {
    uint8_t dec = RP_inv[exact_key[i]];
    printf("%02x", dec);
  }
  puts("");
  return 0;
}

References

Luke O’Connor. On the distribution of characteristics in bijective mappings. In Advances in Cryptology - EUROCRYPT ‘93, volume 765, pages 360-370. Springer, 1993. ↩︎
James A. Muir. A tutorial on white-box AES. In Progress in Cryptology - INDOCRYPT 2013, volume 8250 of Lecture Notes in Computer Science, pages 209-229. Springer, 2013. ↩︎
Olivier Billet, Henri Gilbert, and Charaf Ech-Chatbi. Cryptanalysis of a white box AES implementation. In Selected Areas in Cryptography - SAC 2004, volume 3357 of Lecture Notes in Computer Science, pages 227-240. Springer, 2004. ↩︎ ↩︎