Using Constraint Programming to solve a Cryptanalytic Problem

1 Introduction

Since 2001, AES (Advanced Encryption Standard) is the encryption standard for block ciphers [FIPS 197, 2001]. It guarantees communication confidentiality by using a secret key K to encode an original plaintext X into a ciphertext AESK(X), in such a way that the ciphertext can further be decoded into the original one using the same key, i.e., X = AES-1 K (AESK(X)). Cryptanalysis aims at testing whether confidentiality is actually guaranteed. In particular, differential cryptanalysis [Biham and Shamir, 1991] evaluates whether it is possible to find the key within a reasonable number of trials by considering plaintext pairs (X,X′ ) and studying the propagation of the initial difference X ⊕X′ between X and X′ while going through the ciphering process (where ⊕ is the xor operator). Today, differential cryptanalysis is public knowledge, and block ciphers such as AES have proven bounds against differential attacks. Hence, [Biham, 1993] proposed a new type of attack called related-key attack that allows an attacker to inject differences not only between the plaintexts X and X′ but also between the keys K and K′ (even if the secret key K stays unknown from the attacker).

To mount related-key attacks, the cryptanalist must find optimal related-key differentials, i.e., a plaintext difference δX and a key difference δK that maximize the probability that, for randomly chosen plaintext X and key K, the difference X ⊕ δX in the input plaintexts becomes the difference AESK(X) ⊕AESK⊕δK(X′ ) in the output ciphertexts. Finding the optimal related-key differentials for AES is a highly combinatorial problem that hardly scales. Two main approaches have been proposed to solve this problem: a graph traversal approach [Fouque et al., 2013], and a Branch & Bound approach [Biryukov and Nikolic, 2010]. The approach of [Fouque et al., 2013] requires about 60 GB of memory when the key has 128 bits, and it has not been extended to larger keys. The approach of [Biryukov and Nikolic, 2010] only takes several megabytes of memory, but it requires several days of computation when the key has 128 bits, and several weeks when the key has 192 bits.

∗This research was conducted with the support of the FEDER program of 2014-2020 and the region council of Auvergne

During the process of designing new ciphers, this search generally needs to be performed several times, so it is desirable that it can be done rather quickly. Another point that should not be neglected is the time needed to design and implement these approaches: To ensure that the computation is completed within a "reasonable" amount of time, it is necessary to reduce the branching by introducing clever reasoning. Of course, this hard task is also likely to introduce bugs, and checking the correctness or the optimality of the computed solutions may not be so easy. Finally, reproducibility may also be an issue. Other researchers may want to adapt these algorithms to other problems, with some common features but also some differences, and this may again be very difficult and time-consuming.

In [Gerault et al., 2016], we proposed to use Constraint Programming (CP) to solve this problem. When using CP to solve a problem, one simply has to model the problem by means of constraints. Then, this model is solved by generic solvers usually based on a Branch & Propagate approach: The search space is explored by building a search tree, and constraints are propagated at each tree node in order to prune branches. CP opens new perspectives for this kind of cryptanalysis problems. First, it is very competitive with dedicated approaches: When the key has 128 bits, Chuffed [Chu and Stuckey, 2014] is able to find optimal solutions in less than one hour. Second, the CP model is easier to check or re-use than a full program that not only describes the problem to solve, but also how to solve it. Actually, CP allowed us to prove that a solution claimed to be optimal in [Biryukov and Nikolic, 2010; Fouque et al., 2013] is not optimal by providing a better solution.

2 Problem Statement

AES block cipher. AES ciphers blocks of length n = 128 bits, and each block is a 4 × 4 matrix of bytes. The length of keys is l ∈ {128, 192, 256}. In this paper, we only consider Operations applied at each round i ∈ [0, r - 1]:

keys of length l = 128, and keys are 4 × 4 matrices of bytes. Given a 4×4 matrix of bytes M, we note M[j][k] the byte at row j ∈ [0, 3] and column k ∈ [0, 3].

AES is an iterative process, and we note Xi the ciphertext at the beginning of round i ∈ [0, r]. Each round is composed of the following operations, as displayed in Fig. 1:

• SubBytes S. S is a non-linear permutation which is applied on each byte of Xi separately, i.e., ∀j, k ∈ [0, 3], Xi[j][k] is replaced by S(Xi[j][k]), according to a lookup table. We note SXi = S(Xi).
• ShiftRows SR. SR is a linear mapping that rotates on the left by 1 byte position (resp. 2 and 3 byte positions) the second row (resp. third and fourth rows) of SXi.
• MixColumnsMC. MC is a linear mapping that multiplies each column of SR(SXi) by a 4×4 fixed matrix chosen for its good properties of diffusion [Daemen and Rijmen, 2002]. In particular, it has the Maximum Distance Separable (MDS) property: For each column, the total number of bytes which are different from 0, before and after applying MC, is either equal to 0 or strictly greater than 4. We note Yi = MC(SR(SXi)).
• KeySchedule KS. The subkey at round 0 is the initial key, i.e., K0 = K. For each round i ∈ [0, r - 1], the subkey Ki+1 is generated from Ki by applying KS. It first replaces each byte Ki[j][3] of the last column by S(Ki[j][3]) (where S is the SubBytes operator), and we note SKi[j][3] = S(Ki[j][3]). Then, each column of Ki+1 is obtained by performing a xor operation between bytes coming fromKi, SKi, or Ki+1.
• AddRoundKey ARK. ARK performs a xor operation between bytes of Yi and subkey Ki+1 to obtain Xi+1.

Let us note Bytes the set of all bytes (for all rounds), i.e.,

Optimal related-key differentials. Let us note δB = B ⊕ B′ the difference between two bytes B and B′ , and diffBytes = {δB | B ∈ Bytes} the set of all differential bytes.

Mounting attacks requires finding a related-key differential characteristic, i.e. a plaintext difference δX =X⊕X′ and a key difference δK = K ⊕K′ , such that δX becomes δSXr after r rounds with a probability as hight as possible. This can be achieved by tracking the propagation of the initial differences through the cipher. Once such an optimal differential characteristic is found, the cryptanalist asks for the encryption of messages satisfying the input difference. When this difference propagates as expected, he can infer informations leading to a recovery of the secret key K with less workload than exhaustive search. The difficulty of recovering K decreases as the probability of the used differential characteristic increases. The AES operators SR,MC,ARK, and KS are linear, i.e., they propagate differences in a deterministic way (with probability 1). However, the S operator is not linear: Given a byte difference B ⊕ B′ = δB, the probability that δB becomes S(B) ⊕S(B′ ) = δSB is P r(δB → δSB). It is equal to 1 if δB = 0 (i.e., B = B′ ). However, if δB ̸= 0, then Pr(δB → δSB) ∈ { 2 256 , 4 256 }.

The probability that δX becomes δSXr is equal to the product of all Pr(δB → δSB) such that δB (resp. δSB) is a byte difference before (resp. after) passing through the S operator when ciphering X with K and X′ with K′ , i.e., the product of all Pr(δXi[j][k] → δSXi[j][k]) and all Pr(δKi[j][3] → δSKi[j][3]) with i ∈ [0, r] and j, k ∈ [0, 3]. The goal is to find δX and δK that maximize this probability.

Two step solving process. To find δX and δK, we search for the values of all differential bytes in diffBytes. Both [Biryukov and Nikolic, 2010] and [Fouque et al., 2013] propose to solve this problem in two steps. In Step 1, a Boolean variable ∆B is associated with every differential byte δB ∈ diffBytes such that ∆B = 0 ⇔ δB = 0 and

∆B = 1 ⇔ δB ∈ [1, 255]. The goal of Step 1 is to find a Boolean solution that assigns values to Boolean variables such that the AES transformation rules are satisfied. During this first step, the SubBytes operation S is not considered. Indeed, it does not introduce nor remove differences. Therefore, we have ∆Xi[j][k] = ∆SXi[j][k] and ∆Ki[j][3] = ∆SKi[j][3]. As we search for a solution with maximal probability, the goal of Step 1 is to search for a Boolean solution which minimizes the number of variables ∆Xi[j][k] and ∆Ki[j][3] which are set to 1.

In Step 2, the Boolean solution is transformed into a byte solution: For each differential byte δB ∈ diffBytes, if the corresponding Boolean variable ∆B is assigned to 0, then δB is also assigned to 0; otherwise, we search for a byte value in [1, 255] to be assigned to δB such that the AES transformation rules are satisfied and the probability is maximized. Note that some Boolean solutions may not be transformable into byte solutions. These Boolean solutions are said to be byte-inconsistent.

3 Basic CPModel for Step 1

A first CP model for Step 1 may be derived from the AES transformation rules in a rather straightforward way. A CP model is defined by a set of variables, such that each variable x has a domain D(x), and a set of constraints, i.e., relations that restrict the values that may be simultaneously assigned to the variables.

Variables. For each differential byte δB ∈ diffBytes, we define a Boolean variable ∆B whose domain is D(∆B) = {0, 1}: it is assigned to 0 if δB = 0, and to 1 otherwise.

XOR constraint. We first define a XOR constraint for ARK and KS. Let us consider three differential bytes δA, δB and δC such that δA⊕δB = δC. If δA = δB = 0, then δC = 0. If (δA = 0 and δB ̸= 0) or (δA ̸= 0 and δB = 0) then δC ̸= 0. However, if δA ̸= 0 and δB ̸= 0, then we cannot know if δC is equal to 0 or not: This depends on whether δA = δB or not. When abstracting differential bytes δA, δB and δC with Boolean variables ∆A, ∆B and ∆C (which only model the fact that there is a difference or not), we obtain the following definition of the XOR constraint:

XOR(∆A, ∆B, ∆C) ⇔ ∆A + ∆B + ∆C ̸= 1.

AddRoundKey and KeySchedule constraints. Both ARK and KS are modeled with XOR constraints between ∆Xi, ∆Yi, and ∆Ki variables (for ARK), and between ∆Ki and ∆SKi variables (for KS).

ShiftRows and MixColumns. SR simply shift variables. The MDS property of MC is ensured by posting a constraint on the sum of some variables of ∆Xi and ∆Yi variables, which must belong to the set {0, 5, 6, 7, 8}.

Objective function. We introduce an integer variable objStep1 that must be minimized, and we post an equality constraint between objStep1 and the sum of Boolean variables on which a non linear S operation is performed (all variables ∆Xi[j][k] and ∆Ki[j][3] with i ∈ [0, r] and j, k ∈ [0, 3]).

Let v be the optimal value of objStep1. It may happen that none of the Boolean solutions with objStep1 = v is byte-consistent, or that the maximal probability p of Boolean solutions with objStep1 = v is such that it is possible to have a better probability with a larger value for objStep1 (i.e., p < ( 4 256 )v+1 ). In this case, we need to search for new Boolean solutions, such that objStep1 is minimal while being strictly greater than v. This is done by adding the constraint objStep1 > v before solving again Step 1.

Limitations of the basic CP model. This basic CP model is complete, i.e., for any solution at the byte level (on δ variables), there exists a solution at the Boolean level (on ∆ variables). However, preliminary experiments have shown us that there is a huge number of Boolean solutions which are byte inconsistent. For example, when the number of rounds is r = 4, the optimal cost is objStep1 = 11, and there are more than 90 millions of Boolean solutions with objStep1 = 11. However, none of these solutions is byte-consistent. In this case, most of the Step 1 solving time is spent at generating useless Boolean solutions which are discarded in Step 2.

4 Additional Constraints for Step 1

When reasoning at the Boolean level, many solutions are not byte-consistent because constraints at the byte level have been ignored. To propagate some properties at the byte level, we introduce new variables and constraints that model equality relations between differential bytes.

New equality variables. For each couple of differential bytes (δA, δB) ∈ diffBytes 2 , we introduce a Boolean equality variable EQδA,δB which is equal to 1 if δA = δB, and to 0 otherwise. These variables are constrained to define an equivalence relation by adding a symmetry constraint (EQδA,δB = EQδB,δA ) and a transitivity constraint (if EQδA,δB = EQδB,δC = 1 then EQδA,δC = 1). Also, EQ variables are related to ∆ variables by adding the constraints:

Revisiting the XOR constraint. When defining the constraint XOR(∆A,∆B,∆C), if ∆A = ∆B = 1, then we cannot know if ∆C is equal to 0 or 1. However, whenever ∆C = 0 (resp. ∆C = 1), we know for sure that the corresponding byte δC is equal to 0 (resp. different from 0), meaning that the two bytes δA and δB are equal (resp. different), i.e., that EQδA,δB = 1 (resp. EQδA,δB = 0). The same reasoning may be done for ∆A and ∆B because (δA⊕ δB = δC) ⇔ (δB ⊕ δC = δA) ⇔ (δA⊕ δC = δB). Therefore, we redefine the XOR constraint as follows:

Propagation of MDS at Byte Level. The MDS property ensures that, for each column, the total number of bytes which are different from 0, before and after applying MC, is either equal to 0 or strictly greater than 4. This property also holds for any xor difference between two different columns of X and Y matrices. To propagate this property, for each pair of columns in X and Y matrices, we add a constraint on the sum of equality variables between bytes of these columns.

Constraints derived from KS. The KeySchedule mainly performs xor and S operations. As a consequence, each byte δKi[j][k] may be expressed as a xor between bytes of the original key difference δK0, and bytes of δSKi-1 (which are differences of key bytes that have passed through S during the previous round). Hence, for each byte δKi[j][k], we precompute the set V (i, j, k) such that V (i, j, k) only contains bytes of δK0 and δSKi-1 and δKi[j][k] = ⊕δA∈V (i,j,k) δA. For each set V (i, j, k), we introduce a set variable V1(i, j, k) which is constrained to contain the subset of V (i, j, k) corresponding to the Boolean variables equal to 1. We use these set variables to infer that two differential key bytes that have the same V1 set are equal. Also, if V1(i, j, k) is empty (resp. contains one or two elements), we infer that ∆Ki[j][k] is equal to 0 (resp. a variable, or a xor between 2 variables).

5 CPModel for Step 2

Given a Boolean solution for Step 1, Step 2 aims at searching for the byte-consistent solution with the highest probability (or proving that there is no byte-consistent solution). Hence, for each Boolean variable ∆A of Step 1, we define an integer variable δA whose domain depends on the value of ∆A: If ∆A = 0, then D(δA) = {0} (i.e., δA is also assigned to 0); otherwise, D(δA) = [1, 255]. As we look for a byte-consistent solution with maximal probability, we also add an integer variable PA for each byte A that passes through the S operator, i.e., Xi[j][k] and Ki[j][3]: This variable corresponds to the base 2 logarithm of the probability Pr(δA → δSA) of obtaining the S output difference δSA when the S input difference is δA. The domain of these variables depends on the value of ∆A in the Step 1 solution: If ∆A = 0, then Pr(0 → 0) = 1 and therefore D(PA) = {0}; otherwise, Pr(δA → δSA) ∈ { 2 256 , 4 256 } and D(PA) = {-7, -6}.

The constraints basically follow the AES operations to relate variables, as described for Step 1, but consider the definition of the operations at the byte level, instead of the Boolean level. A main difference is that the SubBytes operation, which has no effect at the Boolean level, must be modeled at the byte level. This is done thanks to a ternary table constraint which extensively lists all triples (X, Y, P) such that there exist two bytes B1 and B2 whose difference before and after passing through S is equal to X and Y , respectively, and such that P is the base 2 logarithm of the probability of this transformation. The goal is to find a byte-consistent solution with maximal differential probability. As we consider logarithms, this amounts to searching for a solution that maximizes the sum of all PA variables.

6 Results and Conclusion

The CP model for Step 1 was implemented with the MiniZinc modeling language [Nethercote et al., 2007], and we compared three CP solvers: Gecode [Gecode Team, 2006], Choco 4 [Prudhomme and Fages, 2013], and Chuffed [Chu and Stuckey, 2014]. The best results were obtained with Chuffed. The basic CP model (described in Section 3) does not scale well because it generates a huge number of Boolean solutions which are not byte-consistent. When adding the additional constraints described in Section 4, most of these byte-inconsistent solutions are filtered out, and Chuffed is able to solve all instances in less than one hour. The CP model for Step 2 was implemented with Choco 3. It solves all instances (i.e., finds an optimal byte-consistent solution for each Boolean solution of Step 1 when it exists) in a few seconds. As a conclusion, CP solvers are much faster than the Branch & Bound approach of [Biryukov and Nikolic, 2010], which needs several days to solve these instances. It is also faster and much less memory consuming than the approach of [Fouque et al., 2013], that needs 60GB and 30 minutes on a 12-core computer to pre-compute the graph.

New results for differential cryptanalysis. For r = 4 rounds, we have found a byte-consistent solution with objStep1 = 12 and a probability equal to 2-79 . This solution is better than the solution claimed to be optimal in [Biryukov and Nikolic, 2010] and [Fouque et al., 2013]: In these papers, the authors say that the best byte-consistent solution has objStep1 = 13, and a probability equal to 2-81 .

We have shown how to extend our CP models to AES with longer keys (such that l ∈ {192, 256}) in [G´erault et al., 2017]. These models allowed us to find optimal solutions for all possible instances of AES with l ∈ {128, 192, 256} in less than 35 hours for Step 1, and in less than 6 minutes for Step 2. This is a clear improvement with respect to existing work as the approach of [Fouque et al., 2013] cannot be extended to l > 128 due to its memory complexity, and the approach of [Biryukov and Nikolic, 2010] needs several weeks to solve instances with l = 192. Furthermore, we have shown that the solution proposed in [Biryukov and Nikolic, 2010] for l = 192 and r = 11 is inconsistent. We have also found better solutions when l = 256, and we have computed the actual optimal solution for AES with l = 256. Its probability is 2-146 (instead of 2-154 for the solution of [Biryukov et al., 2009]). Using this solution, we improved the relatedkey distinguisher and the basic related-key differential attack on the full AES-256 by a factor 26 and the q-multicollisions by a factor 2 (see [G´erault et al., 2017] for more details).

These cryptanalysis problems open new and exciting challenges for the CP community. In particular, these problems are not easy to model. More precisely, naive CP models such as the one described in Section 3 may not scale well. The introduction of equality constraints at the byte level drastically improves the solving process, but these constraints are not straightforward to find and implement. Hence, a challenge is to define new CP frameworks, dedicated to this kind of cryptanalysis problems, in order to ease the development of efficient CP models for these problems.