Improved Division Property Based Cube Attacks Exploiting Algebraic Properties of Superpoly

Cube attack, proposed by Dinur and Shamir [1] in 2009, is one of the general cryptanalytic techniques of analyzing symmetric-key cryptosystems. After its proposal, cube attack has been successfully applied to various ciphers, including stream ciphers [2,3,4,5,6], hash functions [7,8,9], and authenticated encryptions [10, 11]. For a cipher with n secret variables \(\varvec{x} = (x_1,x_2,\ldots ,x_n)\) and m public variables \(\varvec{v} = (v_1,v_2,\ldots ,v_m)\), we can regard the algebraic normal form (ANF) of output bits as a polynomial of \(\varvec{x}\) and \(\varvec{v}\), denoted as \(f(\varvec{x} , \varvec{v})\). For a randomly chosen set \(I= \{i_1 ,i_2 ,...,i_{|I|} \}\subset \{1,\ldots , m\}\), \(f(\varvec{x} , \varvec{v})\) can be represented uniquely as

where \(t_I=v_{i_1}\cdots v_{i_{|I|}}\), \(p(\varvec{x}, \varvec{v})\) only relates to \(v_s\)’s (\(s\notin I\)) and the secret key bits \(\varvec{x}\), and \(q(\varvec{x}, \varvec{v})\) misses at least one variable in \(t_I\). When \(v_s\)’s (\(s\notin I\)) and \(\varvec{x}\) are assigned statically, the value of \(p(\varvec{x}, \varvec{v})\) can be computed by summing the output bit \(f(\varvec{x}, \varvec{v})\) over a structure called cube, denoted as \(C_I\), consisting of \(2^{|I|}\) different \(\varvec{v}\) vectors with \(v_i, i\in I\) being active (traversing all 0-1 combinations) and non-cube indices \(v_s, s\notin I\) being static constants. Traditional cube attacks are mainly concerned about linear or quadratic superpolies. By collecting linear or quadratic equations from the superpoly, the attacker can recover some secret key bits information during the online phase. Aiming to mount distinguishing attack by property testing, cube testers are obtained by evaluating superpolies of carefully selected cubes. In [2], probabilistic tests are applied to detect some algebraic properties such as constantness, low degree and sparse monomial distribution. Moreover, cube attacks and cube testers are acquired experimentally by summing over randomly chosen cubes. So the sizes of the cubes are largely confined. Breakthroughs have been made by Todo et al. in [12] where they introduce the bit-based division property, a tool for conducting integral attacks^{Footnote 1}, to the realm of cube attack. With the help of mixed integer linear programming (MILP) aided division property, they can identify the variables excluded from the superpoly and explore cubes with larger size, e.g., 72 for 832-round Trivium. This enables them to improve the traditional cube attack.

Division property, as a generalization of the integral property, was first proposed at EUROCRYPT 2015 [13]. With division property, the propagation of the integral characteristics can be deduced in a more accurate manner, and one prominent application is the first theoretic key recovery attack on full MISTY1 [14].

The original division property can only be applied to word-oriented primitives. At FSE 2016, bit-based division property [15] was proposed to investigate integral characteristics for bit-based block ciphers. With the help of division property, the propagation of the integral characteristics can be represented by the operations on a set of 0-1 vectors identifying the bit positions with the zero-sum property. Therefore, for the first time, integral characteristics for bit-based block ciphers Simon32 and Simeck32 have been proved. However, the sizes of the 0-1 vector sets are exponential to the block size of the ciphers. Therefore, as has been pointed out by the authors themselves, the deduction of bit-based division property under their framework requires high memory for block ciphers with larger block sizes, which largely limits its applications. Such a problem has been solved by Xiang et al. [16] at ASIACRYPT 2016 by utilizing the MILP model. The operations on 0-1 vector sets are transformed to imposing division property values (0 or 1) to MILP variables, and the corresponding integral characteristics are acquired by solving the models with MILP solvers like Gurobi [17]. With this method, they are able to give integral characteristics for block ciphers with block sizes much larger than 32 bits. Xiang et al.’s method has now been applied to many other ciphers for improved integral attacks [18,19,20,21].

In [12], Todo et al. adapt Xiang et al.’s method by taking key bits into the MILP model. With this technique, a set of key indices \(J=\{j_1,j_2, \ldots , j_{|J|}\}\subset \{1,\ldots , n\}\) is deduced for the cube \(C_I\) s.t. \(p(\varvec{x}, \varvec{v})\) can only be related to the key bits \(x_j\)’s (\(j\in J\)). With the knowledge of I and J, Todo et al. can recover 1-bit of secret-key-related information by executing two phases. In the offline phase, a proper assignment to the non-cube IVs, denoted by \(\varvec{IV}\in \mathbb {F}_2^m\), is determined ensuring \(p(\varvec{x}, \varvec{IV})\) non-constant. Also in this phase, the whole truth table of \(p(\varvec{x}, \varvec{IV})\) is constructed through cube summations. In the online phase, the exact value of \(p(\varvec{x}, \varvec{IV})\) is acquired through a cube summation and the candidate values of \(x_j\)’s (\(j\in J\)) are identified by checking the precomputed truth table. A proportion of wrong keys are filtered as long as \(p(\varvec{x} , \varvec{IV})\) is non-constant.

Due to division property and the power of MILP solver, cubes of larger dimension can now be used for key recoveries. By using a 72-dimensional cube, Todo et al. propose a theoretic cube attack on 832-round Trivium. They also largely improve the previous best attacks on other primitives namely Acorn, Grain-128a and Kreyvium [12, 22]. It is not until recently that the result on Trivium has been improved by Liu et al. [6] mounting to 835 rounds with a new method called the correlation cube attack. The correlation attack is based on the numeric mapping technique first appeared in [23] originally used for constructing zero-sum distinguishers.

1.1 Motivations

Due to [12, 22], the power of cube attacks has been enhanced significantly, however, there are still problems remaining unhandled that we will reveal explicitly.

Finding Proper \(\varvec{IV}\)’s May Require Multiple Trials. As is mentioned above, the superpoly can filter wrong keys only if a proper IV assignment \(\varvec{IV}\in \mathbb {F}_2^m\) in the constant part of IVs is found such that the corresponding superpoly \(p(\varvec{x}, \varvec{IV})\) is non-constant. The MILP model in [12, 22] only proves the existence of the proper \(\varvec{IV}\)’s but finding them may not be easy. According to practical experiments, there are quite some \(\varvec{IV}\)’s making \(p(\varvec{x}, \varvec{IV})\equiv 0\). Therefore, \(t\ge 1\) different \(\varvec{IV}\)’s might be trailed in the precomputation phase before finding a proper one. Since each \(\varvec{IV}\) requires to construct a truth table with complexity \(2^{|I|+|J|}\), the overall complexity of the offline phase can be \(t\times 2^{|I|+|J|}\). When large cubes are used (|I| is big) or many key bits are involved (|J| is large), such a complexity might be at the risk of exceeding the brute-force bound \(2^n\). Therefore, two assumptions are made to validate their cube attacks as follows.

Assumption 1

(Strong Assumption). For a cube \(C_I\), there are many values in the constant part of IV whose corresponding superpoly is balanced.

Assumption 2

(Weak Assumption). For a cube \(C_I\), there are many values in the constant part of IV whose corresponding superpoly is not a constant function.

These assumptions are proposed to guarantee the validity of the attacks as long as \(|I|+|J|<n\), but the rationality of such assumptions is hard to be proved, especially when \(|I|+|J|\) are so close to n in many cases. The best solution is to evaluate different IVs in the MILP model so that the proper \(\varvec{IV}\) of the constant part of IVs and the set J are determined simultaneously before implementing the attack.

Restriction of \(|I|+|J|<n\). The superpoly recovery has always been dominating the complexity of the cube attack, especially in [12], the attacker knows no more information except which secret key bits are involved in the superpoly. Then she/he has to first construct the whole truth table for the superpoly in the offline phase. In general, the truth-table construction requires repeating the cube summation \(2^{|J|}\) times, and makes the complexity of the offline phase about \(2^{|I|+|J|}\). Apparently, such an attack can only be meaningful if \(|I|+|J|<n\), where n is the number of secret variables. The restriction of \(|I|+|J|<n\) barricades the adversary from exploiting cubes of larger dimension or mounting more rounds (where |J| may expand). This restriction can be removed if we can avoid the truth table construction in the offline phase.

1.2 Our Contributions

This paper improves the existing cube attacks by exploiting the algebraic properties of the superpoly, which include the (non-)constantness, low degree and sparse monomial distribution properties. Inspired by the division property based cube attack work of Todo et al. in [12], we formulate all these properties in one framework by developing more precise MILP models, thus we can reduce the complexity of superpoly recovery.

This also enables us to attack more rounds, or employ even larger cubes. Similar to [12], our methods regard the cryptosystem as a non-blackbox polynomial and can be used to evaluate cubes with large dimension compared with traditional cube attack and cube tester. In the following, our contributions are summarized into five aspects.

Flag Technique for Finding Proper IV Assignments. The previous MILP model in [12] has not taken the effect of constant 0/1 bits of the constant part of IVs into account. In their model, the active bits are initialized with division property value 1 and other non-active bits are all initialized to division property value 0. The non-active bits include constant part of IVs, together with some secret key bits and state bits that are assigned statically to 0/1 according to the specification of ciphers. It has been noticed in [22] that constant 0 bits can affect the propagation of division property. But we should pay more attention to constant 1 bits since constant 0 bits can be generated in the updating functions due to the XOR of even number of constant 1’s. Therefore, we propose a formal technique which we refer as the “flag” technique where the constant 0 and constant 1 as well as other non-constant MILP variables are treated properly. With this technique, we are able to find proper assignments to constant IVs (\(\varvec{IV}\)) that makes the corresponding superpoly (\(p(\varvec{x},\varvec{IV})\)) non-constant. With this technique, proper IVs can now be found with MILP model rather than time-consuming trial & summations in the offline phase as in [12, 22]. According to our experiments, the flag technique has a perfect 100% accuracy for finding proper non-cube IV assignments in most cases. Note that our flag technique has partially proved the availability of the two assumptions since we are able to find proper \(\varvec{IV}\)’s in all our attacks.

Degree Evaluation for Going Beyond the \(|I|+|J|<n\) Restriction. To avoid constructing the whole truth table using cube summations, we introduce a new technique that can upper bound the algebraic degree, denoted as d, of the superpoly using the MILP-aided bit-based division property. With the knowledge of its degree d (and key indices J), the superpoly can be represented with its \(\left( {\begin{array}{c}|J|\\ \le d\end{array}}\right) \) coefficients rather than the whole truth table, where \(\left( {\begin{array}{c}|J|\\ \le d\end{array}}\right) \) is defined as

When \(d=|J|\), the complexity by our new method and that by [12] are equal. For \(d<|J|\), we know that the coefficients of the monomials with degree higher than d are constantly 0. The complexity of superpoly recovery can be reduced from \(2^{|I|+|J|}\) to \(2^{|I|}\times \left( {\begin{array}{c}|J|\\ \le d\end{array}}\right) \). In fact, for some lightweight ciphers, the algebraic degrees of their round functions are quite low. Therefore, the degrees d are often much smaller than the number of involved key bits |J|, especially when high-dimensional cubes are used. Since \(d\ll |J|\) for all previous attacks, we can improve the complexities of previous results and use larger cubes mounting to more rounds even if \(|I|+|J|\ge n\).

Precise Term Enumeration for Further Lowering Complexities. Since the superpolies are generated through iterations, the number of higher-degree monomials in the superpoly is usually much smaller than its low-degree counterpart. For example, when the degree of the superpoly is \(d<|J|\), the number of d-degree monomials are usually much smaller than the upper bound \(\left( {\begin{array}{c}|J|\\ d\end{array}}\right) \). We propose a MILP model technique for enumerating all t-degree (\(t=1,\ldots ,d\)) monomials that may appear in the superpoly, so that the complexities of several attacks are further reduced.

Relaxed Term Enumeration. For some primitives (such as 750-round Acorn), our MILP model can only enumerate the d-degree monomials since the number of lower-degree monomials are too large to be exhausted. Alternately, for \(t=1,\ldots , d-1\), we can find a set of key indices \(JR_t\subseteq J\) s.t. all t-degree monomials in the superpoly are composed of \(x_j\), \(j\in JR_t\). As long as \(|JR_t|<|J|\) for some \(t=1,\ldots , d-1\), we can still reduce the complexities of superpoly recovery.

Combining the flag technique and the degree evaluation above, we are able to lower the complexities of the previous best cube attacks in [6, 12, 22]. Particularly, we can further provide key recovery results on 839-round Trivium^{Footnote 2}, 891-round Kreyvium, 184-round Grain-128a, and 750-round Acorn. Furthermore, the precise & relaxed term enumeration techniques allow us to lower the complexities of 833-round Trivium, 849-round Kreyvium, 184-round Grain-128a and 750-round Acorn. Our concrete results are summarized in Table 1.^{Footnote 3} In [26], Todo et al. revisit the fast correlation attack and analyze the key-stream generator (rather than the initialization) of the Grain family (Grain-128a, Grain-128, and Grain-v1). As a result, the key-stream generators of the Grain family are insecure. In other words, they can recover the internal state after initialization more efficiently than by exhaustive search. And the secret key is recovered from the internal state because the initialization is a public permutation. To the best of our knowledge, all our results of Kreyvium, Grain-128a, and Acornare the current best key recovery attacks on the initialization of the targeted ciphers. However, none of our results seems to threaten the security of the ciphers.

Clique View of the Superpoly Recovery. In order to lower the complexity of the superpoly recovery, the term enumeration technique has to execute many MILP instances, which is difficult for some applications. We represent the resultant superpoly as a graph, so that we can utilize the clique concept from the graph theory to upper bound the complexity of the superpoly recovery phase, without requiring MILP solver as highly as the term enumeration technique.

Organizations. Section. 2 provides the background of cube attacks, division property, MILP model etc. Section 3 introduces our flag technique for identifying proper assignments to non-cube IVs. Section 4 details the degree evaluation technique upper bounding the algebraic degree of the superpoly. Combining the flag technique and degree evaluation, we give improved key recovery cube attacks on 4 targeted ciphers in Sect. 5. The precise & relaxed term enumeration as well as their applications are given in Sect. 6. We revisit the term enumeration technique from the clique overview in Sect. 7. Finally, we conclude in Sect. 8.

2.1 Mixed Integer Linear Programming

MILP is an optimization or feasibility program whose variables are restricted to integers. A MILP model \(\mathcal {M}\) consists of variables \(\mathcal {M}.var\), constraints \(\mathcal {M}.con\), and an objective function \(\mathcal {M}.obj\). MILP models can be solved by solvers like Gurobi [17]. If there is no feasible solution at all, the solver simply returns infeasible. If no objective function is assigned, the MILP solver only evaluates the feasibility of the model. The application of MILP model to cryptanalysis dates back to the year 2011 [28], and has been widely used for searching characteristics corresponding to various methods such as differential [29, 30], linear [30], impossible differential [31, 32], zero-correlation linear [31], and integral characteristics with division property [16]. We will detail the MILP model of [16] later in this section.

2.2 Cube Attack

Considering a stream cipher with n secret key bits \(\varvec{x} = (x_1,x_2,\ldots ,x_n)\) and m public initialization vector (IV) bits \(\varvec{v} = (v_1,v_2,\ldots ,v_m)\). Then, the first output keystream bit can be regarded as a polynomial of \(\varvec{x}\) and \(\varvec{v}\) referred as \(f(\varvec{x}, \varvec{v})\). For a set of indices \(I=\{i_1,i_2,\ldots ,i_{|I|}\} \subset \{1,2,\ldots ,n\}\), which is referred as cube indices and denote by \(t_I\) the monomial as \(t_I = v_{i_1} \cdots v_{i_{|I|}}\), the algebraic normal form (ANF) of \(f(\varvec{x}, \varvec{v})\) can be uniquely decomposed as

where the monomials of \(q(\varvec{x}, \varvec{v})\) miss at least one variable from \(\{v_{i_1},v_{i_2},\ldots ,v_{i_{|I|}}\}\). Furthermore, \(p(\varvec{x}, \varvec{v})\), referred as the superpoly in [1], is irrelevant to \(\{v_{i_1},v_{i_2},\ldots ,v_{i_{|I|}}\}\). The value of \(p(\varvec{x}, \varvec{v})\) can only be affected by the secret key bits \(\varvec{x}\) and the assignment to the non-cube IV bits \(v_s\) (\(s\notin I\)). For a secret key \(\varvec{x}\) and an assignment to the non-cube IVs \(\varvec{IV}\in \mathbb {F}_2^m\), we can define a structure called cube, denoted as \(C_I(\varvec{IV})\), consisting of \(2^{|I|}\) 0-1 vectors as follows:

It has been proved by Dinur and Shamir [1] that the value of superpoly p corresponding to the key \(\varvec{x}\) and the non-cube IV assignment \(\varvec{IV}\) can be computed by summing over the cube \(C_I(\varvec{IV})\) as follows:

In the remainder of this paper, we refer to the value of the superpoly corresponding to the assignment \(\varvec{IV}\) in Eq. (3) as \(p_{\varvec{IV}}(\varvec{x})\) for short. We use \(C_I\) as the cube corresponding to arbitrary \(\varvec{IV}\) setting in Eq. (2). Since \(C_I\) is defined according to I, we may also refer I as the “cube” without causing ambiguities. The size of I, denoted as |I|, is also referred as the dimension of the cube.

Note: since the superpoly p is irrelevant to cube IVs \(v_i, i\in I\), the value of \(\varvec{IV}[i], i\in I\) cannot affect the result of the summation in Eq. (3) at all. Therefore in Sect. 5, our \(\varvec{IV}[i]\)’s (\(i\in I\)) are just assigned randomly to 0-1 values.

2.3 Bit-Based Division Property and Its MILP Representation

At 2015, the division property, a generalization of the integral property, was proposed in [13] with which better integral characteristics for word-oriented cryptographic primitives have been detected. Later, the bit-based division property was introduced in [15] so that the propagation of integral characteristics can be described in a more precise manner. The definition of the bit-based division property is as follows:

Definition 1

((Bit-Based) Division Property). Let \(\mathbb {X}\) be a multiset whose elements take a value of \(\mathbb {F}_2^n\). Let \(\mathbb {K}\) be a set whose elements take an n-dimensional bit vector. When the multiset \(\mathbb {X}\) has the division property \(\mathcal{D}_\mathbb {K}^{1^n}\), it fulfills the following conditions:

where \(\varvec{u} \succeq \varvec{k}\) if \(u_i \ge k_i\) for all i, and \(\varvec{x}^{\varvec{u}}=\prod _{i=1}^{n} x_i^{u_i}\).

When the basic bitwise operations COPY, XOR, AND are applied to the elements in \(\mathbb {X}\), transformations of the division property should also be made following the propagation corresponding rules copy, xor, and proved in [13, 15]. Since round functions of cryptographic primitives are combinations of bitwise operations, we only need to determine the division property of the chosen plaintexts, denoted by \(\mathcal{D}_{\mathbb {K}_0}^{1^n}\). Then, after r-round encryption, the division property of the output ciphertexts, denoted by \(\mathcal{D}_{\mathbb {K}_r}^{1^n}\), can be deduced according to the round function and the propagation rules. More specifically, when the plaintext bits at index positions \(I=\{i_1, i_2, \ldots , i_{|I|}\} \subset \{1,2,\ldots ,n\}\) are active (the active bits traverse all \(2^{|I|}\) possible combinations while other bits are assigned to static 0/1 values), the division property of such chosen plaintexts is \(\mathcal{D}_{\varvec{k}}^{1^n}\), where \(k_i=1\) if \(i \in I\) and \(k_i=0\) otherwise. Then, the propagation of the division property from \(\mathcal{D}_{\varvec{k}}^{1^n}\) is evaluated as

where \(\mathcal{D}_{\mathbb {K}_i}\) is the division property after i-round propagation. If the division property \(\mathbb {K}_r\) does not have an unit vector \(\varvec{e}_i\) whose only ith element is 1, the ith bit of r-round ciphertexts is balanced.

However, when round r gets bigger, the size of \(\mathbb {K}_r\) expands exponentially towards \(O(2^n)\) requiring huge memory resources. So the bit-based division property has only been applied to block ciphers with tiny block sizes, such as Simon32 and Simeck32 [15]. This memory-crisis has been solved by Xiang et al. using the MILP modeling method.

Propagation of Division Property with MILP. At ASIACRYPT 2016, Xiang et al. first introduced a new concept division trail defined as follows:

Definition 2

(Division Trail [16]). Let us consider the propagation of the division property \(\{\varvec{k}\} \overset{\underset{\mathrm {def}}{}}{=} \mathbb {K}_0 \rightarrow \mathbb {K}_1 \rightarrow \mathbb {K}_2 \rightarrow \cdots \rightarrow \mathbb {K}_r\). Moreover, for any vector \(\varvec{k}^*_{i+1} \in \mathbb {K}_{i+1}\), there must exist a vector \(\varvec{k}^*_{i} \in \mathbb {K}_{i}\) such that \(\varvec{k}^*_{i}\) can propagate to \(\varvec{k}^*_{i+1}\) by the propagation rule of the division property. Furthermore, for \((\varvec{k}_0, \varvec{k}_1,\ldots , \varvec{k}_r) \in (\mathbb {K}_0 \times \mathbb {K}_1 \times \cdots \times \mathbb {K}_r)\) if \(\varvec{k}_{i}\) can propagate to \(\varvec{k}_{i+1}\) for all \(i \in \{0,1,\ldots ,r-1\}\), we call \((\varvec{k}_0 \rightarrow \varvec{k}_1 \rightarrow \cdots \rightarrow \varvec{k}_r)\) an r-round division trail.

Let \(E_k\) be the target r-round iterated cipher. Then, if there is a division trail \(\varvec{k}_0 \xrightarrow {E_k} \varvec{k}_r=\varvec{e}_j\) (\(j=1,...,n\)), the summation of jth bit of the ciphertexts is unknown; otherwise, if there is no division trial s.t. \(\varvec{k}_0 \xrightarrow {E_k} \varvec{k}_r=\varvec{e}_j\), we know the ith bit of the ciphertext is balanced (the summation of the ith bit is constant 0). Therefore, we have to evaluate all possible division trails to verify whether each bit of ciphertexts is balanced or not. Xiang et al. proved that the basic propagation rules copy, xor, and of the division property can be translated as some variables and constraints of an MILP model. With this method, all possible division trials can be covered with an MILP model \(\mathcal {M}\) and the division property of particular output bits can be acquired by analyzing the solutions of the \(\mathcal {M}\). After Xiang et al.’s work, some simplifications have been made to the MILP descriptions of copy, xor, and in [12, 18]. We present the current simplest MILP-based copy, xor, and as follows:

Proposition 1

(MILP Model for COPY [18]). Let \(a \xrightarrow {COPY} (b_1,b_2,\ldots ,b_m)\) be a division trail of COPY. The following inequalities are sufficient to describe the propagation of the division property for copy.

Proposition 2

(MILP Model for XOR [18]). Let \((a_1, a_2, \ldots , a_m) \xrightarrow {XOR} b\) be a division trail of XOR. The following inequalities are sufficient to describe the propagation of the division property for xor.

Proposition 3

(MILP Model for AND [12]). Let \((a_1, a_2, \ldots , a_m) \xrightarrow {AND} b\) be a division trail of AND. The following inequalities are sufficient to describe the propagation of the division property for and.

Note: Proposition 3 includes redundant propagations of the division property, but they do not affect preciseness of the obtained characteristics [12].

2.4 The Bit-Based Division Property for Cube Attack

When the number of initialization rounds is not large enough for a thorough diffusion, the superpoly \(p(\varvec{x}, \varvec{v})\) defined in Eq. (2) may not be related to all key bits \(x_1,\ldots , x_n\) corresponding to some high-dimensional cube I. Instead, there is a set of key indices \(J\subseteq \{1,\ldots ,n\}\) s.t. for arbitrary \(\varvec{v}\in \mathbb {F}_2^m\), \(p(\varvec{x}, \varvec{v})\) can only be related to \(x_j\)’s (\(j\in J\)). In CRYPTO 2017, Todo et al. proposed a method for determining such a set J using the bit-based division property [12]. They further showed that, with the knowledge of such J, cube attacks can be launched to recover some information related to the secret key bits. More specifically, they proved the following Lemma 1 and Proposition 4.

Lemma 1

Let \(f(\varvec{x})\) be a polynomial from \(\mathbb {F}_2^n\) to \(\mathbb {F}_2\) and \(a_{\varvec{u}}^f\in \mathbb {F}_2\,(\varvec{u} \in \mathbb {F}_2^n)\) be the ANF coefficients of f(x). Let \(\varvec{k}\) be an n-dimensional bit vector. Assuming there is no division trail such that \(\varvec{k} \xrightarrow {f} 1\), then \(a_{\varvec{u}}^f\) is always 0 for \(\varvec{u} \succeq \varvec{k}\).

Proposition 4

Let \(f(\varvec{x}, \varvec{v})\) be a polynomial, where \(\varvec{x}\) and \(\varvec{v}\) denote the secret and public variables, respectively. For a set of indices \(I = \{i_1 ,i_2 ,\ldots ,i_{|I|} \} \subset \{1,2,\ldots ,m\}\), let \(C_I\) be a set of \(2^{|I|}\) values where the variables in \(\{v_{i_1} ,v_{i_2} ,\ldots , v_{i_{|I|}} \}\) are taking all possible combinations of values. Let \(\varvec{k}_I\) be an m-dimensional bit vector such that \(\varvec{v}^{\varvec{k}_I} = t_I = v_{i_1} v_{i_2} \cdots v_{i_{|I|}}\), i.e. \(k_i=1\) if \(i \in I\) and \(k_i=0\) otherwise. Assuming there is no division trail such that \((\varvec{e}_\lambda , \varvec{k}_I) \xrightarrow {f} 1\), \(x_\lambda \) is not involved in the superpoly of the cube \(C_I\).

When f represents the first output bit after the initialization iterations, we can identify J by checking whether there is a division trial \((\varvec{e}_\lambda , \varvec{k}_I) \xrightarrow {f} 1\) for \(\lambda =1,\ldots , n\) using the MILP modeling method introduced in Sect. 2.3. If the division trial \((\varvec{e}_\lambda , \varvec{k}_I) \xrightarrow {f} 1\) exists, we have \(\lambda \in J\); otherwise, \(\lambda \notin J\).

When J is determined, we know that for some assignment to the non-cube \(\varvec{IV}\in \mathbb {F}_2^m\), the corresponding superpoly \(p_{\varvec{IV}}(\varvec{x})\) is not constant 0, and it is a polynomial of \(x_j, j\in J\). With the knowledge of J, we recover offline the superpoly \(p_{\varvec{IV}}(\varvec{x})\) by constructing its truth table using cube summations defined as Eq. (3). As long as \(p_{\varvec{IV}}(\varvec{x})\) is not constant, we can go to the online phase where we sum over the cube \(C_I(\varvec{IV})\) to get the exact value of \(p_{\varvec{IV}}(\varvec{x})\) and refer to the precomputed truth table to identify the \(x_j, j\in J\) assignment candidates. We summarize the whole process as follows:

1. 1.

   Offline Phase: Superpoly Recovery. Randomly pick an \(\varvec{IV}\in \mathbb {F}_2^m\) and prepare the cube \(C_I(\varvec{IV})\) defined as Eq. (2). For \(\varvec{x}\in \mathbb {F}_2^n\) whose \(x_j,j\in J\) traverse all \(2^{|J|}\) 0-1 combinations, we compute and store the value of the superpoly \(p_{\varvec{IV}}(\varvec{x})\) as Eq. (3). The \(2^{|J|}\) values compose the truth table of \(p_{\varvec{IV}}(\varvec{x})\) and the ANF of the superpoly is determined accordingly. If \(p_{\varvec{IV}}(\varvec{x})\) is constant, we pick another \(\varvec{IV}\) and repeat the steps above until we find an appropriate one s.t. \(p_{\varvec{v}}(\varvec{x})\) is not constant.

2. 2.

   Online Phase: Partial Key Recovery. Query the cube \(C_I(\varvec{IV})\) to encryption oracle and get the summation of the \(2^{|I|}\) output bits. We denoted the summation by \(\lambda \in \mathbb {F}_2\) and we know \(p_{\varvec{IV}}(x)=\lambda \) according to Eq. (3). So we look up the truth table of the superpoly and only reserve the \(x_j, j\in J\) s.t. \(p_{\varvec{IV}}(x)=\lambda \).

3. 3.

   Brute-Force Search. Guess the remaining secret variables to recover the entire value in secret variables.

Phase 1 dominates the time complexity since it takes \(2^{|I|+|J|}\) encryptions to construct the truth table of size \(2^{|J|}\). It is also possible that \(p_{\varvec{IV}}(\varvec{x})\) is constant so we have to run several different \(\varvec{IV}\)’s to find the one we need. The attack can only be meaningful when (1) \(|I|+|J|<n\); (2) appropriate \(\varvec{IV}\)’s are easy to be found. The former requires the adversary to use “good” cube I’s with small J while the latter is the exact reason why Assumptions 1 and 2 are proposed [12, 22].

In the initial state of stream ciphers, there are secret key bits, public modifiable IV bits and constant 0/1 bits. In the previous MILP model, the initial bit-based division properties of the cube IVs are set to 1, while those of the non-cube IVs, constant state bits or even secret key bits are all set to 0.

Obviously, when constant 0 bits are involved in multiplication operations, it always results in an constant 0 output. But, as is pointed out in [22], such a phenomenon cannot be reflected in previous MILP model method. In the previous MILP model, the widely used COPY+AND operation:

can result in division trials \((x_1, x_2)\xrightarrow {COPY+AND}(y_1, y_2, a)\) as follows:

Assuming that either \(s_1\) or \(s_2\) of Eq. (4) is a constant 0 bit, \((s_1 \wedge s_2)\) is always 0. In this occasion, the division property of \((s_1 \wedge s_2)\) must be 0 which is overlooked by the previous MILP model. To prohibit the propagation above, an additional constraint \(\mathcal {M}.con \leftarrow a = 0\) should be added when either \(s_1\) or \(s_2\) is constant 0.

In [22], the authors only consider the constant 0 bits. They thought the model can be precise enough when all the state bits initialized to constant 0 bits are handled. But in fact, although constant 1 bits do not affect the division property propagation, we should still be aware because 0 bits might be generated when even number of constant 1 bits are XORed during the updating process. This is shown in Example 2 for Kreyvium in Appendix A [25].

Therefore, for all variables in the MILP \(v\in \mathcal {M}.var\), we give them an additional flag \(v.F\in \{1_c, 0_c, \delta \}\) where \(1_c\) means the bit is constant 1, \(0_c\) means constant 0 and \(\delta \) means variable. Apparently, when \(v.F=0_c/1_c\), there is always a constraint \(v=0\in \mathcal {M}.con\). We define \(=\), \(\oplus \) and \(\times \) operations for the elements of set \(\{1_c, 0_c, \delta \}\). The \(=\) operation tests whether two elements are equal(naturally \(1_c=1_c\), \(0_c=0_c\) and \(\delta =\delta \) ). The \(\oplus \) operation follows the rules:

The \(\times \) operation follows the rules:

Therefore, in the remainder of this paper, the MILP models for COPY, XOR and AND should also consider the effects of flags. So the previous copy, xor, and and should now add the assignment to flags. We denote the modified versions as copyf, xorf, and andf and define them as Propositions 5, 6 and 7 as follows.

Proposition 5

(MILP Model for COPY with Flag). Let \(a \xrightarrow {COPY} (b_1,b_2,\ldots ,b_m)\) be a division trail of COPY. The following inequalities are sufficient to describe the propagation of the division property for copyf.

We denote this process as \((\mathcal {M}, b_1,\ldots , b_m)\leftarrow \mathtt{copyf}(\mathcal {M}, a,m)\).

Proposition 6

(MILP Model for XOR with Flag). Let \((a_1, a_2, \ldots , a_m) \xrightarrow {XOR} b\) be a division trail of XOR. The following inequalities are sufficient to describe the propagation of the division property for xorf.

We denote this process as \((\mathcal {M}, b)\leftarrow \mathtt{xorf}(\mathcal {M}, a_1,\ldots , a_m)\).

Proposition 7

(MILP Model for AND with Flag). Let \((a_1, a_2, \ldots , a_m) \xrightarrow {AND} b\) be a division trail of AND. The following inequalities are sufficient to describe the propagation of the division property for andf.

We denote this process as \((\mathcal {M}, b)\leftarrow \mathtt{andf}\) \((\mathcal {M}, a_1,\ldots , a_m)\).

With these modifications, we are able to improve the preciseness of the MILP model. The improved attack framework can be written as Algorithm 1. It enables us to identify the involved keys when the non-cube IVs are set to specific constant 0/1 values by imposing corresponding flags to the non-cube MILP binary variables. With this method, we can determine an \(\varvec{IV}\in \mathbb {F}_2^m\) s.t. the corresponding superpoly \(p_{\varvec{IV}}(\varvec{x})\ne 0\).

For an \(\varvec{IV}\in \mathbb {F}_2^m\) s.t. \(p_{\varvec{IV}}(\varvec{x})\ne 0\), the ANF of \(p_{\varvec{IV}}(\varvec{x})\) can be represented as

where \(a_{\varvec{u}}\) is determined by the values of the non-cube IVs. If the degree of the superpoly is upper bounded by d, then for all \(\varvec{u}\)’s with Hamming weight satisfying \(hw(\varvec{u})>d\), we constantly have \(a_{\varvec{u}} =0\). In this case, we no longer have to build the whole truth table to recover the superpoly. Instead, we only need to determine the coefficients \(a_{\varvec{u}} \) for \(hw(\varvec{u})\le d\). Therefore, we select \(\sum _{i=0}^{d}\left( {\begin{array}{c}|J|\\ i\end{array}}\right) \) different \(\varvec{x}\)’s and construct a linear system with \(\left( \sum _{i=0}^{d}\left( {\begin{array}{c}|J|\\ i\end{array}}\right) \right) \) variables and the coefficients as well as the whole ANF of \(p_{\varvec{IV}}(\varvec{x} )\) can be recovered by solving such a linear system. So the complexity of Phase 1 can be reduced from \(2^{|I|+|J|}\) to \(2^{|I|}\times \sum _{i=0}^{d}\left( {\begin{array}{c}|J|\\ i\end{array}}\right) \). For the simplicity of notations, we denote the summation \(\sum _{i=0}^{d}\left( {\begin{array}{c}|J|\\ i\end{array}}\right) \) as \(\left( {\begin{array}{c}|J|\\ \le d\end{array}}\right) \) in the remainder of this paper. With the knowledge of the involved key indices \(J=\{ {j_1}, {j_2}, \ldots , {j_{|J|}}\}\) and the degree of the superpoly \(d=\deg p_{\varvec{IV}}(\varvec{x})\), the attack procedure can be adapted as follows:

1. 1.

   Offline Phase: Superpoly Recovery. For all \(\left( {\begin{array}{c}|J|\\ \le d\end{array}}\right) \) \(\varvec{x}\)’s satisfying \(hw(\varvec{x})\le d\) and \(\bigoplus _{j\in J}\varvec{e}_j\succeq \varvec{x}\), compute the values of the superpolys as \(p_{\varvec{IV}}(\varvec{x})\) by summing over the cube \(C_I(\varvec{IV})\) as Eq. (3) and generate a linear system of the \(\left( {\begin{array}{c}|J|\\ \le d\end{array}}\right) \) coefficients \(a_{\varvec{u}}\) (\(hw(\varvec{u})\le d\)). Solve the linear system, determine the coefficient \(a_{\varvec{u}}\) of the \(\left( {\begin{array}{c}|J|\\ \le d\end{array}}\right) \) terms and store them in a lookup table T. The ANF of the \(p_{\varvec{IV}}(\varvec{x})\) can be determined with the lookup table.

2. 2.

   Online Phase: Partial Key Recovery. Query the encryption oracle and sum over the cube \(C_I(\varvec{IV})\) as Eq. (3) and acquire the exact value of \(p_{\varvec{IV}}(\varvec{x})\). For each of the \(2^{|J|}\) possible values of \(\{x_{j_1}, \ldots , x_{j_{|J|}}\}\), compute the values of the superpoly as Eq. (7) (the coefficient \(a_{\varvec{u}}\) are acquired by looking up the precomputed table T) and identify the correct key candidates.

3. 3.

   Brute-force search phase. Attackers guess the remaining secret variables to recover the entire value in secret variables.

The complexity of Phase 1 becomes \(2^{|I|}\times \left( {\begin{array}{c}|J|\\ \le d\end{array}}\right) \). Phase 2 now requires \(2^{|I|}\) encryptions and \(2^{|J|}\times \left( {\begin{array}{c}|J|\\ \le d\end{array}}\right) \) table lookups, so the complexity can be regarded as \(2^{|I|}+ 2^{|J|}\times \left( {\begin{array}{c}|J|\\ \le d\end{array}}\right) \). The complexity of Phase 3 remains \(2^{n-1}\). Therefore, the number of encryptions a feasible attack requires is

The previous limitation of \(|I|+|J|<n\) is removed.

The knowledge of the algebraic degree of superpolys can largely benefit the efficiency of the cube attack. Therefore, we show how to estimate the algebraic degree of superpolys using the division property. Before the introduction of the method, we generalize Proposition 4 as follows.

Proposition 8

Let \(f(\varvec{x}, \varvec{v})\) be a polynomial, where \(\varvec{x}\) and \(\varvec{v}\) denote the secret and public variables, respectively. For a set of indices \(I = \{i_1 ,i_2 ,\ldots ,i_{|I|} \} \subset \{1,2,\ldots ,m\}\), let \(C_I\) be a set of \(2^{|I|}\) values where the variables in \(\{v_{i_1} ,v_{i_2} ,\ldots , v_{i_{|I|}} \}\) are taking all possible combinations of values. Let \(\varvec{k}_I\) be an m-dimensional bit vector such that \(\varvec{v}^{\varvec{k}_I} = t_I = v_{i_1} v_{i_2} \cdots v_{i_{|I|}}\). Let \(\varvec{k}_\varLambda \) be an n-dimensional bit vector. Assuming there is no division trail such that \((\varvec{k}_\varLambda || \varvec{k}_I) \xrightarrow {f} 1\), the monomial \(\varvec{x}^{\varvec{k}_\varLambda }\) is not involved in the superpoly of the cube \(C_I\).

Proof

The ANF of \(f(\varvec{x}, \varvec{v})\) is represented as follows

where \(a_{\varvec{u}}^f \in \mathbb {F}_2\) denotes the ANF coefficients. The polynomial \(f(\varvec{x}, \varvec{v})\) is decomposed into

Therefore, the superpoly \(p(\varvec{x}, \varvec{v})\) is represented as

Since there is no division trail \(( \varvec{k}_\varLambda \Vert \varvec{k}_I) \xrightarrow {f} 1\), \(a_{\varvec{u}}^f = 0\) for \(\varvec{u} \succeq (\varvec{k}_\varLambda \Vert \varvec{k}_I)\) because of Lemma 1. Therefore,

This superpoly is independent of the monomial \(\varvec{x}^{\varvec{k}_\varLambda }\) since \(\varvec{u}^{\varvec{k}_\varLambda \Vert \varvec{0}}\) is always 0.    \(\square \)

According to Proposition 8, the existence of the division trial \((\varvec{k}_\varLambda || \varvec{k}_I) \xrightarrow {f} 1\) is in accordance with the existence of the monomial \(x^{\varvec{k}_{\varLambda }}\) in the superpoly of the cube \(C_I\).

If there is \(d\ge 0\) s.t. for all \(\varvec{k}_\varLambda \) of hamming weight \(hw(\varvec{k}_\varLambda )>d\), the division trail \(x^{\varvec{k}_{\varLambda }}\) does not exist, then we know that the algebraic degree of the superpoly is bounded by d. Using MILP, this d can be naturally modeled as the maximum of the objective function \(\sum _{j=1}^n x_j\). With the MILP model \(\mathcal {M}\) and the cube indices I, we can bound the degree of the superpoly using Algorithm 2. Same with Algorithm 1, we can also consider the degree of the superpoly for specific assignment to the non-cube IVs. So we also add the input \(\varvec{IV}\) that can either be a specific assignment or a NULL referring to arbitrary assignment. The solution \(\mathcal {M}.obj=d\) is the upper bound of the superpoly’s algebraic degree. Furthermore, corresponding to \(\mathcal {M}.obj=d\) and according to the definition of \(\mathcal {M}.obj\), there should also be a set of indices \(\{l_1,\ldots , l_d\}\) s.t. the variables representing the initially declared \(\varvec{x}\) (representing the division property of the key bits) satisfy the constraints \(x_{l_1}=\ldots =x_{l_d}=1\). We can also enumerate all t-degree (\(1\le t\le d\)) monomials involved in the superpoly using a similar technique which we will detail later in Sect. 6.

We apply our method to 4 NLFSR-based ciphers namely Trivium, Kreyvium, Grain-128a and Acorn. Among them, Trivium, Grain-128a and Acornare also targets of [12]. Using our new techniques, we can both lower the complexities of previous attacks and give new cubes that mount to more rounds. We give details of the application to Trivium in this section, and the applications to Kreyvium, Grain-128a and Acornin our full version [25].

Structure of Trivium

Full size image

5.1 Specification of Trivium

Trivium is an NLFSR-based stream cipher, and the internal state is represented by 288-bit state \((s_1,s_2,\ldots ,s_{288})\). Figure 1 shows the state update function of Trivium. The 80-bit key is loaded to the first register, and the 80-bit IV is loaded to the second register. The other state bits are set to 0 except the least three bits in the third register. Namely, the initial state bits are represented as

The pseudo code of the update function is given as follows.

Here z denotes the 1-bit key stream. First, in the key initialization, the state is updated \(4 \times 288 = 1152\) times without producing an output. After the key initialization, one bit key stream is produced by every update function.

5.2 MILP Model of Trivium

The only non-linear component of Trivium is a 2-degree core function denoted as \(f_{core}\) that takes as input a 288-bit state \(\varvec{s}\) and 5 indices \(i_1,\ldots , i_5\), and outputs a new 288-bit state \(\varvec{s}'\leftarrow f_{core}(\varvec{s}, i_1,\ldots , i_5)\) where

The division property propagation for the core function can be represented as Algorithm 3. The input of Algorithm 3 consists of \(\mathcal {M}\) as the current MILP model, a vector of 288 binary variables \(\varvec{x}\) describing the current division property of the 288-bit NFSR state, and 5 indices \(i_1, i_2, i_3, i_4, i_5\) corresponding to the input bits. Then Algorithm 3 outputs the updated model \(\mathcal {M}\), and a 288-entry vector \(\varvec{y}\) describing the division property after \(f_{core}\).

With the definition of Core, the MILP model of R-round Trivium can be described as Algorithm 4. This algorithm is a subroutine of Algorithm 1 for generating the MILP model \(\mathcal {M}\), and the model \(\mathcal {M}\) can evaluate all division trails for Trivium whose initialization rounds are reduced to R. Note that constraints to the input division property are imposed by Algorithm 1.

5.3 Experimental Verification

Identical to [12], we use the cube \(I = \{ 1, 11, 21, 31, 41, 51, 61, 71 \}\) to verify our attack and implementation. The experimental verification includes: the degree evaluation using Algorithm 2, specifying involved key bits using Algorithm 1 with \(\varvec{IV}=\mathtt {NULL}\) or specific non-cube IV settings.

Example 1

(Verification of Our Attack against 591-round Trivium ). With \(\varvec{IV}=\mathtt {NULL}\) using Algorithm 1, we are able to identify \(J=\{23,24,25,66,67\}\). We know that with some assignment to the non-cube IV bits, the superpoly can be a polynomial of secret key bits \(x_{23}, x_{24}, x_{25}, x_{66}, x_{67}\). These are the same with [12]. Then, we set \(\varvec{IV}\) to random values and acquire the degree through Algorithm 2, and verify the correctness of the degree by practically recovering the corresponding superpoly.

• When we set \(\varvec{IV}=\mathtt{0xcc2e487b,0x78f99a93,0xbeae}\), and run Algorithm 2, we get the degree 3. The practically recovered superpoly is also of degree 3:

  $$\begin{aligned} p_{\varvec{v}}(\varvec{x}) = x_{66}x_{23}x_{24}+x_{66}x_{25}+x_{66}x_{67} +x_{66}, \end{aligned}$$

  which is in accordance with the deduction by Algorithm 2 through MILP model.

• When we set \(\varvec{IV}=\mathtt{0x61fbe5da, 0x19f5972c, 0x65c1}\), the degree evaluation of Algorithm 2 is 2. The practically recovered superpoly is also of degree 2:

  $$\begin{aligned} p_{\varvec{v}}(\varvec{x}) = x_{23}x_{24} + x_{25} + x_{67}+ 1. \end{aligned}$$

• When we set \(\varvec{IV}=\mathtt{0x5b942db1,0x83ce1016,0x6ce}\), the degree is 0 and the superpoly recovered is also constant 0.

On the Accuracy of MILP Model with Flag Technique. As a comparison, we use the cube above and conduct practical experiments on different rounds namely 576, 577, 587, 590, 591 (selected from Table 2 of [22]). We try 10000 randomly chosen \(\varvec{IV}\)’s. For each of them, we use the MILP method to evaluate the degree d, in comparison with the practically recovered ANF of the superpoly \(p_{\varvec{IV}}(\varvec{x})\). For 576, 577, 587 and 590 rounds, the accuracy is 100%. In fact, such 100% accuracy is testified for most of our applied ciphers, which is shown in [25]. For 591-round, the accuracies are distributed as:

1. 1.

   When the MILP model gives degree evaluation \(d=0\), the accuracy is 100% that the superpoly is constant 0.

2. 2.

   When the MILP model gives degree evaluation \(d=3\), there is an accuracy 49% that the superpoly is a 3-degree polynomial. For the rest, the superpoly is constant 0.

3. 3.

   When the MILP model gives degree evaluation \(d=2\), there is accuracy 43% that the superpoly is a 2-degree polynomial. For the rest, the superpoly is constant 0.

The ratios of error can easily be understood: for example, in some case, one key bit may multiply with constant 1 in one step \(x_i\cdot 1\) and be canceled by XORing with itself in the next round, this results in a newly generated constant 0 bit (\((x_i\cdot 1)\oplus x_i=0\)). However, by the flag technique, this newly generated bit has flag value \(\delta =(\delta \times 1_c)+\delta \). In our attacks, the size of cubes tends to be large, which means most of the IV bits become active, the above situation of \((x_i\cdot 1)\oplus x_i=0\) will now become \((x_i\cdot v_j)\oplus x_i\). Therefore, when larger cubes are used, fewer constant 0/1 flags are employed, and the MILP models are becoming closer to those of \(\varvec{IV}=NULL\). It is predictable that the accuracy of the flag technique tends to increase when larger cubes are used. To verify this statement, we construct a 10-dimensional cube \(I=\{5,13,18,22,30,57,60,65,72,79\}\) for 591-round Trivium. When \(\varvec{IV}=NULL\), we acquire the same upper bound of the degree \(d=3\). Then, we tried thousands of random IVs, and get an overall accuracy 80.9%. From above, we can conclude that the flag technique has high preciseness and can definitely improve the efficiency of the division property based cube attacks.

5.4 Theoretical Results

The best result in [12] mounts to 832-round Trivium with cube dimension \(|I|=72\) and the superpoly involves \(|J|=5\) key bits. The complexity is \(2^{77}\) in [12]. Using Algorithm 2, we further acquire that the degree of such a superpoly is 3. So the complexity for superpoly recovery is \(2^{72}\times \left( {\begin{array}{c}5\\ \le 3\end{array}}\right) =2^{76.7}\) and the complexity for recovering the partial key is \(2^{72}+2^3\times \left( {\begin{array}{c}5\\ 3\end{array}}\right) \). Therefore, according to Eq. (8), the complexity of this attack is \(2^{76.7}\).

We further construct a 77-dimensional cube, \(I=\{1,\ldots ,80\}\setminus \{5,51,65\}\). Its superpoly after 835 rounds of initialization only involves 1 key bit \(J=\{57\}\). So the complexity of the attack is \(2^{78}\). Since there are only 3 non-cube IVs, we let \(\varvec{IV}\) be all \(2^3\) possible non-cube IV assignments and run Algorithm 1. We find that \(x_{57}\) is involved in all of the \(2^3\) superpolys. So the attack is available for any of the \(2^3\) non-cube IV assignments. This can also be regarded as a support to the rationality of Assumption 1.

According previous results, Trivium has many cubes whose superpolys only contain 1 key bit. These cubes are of great value for our key recovery attacks. Firstly, the truth table of such superpoly is balanced and the Partial Key Recovery phase can definitely recover 1 bit of secret information. Secondly, the Superpoly Recovery phase only requires \(2^{|I|+1}\) and the online Partial Key Recovery only requires \(2^{|I|}\) encryptions. Such an attack can be meaningful as long as \(|I|+1<80\), so we can try cubes having dimension as large as 78. Therefore, we investigate 78-dimensional cubes and find the best cube attack on Trivium is 839 rounds. By running Algorithm 1 with \(2^2=4\) different assignments to non-cube IVs, we know that the key bit \(x_{61}\) is involved in the superpoly for \(\varvec{IV}=\mathtt {0x0, 0x4000, 0x0}\) or \(\varvec{IV}=\mathtt {0x0, 0x4002, 0x0}\). In other words, the 47-th IV bit must be assigned to constant 1. The summary of our new results about Trivium is in Table 2.

In this section, we show how to further lower the complexity of recovering the superpoly (Phase 1) in Sect. 4.

With cube indices I, key bits J and degree d, the complexity of the current superpoly recovery is \(2^{I}\times \left( {\begin{array}{c}|J|\\ \le d\end{array}}\right) \), where \(\left( {\begin{array}{c}|J|\\ \le d\end{array}}\right) \) corresponds to all 0-, 1-\(\ldots \), d-degree monomials. When \(d\le |J|/2\) (which is true in most of our applications), we constantly have \(\left( {\begin{array}{c}|J|\\ 0\end{array}}\right) \le \ldots \le \left( {\begin{array}{c}|J|\\ d\end{array}}\right) \). But in practice, high-degree terms are generated in later iterations and the high-degree monomials should be fewer than their low-degree counterparts. Therefore, for all \(\left( {\begin{array}{c}|J|\\ i\end{array}}\right) \) monomials, only very few of them may appear in the superpoly. Similar to Algorithm 1 that decides all key bits appear in the superpoly, we propose Algorithm 5 that enumerates all t-degree monomials that may appear in the superpoly. Apparently, when we use \(t=1\), we can get \(J_1=J\), the same output as Algorithm 1 containing all involved keys. If we use \(t=2,3,\ldots , d\), we get \(J_2,\ldots , J_d\) that contains all possible monomials of degrees \(2,3,\ldots , d\). Therefore, we only need to determine \(1+|J_1|+|J_2|+\ldots +|J_d|\) coefficients in order to recover the superpoly and apparently, \(|J_t|\le \left( {\begin{array}{c}|J|\\ t\end{array}}\right) \) for \(t=1,\ldots d\). With the knowledge of \(J_t, t=1,\ldots , d\), the complexity for Superpoly Recovery (Phase 1) has now become

And the size of the lookup table has also reduced to \((1+\sum _{t=1}^d|J_t|)\). So the complexity of the attack is now

Furthermore, since high-degree monomials are harder to be generated through iterations than low-degree ones, we can often find \(|J_i|<\left( {\begin{array}{c}|J|\\ i\end{array}}\right) \) when i approaches d. So the complexity for superpoly recovery has been reduced.

Note: \(J_t\)’s (\(t=1,\ldots , d\)) can be generated by TermEnum of Algorithm 5 and they satisfy the following Property 1. This property is equivalent to the “Embed Property” given in [19].

Property 1

For \(t=2,\ldots , d\), if there is \(T=(i_1, i_2,\ldots , i_t)\in J_t\) and \(T'=(i_{s_1},\ldots , i_{s_l})\) (\(l<t\)) is a subsequence of T (\(1\le s_1< \ldots < s_l \le t\)). Then, we constantly have \(T'\in J_l\).

Before proving Property 1, we first prove the following Lemma 2.

Lemma 2

If \(\varvec{k}\succeq \varvec{k}'\) and there is division trial \(\varvec{k}\xrightarrow {f} \varvec{l}\), then there is also division trial \(\varvec{k}'\xrightarrow {f} \varvec{l}'\) s.t. \(\varvec{l}\succeq \varvec{l}'\).

Proof

Since f is a combination of COPY, AND and XOR operations, and the proofs when f equals to each of them are similar, we only give a proof of the case when f equals to COPY. Let \(f: (*,\ldots , *, x)\xrightarrow {COPY} (*,\ldots , *, x,x)\).

First assume the input division property be \(\varvec{k}=(\varvec{k}_1,0)\), since \(\varvec{k}\succeq \varvec{k}'\), there must be \(\varvec{k}'=(\varvec{k}_1',0)\) and \(\varvec{k}_1\succeq \varvec{k}_1'\). We have \(l=k\), \(l'=k'\), thus the property holds.

When the input division property is \(\varvec{k}=(\varvec{k}_1,1)\), we know that the output division property can be \(\varvec{l}\in \{(\varvec{k}_1, 0,1),(\varvec{k}_1, 1,0)\}\). Since \(\varvec{k}\succeq \varvec{k}'\), we know \(\varvec{k}'=(\varvec{k}_1',1)\) or \(\varvec{k}'=(\varvec{k}_1',0)\), and \(\varvec{k}_1\succeq \varvec{k}_1'\). When \(\varvec{k}'=(\varvec{k}_1',0)\), then \(l'=k'=(k_1',0)\), the relation holds. When \(\varvec{k}'=(\varvec{k}_1',1)\), we know \(\varvec{l}'\in \{(\varvec{k}_1', 0,1),(\varvec{k}_1', 1,0)\}\), the relation still holds.   \(\square \)

Now we are ready to prove Property 1.

Proof

Let \(\varvec{k},\varvec{k}\in \mathbb {F}_2^n\) satisfy \(k_i=1\) for \(i\in T\) and \(k_i=0\) otherwise; \(k_i'=1\) for \(i\in T'\) and \(k_i'=0\) otherwise. Since \(T\in J_t\), we know that there is division trial \((\varvec{k}, \varvec{k}_I)\xrightarrow {R-Rounds}(\varvec{0} ,1)\) Since \(k\succeq k'\), we have \((\varvec{k}, \varvec{k}_I)\succeq (\varvec{k}', \varvec{k}_I)\) and according to Lemma 2, there is division trial s.t. \((\varvec{k}', \varvec{k}_I)\xrightarrow {R-Rounds}(\varvec{0}^{m+n},s)\) where \((\varvec{0}^{m+n},1)\succeq (\varvec{0}^{m+n},s)\). Since the hamming weight of \((\varvec{k}', \varvec{k}_I)\) is larger than 0 and there is no combination of COPY, AND and XOR that makes non-zero division property to all-zero division property. So we have \(s=1\) and there exist division trial \((\varvec{k}', \varvec{k}_I)\xrightarrow {R-Rounds}(\varvec{0} ,1)\).    \(\square \)

Property 1 reveals a limitation of Algorithm 5. Assume the superpoly is

We can acquire \(J_3= \{(1,2,3)\}\) by running TermEnum of Algorithm 5. But, if we run TermEnum with \(t=2\), we will not acquire just \(J_2=\{(1,4)\}\) but \(J_2=\{(1,4), (1,2), (1,3), (2,3)\}\) due to \((1,2,3)\in J_3\) and (1, 2), (1, 3), (2, 3) are its subsequences. Although there are still redundant terms, the reduction from \(\left( {\begin{array}{c}|J|\\ d\end{array}}\right) \) to \(|J_d|\) is usually huge enough to improve the existing cube attack results.

Applying such term enumeration technique, we are able to lower complexities of many existing attacks namely: 832-, 833-round Trivium, 849-round Kreyvium, 184-round Grain-128a and 704-round Acorn. The attack on 750-round Acorncan also be improved using a relaxed version of TermEnum which is presented as RTermEnum on the righthand side of Algorithm 5. In the relaxed algorithm, RTermEnum is acquired from TermEnum by replacing some states which are marked in red in Algorithm 5, and we state details later in Sect. 6.4.

6.1 Application to Trivium

As can be seen in Table 2, the attack on 832-round Trivium has \(J=J_1=5\) and degree \(d=3\), so we have \(\left( {\begin{array}{c}5\\ \le 3\end{array}}\right) =26\) using previous technique. But by running Algorithm 5, we find that \(|J_2|=5\), \(|J_3|=1\), so we have \(1+\sum _{t=1}^3|J_t|=12<\left( {\begin{array}{c}5\\ \le 3\end{array}}\right) =26\). Therefore, the complexity has now been reduced from \(2^{76.7}\) to \(2^{75.8}\). Similar technique can also be applied to the 73 dimensional cube of Table 2. Details are shown in Table 3.

6.2 Applications to Kreyvium

We revisit the 61-dimensional cube first given in [23] and transformed to a key recovery attack on 849-round Kreyvium in [22]. The degree of the superpoly is 9, so the complexity is given as \(2^{81.7}\) in Appendix A of [25]. Since \(J=J_1\) is of size 23, we enumerate all the terms of degree 2–9 and acquire the sets \(J_2,\ldots , J_9\). \(1+\sum _{t=1}^d|J_t|=5452\approx 2^{12.41}\). So the complexity is now lowered to \(2^{73.41}\). The details are listed in Table 4.

6.3 Applications to Grain-128a

For the attack on 184-round Grain-128a, the superpoly has degree \(d=14\), the number of involved key bits is \(|J|=|J_1|=21\) and we are able to enumerate all terms of degree 1–14 as Table 5.

6.4 Applications to Acorn

For the attack on 704-round Acorn, with the cube dimension 64, the number of involved key bits in the superpoly is 72, and the degree is 7. We enumerate all the terms of degree from 2 to 7 as in Table 6, therefore we manage to improve the complexity of our cube attack in the previous section.

Relaxed Algorithm 5. For the attack on 750-round Acorn(the superpoly is of degree \(d=5\)), The left part of Algorithm 5 can only be carried out for the 5-degree terms \(|J_5|=46\). For \(t=2,3,4\), the sizes of \(J_t\) are too large to be enumerated. We settle for the index set \(JR_t\) containing the key indices that composing all the t-degree terms. For example, when \(J_3=\{(1,2,3), (1,2,4)\}\), we have \(JR_3=\{1,2,3,4\}\). The relationship between \(J_t\) and \(JR_t\) is \(|J_t|\le \left( {\begin{array}{c}|JR_t|\\ t\end{array}}\right) \) and \(J_1=JR_1\). The searching space for \(J_t\) in Algorithm 5 is \(\left( {\begin{array}{c}|J_1|\\ t\end{array}}\right) \) while that of the relaxed algorithm is only \(\left( {\begin{array}{c}|JR_t|\\ t\end{array}}\right) \). So it is much easier to enumerate \(JR_t\), therefore the complexity can still be improved (in comparison with Eq. (8)) as long as \(|JR_t|<|J_1|\). The complexity of this relaxed version can be written as

For 750-round Acorn, we enumerate \(J_5\) and \(JR_1,\ldots , JR_4\) whose sizes are listed in Table 7. The improved complexity, according to Eq. (12), is \(2^{120.92}\), lower than the original \(2^{125.71}\) given in Appendix A in [25].

The precise & relaxed term enumeration technique introduced in Sect. 6 have to execute many MILP instances, which is difficult for some applications. In this section, we represent the resultant superpoly as a graph, which is called superpoly graph, so that we can utilize the clique concept from the graph theory to upper bound the complexity of the superpoly recovery phase in our attacks, without requiring MILP solver as highly as the term enumeration technique.

Definition 3

(Clique [33]). In a graph \(G =(V,E)\), where V is the set of vertices and E is the set of edges, a subset \(C \subseteq V\), s.t. each pair of vertices in C is connected by an edge is called a clique.

A i-clique is defined as a clique consists of i vertices, and i is called the clique number. A 1-clique is a vertex, a 2-clique is just an edge, and a 3-clique is called a triangle.

Given a cube \(C_I\), by running Algorithm 5 for degree i, we determine \(J_i\), which is the set of all the degree-i terms that might appear in the superpoly \(p(\varvec{x},\varvec{v})\) (see Sect. 6). Then we represent \(p(\varvec{x}, \varvec{v})\) as a graph \(G=(J_1, J_2)\), where the vertices in \(J_1\) correspond to the involved secret key bits in \(p(\varvec{x},\varvec{v})\), the edges between any pairs of the vertices reveal the quadratic terms involved in \(p(\varvec{x},\varvec{v})\), We call the graph \(G=(J_1, J_2)\) the superpoly graph of the cube \(C_I\). The set of i-cliques in the superpoly graph is denoted as \(\mathcal {K}_i\). Note that there is a natural one-to-one correspondence between the sets \(J_i\) and \(\mathcal {K}_i\) for \(i=1,2\).

It follows from the definition of a clique that any i-clique in \(\mathcal {K}_i\) (\(i\ge 2\)) represents a monomial of degree i whose all divisors of degree 2 belong to \(J_{2}\). On the other hand, due to the “embed” Property 1 in Sect. 6, we have that all its quadratic divisors must be in \(J_{2}\). Then any monomial in \(J_i\) can be represented by an i-clique in \(\mathcal {K}_i\). Hence for all \(i\ge 2\), \(J_i\) corresponds to a subset of \(\mathcal {K}_i\). Denote the number of i-cliques as \(|\mathcal {K}_i|\), then \(|J_{i}|\le |\mathcal {K}_i|\). Apparently, \(|\mathcal {K}_i|\le \left( {\begin{array}{c}|J|\\ i\end{array}}\right) \) for all \(1\le i\le d\).

Now we show a simple algorithm for constructing \(\mathcal {K}_i\) from \(J_1\) and \(J_2\) for \(i\ge 3\). For instance, when constructing \(\mathcal {K}_{3}\), we take the union operation of all possible combinations of three elements from \(J_{2}\), and only keep the elements of degree 3. Similarly, we construct \(\mathcal {K}_{i}\) for \(3 < i\le d\), where d is the degree of the superpoly. Therefore, all the i-cliques (\(3\le i\le d\)) are found by the simple algorithm, i.e. the number of i-cliques \(|\mathcal {K}_i|\) in \(G(J_1,J_2)\) is determined. We therefore can upper bound the complexity of the offline phase as

Note that we have \(|J_i| \le \ |\mathcal {K}_i|\le \left( {\begin{array}{c}|J_1|\\ i\end{array}}\right) \). It indicates that the upper bound of the superpoly recovery given by clique theory in Eq. (13) is better than the one provided by our degree evaluation in Eq. (8), while it is weaker than the one presented by our term enumeration techniques in Eq. (10). However, it is unclear if there exists a specific relation between \(|\mathcal {K}_i|\) and \(\left( {\begin{array}{c}|JR_i|\\ i\end{array}}\right) \) in the relaxed terms enumeration technique.

Advantage over the Terms Enumeration Techniques. In Sect. 6 when calculating \(J_i\) \((i\ge 3)\) by Algorithm 5, we set the target degree as i and solve the newly generated MILP to obtain \(J_i\), regardless of the knowledge of \(J_{i-1}\) we already hold. On the other hand, as is known in some cases, the MILP solver might take long time before providing \(J_i\) as desired. However, by using clique theory, we first acquire \(J_1\) and \(J_2\), which are essential for the term enumeration method as well. According to the “embed” property, we then make full use of the knowledge of \(J_1\) and \(J_2\), to construct \(\mathcal {K}_i\) for \(i\ge 3\) by an algorithm which is actually just performing simple operations (like union operations among elements, or removal of repeated elements, etc) in sets. So hardly any cost is required to find all the \(\mathcal {K}_i\) (\(3\le i\le d\)) we want. This significantly saves the computation costs since solving MILP is usually very time-consuming.

Algebraic properties of the resultant superpoly of the cube attacks were further studied. We developed a division property based framework of cube attacks enhanced by the flag technique for identifying proper non-cube IV assignments. The relevance of our framework is three-fold: For the first time, it can identify proper non-cube IV assignments of a cube leading to a non-constant superpoly, rather than randomizing trails & summations in the offline phase. Moreover, our model derived the upper bound of the superpoly degree, which can break the \(|I|+|J|<n\) barrier and enable us to explore even larger cubes or mount to attacks on more rounds. Furthermore, our accurate term enumeration techniques further reduced the complexities of the superpoly recovery, which brought us the current best key recovery attacks on ciphers namely Trivium, Kreyvium, Grain-128a and Acorn.

Besides, when term enumeration cannot be carried out, we represent the resultant superpoly as a graph. By constructing all the cliques of our superpoly graph, an upper bound of the complexity of the superpoly recovery can be obtained.