SoK: Modeling for Large S-boxes Oriented to Differential Probabilities and Linear Correlations (Long Paper)

. Automatic methods for differential and linear characteristic search are well-established at the moment. Typically, the designers of novel ciphers also give preliminary analytical findings for analysing the differential and linear properties using automatic techniques. However, neither MILP-based nor SAT/SMT-based approaches have fully resolved the problem of searching for actual differential and linear characteristics of ciphers with large S-boxes. To tackle the issue, we present three strategies for developing SAT models for 8-bit S-boxes that are geared toward differential probabilities and linear correlations. While these approaches cannot guarantee a minimum model size, the time needed to obtain models is drastically reduced. The newly proposed SAT model for large S-boxes enables us to establish that the upper bound on the differential probability for 14 rounds of SKINNY-128 is 2 − 131 , thereby completing the unsuccessful work of Abdelkhalek et al. We also analyse the seven AES -based constructions C1 - C7 designed by Jean and Nikolić and compute the minimum number of active S-boxes necessary to cause an internal collision using the SAT method. For two constructions C3 and C5 , the current lower bound on the number of active S-boxes is increased, resulting in a more precise security analysis for these two structures.

Automatic methods for differential and linear characteristic search are currently wellestablished. Typically, the designers of novel ciphers also provide preliminary analytical results for analysing the differential and linear properties via automatic approaches. 1 However, neither MILP-based nor SAT/SMT-based approaches have yet completely solved the challenge of searching for true differential and linear characteristics of ciphers with large S-boxes.
[AST + 17] investigated the differential distribution table (DDT) of large S-boxes and provided the first method for constructing MILP models for large S-boxes in the development of MILP modelling. The new approach was implemented on SKINNY-128 [BJK + 16] and two AES-based constructions in [JN16]. The upper bound on differential probability for 14 rounds of SKINNY-128 is not completely resolved due to the extremely lengthy runtime of the MILP optimiser. While the minimum number of differential active S-boxes for one AES-based construction is raised, the precise lower bound on the number of active S-boxes is not entirely known. Despite the fact that Abdelkhalek et al. made it possible to describe large S-boxes, the model for the * -DDT, which is a reduced version of the DDT in which all non-zero entries are replaced by one, of the S-box in AES consists of 8302 linear inequalities, which is sometimes too complicated for practical applications. To tackle this issue, Boura and Coggia [BC20] proposed two strategies for constructing efficient algorithms for modelling 8-bit S-boxes with the minimum amount of linear inequalities. Notably, the enhanced MILP models for the S-boxes of AES and SKINNY-128 are not used to search for differential characteristics of the two ciphers.
In the SAT/SMT modelling development, Ankele and Kölbl [AK18] stated that their S-box modelling technique was applicable to 8-bit S-boxes; yet, all of the ciphers analysed in the study use 4-bit S-boxes. The SMT model suggested by Liu et al. [LLL + 19] can manage 8-bit S-boxes. The authors used their model to investigate the security of AES and SKINNY-128 against the impossible differential attack; once again, the differential probability of the two ciphers is not taken into account.
Motivated by resolving the aforementioned unsolved issues pertaining to SKINNY-128 and AES-based constructions, we want to identify a way for efficiently creating SAT models of large S-boxes.
Our contributions. After examining the evolution of MILP and SAT/SMT modelling techniques, we see that the search for actual differential and linear characteristics for ciphers with large S-boxes is not conclusive. We intend to provide effective ways for generating SAT models for large S-boxes.
Three strategies for creating SAT models for large S-boxes oriented to differential probabilities and linear correlations are proposed. Although these methods cannot ensure a minimum model size, the time required to obtain models is greatly decreased. The three tactics are summarised in the following.
1. Since developing SAT models for S-boxes is equal to finding a simplification of a given Boolean function, we attempt to strike a balance between the degree of simplification and the runtime. The first strategy utilises the option of the ESPRESSO logic minimizer [BHH + 82, BHMS84], which is an algorithm used to simplify Boolean functions.
2. Noting that the complexity of simplification rises exponentially with the number of input variables of the given Boolean function, we attempt to lower the size of the function that ESPRESSO delivers. The primary concept of the second technique is to divide the description of a big S-box into two steps, so transforming the simplification of large-scale Boolean functions into the simplification of two relatively small-scale functions.
3. In the test, we observe that simplification of a large-scale function is not difficult if the number of terms is not extremely enormous. This discovery motivates the partitioning mechanism of third strategy. Specifically, the target n-bit Boolean function f is subdivided into many n-bit functions f 0 , f 1 , . . ., and f ℓ−1 , where ℓ is an integer, so that the conjunction ℓ−1 i=0 f i is identical to f . Then, the simplification of f is converted into the simplification of f 0 , f 1 , . . ., and f ℓ−1 , hence accelerating the production of the SAT model.
SKINNY-128 is used to evaluate the three strategies. Numerous SAT models are applied to calculate the upper bound of differential probability for SKINNY-128 from 1 to 14 rounds. After analysing all the test results, we believe that the first encoding technique achieves a good compromise between the runtime of ESPRESSO and the execution time of the SAT solver. The test results also reveal that lowering the number of constraints in SAT problems does not necessarily reduce the execution time of the SAT solver, echoing the observation made by Sasaki and Todo [ST17] in the MILP approach. In addition, we empirically demonstrate that reducing the number of variables in SAT problems does not always reduce the execution time of the SAT solver.
The newly presented SAT model for large S-boxes lets us to demonstrate that the upper bound on the differential probability for 14 rounds of SKINNY-128 is 2 −131 , thereby completing the failed work of Abdelkhalek et al. Additionally, the related-key differential properties of both versions of PIPO [KJK + 20] are investigated. For PIPO-128, we identify 1792 full-round characteristics with a probability of 2 −24 , extending the findings of Yadav and Kumar [YK21]. For PIPO-256, we are the first to report full-round differential characteristics.
The encoding approach for large S-boxes is applied to seven AES-based constructions C1 -C7 devised by Jean and Nikolić [JN16]. The SAT method is used to calculate the minimum number of active S-boxes required to induce an internal collision. The test results for all seven constructions, ranging from two to eight stages, are shown in Table 9. For two constructions C3 and C5, the current lower bound on the number of active S-boxes is increased, resulting in a more precise security analysis for these two structures. For C1, C2, C4, and C6, we validate the previously determined lower bound for the number of active S-boxes. In addition, all differential patterns for C1 -C5 that have a minimum number of active S-boxes are supplied.
Outline. The remainder of the paper is structured as follows. In Section 2, notations and knowledge related to the automatic search method are presented. Section 3 is devoted to discussing the MILP modelling work for large S-boxes, while Section 4 discusses the SAT modelling progress for large S-boxes. In Section 5, we propose three efficient methods for building SAT models for large S-boxes. In Section 6, the newly suggested encoding techniques are implemented on SKINNY-128, and their performance is compared. Based on the outcome of the test around SKINNY-128, we recognise that the first modelling technique is superior, and it is subsequently implemented in the studies of PIPO and seven AES-based constructions in Sections 7 and 8, respectively. Section 9 is the conclusion of the paper. Supplementary Material of this paper can be found at https://github.com/ SunLing134340/Large-Sbox-Supplementary-Material.

Preliminary
In this part, after presenting the notations used throughout the work, we will review two standard ways for simplifying Boolean functions.

Notations
The Hamming weight of a vector x = (x 0 , x 1 , . . . , x n−1 ) ∈ F n 2 , indicated by wt(x), is the number of non-zero elements in x, i.e., wt(x) = n−1 i=0 x i . A vector may be represented in binary representation, decimal representation, and hexadecimal representation; we utilise three fonts to differentiate between these representations: sans serif font for binary representation, roman typeface for decimal representation, and typewriter font for hexadecimal representation. Consequently, 1011, 11, and 0xb represent the same vector in this study.
The support of x, indicated by supp(x), is the set of coordinates of elements in x that are non-zero, i.e., supp(x) = {i ∈ [0, n − 1] | x i = 1 }. Given a Boolean function f over F n 2 , the support of f , represented as supp(f ), is the set of vectors in F n 2 , for which f is non-zero, i.e., supp(f ) = {x ∈ F n 2 | f (x) = 1 }.

Sum of Products and Product of Sums
For an n-variable Boolean function over x = (x 0 , x 1 , . . . , x n−1 ), a conjunction in which each of the n variables appears once (either in its complemented or uncomplemented form) is known as a minterm, i.e., it is a logical expression of n variables using just the complement operator (·) 1 and the conjunction operator (∧). For n variables, there are in total 2 n minterms. For u = (u 0 , u 1 , . . . , u n−1 ) ∈ F n 2 , the minterm Based on the notion of minterm, the canonical sum of products (SOP) form of f is When the output of f is true, minterms are included in the canonical SOP, which is also referred to as the sum of minterms.
A maxterm over n Boolean variables is defined as a disjunction (∨) where each of the n variables occurs once (either in its complemented or uncomplemented form). There are 2 n maxterms over n variables. Denote the maxterm where supp(f ) = {x ∈ F n 2 | f (x) = 0 } represents the complement of the support for the n-bit Boolean function f . It can be seen that M u (x) = 0 if and only if x = u. Therefore, the creation of the canonical POS of f may be seen as the elimination of any vector over which f is false.
For a given function in canonical SOP or POS form, computing the minimum SOP or POS form, which implies the minimum number of disjunctions or conjunctions, is a wellresearched issue in the design of digital logic gate circuits. Two conventional simplification techniques will be discussed in Sections 2.4 and 2.5, as they are closely related to our newly suggested encoding methods for large S-boxes in Section 5.

SAT Problem and POS Form Simplification
Mouha and Preneel [MP13] initiated the use of Boolean satisfiability problem (SAT) and satisfiability modulo theories (SMT) as a supporting tool in symmetric-key cryptography. The fundamental concept is to model the distinguisher searching issue in cryptanalysis as an SAT/SMT problem that can be addressed by SAT/SMT solvers such as STP [GD07] and Cryptominisat [SNC09].
All variables taken into account in the SAT problem are binary. A Boolean equation is a clause if it is a disjunction (∨) of one or more (possibly negated) Boolean variables; a Boolean formula is in conjunctive normal form (CNF) if it is a conjunction (∧) of one or more clauses. The CNF format is the input format that practically all SAT solvers claim to support. The essential step for the automatic search is thus to express the distinguisher searching issue using the CNF formula.
Since the SAT problem is NP-complete, all available solvers are only able to address problems of practical magnitude. The performance of the SAT solver may be affected if the SAT issue has an excessive number of clauses. As a result, reducing the number of clauses in the CNF format is a common issue when creating SAT problems regarding distinguisher search issues, particularly when constructing SAT models for S-boxes. It is observable that the CNF format of the SAT problem resembles the canonical POS form of a Boolean function. Prior research [AST + 17, SWW18] has demonstrated that the POS simplification approach may be used to simplify the CNF format of the SAT problem.

Quine-McCluskey Algorithm
For a given function, the Quine-McCluskey algorithm is used to calculate the minimum SOP or POS. Quine [Qui52,Qui55] devised the method, while McCluskey [McC56] expanded upon it. The Quine-McCluskey method, whether dealing with SOP or POS, acts on the vector u in the minterm m u (x) or the maxterm M u (x) and is independent of the operator between variables. In the description that follows, we use the simplification of the SOP form as an illustration and remark that the simplification of the POS form may be accomplished in a similar manner. The Quine-McCluskey algorithm consists of two phases.
Phase 1: finding all prime implicants. We use the vector u to represent the minterm m u (x) in this subsection for simplicity. In this phase, we attempt to merge two or more minterms. If two minterms only vary by one bit, that bit may be substituted with an asterisk ( * ) to indicate that it is arbitrary. This allows two minterms to be represented by a single string. As an example, two 4-bit minterms m 0x1 and m 0x5 may be combined to get 0 * 01, which is referred to as an implicant of size two since it represents two minterms with no constraint on the second bit. Two implicants of size two may be merged further to yield an implicant of size four. In this process, matching asterisks comes first. Then, two implicants that vary by one bit in the remaining positions may be merged.
To locate all prime implicants, it is necessary to compare and combine all feasible pairings of minterms wherever possible. Grouping minterms according to their Hamming weight decreases the number of comparisons. Suppose the minterms of a given function are split into ϱ groups G 0 , G 1 , . . . , G ϱ−1 , and the Hamming weight of minterms in group G i is i for 0 ⩽ i ⩽ ϱ − 1. Consequently, there is no need to compare minterms of nonadjacent groups, since they will always vary in more than one variable. Similarly, comparing minterms within one group is unnecessary as each minterm has the same Hamming weight and any two minterms vary in at least two variables. Therefore, only the minterms of neighbouring groups need to be compared.
A checkmark (✓) is placed next to a minterm after it has been coupled with another minterm to signify that it has been included into an implicant. Thereafter, by merging implicants of size two, implicants of size four are created. Once again, a checkmark is put next to an implicant that may be paired with another implicant. In general, this procedure should be repeated until no more implicant combinations can be made. After that, the unchecked minterms and implicants are known as prime implicants. Since each minterm is contained in at least one prime implicant, the original function equals the disjunction of its prime implicants. Despite the fact that the canonical SOP has been reduced to some degree, the outcome is often not the minimum SOP. To do this, we should complete the second phase. Example 1 provides a tiny illustration of this phase.
Phase 2: determining the minimum SOP. First, a prime implicant chart is constructed in this phase. Along the side are the newly created prime implicants, and along the top are the minterms of the supplied function. If the prime implicant in the row covers the minterm in the column, the symbol ( ) is inserted into the table cell. Then, we search for columns containing a single symbol ( ). If a column contains just one symbol, the minterm can only be covered by a single prime implicant, which is termed as an essential prime implicant. The minimum SOP expression must include all essential prime implicants. In certain instances, the essential prime implicants do not cover all minterms, necessitating the use of extra chart reduction processes. Petrick's method [Pet56] should be a more methodical approach, as opposed to trial and error. Using Petrick's approach, all minimum SOP solutions are derived from a prime implicant chart. This methodology is described algorithmically in Algorithm 1.
Remark 1. Boura and Coggia [BC20] offered a novel formulation for the first phase of the Quine-McCluskey algorithm. By substituting the second phase of the Quine-McCluskey algorithm with a minimisation problem in MILP, an improved algorithm for creating MILP models for S-boxes was created. Although the MILP optimiser does not always find the optimal solution, the output provided by the optimiser might be satisfactory. See Section 3.3.1 and [BC20] for more information.
We present a tiny example to easily comprehend the Quine-McCluskey algorithm.
Example 1. The target is a 4-bit Boolean function f t , and the support is supp(f t ) = {0x0, 0x1, 0x2, 0x5, 0x6, 0x7, 0x8, 0x9, 0xa, 0xe}. Therefore, the canonical SOP form of f t contains ten minterms. As shown in Table 1, the Hamming weight is used to classify the Algorithm 1 Petrick's method 1: Delete the rows for essential prime implicants and their associated columns with the symbol ( ) to reduce the size of the prime implicant chart. 2: Mark the rows of the reduced prime implicant chart as P 0 , P 1 , P 2 , etc. 3: Create a logical function P consisting of the conjunction of disjunctions. Each disjunction relates to a column, and the i-th disjunction is (P i0 ∨ P i1 ∨ · · · ∨ P i k ), where each P ij represents a row covering the i-th column. 4: Multiply out the function P and use the absorption rule P i ∨ (P i ∧ P j ) = P i to reduce it to a minimum SOP. Each term in the reduced expression corresponds to a set of prime implicants including all the minterms in the reduced prime implicant chart. 5: Find each term of the reduced expression of P that has the smallest number of prime implicants. For each term containing a minimum number of prime implicants, count the number of variables in each prime implicant and calculate the overall number of variables. 6: Choose the term or terms that include the minimum total variables and write down the conjunction of prime implicants that corresponds.
minterms into four groups G 0 , G 1 , G 2 , and G 3 . To generate implicants of size two, we begin with group G 0 and compare all of its minterms to those of group G 1 . m 0x0 and m 0x1 may be joined to generate 000 * by replacing the final bit with an asterisk. A checkmark (✓) is then put next to the two minterms to indicate that they have been included into an implicant. m 0x0 and m 0x2 are then joined to make 00 * 0, and m 0x0 and m 0x8 form * 000.
After repeating this process for (G 1 , G 2 ) and (G 2 , G 3 ), we get a list of implicants with a size of two. Subsequently, implicants of size four are derived by merging implicants of size two, and the result can be found in Table 1. Again, a checkmark (✓) is put next to an implicant when it may be paired with another implicant. In this example, no more combinations of implicants of size four are possible. Note that three implicants with a size of four are duplicates generated by merging the same set of four minterms in different orders. The following analysis will not include duplicates. From Table 1, we observe that f t has six prime implicants: 0 * 01, 01 * 1, 011 * , * 00 * , * 0 * 0, and * * 10.  With the prime implicants, the prime implicant chart of f t can be constructed (cf. Table 2). Observe that the columns m 0x9 and m 0xe only contain a single symbol ( ).
* 00 * and * * 10 are hence two essential prime implicants. Table 3 displays the reduced prime implicant chart for f t , which is derived by deleting the two essential prime implicant rows and the corresponding columns with the symbol ( ). The three rows in Table 3 are designated P 0 , P 1 , and P 2 respectively. According to Algorithm 1, the logical function P is generated as Consequently, P 1 (01 * 1), which corresponds to the expression x 0 ∧ x 1 ∧ x 3 , is the term of P having the fewest prime implicants. Together with the two essential prime implicants * 00 * and * * 10, the minimum SOP for f t should be

ESPRESSO Logic Minimizer
The Quine-McCluskey method makes the creation of all prime implicants (Phase 1) more efficient. The number of prime implicants of an n-bit logic function, however, may be proven to equal O(3 n /n). In addition, the second phase attempts to solve a minimum covering problem that is known to be NP-complete. Since the number of components in the covering issue may be proportional to the exponential of the number of input variables of the logic function, the Quine-McCluskey approach becomes unworkable for even moderately large situations. For instance, Abdelkhalek et al.
[AST + 17] failed to implement the Quine-McCluskey method to reduce a 16-bit function pertaining to AES [DR02]. Various heuristic approaches to the issue were prompted by the need to simplify large-scale functions. These heuristic procedures use two routes. One method utilises the structure of the traditional method by first producing all prime implicants. However, rather than creating a minimum cover, a near minimum cover is heuristically determined. This approach still has the potential to produce a very large quantity of prime implicants. The second method attempts to simultaneously detect and choose implicants (who are not necessarily prime) for the cover. A lot of follow-up work centered on the second approach, and the first and most successful was the IBM programme MINI [HCO74]. Later, Brown [Bro81] proposed a method called PRESTO. The core of MINI and PRESTO is expanding each implicant of the supplied function and eliminating implicants that are covered by the expanded implicant. MINI expands each implicant to its maximum size both in the input and output parts of the function. PRESTO expands the input part of each implicant to its maximum size, but then reduces the output parts of the implicants maximally by removing the implicants from as many output spaces as possible. In both methods, implicants are expanded and covered, and thus handling the problem of creation and storage of all minterms.
Nonetheless, this primary technique results in a local minimum that may be far off from the global minimum. In order to prevent this issue, following the first expansion and removal of covered implicants, MINI reduces the remaining implicants to their smallest size while preserving the cover of the function. The implicants are then assessed in pairs for reshaping, with one implicant being enlarged and the other being reduced by the same set of minterms. The expansion process is then resumed, and the whole method is repeated until there is no further decrease in the number of product terms. This metric enables the investigation of wider portions of the optimisation space, with the goal of achieving better results at the cost of increased computation time. PRESTO differs from MINI in the way the expansion process is carried out. MINI generates the complement of the logic function to check is the expansion of an implicant does not change the cover of the function. PRESTO avoids the initial cost of computing the complement, but then the input expansion process requires a check on whether all the minterms covered by the expanded implicant are covered by someother implicant of the cover, which, in general, costs more computation time.
Comparing the different tactics adopted by MINI and PRESTO in a specific context inspired the development of ESPRESSO. The first version ESPRESSO-I was created by Brayton

MILP Modelling Progress for Large S-boxes
In this part, we review the progress made in the MILP modelling of large S-boxes.

First Bit-Oriented Model for S-boxes
Sun et al.
[SHW + 14b] proposed the first bit-oriented model in the search with MILP approach. Each bit of the inner state should have a binary variable allocated to it. The simplest model for an S-box is to describe its * -DDT, which is a reduced form of the DDT in which all non-zero entries are substituted with 1. In this model, we are solely concerned with the feasibility of differential propagation, and a propagation is deemed feasible if its corresponding item in * -DDT equals 1. With this notation, two sets F = {x∥y | x → y is a possible propagation, x, y ∈ F n 2 } , I = {x∥y | x → y is an impossible propagation, x, y ∈ F n 2 } (2) can be constructed for an n-bit S-box. * -DDT can be modelled by employing linear inequalities to either characterise all vectors in F or exclude all vectors in I. These two ideas correspond to the two approaches proposed in [SHW + 14b].

H-representation of the Convex Hull
Note that F is a subset of R 2n , where R is the real number field, and that its convex hull is defined as the smallest convex set that contains F. A convex hull can be represented as the common solutions of a finite number of linear inequalities; this is known as the H-representation of a convex hull. Using the inequality_generator() function in the sage.geometry.polyhedron class of the open-source mathematical software SageMath 2 [The22], it is possible to derive the convex hull. Given that the number of inequalities in the H-representation of a convex hull is typically quite high, Sun et al. developed a greedy technique to reduce the number of inequalities. The greedy algorithm builds a collection of valid cutting-off inequalities by picking at each step an inequality from the convex hull that maximises the number of impossible differential propagations removed from the current feasible region.

Logical Condition Modelling
The main idea of logical condition modelling is that CNF/POS formulas are directly expressible as inequalities in MILP. This method may be thought of as focusing on the set I and generating linear inequalities to eliminate impossible differential propagations one by one. Specifically, for every vector a∥b = a 0 ∥ · · · ∥a n−1 ∥b 0 ∥ · · · ∥b n−1 in the set I, an inequality about 2n Boolean variables x∥y = x 0 ∥ · · · ∥x n−1 ∥y 0 ∥ · · · ∥y n−1 is generated Only when x∥y = a∥b does the inequality become invalid. Applying this method directly will yield an excessive amount of inequalities. To overcome this issue, Sun et al. found that certain S-boxes contain conditional differential features, allowing us to rule out many impossible propagations with a single inequality. If the input difference of the 4-bit S-box of PRESENT [BKL + 07] is 0x9, for example, the least significant bit of the output difference must be 0. Then, eight impossible propagations 1001 ↛ * * * 1 are eliminated by the linear inequality −x 0 + x 1 + x 2 − x 3 − y 3 + 2 ⩾ 0. However, such conditional differential properties do not present in all circumstances. When transitioning from * -DDT to DDT in a model, additional variables must be included to encode differential probabilities. [SHW + 14a] provides an example of how to describe the differential probability of the S-box of PRESENT. Observe that all applications in [SHW + 14b, SHW + 14a] are ciphers with at most 6-bit S-boxes. [SHW + 14a] also stated that the MILP model is mostly applicable to lightweight ciphers. In [SGL + 17], the authors claim unequivocally that the MILP method is incapable of searching for true differential characteristics of 8-bit S-box ciphers.

Modelling for Large S-boxes
In response to the limitations of previous MILP models for large S-boxes, Abdelkhalek et al.
[AST + 17] suggested a novel modelling method for large S-boxes oriented to * -DDT and differential probabilities. The new method is founded on logical condition modelling introduced in Section 3.1.2, in which the number of linear inequalities is minimised or kept as small as possible.

Modelling * -DDT of Large S-boxes
The technique commences with a description of * -DDT. A 2n-bit Boolean function f is built using the * -DDT, and f (x∥y) equals one if and only if x∥y ∈ F = F 2n 2 \ I. The canonical POS form of f is Observe that f (x∥y) is non-zero only when for all a∥b ∈ I, the following expression about x∥y holds (4) is identical to the expression in Equation (3). In contrast to the technique employed in [SHW + 14b] to reduce the number of inequalities by monitoring conditional differential characteristics, Abdelkhalek et al. asserted that lowering the number of inequalities is akin to simplifying the POS form of f , a well-studied subject. The Quine-McCluskey algorithm can be used to determine the minimum POS, which ensures the minimum number of linear inequalities under logical condition modelling. Since the algorithm becomes inapplicable for even moderately large problems, Abdelkhalek et al. were unable to apply it for reducing a 16-bit function corresponding to the * -DDT of AES. In this instance, the heuristic algorithm ESPRESSO was suggested for simplifying large-scale functions so that the number of linear inequalities is maintained as small as possible.

Modelling DDT of Large S-boxes
Abdelkhalek et al. split the entries of a DDT into numerous pieces in order to encode the differential probability and control the scale of the Boolean function, from which the concept of p-DDT is generated.
. For a given S-box and its DDT, the corresponding entry of the p-DDT is one if the probability of entry in the DDT is p; otherwise, it is zero.
By applying the approach to model * -DDT, one set of linear inequalities L p is obtained for each p-DDT, where and ℓ p represents the number of linear inequalities for describing the p-DDT. In addition, a binary variable Q p is allocated to each p-DDT to determine if the linear inequalities in L p should be included in the solution phase. Specifically, the linear inequalities in L p are only applicable when Q p = 1, and they are disregarded otherwise. This technique is realised by adding a suitably large integer 3 M and transforming each linear inequality in L p into a conditional inequality as When Q p = 0, x∥y can take any value; otherwise, the conditional constraint reverts to the original linear inequality. To ensure that the automatic searching model imposes no more than one set of linear inequalities for a certain p-DDT, no more than one Q p may take on a non-zero value. Consequently, an additional binary variable Q is allocated to each S-box, and Q = 1 if and only if the S-box is active. The constraint Q = p Q p establishes the restriction for Q p . Consequently, the base-2 logarithm of the differential probability for a single S-box is calculated to be

Efficient Modelling for Large S-boxes
Despite the fact that Abdelkhalek et al. enabled the description of large S-boxes, the model for the * -DDT of the S-box in AES consists of 8302 linear inequalities, which is sometimes too cumbersome for practical usage. Boura and Coggia [BC20] offered two methods for developing efficient algorithms for modelling 8-bit S-boxes with the fewest number of linear inequalities. These approaches are conceptually related to logical condition modelling (cf. Section 3.1.2), in which an inequality is generated to eliminate an impossible vector. The contrast lies in the fact that the new approaches aim to remove more impossible vectors, which often possess algebraic structures, using a single inequality.

Generalised Logical Condition Modelling
The first method derives simple inequalities to eliminate vector spaces of the type a ⊕ Prec(u) using the following proposition.

Proposition 1 ([BC20]). For two vectors a and u in
Algorithm 2 in [BC20] determines, given the set of impossible propagations I (cf. Equation (2)), all subsets of the type a ⊕ Prec(u) that are not subsets of others. For each a ∈ I, spaces a ⊕ Prec(u) ⊂ I are generated by incrementally increasing the weight of u, where u should satisfy a ⊕ u ∈ I and supp(a) ∩ supp(u) = ∅, and determining, for all v ⪯ u and wt(v) = wt(u) − 1, if a ⊕ Prec(v) has previously been designated as a subset of I. This algorithm, according to Boura and Coggia, is closely connected to the first step of the Quine-McCluskey algorithm. By applying this algorithm, an initial set of inequalities is constructed, with a total of 70336 for the * -DDT of AES. To select a representative set from the initial inequalities, Boura and Coggia formulated a minimisation problem in MILP according to the method described in [ST17] (cf. Section 3.4.1) and solved it using the Groubi optimiser [Gur22]. Although the optimiser did not achieve the minimum for the S-box of AES even after more than twenty days, the solution supplied by the optimiser is much superior than that returned by ESPRESSO in [AST + 17]. The number of inequities falls from 8302 to 7461.

Modelling S-boxes with Inequalities Issued from Balls
The second method creates inequalities to exclude points from a ball B(r, c) whose definition is shown below.
Definition 2 ([BC20]). A ball of F m 2 of radius r centred at c ∈ F m 2 is the set of all points whose Hamming distance from the center c is no more than r, i.e., B(r, With the following proposition, all points in the ball B(r, c) may be eliminated.
Since the set of possible propagations F (cf. Equation (2)) is not sparse for S-boxes with a low differential uniformity, the number of balls consisting purely of impossible propagations is normally rather modest, therefore deleting complete balls is frequently ineffective. In this instance, it is feasible to create distorted balls in which any possible propagations at the sphere are eliminated. The following proposition provides instructions for making distorted balls.
2 be a ball of radius r from which we remove the set of points Q = (c ⊕ Prec(q)) ∩ S(r, c) for some q ∈ F m 2 . The vectors in Q represent possible propagations towards the sphere of the ball. We define a ∈ Q m , where Q is the rational number field, so that a i = r+1 r if q i = 1 and a i = 1 otherwise. Then, Boura and Coggia proposed a way to remove all points from three distorted balls with a radius of one using a single inequality in order to obtain better candidates for the initial set of inequalities; this innovative notion is used in Algorithm 3 of [BC20]. After constructing an initial set of inequalities with Algorithm 3, a minimisation problem in MILP is created which is then resolved by the Groubi optimiser. The number of inequalities reaches 2882, despite the fact that the optimisation for the S-box of AES did not terminate after a few days. In terms of the amount of inequalities, the proposed solution is considerably superior to the ESPRESSO-simplified version [AST + 17]. Notably, despite the creation of improved MILP models for the S-boxes of AES and SKINNY-128, these models are not employed to analyse the differential characteristics of the two ciphers. The authors used the new models to investigate the resistances of the two ciphers against the impossible differential attack [BBS99].

Modelling by Minimum Number of Inequalities
Whether reducing the number of inequalities for the S-box can reduce the runtime to solve the MILP problem is a debate from the pioneer work [SHW + 14b, SHW + 14a]. As a result, modelling S-boxes with the fewest inequalities is a hot topic for research in the area of automatic search. Numerous studies have been conducted to investigate this matter.

Reduction Algorithm for Small S-boxes Based on MILP
The research of Sasaki and Todo [ST17] is able to overcome the limitations of greedy algorithm in terms of locating the fewest inequalities. This technique is derived from the Hrepresentation supplied by SageMath. Assume that the H-representation comprises ℏ linear inequalities designated L 0 , L 1 , . . ., L ℏ−1 . Selecting inequalities from the H-representation is transformed into a MILP problem by the reduction algorithm.
For each inequality L i , the reduction method introduces one binary variable ℓ i , i = 0, 1, . . . ℏ − 1. ℓ i = 1 indicates that the inequality L i is employed in the reduced system, whereas ℓ i = 0 implies that L i is not used. Then, as in Equation (2), the set I of impossible propagation is constructed. Assuming I consists of Λ vectors labelled I 0 , I 1 , . . ., I Λ−1 . For each vector I j in I, the inequalities that can exclude I j are computed, and the indices of the inequalities that satisfy this requirement are arranged in the set Index j ⊂ {0, 1, . . . , ℏ − 1}. To guarantee that each impossible propagation is eliminated by at least one inequality, the following Λ constraints Sasaki and Todo demonstrated empirically, via the reduction algorithm, that minimising the number of inequalities does not reduce the execution time of the MILP optimiser.

Finding Smallest MILP Models
Although minimising the number of inequalities does not always reduce the runtime of the MILP optimiser, many subsequent studies continue to seek a breakthrough on the number of inequalities. The task of discovering the smallest MILP models might be of theoretical importance in its own right. For further detail, interested readers might see [LS22,Udo21].

SAT/SMT Modelling Progress for S-boxes
Mouha and Preneel [MP13] introduced the use of SAT/SMT in symmetric-key cryptography for the search of distinguisher. Initially, this technique was used to ARX ciphers [KLT15,LWR16,SHY16]. Two studies utilising the SAT/SMT approach to search for differential characteristics appeared around the same time [AK18,SWW18]. [SWW21] is also associated with the SAT-based search for differential characteristics. However, the focus of the study was on modifying the objective function to achieve acceleration for the SAT method, and the modelling for S-boxes adheres to the method described in [SWW18].

Logical Condition Modelling
The S-box modelling approach described by Ankele and Kölbl [AK18] is comparable to the logical condition modelling using MILP method (cf. Section 3.1.2). The clause is built for each vector a∥b = a 0 ∥ · · · ∥a n−1 ∥b 0 ∥ · · · ∥b n−1 in I and then added to the SAT problem. This expression is invalid if and only if x∥y = a∥b. Adding all clauses thereby removes all impossible propagations. Although the authors claim that their approach is compatible with 8-bit S-boxes, all of the ciphers analysed in [AK18] employ 4-bit S-boxes.

Logical Condition Modelling Using Simplification
Sun et al. [SWW18] suggested a modelling technique with simplification, which was inspired by the work of Abdelkhalek et al. [AST + 17] in MILP approach. We first review the model oriented to active S-boxes.

Modelling Oriented to Active S-boxes
The model for linear active S-boxes is a natural extension of the model for differential active S-boxes, which we apply as an example. Apart from 2n Boolean variables x = (x 0 , x 1 , . . . , x n−1 ) and y = (y 0 , y 1 , . . . , y n−1 ) to represent the input and output differences of an n-bit S-box, there should be an auxiliary Boolean variable w. The value of w for active S-boxes is one and zero for inactive S-boxes under the premise that x → y is a possible propagation. Limiting the total number of active S-boxes for the differential characteristic is equal to setting a maximum for the total sum of auxiliary variables for all S-boxes in the characteristic. Pr(x → y) represents the differential probability of propagation x → y. According to this criteria, possible values for x∥y∥w fall into the set where the three components in the subscript of F are, in order, the number of input bits of the S-box, the number of output bits of the S-box, and the number of auxiliary variables involved. The direct method removes impossible vectors one by one with clauses The number of clauses in the resulting SAT problem will be extremely high if the scale of the set F 2n+1 2 \F ⟨n,n,1⟩ is enormous, which could negatively affect the performance of the SAT solver. Finding a solution to reduce the clauses is thus necessary.

Modelling Oriented to Differential Probabilities and Linear Correlations
We only introduce modelling for differential probabilities because the methods for constructing SAT models for linear correlations and differential probabilities are comparable. Note that additional auxiliary variables are required to encode probability information.
Given an n-bit S-box, we first construct a set P consisting of all probabilities for possible differential propagation x → y. The weight of a possible propagation with a probability of p is − log 2 (p) and can take on non-integer values. For non-integral weights, the set D = {⌈log 2 (p)⌉ − log 2 (p) | log 2 (p) is not an integer, and p ∈ P } contains all feasible decimal values. Assume there are ν items in the set D, which are denoted by d 0 , d 1 , . . ., d ν−1 . The plan is to introduce two groups of auxiliary variables that are used to encode the integral and decimal components of the weight, respectively. To be precise, max{⌊− log 2 (p)⌋ | p ∈ P } ≜ µ variables (u 0 , u 1 , . . . , u µ−1 ) ≜ u are introduced for the integral portion and ν variables (v 0 , v 1 , . . . , v ν−1 ) ≜ v are introduced for the decimal portion.
Similarly to the scenario concerning active S-boxes, we subsequently determine the possible values for the vector x∥y∥u∥v. If x → y is a possible propagation with a probability of p, then x∥y∥u∥v assignments must fulfil the equation where e l is a unit vector in F ν 2 and the l-th bit is the only non-zero bit. Again, F ⟨n,n,µ+ν⟩ can be considered as the support of the function f ⟨n,n,µ+ν⟩ (x∥y∥u∥v) = 1, if x∥y∥u∥v ∈ F ⟨n,n,µ+ν⟩ 0, otherwise .
The SAT model for differential probabilities can be derived by simplifying the canonical POS form of f ⟨n,n,µ+ν⟩ .
Example 2. Consider a 4-bit S-box whose DDT contains the values 0, 2, 4, 6, and 16. The corresponding set of probability for all feasible differential propagations is {2 −3 , 2 −2 , 2 −1.415 , 1}. Thus, we claim three variables (u 0 , u 1 , u 2 ) ≜ u to represent the integral portion of the weight and one variable v 0 to represent the decimal portion of the weight. The optional set of possible values for x∥y∥u∥v 0 is The vector in the set confirms the equation u 0 + u 1 + u 2 + 0.415 · v 0 = − log 2 (p) for a possible propagation x → y with a probability of p.

SMT Modelling for Large S-boxes
Note that the modelling techniques described in [AK18,SWW18] mostly apply to ciphers with 4-bit S-boxes. Liu et al. [LLL + 19] were the first to introduce the SMT modelling technique for 8-bit S-boxes. The effectiveness of SMT modelling is dependent on the reasoning statement supplied by the SMT solver STP [GD07].
In addition to two n-bit vectors x and y representing the input and output differences of the S-box, the SMT model requires a 1-bit vector w for the search of valid differential characteristics. The reasoning statement is added to the SMT model for each vector a∥b in F where a → b represents a possible propagation. This phrase specifies that w is assigned the value 1 if x∥y is a member of the set F. Then, for each vector a∥b in I, the reasoning statement is added to the SMT model, resulting in w = 0 if x∥y takes values from the set I. The SMT model also asserts that w = 1 so that any impossible propagations for the S-box are ignored by the solver.
SMT modelling for S-boxes oriented to differential probabilities and linear correlations is comparable to the previously described fundamental model. The SMT model employs a direct way to eliminate all impossible propagations and does not involve simplification.

Fast SAT Models for Large S-boxes
As the reduction of canonical POS is NP-complete, it is typically inefficient to generate reduced sets of clauses for big S-boxes (e.g., 8-bit S-boxes) oriented to differential probabilities and linear correlations. Thus, there has been no SAT model 4 for large S-boxes till now. In this section, three fast SAT modellings for large S-boxes oriented to differential probabilities and linear correlations are proposed. Although these approaches cannot guarantee a minimum number of clauses, the time required to acquire clauses is significantly reduced.
The 8-bit S-box S 8 (cf. Table 5) of SKINNY-128 [BJK + 16] will serve as an illustration for the presentation of the subsequent three methods. The set consists of the probabilities of possible differential propagations x → y for S 8 . Following the method in [SWW18], seven variables (u 0 , u 1 , . . . , u 6 ) ≜ u are then introduced to describe the integral portion of the weight. To encode the decimal portion of the weight, three variables (v 0 , v 1 , v 2 ) ≜ v are used, with v 0 = 1 (resp., v 1 = 1, v 2 = 1) if and only if the decimal portion of the weight equals 0.678 (resp., 0.415, 0.193). The option for the set of possible values for x∥y∥u∥v is The vector in the set verifies the equation for a propagation x → y that is possible with probability p. F ⟨8,8,10⟩ can be viewed as the support for a 26-bit Boolean function f ⟨8,8,10⟩ . We attempted to use ESPRESSO to identify a simplification for the canonical POS form of f ⟨8,8,10⟩ , but after more than a hundred days, nothing was returned. Therefore, we wish to discover an alternate method for constructing a SAT model for S 8 that demonstrates differential probability. All of the experiments in this study are done on a 16-core AMD EPYC™ 7302 processor , of which only one core is utilised. Although the subsequent approaches are demonstrated using models for differential probabilities, the creation of SAT models for linear correlation may be accomplished in a similar fashion.

Trade-off Between Level of Simplification and Time
The initial attempt is centred on ESPRESSO itself. We notice that ESPRESSO provides multiple options and commands for minimisation and creates a trade-off between simplification level and execution time. Listed below are several options that may reduce the runtime.

-efast This option stops ESPRESSO after the first [Expand] and [Irredundant Cover]
operations (cf. Algorithm 2) and prevents it from iterating over the solution.
-eness With this setting, essential prime implicants will not be identified.
-enirr With this option, the result will not necessarily be made irredundant in the last step which removes redundant literals.
-eonset This option recalculates the support of the input function prior to minimisation, which is advantageous when the canonical POS form includes a large number of maxterms.
Due to the heuristic nature of these options, the number of clauses in the result may not be the minimal. However, in the vast majority of applications, the number of clauses is acceptable for modern SAT solvers. ESPRESSO is used to conduct the simplification of f ⟨8,8,10⟩ with the four newly introduced options. The simplification is only accomplished with the option -eonset. There are 820 clauses in the output, and the total execution time is 3521.42 seconds. Although the remaining three options claim that certain procedures are skipped and the simplification may be faster, the programmes do not terminate after fifty days. Consequently, the first fast SAT model for large S-boxes leverages the -eonset option of ESPRESSO. The applicability of this strategy is also examined for the S-box of PIPO [KJK + 20], which is a simplification issue for a 25-bit Boolean function. Again, only the -eonset option delivers the output after 3769.86 seconds, and the simplified result has 6066 clauses.

Two-Step Encoding Method
Noting that the complexity of simplification increases exponentially with the number of input variables, we wondered whether the size of the function may be decreased to alleviate the difficulty of developing the SAT model. Consequently, the encoding mechanism for auxiliary variables must be modified.
Initially, we examine the set F ⟨8,8,10⟩ and notice that u 5 and u 6 always have the same value. Therefore, it is logical to reduce the number of auxiliary variables for the integral portion to six, i.e., u ∈ F 6 2 , and reconstruct the set of possible values for x∥y∥u∥v as As a result, the evaluation of the weight for a feasible propagation is modified using , the time required to achieve the solution is decreased although the number of clauses is almost the same. Then, we evaluate the feasibility of reducing the size of the function received by ESPRESSO further. The main idea is dividing the encoding phase for an n-bit S-box into two steps. In addition to the auxiliary variables u and v introduced regarding the integral and decimal parts of the weight, we claim a set of chaining variables z.
In the first step, we map the various probabilities to distinct values of the vector z. We recall that P represents the set containing all probabilities for possible differential propagations and denote the elements in P as {p 0 , p 1 , . . . , p |P|−1 }. Thus, the number of required chaining variables is ⌈log 2 (|P|)⌉ ≜ ζ. To be specific, we introduce ζ variables (z 0 , z 1 , . . . , z ζ−1 ) = z. An option for the set of possible values for the vector x∥y∥z is In the second step, we create a connection between chaining variables z and auxiliary variables u∥v so that d i ·v i can correctly compute the weight of the propagation.
An option for the set of possible values for the vector z∥u∥v is where the two components in the subscript of F are the number of chaining variables and the number of auxiliary variables for the integral and decimal parts of the weight.
Thus, with these two steps, we convert the simplification of a large-scale function into the simplification of two relatively small-scale functions f (1) ⟨n,n,ζ⟩ and f (2) ⟨ζ,µ+ν⟩ , whose supports are F (1) ⟨n,n,ζ⟩ and F (2) ⟨ζ,µ+ν⟩ , respectively. We return to the S-box S 8 of SKINNY-128. In first step, we introduce four chaining variables (z 0 , z 1 , z 2 , z 3 ) ≜ z and create an option for the set of possible values for x∥y∥z as In the second step, we create the connection between z and u∥v. An option for the set of possible values for z∥u∥v is (cf. Algorithm 2) , which is more expensive but sometimes yields superior results.
-Dexact This command executes the precise minimisation procedure, which guarantees the smallest possible amount of clauses and heuristically minimises the number of literals. It might potentially be costly.
We apply different options in ESPRESSO to simplify f (1)

Simplifying by Partitioning Method
In the test regarding ESPRESSO, we find that the simplification of a large-scale function is not difficult if the number of maxterms in the function is not excessively huge. This observation is the basis for the third method.
Definition 3. Given a set X , a family of sets Ψ is a partition of X if and only if the following conditions are met.
• The family Ψ does not contain the empty set.
• X is equal to the union of the sets contained in Ψ .
• In Ψ, the intersection of any two different sets is empty set.
Assume that Ψ = {ψ 0 , ψ 1 , . . . , ψ χ−1 } is a partition of the set supp(f ), where supp(f ) is the support of the n-bit function f in Equation (1). For each set ψ i in Ψ, we define a function f i with ψ i as its support. With this definition, f can be expressed as    . This process is done until the clauses of two functionsf | 1 0 andf | 1 1 are merged into a single function f | 0 ⊘ . By using ESPRESSO to simplify the function f | 0 ⊘ , f | 0 ⊘ is generated, which is a simplification of the original function f . It can be imaged that the level of simplification of the final output and the runtime are affected by the number of components in the initial partition Ψ n ⟨s⟩ f . To discover the feasibility of the iterative approach, we implement tests for the 26-bit Boolean function f ⟨8,8,10⟩ regarding S 8 . In the test, the value of s is traversed from 1 to 15, and the outcome of simplification using the -eonset option is depicted in Figure 2(a). From Figure 2(a), we find that the runtime shows a decrease followed by an increase with the increasing of s. When s is set to 10, the time required to obtain the simplification is 1076s, which is significantly shorter than the time required without the partition approach. Additionally, some values of s make it possible to achieve a simplification with fewer clauses. For instance, when s is set to 2, 4, 5, 12, or 13, 807 clauses are returned. Since the number of maxterms inf | s x1 is significantly fewer than that of f , the -estrong option can be used to simplify f ⟨8,8,10⟩ for certain values of s. As illustrated in Figure 2(b), it is possible to reduce the number of clauses even if the runtime may be exceptionally long. Notably, when s is set to 12, the result of employing the -estrong option is reduced to 792 clauses. To illustrate that simplification by partition is a widely used technique, it is also applied to the 25-bit Boolean function f ⟨8,8,9⟩ and the 20-bit Boolean function f (1) ⟨8,8,4⟩ with respect to SKINNY-128. As seen in Figure 3, the partition approach reduces the time required to get simplified clauses with the -eonset option, and when s is set to 12, the number of clauses reaches 805, which is fewer than the number of clauses that would be generated without partitioning. In addition, when s ⩾ 4, the -estrong option can be used, and the runtime is less than 1883.51 seconds when s is 11 or 12. In general, the number of clauses returned by the -estrong option is always smaller than the number returned by the -eonset option. The number of clauses reaches 792 when s equals 12. Figure 4 illustrates the outcome of simplifying f (1) ⟨8,8,4⟩ . Using the partition technique for simplification with the -eonset option might marginally reduce the time required to obtain the simplified clauses, as the runtime without partitioning is already quite quick. A benefit is that the number of clauses can be somewhat decreased. For simplifications utilising the -estrong option, the time required to achieve a simplified result is drastically decreased when compared to the simplification that do not employ the partition approach. When s is set to 7, the execution time is lowered from 101051 seconds to 168 seconds, and the number of clauses is decreased from 730 to 727. These examples demonstrate that the partition approach is effective not just for simplifying large-scale functions, but also for simplifying small-scale functions with the -estrong option.

Tight Probability Bound for 14 Rounds of SKINNY-128
In this section, we will first review the specifications of SKINNY-128 and then present the probability bound for 14-round encryption. In the latter portion of this section, we take SKINNY-128 as an example to analyse the runtime for solving SAT problems constructed using different encoding methods in Section 5.

Specification of SKINNY-128
SKINNY-128 is one variant of the SKINNY family of lightweight block ciphers [BJK + 16]. The block size is n = 128 bits, and during the encryption phase, the internal state is viewed as a 4 × 4 square array of bytes. SKINNY-128 is based on the TWEAKEY framework in [JNP14] and accepts tweakey input. The user can select which portion of the tweakey input serves as the key and/or tweak material. There are three variations of SKINNY-128 based on tweakey size t, with tweakey sizes of 128, 256, and 384, respectively. The tweakey state is seen as a collection of t/n 4 × 4 square arrays of bytes.
Initialisation After receiving a plaintext m = m 0 ∥m 1 ∥ · · · ∥m 15 , where m i ∈ F 8 2 , 0 ⩽ i ⩽ 15, the internal state is created by setting the internal state IS to As illustrated in Figure 5, a single encryption round consists of the following five operations: SubCells (SC), AddConstants (AC), AddRoundTweakey (ART), ShiftRows (SR), and MixColumns (MC). The number of rounds r varies according to the block and tweakey sizes. The minimal number of rounds required for SKINNY-128 is 40 rounds.

SubCells (SC)
The 8-bit S-box S 8 (cf. Table 5) is applied to each byte of the state.
The six bits of the LFSR are initialised to zero, and the LFSR is updated before being used in the current round. The bits of the LFSR are organised in a 4 × 4 array where c 0 = 0∥0∥0∥0∥rc 3 ∥rc 2 ∥rc 1 ∥rc 0 , c 1 = 0∥0∥0∥0∥0∥0∥rc 5 ∥rc 4 , and c 2 = 0x02. The round constants are mixed with the internal state via bitwise XOR operations, while array placement is respected.
AddRoundTweakey (ART) The first and second rows of each tweakey array are extracted and XORed with the internal state. Since the value of tweakey has no effect on the search for differential characteristics, we do not describe its construction.

ShiftRows (SR)
The second, third, and fourth rows are rotated to the right by one, two, and three bytes, respectively.

Tight Probability Bound for 14-Round of SKINNY-128
Due to the large size of the state, the designers of SKINNY-128 gave only lower bounds for the number of differential active S-boxes. Given that the number of active S-boxes for 14 rounds is 61, the probability for 14-round differential characteristics has an upper bound of 2 −122 . Abdelkhalek et al. attempted to propose tight upper bounds for the probability of SKINNY-128 utilising a MILP model for large S-boxes. The task was completed up to 13 rounds, and the search on 13 rounds took 16 days, according to [AST + 17]. For 14-round of SKINNY-128, they merely demonstrated that no differential characteristic had a probability greater than 2 −128 . Using the fast SAT models for big S-boxes, we find in this paper that the upper bound on the probability for 14 rounds of SKINNY-128 is 2 −131 , thereby finishing the remaining task of Abdelkhalek et al. The outcome of searching for the upper bound on the differential probability for up to 14 rounds of SKINNY-128 using the SAT approach is shown in Table 6. The probability bounds for 10-round and 13-round differential characteristics conform to the result in [AST + 17]. Figure 6 illustrates a 14-round differential characteristic of SKINNY-128 with probability 2 −131 .

Runtime with Different Encoding Methods for S-boxes
The time required to solve SAT problems generated with different S-box encoding techniques varies. The numerous sets of clauses in Section 5 produced by ESPRESSO under various settings and encoding techniques are utilised to determine the upper bound of differential probability for SKINNY-128 from 1 to 14 rounds. The solver Cryptominisat is used to solve all of the SAT problems in this study. In the test, we use the models returned by ESPRESSO without modifying the order of the clauses in the model. We should remind the reader that the order of clauses and/or variables may also influence the execution time of the SAT solver. In addition, the SAT solver may also have a considerable variance in the execution time, and we believe that a comparison with other SAT solvers would be an interesting research project.
As seen in Figure 2, for the 26-bit encoding of S 8 , the partition approach can minimise the time required for ESPRESSO to create reduced clauses, but the time required for the SAT solver to solve the SAT problem rises. The good news is that for some values of s, the overall runtime, which consists of the runtimes for ESPRESSO and SAT solver, can be decreased. For example, when s is set to 8, the overall runtime is 10754.7 seconds, however it is 11123.1 seconds when using the clauses provided in Section 5.1. Figure 3 depicts the test outcomes relevant to the 25-bit encoding of S 8 . Similar to the scenario presented in Figure 2, the partition approach shortens the execution time for ESPRESSO but cannot reduce the execution time for the SAT solver. In addition, a comparison of Figure 2 and Figure 3 reveals that lowering the number of variables may not reduce the execution time of SAT solver. If the partition technique is not employed, for instance, it takes 7601.7 seconds to solve a SAT problem with S-boxes containing 820 clauses and 26 variables; whereas it takes 7952.3 seconds to solve a SAT problem with S-boxes containing 818 clauses and 25 variables. Recall that Sasaki and Todo [ST17] Round 0 -1 0x20 0x20 0x20 0x20 SC 0x20 0x20 0x20 0x20 0x20 0x20 0x20 0x20 0x20 0x20 SR • ART • AC 0x20 0x20 0x20 0x20 0x20 0x20 0x20 0x20 0x20 0x20 0x20 0x20 0x20 0x20 0x20 0x20 0x20 0x20 An active S-box with probability 2 −2 An active S-box with probability 2 −3 A byte with nonzero difference A byte with zero difference Figure 6: 14-round differential characteristic of SKINNY-128 with probability 2 −131 .
shown that, for 4-bit S-boxes, lowering the number of inequalities in MILP models does not necessarily reduce the runtime of the MILP optimiser. Figure 2 and Figure 3 illustrate that, for 8-bit S-boxes, reducing the number of clauses in SAT problems does not always decrease the execution time of SAT solver. Moreover, a comparison of Figure 2 and Figure 3 demonstrates that lowering the number of variables in SAT problems may not reduce the execution time of SAT solver. The sets of clauses (see Table 4) created by the two-step encoding approach are also incorporated into SAT problems so that the runtime for searching for 1-round to 14-round upper bound on differential probability for SKINNY-128 may be evaluated. The test result is given in Figure 7. Figure 7 demonstrates that for two sets of clauses with the same amount of clauses, the runtime for the SAT solver to solve their respective SAT problems may change. When comparing the results of Figure 7 to those of Figure 2 and Figure 3, it is evident that the time required for ESPRESSO to locate simplified clauses has been drastically decreased. If the ESPRESSO option -eonset is fixed, the runtime for the simplification phase of the two-step technique is 99.0 seconds, whereas the runtimes for the 26-bit and 25-bit encodings are 3521.4 and 1883.5 seconds, respectively. The bad news is that the total runtime increases; the minimum total runtime for the two-step technique is 12659.4 seconds, whereas the minimum total runtimes for the 26-bit and 25-bit encodings are 11123.1 seconds and 9835.8 seconds, respectively. Figure 4 depicts the cumulative impact of the two-step technique and partition method on the execution time of SAT solver. As shown in Figure 4(a), incorporating the partition into the two-step technique may have a beneficial impact on the execution time of SAT solver when ESPRESSO is implemented with the -eonset option. For instance, when s is set to 1, the runtime of the SAT solver is 11087.2 seconds, which is shorter than the 12639.4 seconds required by the technique without the partition method. In addition, the overall execution time for s = 1 is 11126.3 seconds, which is faster than the total execution time when the partition technique is not used. Comparing the findings in Figure 4(a) and Figure 4(b) demonstrates the previously indicated fact, namely that lowering the number of clauses in SAT problems does not necessarily reduce the execution time of the SAT solver.
After analysing all the test results, we conclude that the encoding strategy described in Section 5.1, which relies on the -eonset option given by ESPRESSO, is a good solution for balancing the runtime for ESPRESSO and the execution time for the SAT solver. In the remaining applications for PIPO and AES-based constructions, we employ the procedure described in Section 5.1 to create S-box clauses.

Related-Key Differential Properties of PIPO
In this section, the modelling approach for large S-boxes is applied to PIPO, and the related-key differential properties of the two variants of PIPO are studied.

Description of PIPO
PIPO is a lightweight block cipher proposed at ICISC 2020 by Kim et al. [KJK + 20]. It is a 64-bit block cipher with two instances that accept 128-bit and 256-bit keys, and we distinguish between them using the notations PIPO-128 and PIPO-256. The number of rounds r during encryption is dependent on the size of the key: r = 13 for PIPO-128 and r = 17 otherwise.
The internal state is regarded as an 8 × 8 square array of bits during the encryption and decryption processes. We refer to X as the internal state and X[i] as the i-th row of X for 0 ⩽ i ⩽ 7. Plaintext m = m 63 ∥m 62 ∥ · · · ∥m 0 is supplied to the cipher. The internal state is initialised by setting X in the following manner The key schedule is quite simple.
▷ The 128-bit master key K is divided into two 64-bit subkeys K 0 and K 1 as K = K 1 ∥K 0 for PIPO-128. The whitening and round keys are defined as RK i = K i mod 2 , where i = 0, 1, . . . , 13.
Before running the round function in both variants, the whitening key RK 0 is XORed with the internal state. One encryption round in PIPO consists of the following three operations: S-layer, R-layer, and round key and constant XOR. See Figure 8 for a demonstration of the round function.
S-layer An 8-bit S-box denoted as S 8 is applied to each column of the 8 × 8 square array of bits, with the uppermost bit being the least significant. The definition of S 8 is available in [KJK + 20].

R-layer
The rows of the array of cipher state bits are rotated in this layer, as depicted in Figure 8.
Round key and constant XOR Internal state is XORed with the i-th round key RK i and the i-th constant c i = i, where 1 ⩽ i ⩽ r.

Application to AES-Based Constructions
Jean and Nikolić [JN16] suggested a number of AES-based constructions that can be utilised as building blocks for Message Authentication Codes (MACs) and Authenticated Encryption (AE) schemes. They intended to explore the efficiency boundaries of constructions that may be realised with AES-NI instruction aesenc, which executes one round of AES. The internal states of the designs are made up of many 128-bit words, known as blocks, and the step functions are based only on aesenc and bitwise XOR operations. The state size, the amount of aesenc calls each step, and the selection of state words to which aesenc is applied assure the high efficiency of the designs. In addition to the efficiency, Jean and Nikolić were particularly concerned with the security of the designs. Since the most common attacks against MACs and AE are internal collisions based on high probability differential characteristics that begin and end with zero state differences, the security of these constructions is determined by counting the number of active S-boxes required to produce an internal collision. The minimum number of active S-boxes should be 22, given that the key size is 128 bits and the maximum differential probability of the S-box in AES is 2 −6 . Jean and Nikolić used a MILP-based search to determine the number of active S-boxes and provided seven secure constructions C1 -C7 (cf. Figure 11(a) - Figure 11(g)) with excellent state size and efficiency trade-offs. The lower bounds for the number of active S-boxes for C1 -C7 using the MILP model are listed in Table 9, however the number of step functions necessary to attain the lower bound is not specified. As noted by the authors, the automatic search strategy they employ targets truncated differentials. There may not be a true differential characteristic that corresponds to the search outcomes. In other words, it is possible that the security of these structures was underestimated. (g) Step function of C7. Abdelkhalek et al.
[AST + 17] were aware of this issue and analysed the security of two arbitrarily selected constructions, i.e., C5 and C1, using the MILP model for large S-boxes (cf. Table 9). They demonstrated that it is impossible for C5 to have any 4-step truncated differential characteristics with no more than 23 active S-boxes. Therefore, they concluded that the minimum number of active S-boxes for C5 is 24, an increase from the initial estimate of 22 by the designers. For C1, they confirmed that the minimal number of active S-boxes required to cause a collision is 22, while the best differential characteristic has a probability of 2 −134 , as opposed to 2 −132 .
In this study, the encoding approach for large S-boxes proposed in Section 5.1 is applied to the constructions C1 -C7. The SAT method is used to determine the minimum number of active S-boxes necessary to cause an internal collision. All seven constructions, ranging from two to eight steps, are analysed. The test outcomes are shown in Table 9.

Results on C1
For the construction C1, we first verify that no differential characteristic may lead to an internal collision if the number of steps n s is fixed at 2. When n s ⩾ 3, the lower bound on the number of active S-boxes reaches 22, in agreement with the analyses in [JN16] and : There is no differential characteristic with the specified number of step functions.
[AST + 17]. According to our analysis of all 3-step differential characteristics with 22 active S-boxes, there are a total of four differential patterns (cf. Supplementary Material), which match to the four patterns described in [AST + 17].
In addition to examining the lower bound on the number of active S-boxes, Abdelkhalek et al. [AST + 17] evaluate the upper bound on the differential probability for 3-step differential characteristics. The probability upper bound is demonstrated by analysing each of the four differential patterns and determining that not all of the active S-boxes in the characteristic attain the optimal probability, say 2 −6 . We also study the upper bound on the differential probability for C1 using the automatic technique. The SAT modelling approach described in Section 5.1 for large S-boxes oriented to differential probability is implemented. The test results validate the conclusion presented by Abdelkhalek et al., namely that the maximum differential probability for 3-step characteristics that facilitate an internal collision is 2 −134 . Figure 12 demonstrates a newly identified 3-step differential characteristic with a probability of 2 −134 .

Results on C2
For C2, we first determine that there is no differential characteristic that would support an internal collision if the number of steps n s is set at 2. When n s is set to 3, the minimum number of active S-boxes is 50. When n s ⩾ 4, the minimum number of active S-boxes for differential characteristics guaranteeing an internal collision reaches 25, as specified in [JN16]. The differential patterns of all 4-step differential characteristics of 25 active S-boxes are analysed, resulting in a total of 16 differential patterns. The 16 differential patterns are supplied in the Supplementary Material.

Results on C3
For C3, we begin by ensuring that, while the number of steps n s is set at 2, there are no differential characteristics that might cause an internal collision. The minimum number of active S-boxes for differential characteristics is 47 when n s is set to 3 or 4. For n s ⩾ 5, we discover that the minimum number of active S-boxes for differential characteristics is 36, and this bound holds for up to eight steps. Given the 34 active S-boxes established by [JN16], our new result increases the lower bound on the number of active S-boxes for C3. There are 32 differential patterns after analysing the differential patterns of all 5-step differential characteristics with 36 active S-boxes. The 32 differential patterns are presented in the Supplementary Material.

Results on C4
For the C4 construction, we first discover that if the number of steps n s is set at 2, no differential characteristic can sustain an internal collision. When n s is increased to 3, the number of active S-boxes for differential characteristics reaches a minimum of 33. If n s is set to 4, the minimum number of active S-boxes is reduced to 25; this bound is consistent with the one provided by [JN16]. The differential patterns of all 4-step characteristics with 25 active S-boxes are then analysed, and a total of 12 differential patterns are discovered. The 12 differential patterns can be seen in the Supplementary Material. An active S-box with probability 2 −6 An active S-box with probability 2 −7 A byte with nonzero difference An inactive S-box Figure 12: 3-step differential characteristic of probability 2 −134 for C1.

Results on C5
For the construction C5, when the number of step functions n s is fixed at 2, there is no differential characteristic that may cause an internal collision, and an internal collision is only feasible if n s ⩾ 3. When n s is set to 3, there must be at least 40 active S-boxes for differential characteristics to cause a collision. The minimum number of active S-boxes is reduced to 25 when n s = 4. Table 9 demonstrates that our results raise the lower bound on the number of active S-boxes for C5 when compared to those of [AST + 17] and [JN16]. We examine the 4-step differential characteristics of 25 active S-boxes and discover that all of these characteristics can be classified into 12 differential patterns. The 12 differential patterns are available in the Supplementary Material.

Results on C6
First, we verify that 2-step differential characteristics cannot lead to an internal collision for C6. When the number of steps n s is set to 3, the minimum number of active S-boxes is determined to be 48. Due to the extremely lengthy runtime of the SAT solver, when n s is fixed to 4, we only validate that all 4-step differential characteristics have more than 41 active S-boxes, but we cannot guarantee an exact lower bound for the number of active S-boxes. For n s ⩾ 5, the minimum number of active S-boxes is confirmed to be 23; this value corresponds to the one specified in [JN16].

Results on C7
For the C7 construction, we find that 2-step differential characteristics can lead to an internal collision, and that the minimal number of active S-boxes is 48. When n s is set to 3, the minimum number of active S-boxes remains unchanged. For 4-step differential characteristics, we only validate that the minimal number of active S-boxes is more than 37, and the concrete bound cannot be determined until the SAT solver is run for an extremely long period of time. The minimum number of active S-boxes is 28 when n s = 5 or 6. Due to the restricted problem-solving capacity of the SAT solver, the lower bound on the number of active S-boxes for n s ⩾ 7 is uncertain. We only validate that the minimum number of active S-boxes exceeds 24, tying the result with [JN16].

Conclusion
This study begins with a summary of the MILP and SAT modelling developments for large S-boxes. Then, we provide three techniques for quickly constructing SAT models for large S-boxes oriented to differential probabilities and linear correlations. The newly suggested encoding techniques are initially used to the analysis of SKINNY-128. Comparing the runtime of alternative encoding techniques for S-boxes reveals that the first encoding approach achieves a good compromise between the runtime of ESPRESSO and the execution time of the SAT solver. On the other hand, we determine that the upper bound on the differential probability for 14 rounds of SKINNY-128 is 2 −131 , completing the remaining work of Abdelkhalek et al. The first approach of encoding large S-boxes is also implemented on PIPO and the seven AES-based constructions C1 -C7. For two constructions C3 and C5, the current lower bound on the number of active S-boxes is lifted, resulting in a more precise security analysis for these two structures. Additionally, all differential patterns for C1 -C5 that achieve a minimum number of active S-boxes are reported.