Multidimensional Linear Cryptanalysis

Linear cryptanalysis introduced by Matsui is a statistical attack which exploits a binary linear relation between plaintext, ciphertext and key, either in Algorithm 1 for recovering one bit of information of the secret key of a block cipher, or in Algorithm 2 for ranking candidate values for a part of the key. The statistical model is based on the expected and observed bias of a single binary value. Multiple linear approximations have been used with the goal to make the linear attack more efficient. More bits of information of the key can potentially be recovered possibly using less data. But then also more elaborated statistical models are needed to capture the joint behaviour of several not necessarily independent binary variables. Also more options are available for generalising the statistics of a single variable to several variables. The multidimensional extension of linear cryptanalysis to be introduced in this paper considers using multiple linear approximations that form a linear subspace. Different extensions of Algorithm 1 and Algorithm 2 will be presented and studied. The methods will be based on known statistical tools such as goodness-of-fit test and log-likelihood ratio. The efficiency of the different methods will be measured and compared in theory and experiments using the concept of advantage introduced by Selçuk. The block cipher Serpent with a reduced number of rounds will be used as test bed. The multidimensional linear cryptanalysis will also be compared with previous methods that use biasedness of multiple linear approximations. It will be shown in the simulations that the multidimensional method is potentially more powerful. Its main theoretical advantage is that the statistical model can be given without the assumption about statistical independence of the linear approximations.


General
The purpose of this paper is to examine how the efficiency of key recovery attacks using linear cryptanalysis can be improved by extending the attack to multiple dimensions. One-dimensional linear cryptanalysis was introduced by Matsui [27] in 1993. The method is based on one binary linear relation involving some bits of the plaintext, ciphertext and the secret key of a block cipher. The linear relation holds for a fraction of plaintexts and therefore is called a linear approximation of the cipher. If the fraction of plaintexts that satisfy the relation deviates significantly from one half, the approximation has a large bias. In this case, the approximation is called strong, and it can be used for recovering information about the key, provided that the attacker has enough data, i.e. plaintext-ciphertext pairs.
Matsui presented two algorithms, Algorithm 1 (Alg. 1) and Algorithm 2 (Alg. 2) for iterated block ciphers. While Alg. 1 extracts one bit of information about the secret key, Alg. 2 can be used in finding several bits of the last round key of the cipher. The candidate values for the part of the last round key to be recovered are ranked according to a test statistic and the right value is expected to have the highest rank. Using the recovered part of the last round key, it is then possible to extract one more bit of information about the secret key using Alg. 1.

Related Work
In practice, it is often possible to find several strong linear approximations. Hence, the obvious enhancement to linear cryptanalysis is to use multiple linear approximations simultaneously. Matsui considered using two approximations already in [26]. Junod and Vaudenay analysed this approach in [22]. In an attempt to reduce the data complexities of Matsui's algorithms, Kaliski and Robshaw used several linear approximations involving the same key bits [25]. Multiple linear approximations were also used by Biryukov et al. [4] for extracting several bits of information about the key in an Alg. 1-type attack. They also extended this basic method to a combination of Alg. 1-and Alg. 2-type attacks. However, the results in both [4,25] depend on theoretical assumptions about the statistical properties of the one-dimensional linear approximations. In particular, it is assumed that they are statistically independent. Murphy noted that this assumption does not hold in general [29]. Moreover, practical experiments by Hermelin et al. [20] showed that the assumption does not always hold in the case of the block cipher Serpent.
The statistical linear distinguisher presented by Baignères et al. does not suffer from this limitation, since it works with distributions of data values instead of biases of single linear approximations [1]. This approach is applicable in case the linear approximations form a linear subspace. A linear subspace of linear approximations is called a multidimensional linear approximation. The solution presented in [1] has also another advantage over the previous approaches in [4,25]: it is based on a well-established statistical theory of log-likelihood ratio, LLR; see also [24]. Early work to this direction by Vaudenay pro-poses to use χ 2 -cryptanalysis [34]. More recently, Baignères and Vaudenay [2] studied hypothesis testing related to different distinguishing scenarios.
To realise a multidimensional linear distinguishing attack, it is necessary to calculate the probability distribution for the multidimensional approximation of the cipher. Englund and Maximov [16] and Maximov and Johansson [28] studied algorithms for computing large distributions and applied them to compute probability distributions over output domains of functions and their compositions used in building the cipher. However, this approach gets soon infeasible when the input and output sizes of the functions and their compositions exceed 32 bits. Multidimensional linear cryptanalysis allows to focus on the most essential information and control the sizes of domains over which the probability distributions are computed. Given a suitable number of strong one-dimensional linear approximations, the linear space spanned by them forms the multidimensional linear approximation for the attack. From the estimated values of biases of all the linear approximations in this linear space, one can estimate the probability distribution of the multidimensional values. The problem is then reduced to the problem of finding strong one-dimensional approximations.

Contributions
In this paper, the statistical theory of multidimensional linear distinguisher will be developed for extending Matsui's Alg. 1 and Alg. 2 to multiple dimensions. While in dimension one there is essentially only one statistic, the bias of a linear approximation, for realising Alg. 1 and Alg. 2 in multiple dimensions, there are several possible statistical interpretations for the key recovery problem and different statistical tests for solving them. The purpose of this work is to compare different key recovery methods for both Alg. 1 and Alg. 2.
The different methods will be compared by using the concept of advantage proposed by Selçuk [33]. Originally, the advantage was proposed to be used in measuring the success of key ranking in the one-dimensional Alg. 2. In this paper, the theory of advantage is extended to multiple dimensions and applied to study the key ranking in both Alg. 1 and Alg. 2. The advantage for different methods is determined in theory and evaluated in experiments for block cipher Serpent with a reduced number of rounds.
Both Alg. 1 and Alg. 2 can be interpreted as goodness-of-fit problems that can be solved using χ 2 -test. A method based on LLR will also be studied for both algorithms. While for Alg. 1 these two methods seem to be equally efficient in practice, for Alg. 2 both theory and practice demonstrate superiority of the LLR-method over the χ 2 -based method. The multidimensional methods will also be compared to the classical onedimensional method. We also define a combination of the χ 2 -based Alg. 1 and Alg. 2 and show that it is essentially identical to the combined method proposed by Biryukov et al. if in the latter, all (nonzero) linear combinations of the base approximations are included in the attack.

Outline
The structure of this paper is as follows: in Sect. 2, the basic statistical theory and notation is given. Results about multidimensional linear distinguishers are presented as well as the construction of the multidimensional probability distribution using the one-dimensional linear approximations of the cipher. The advantage and the generalisation of Selçuk's theory to multiple dimensions are presented in Sect. 4. In Sect. 3, the multidimensional linear approximation of a block cipher is discussed. Different extensions of Alg. 1 to multiple dimensions are introduced in Sect. 5 with the theoretical and experimental results. The time, data and memory complexities of Alg. 1 are also considered. Similarly, Alg. 2 is studied in Sect. 6.

Boolean Function and Probability Distribution
The space of n-dimensional binary vectors is denoted by F n 2 . The sum modulo 2 is denoted by ⊕. The inner product for a = (a 1 , . . . , a n ), b = (b 1 , . . . , b n ) ∈ F n 2 is defined as a · b = a 1 b 1 ⊕ · · · ⊕ a n b n . Then the vector a is called the (linear) mask of b.
A function f : F n 2 → F 2 is called a Boolean function. A linear Boolean function is a mapping . . , f m ), where f i are Boolean functions, is called a vector Boolean function of dimension m. A linear Boolean function from F n 2 to F m 2 is represented by an m × n binary matrix U . The m rows of U are denoted by u 1 , . . . , u m , where each u i is a linear mask.
The correlation between a Boolean function f : F n 2 → F 2 and zero is and it is also called the correlation of f. We say that the vector p = ( p 0 , . . . , p M ) is a probability distribution (p.d.) of a random variable X taking on values in {0, 1, . . . , M} and denote X ∼ p, if Pr(X = η) = p η , for all η = 0, . . . , M. We will denote the uniform p.d. by θ . Let f : F n 2 → F m 2 and X ∼ θ. We call the p.d. p of the random variable f (X ) the p.d. of f . Let us recall some general properties of p.d.'s. Let p = ( p 0 , . . . , p M ) and q = (q 0 , . . . , q M ) be p.d.'s of random variables taking on values in a set with M +1 elements. The Kullback-Leibler distance between p and q is defined as follows: Definition 2.1. The relative entropy or Kullback-Leibler distance between p and q is with the conventions 0 log 0/b = 0, for all b ≥ 0 and b log b/0 = ∞ for b > 0.
The following property usually holds for p.d.'s related to any real ciphers, so it will be frequently used throughout this work.
This is a natural property of modern ciphers with q = θ as one of their design criteria is to resist one-dimensional linear cryptanalysis, and therefore, they must have as uniform p.d.'s as possible.
If p is close to q, we can approximate their Kullback-Leibler distance using the Taylor series [1] such that where ε = max η∈{0,1,...,M} | p η − q η | and the capacity C( p, q) of p and q is defined as follows: Definition 2.2. The capacity between two p.d.'s p and q is defined by If q is the uniform distribution, then C( p, q) will be denoted by C( p) and called the capacity of p.
The normed normal distribution with mean 0 and variance 1 is denoted by N (0, 1). Its probability density function (p.d.f.) is The normal distribution with mean μ and variance σ 2 is denoted by N (μ, σ 2 ), and its p.d.f. and c.d.f. are φ μ,σ 2 and Φ μ,σ 2 , respectively. The following standard approximation of c.d.f. Φ(x) can be found in common statistical reference books, e.g. [32]. The approximation holds for large values of x and will be used in the derivations of this paper. Let X 1 , . . . , X M be independent random variables, with . . , M, are equal to zero, then the sum of the squares is χ 2 M -distributed with M degrees of freedom and has mean M and variance 2M. Otherwise the sum (2) follows the non-central χ 2 M (λ) distribution which has mean λ + M and variance 2(M + 2λ). If M > 30, we may approximate χ 2 M (λ) ∼ N (λ + M, 2(M + 2λ)) [13].
Let X 1 , . . . , X N be a sequence of independent and identically distributed (i.i.d.) random variables where either X 1 , . . . , X N ∼ p (corresponding to null hypothesis H 0 ) or X 1 , . . . , X N ∼ q = p (corresponding to alternate hypothesis H 1 ). The hypothesis testing problem is then to determine whether to accept or reject H 0 .
We can make two types of error in the test. In type I error, we reject H 0 when it is true. The level α of the test measures how probably this will happen: α = Pr(H 1 | H 0 ). In type II error, we accept H 0 , when it is not true. This is measured by the power 1 − β of the test, defined as β = Pr(H 0 | H 1 ).
According to the Neyman-Pearson lemma [11], given empirical datax 1 , . . . ,x N , the optimal statistic for solving this problem, that is, distinguishing between p and q, is the log-likelihood ratio (LLR) defined by The distinguisher accepts H 0 and outputs p (or rejects H 0 and outputs q, ) if LLR(q, p, q) ≥ τ (or LLR(q, p, q) < τ, respectively) where τ is the threshold. The threshold depends on the level and the power of the test. Usually τ = 0.
The proof of the following result can be found in [11], see also [1].  (3) is asymptotically normal with mean and variance N μ 0 and N σ 2 0 (N μ 1 and N σ 2 1 , resp.) if the data are drawn from p (q, resp.). The means and variances are given by Moreover, if p is close to q, the following estimates hold The data complexity of the test is defined as the amount of data N needed for successfully solving the hypothesis testing problem between p and q with given power and level of test. Assuming that p is close to q, the following corollary is obtained from Proposition 2.1 [1].

Multidimensional Approximation of Boolean Functions
Let f : F n 2 → F n 2 be a vector Boolean function. Its one-dimensional linear approximation with input mask u ∈ F n 2 and output mask w ∈ F n 2 is the Boolean function with some (non-negligible) correlation c. For many ciphers, it is possible to find several such approximations (4) with non-negligible correlations. The aim of multidimensional linear cryptanalysis is to efficiently exploit given one-dimensional approximations with non-negligible correlations to obtain information about the cipher. Linear approximations are said to be linearly independent if the mask pairs (u i , w i ), i = 1, . . . , m, considered as concatenated vectors of length n + n are linearly independent. Given a set of one-dimensional approximations of f , let m be the dimension of the linear space spanned by them. We call this set a multidimensional linear approximation of f of dimension m and it can be given by where ⊕ is the component-wise modulo 2 sum and U = (u 1 , . . . , u m ) and W = (w 1 , . . . , w m ) are linear matrices from F n 2 → F m 2 and F n 2 → F m 2 , respectively. The linear approximations u i · x ⊕ w i · f (x), i = 1, . . . , m, are called base approximations.
Then we need to calculate the multidimensional p.d. p of the m-dimensional approximation. We observe that it defines a vector Boolean function and recall that by the p.d. of a vector Boolean function g : F n 2 → F m 2 we mean the p.d. of the random variable g(X ), where X ∼ θ . The following result connects the p.d. of a vector Boolean function g and its one-dimensional correlations c(a · g); see, for example, [20] or [18]. Hence, for determining the p.d. p of approximation (6), we have to determine the correlations c(a · (U x ⊕ W f (x))) of all the one-dimensional approximations of f. That is, p is determined based on one-dimensional projections of f , which is a known principle in statistics due to Cramér and Wold [12].
The following corollary is obtained from Proposition 2.2 using Parseval's theorem. An equivalent form of it can be found in [1], where the proof was based on the inverse Walsh-Hadamard transform of the deviations ε η from the uniform distribution, ε η = p η − 2 −m . Corollary 2.2. Let g : F n 2 → F m 2 be a Boolean function with p.d. p and one-dimensional correlations c(a · g). Then c(a · g) 2 .
We will need this equality in the next section where we study how linear distinguishing is done in multiple dimensions.

Multidimensional Linear Distinguishers
A linear distinguisher determines whether given datax t , t = 1, . . . , N , are drawn from a cipher with p.d. p = θ or a random source with p.d. θ . In the one-dimensional case, the attacker uses one linear approximation such as (4) with correlation c. The data complexity is inversely proportional to c 2 [23,27]. When using multiple linear approximations with non-negligible correlations c i = c(u i · x ⊕ w i · f (x)) and drawing datax t , t = 1, . . . , N , from the cipher, the empirical correlationsĉ i , i = 1, . . . , m, of the approximations (7) are calculated aŝ The distinguisher based on the method of Biryukov et al. [4] is given as the 2 -distance between the vectors c = (c 1 , . . . , c m ) andĉ = (ĉ 1 , . . . ,ĉ m ): Under the assumption that the approximations (7) are statistically independent, it was proved in [4] that the data complexity of the distinguisher (8) is inversely proportional toc This notion was defined in [4] and called as capacity of the set of linear approximations. This result means a significant improvement in data complexity when compared to one-dimensional method, but relies on the assumption that the approximations are statistically independent. Later in [29], this assumption was questioned and shown that it does not hold in general. The experiments in [20] on reduced round Serpent confirmed this. Moreover, verifying statistical independence of a set of linear approximations is essentially equally hard as computing correlations of all their linear combinations. Indeed, given independence, all correlations can be computed given the correlations of a set of base approximations using the piling-up lemma. And given the correlations, statistical independence can be proved by the inverse of piling-up lemma. A truly multidimensional approach to the distinguishing problem was given in [1] based on the LLR-statistic (3). By Corollary 2.1, the data complexity of a multidimensional linear distinguisher with p.d. p is inversely proportional to C( p). By Corollary 2.2, we see that C( p) ≥c 2 . This indicates that data complexity for the distinguisher (8) with statistically independent linear approximations is larger than for the multidimensional LLR distinguisher. Moreover, using LLR and multidimensional p.d., we do not need to assume that the used approximations are statistically independent.

Linear Approximation of a Block Cipher
We will be applying the statistical methods of Sect. 2.1 in the case where M = 2 m − 1 from now on. Biryukov et al. proposed in Sect. 3.4. in [4] to extend the set of m linearly independent approximations to m , m ≤ m ≤ M, approximations by including linear combinations of the base approximations with non-negligible correlations. It was argued that the same rule that the data complexity is inversely proportional applies also for the larger set of m linearly and statistically dependent approximations. If this holds, then by Corollary 2.2 Biryukov's distinguisher should converge to the multidimensional distinguisher by adding all linear combinations of the base approximations such that m = M. However, since the statistic (8) is based on the assumption of statistical independence, its data complexity cannot be determined accurately.
An optimal case would be where all the M linear approximations used for a multidimensional distinguisher have equally large correlations (in absolute value). Then using all M of them might allow reducing data complexity. On the other hand, if only a single one-dimensional approximation from a set of M approximations has a large correlation, then it is not useful to include the others in the distinguisher. In reality, all cases between these too extremes occur.
In this paper, we will present a number of methods for generalising the key recovery attacks based on Matsui's Alg. 1 and Alg. 2 to multiple dimensions using the statistical framework of multidimensional distinguishers. For each method, we will derive an explicit estimate of data complexity and show that while the data complexity decreases as the capacity of the linear approximations increases, the mere number of the linear approximations can have an opposite effect which varies depending on the method.
We have chosen the block cipher Serpent [3] as the test bed for our methods. Unlike the block cipher DES, it has been designed to resist linear cryptanalysis and therefore does not have individual strong linear approximations. Still it has many linear approx-imations with correlations significantly stronger than the average that could potentially be combined and used in a multidimensional linear attack. Such linear approximations for Serpent have been previously found and used in experiments for Biryukov's multiple linear cryptanalysis by Collard et al. [10].
Let us study an iterated block cipher with block size n. Let x be the plaintext and y the output of the cipher after R rounds. The cipher key is expanded using a key-scheduling algorithm to a sequence of round keys. We denote by K the expanded key, that is, the vector consisting of all round key bits used in R rounds, and by h the length of K . Then a block cipher is a vector Boolean function with input (x, K ) ∈ F n 2 × F h 2 . By (6), an m-dimensional linear approximation of the block cipher can be considered as a vector Boolean function where U and W are m × n binary matrices. The matrix V has also m rows and it divides the expanded keys into 2 m equivalence classes g = V K , g ∈ F m 2 . Let us denote by p(η|K ) the probabilities of the m-bit values η of (10) for a fixed key K . In the analysis of this paper, it is assumed that for each η ∈ F m 2 the probabilities p(η|K ) are (about) equal for all K . We denote by p be the p.d. of (10) for a fixed key K , for all K ∈ F h 2 . Then for each g ∈ F m 2 , the data Ux t ⊕ Wŷ t , t = 1, . . . , N , are drawn from p.d. p g , a fixed permutation of p determined by g, and all p.d. p g , g ∈ F m 2 , are each other's permutations. In particular, p g η⊕k = p g⊕k η , for all g, η, k ∈ F m 2 . By symmetry, the following results apply: and min g,g =g where the minimum Kullback-Leibler distance and capacity will be denoted by D min ( p) and C min ( p), respectively. Moreover, if D min ( p) = 0, we need to join the corresponding key classes to one class such that we may assume D min ( p) = 0 and C min ( p) = 0.

Advantage in Key Ranking
In a key recovery attack, a set of candidate keys is given, and the problem is to determine which candidate is the right one. Let the keys be searched from a set F l 2 of all 2 l strings of l bits. The algorithm consists of four phases: the counting phase, analysis phase, sorting phase and searching phase [34]. In the counting phase, data, for example plaintextciphertext pairs, are collected from the cipher. In the analysis phase, a real-valued statistic T is used in calculating a mark T (κ) for all candidates κ ∈ F l 2 . In the sorting phase, the candidates κ are sorted, i.e. ranked, according to their marks T (κ). Optimally, the right key, denoted by κ 0 , should be at the top of the list. If this is not the case, then in the search phase the candidates in the list are tested until κ 0 is found. The goal of this paper is to find a ranking statistic T that is easy to compute and that is also reliable and efficient in finding the right key.
The time complexity of the search phase, given amount N of data, was measured using a special purpose quantity "gain" in [4]. A similar but more generally applicable concept of "advantage" was introduced by Selçuk [33], where it was defined as follows: Given a data sample obtained using an unknown fixed key, a key recovery attack for an l-bit key is said to achieve an advantage of a bits over exhaustive key search, if the correct key is ranked among the top r = 2 l−a out of all 2 l key candidates.
Statistical tests for key recovery attacks are based on the wrong key hypothesis [17]. We state it as follows: Assumption 1. (Wrong key hypothesis) There are two p.d.'s q and q , q = q , such that for the right key κ 0 , the data are drawn from q and for any wrong candidate key κ = κ 0 the data are drawn from q = q.
A real-valued statistic T is computed from q and q , where one of these p.d.'s may be unknown, and the purpose of a statistic T is to distinguish between q and q . We , with parameters μ R and σ R , as this is the case with all statistics in this paper. Then μ R and σ R are determined with the help of linear cryptanalysis. We denote by D W the p.d. known based on the wrong key hypothesis such that Ranking the candidates κ according to T means rearranging the 2 l random variables T (κ), κ ∈ F l 2 , in decreasing (or sometimes increasing) order of magnitude. Writing the ordered random variables as T (0) ≥ T (1) ≥ · · · ≥ T (i) ≥, we call T (i) the ith order statistic. Let us fix the advantage a such that the right key should be among the r = 2 l−a highest ranking key candidates. Hence, the right key should be at least as high as the r th wrong key with mark T (r ) . If the random variables T (κ), κ = κ 0 , are statistically independent, then by Theorem 1 in [33] the random variable T (r ) is distributed as Let us define the success probability P S of having κ 0 among the r highest ranking candidates. If all the random variables T (κ), κ ∈ F n 2 are statistically independent, we have since As the data complexity N depends on the parameters μ R − μ a and σ 2 R + σ 2 a , we can solve N from (14) as a function of a and vice versa. Hence, (14) describes the trade-off between the data complexity N and the time complexity of the search phase.

Finding the Right Key Class
Assume that we have found an m-dimensional approximation (10) with non-uniform p.d. p. The output y is the ciphertext obtained from the cipher by encrypting plaintext x over R rounds using the key K . We denote by g 0 the right key class to which K belongs. Our goal is to find g 0 .
To recall how Alg. 1 works for m = 1, let us denote by c (ĉ) the theoretical (empirical) correlation of u · x ⊕ w · y. The decision in Alg. 1 in one dimension is based on the following test: the key class candidate bit g is chosen to be 0 if cĉ > 0. Otherwise, g = 1. In other words, the statistical decision problem is to determine which of the two gives the best fit with the data. In multiple dimensions, the same data will be used for ranking the different candidate classes g ∈ F m 2 . Hence, the corresponding marks T (g) will be statistically dependent. We are not aware of any general method of calculating the c.d.f. of the order statistic T (r ) of statistically dependent random variables. The asymptotic c.d.f. of the maximum of normal, identically distributed but dependent random variables for large M = 2 m − 1, m ≥ 7 is derived in Sect. 9.3. in [14]. However, the problem still remains as the random variables T (g 0 ) and max g =g 0 T (g) are statistically dependent.
Denote by N (g) the data complexity of ranking g with advantage a, if g is the right key. We will assume T (g)'s to be statistically independent to determine N (g).
The assumption of statistical independence of T (g)'s could be avoided by drawing g∈V m N (g) ≈ 2 m max g N (g) words of data, as then the right key class g 0 would be ranked with advantage a and each mark T (g), g ∈ F m 2 could be calculated from different data. However, the resulting complexity estimate would be far too large to be of practical value.
We will study three different ways to generalise the one-dimensional Alg. 1 to multiple dimensions. Since the data are drawn i.i.d. from the p.d. p g 0 and not from any other p.d. p g , g = g 0 , we can interpret the problem of finding g 0 as a generalisation of the goodness-of-fit test where we determine whether given data are drawn from p.d. p g or not. The candidate key class g ∈ F m 2 which is most strongly indicated by this test to fit the data is chosen to be the right key class. The classical goodness-of-fit tests are the χ 2 -test and the G test based on the Kullback-Leibler distance. The first two methods to be presented in this paper are generalisations of these tests into the case of multiple distributions, i.e. finding one distribution from a set of distributions. The χ 2 -method based on the χ 2 -test and the log-likelihood method based on the G test are studied in Sects. 5.2 and 5.3, respectively.
Our third method is the LLR-method to be studied in Sect. 5.4. In [2], the problem of distinguishing one known p.d. from a set of other p.d.'s was studied. It was then possible to use the optimal distinguisher, the LLR-statistic, in solving the problem. However, since g 0 is unknown, we cannot apply the results of [2] in our work directly. Instead, we will use the following heuristic: since the data corresponding toq are drawn from the unknown p.d. p g 0 = θ , it should be easiest to distinguish the right p.d. p g 0 rather than any other p.d. p g , g = g 0 from the uniform distribution using the LLR-statistic. Hence, the candidate key class g ∈ F m 2 that gives the strongest distinguisher between the corresponding p.d p g and θ is chosen to be the right key. This can be seen as a variant of Assumption 1.
In all our analysis, it is assumed that p g and p g , g = g , are close to each other, and all these distributions p g are close to θ in the sense of Property 1. The assumption about closeness must be verified with practical experiments.

χ 2 -Method
The mark for each candidate class g ∈ F m 2 based on χ 2 -statistic is defined as follows: where N is the amount of data used in the attack, with M = 2 m − 1 degrees of freedom. The empirical distributionq should be near to the correct p.d. p g 0 while being further away from all the other p.d.'s p g , g = g 0 . Hence, the candidate class corresponding to the smallest s(g) is chosen to be the right key class. By [15], the distribution of s(g) can be approximated by χ 2 M (N C( p g , p g 0 )) which can further be approximated by s(g) ∼ N (μ g , σ 2 g ), with mean μ g = M + N C( p g , p g 0 ) and variance σ 2 g = 2(M + 2N C( p g , p g 0 )), provided that μ g > 30 [13], i.e. m should be at least 5.
For simplicity, we have only determined the data complexity of full advantage of a = m bits and assumed that s(g)'s are statistically independent. That is, we have determined the data complexity of ranking g 0 on the top of the list with s(g). Let P S be the probability of finding g 0 such that Using Property (11), we can do calculations and approximations similar to those done in Sect. 4 in [1] or in the proof of Theorem 2 in [20] and get the following estimate of the data complexity N that would be sufficient for finding g 0 with success probability where β S = ln( √ 2π ln P −1 S ). Note the exponential dependence of N on m as M = 2 m − 1. Our experiments showed, however, that much less data are needed than what is predicted by (16). The reason may be that the statistical dependence of the marks s(g) strengthens the method. However, as noted previously, we could not find a way to do the derivations without the assumption of statistical independence. The formula for the advantage of the χ 2 -method could also be derived but we have omitted it here, since it will be given for the LLR-method studied in Sect. 5.4 which performs better in terms of data complexity.

The Log-Likelihood Method
Another popular goodness-of-fit test is the log-likelihood test, also known as G test. The experiments on Alg. 1 done in [20] used this test. The mark based on G test uses the Kullback-Leibler distance between the empirical p.d.q and the theoretical p.d. p g . In [15], it is shown that for each candidate class g ∈ F m 2 the random variable G(g) can be approximated to be distributed as Since p g are near to p g 0 , the parameters δ g ≈ N C( p g , p g 0 ) and ξ g ≈ 0 and the G test performs similarly as the χ 2 -test [15].

Log-Likelihood Ratio Method
The log-likelihood ratio is the optimal statistic for distinguishing two distributions [11]. It is also asymptotically normal as stated in Proposition 2.1. Hence, we would like to use it for key ranking. The idea is that the empirical distribution can be used for distinguishing the p.d. p g 0 related to the correct key class from the uniform p.d. with large LLR value, while any wrong p.d. p g , g = g 0 is less distinguishable from θ . For each g ∈ F m 2 , we compute the mark (g) = LLR(q, p g , θ).
We choose the candidate class g with largest (g) to be the right key class. We cannot apply [2] here to determine the data complexity of finding g 0 as the result would be too optimistic. The task is to distinguish an unknown p g 0 from a set of p.d.'s { p g | g ∈ F m 2 }. In [2], one distinguishes only p g 0 from p g 0 , g 0 = g 0 , the p.d. closest to p g 0 in Kullback-Leibler distance. We have to consider all the other candidate classes as well, which increases the data complexity. Applying the theory of key ranking described in Sect. 4, we derive the following result.
Theorem 5.1. Assume that the random variables (g) are statistically independent and that (g) ∼ N (μ W , σ 2 W ), g = g 0 where μ W = 0 and σ 2 W ≈ σ 2 R , where σ 2 R is the variance of the random variable (g 0 ). If the p.d.'s p g , g ∈ F m 2 and θ are close to each other in the sense of Property 1, the advantage a of the LLR-method with marks calculated by (17) can be approximated by where P S ( ≥ 0.5) is the probability of success, N is the amount of data used in the attack and C( p) and m are the capacity and the dimension of the linear approximation (10), respectively.
The assumption of statistical independence of (g)'s was discussed in Sect. 5.1. The assumption about normal distribution of the wrong keys (g), g = g 0 is based on the law of large numbers [11]. The approximation of variance σ 2 W = σ 2 R is commonly used, for example, in [34], and the approximation of mean is based on the idea that the empirical data are not closer to any p g , g = g 0 than θ implying that μ W ≤ 0. In the worst case, with largest data complexity, μ W = 0.
Proof. Let us proceed first by finding the p.d.'s for the random variables (g), g ∈ F m 2 . By Proposition 2.1 and property (11) . By the assumptions, we may use (13) .
By approximation (1), On the other hand and by (18) we have Since a ≤ m, and b > 1, we have showed that σ 2 a /σ 2 W 1 and also σ 2 a σ 2 R . Then from which we can solve an estimate of N as a function of a to be where the last approximation is obtained using approximation (1). It By inversion, we get an estimate of a as a function of N as desired.
The experimental advantages for the different methods are studied in the next section. A common choice for P S is 0.5 ≤ P S ≤ 0.99. Hence, the value of Φ −1 (P S ) is typically a small positive number less than 3, and if m ≥ Φ −1 (P S ) 2 , the numerator of (19) is bounded above by 16m. This shows that the dependence of the data complexity on the dimension m of the multidimensional linear approximation is linear for the LLR-method, while it is exponential in (16) for the χ 2 -method. Since in practice C( p) ≈ C min ( p), the comparison of (19) and (16) indicates that the LLR-method is more efficient than the χ 2 -method or the log-likelihood method.

Algorithms and Complexities
For comparing the two methods, LLR and χ 2 , we are interested in the complexities of the first two phases of Alg. 1 since the sorting and searching phase with fixed advantage a do not depend on the chosen method. The counting phase is done online and all the other phases can be done offline. However, we have not followed this division [34] in our implementation, as we do part of the analysis phase online. We will divide the algorithm into two phases as follows: In the online phase, depicted in Fig. 1, we calculate the empirical p.d.q. The marks s(g) for the χ 2 -method and (g) for the LLR-method are then assigned to the keys in the offline phase. The offline phases for χ 2 -method and LLR-method are depicted in Figs. 2 and 3, respectively.  The data complexities of the online phase for χ 2 and LLR are given by (16) and (19), respectively. The dependence of the data complexity of the χ 2 -method on the advantage is similar to the LLR-method. The main difference is in the predicted dependence on m. The time complexity of the online phase is N m, where N is the data complexity. The memory usage is 2 m , the size of the empirical distribution.
In the offline phase, the time and memory complexities for both methods are 2 m and 2 2m , respectively.

Experiments on Four-Round Serpent
Similarly as in previous experiments on multiple linear cryptanalysis, see [10], the Serpent block cipher was used as a test bed. Description of Serpent can be found in [3]. We tested the different methods for multidimensional Alg. 1 described in this paper on four-round Serpent, 4th-7th rounds, by selecting linearly independent one-dimensional base approximations u i · x ⊕ w · y, i = 1, . . . , m, to construct a linear approximation of the form (10) with m = 7 and m = 10. The used masks w and u i , i = 1, . . . , m, can be found in [20].
We checked the assumption about closeness of the hypothetical distributions p g and θ and saw that it holds as | p g η − p g η | < 1 150 p g η , for all g, g and η ∈ F m 2 . We also checked that C min ( p) = 0 and actually, C min ( p) ≈ C( p).
The experiments showed that the empirical advantage when ranking the key classes was exactly the same for all multidimensional methods. Hence, we only depict the LLRmethod. In particular, all methods were equally efficient in determining the correct key class. Equations (19) and (16) predict that the LLR-method should be the most efficient: when m increases, the data requirement of χ 2 -based tests increases exponentially with m, whereas the increase is linear for the LLR-method. It is possible that the variance  of the χ 2 -method is not as large as the theory predicts, or the statistical dependence of random variables s(g) strengthens the χ 2 -method more than expected.
The statistical model of the relationship between the advantage a and data complexity N derived in this paper was tested in experiments. The results are given in Figs. 4 and 5. The empirical advantage is determined by averaging the advantages obtained for 16 randomly selected keys using the LLR-method. The theoretical advantage given by Theorem 5.1 is depicted for three different values of P S . The experiments show that the data complexity predicted this way is too high, but we believe it gives a good upper bound.
As discussed after the proof of Theorem 5.1, for full advantage a = m, we can in usual cases approximate N ≈ 16m/C( p). By increasing m, that is, using more linear approximations, the data complexity N decreases as long as the ratio C( p)/m increases. This sets an upper bound for m to be used in practice. In case of four-round Serpent, the practical upper bound around m = 12 was found in experiments.
In both cases m = 7 and m = 10, we also show how much better the m-dimensional LLR-method is compared to the binomial method where the same set of m one-dimensional approximations and Matsui's Alg. 1 is used to determine each key class candidate bit separately and independently.  [10,20] for Serpent already confirmed that this is a favourable thing to do in spite of the lack of theoretical justification. In [20], also the enhanced method of Biryukov et al. and the full multidimensional log-likelihood method were compared in experiments on Serpent and the latter was shown to be more powerful.

Algorithm 2
Let us study a cipher with R + 1 rounds. Let x be the plaintext and z be the cipher text after R + 1 rounds. Let the (R + 1)th round function and round key be f and k ∈ F l 2 , respectively. Then the output after R rounds is y = f −1 (z, k). Algorithm 2 uses a multidimensional linear approximation over R rounds given by (10) with p.d. p. The task is to find the right last round key k 0 and possibly, in addition, the right key class g 0 for the key used in the first R rounds.

Statistical Setting for Alg. 2
In the counting phase, we draw N data pairs (x t ,ẑ t ), t = 1, . . . , N . In the analysis phase, for each last round key candidate k, we first calculateŷ k t = f −1 (ẑ t , k), t = 1, . . . , N . Then, for each k, we calculate the empirical p.d.q k = (q k 0 , . . . ,q k M ), wherê The keys are then given mark T (k) by using some statistic T that is calculated from different data Ux t ⊕ Wŷ k t , t = 1, . . . , N , for each key candidate k ∈ F l 2 . Hence, the random variables T (k) are statistically independent.
If we use the wrong key k = k 0 to decrypt the ciphertext, it means we essentially encrypt over one more round and the resulting data will be more uniformly distributed. This heuristics is behind the original wrong key randomisation hypothesis [22,27], which in our case means that the data Ux t ⊕ Wŷ k t , t = 1, . . . , N , k = k 0 are drawn i.i.d. from the uniform distribution. On the other hand, when decrypting with the correct key k 0 the data Ux t ⊕ Wŷ k 0 t , t = 1, . . . , N are drawn i.i.d. from p g 0 , where g 0 ∈ F m 2 is unknown.

Key Ranking in One-dimensional Alg. 2
Key ranking and advantage in the one-dimensional case, m = 1, of Alg. 2 was studied in [33]. We will present it here briefly for completeness. Let c > 0 be the correlation of (10) (the calculations are similar if c < 0) and letĉ k be the empirical correlation calculated from the data. The mark used in ranking the keys is then given by s (k) = |ĉ k |.
The random variableĉ k 0 is binomially distributed with mean μ R = c and variance The wrong key random variablesĉ k , k = k 0 , are binomially distributed with mean μ W = 0 (following Assumption 1) and variance σ 2 where FN is the folded normal distribution; see Appendix A in [33]. Then for finding g 0 with given success probability P S and advantage a the data complexity N is estimated as follows

Different Scenarios in Multiple Dimensions
When considering generalisation of Alg. 2 to the case, where multidimensional linear approximation (10) is used, basically two different standard statistical methods can be used: • Goodness of fit (usually based on χ 2 -statistic, see also [29,34]) and • Distinguishing of an unknown p.d. from a given set of p.d.'s (usually based on LLR-statistic) The goodness-of-fit approach is a straightforward generalisation of one-dimensional Alg. 2. It can be used in searching for the last round key. It measures whether the data are drawn from the uniform distribution, or not, by measuring the deviation from the uniform distribution. It ranks highest the key candidate whose empirical distribution is farthest away from the uniform distribution. The method does not depend on the key class candidate for the key used before the last round. Information about p.d. p is required only for measuring the strength of the test. We will study this method in Sect. 6.4. After the right last round key k 0 is found, the data derived by Alg. 2 can then be used in any form of Alg. 1 for finding the right key class g 0 of the key used in the R first rounds. In this manner, the χ 2 -approach allows separating between Alg. 1 and Alg. 2.
The idea behind the goodness-of-fit approach is that the "wrong key" distributionŝ q k , k = k 0 are close to the uniform distribution, whereasq k 0 is further away from it. Moreover, the expected value of C(q k 0 , θ) = C(q k 0 ) should be approximately C( p).
The LLR-method uses the information about the p.d. related to the key class candidate g also in Alg. 2. In this sense, it is similar to the method presented in [4], where both g 0 and k 0 were searched at the same time. As we noted in Sect. 2.1, the LLR-statistic is the optimal distinguisher between two known p.d.'s. If we knew the right key class g 0 , we could simply use the empirical p.d.'sq k for distinguishing p g 0 and the uniform distribution and then choose the last round candidate k for which this distinguisher is strongest [1]. In practice, the right key class g 0 is unknown when running Alg. 2 for finding the last round key.
Our approach is the following. In [2], it was described how one can use LLR to distinguish one known p.d. from a set of p.d.'s. We will use this distinguisher for distinguishing θ from the given set p g , g ∈ F m 2 . In the setting of Alg. 2, we can expect that for the right k 0 , it should be possible to clearly conclude that the data (x t ,ẑ t ), t = 1, . . . , N , yield data Ux t ⊕ Wŷ k 0 t , t = 1, . . . , N , which follow a p.d. p g , for some g ∈ F m 2 , rather than the uniform distribution. On the other hand, for the wrong k = k 0 , the data follow the uniform distribution, rather than any p g , g ∈ F m 2 . To distinguish k 0 from the wrong key candidates, we determine, for each round key candidate k, the key class candidate g, for which the LLR-statistic is the largest with the given data. The right key k 0 is expected to have g 0 such that the LLR-statistic with this pair (k 0 , g 0 ) is larger than for any other pair (k, g) = (k 0 , g 0 ). In this manner, we also recover g 0 in addition to k 0 . The LLR-method is studied in Sect. 6.5.

The χ 2 -Method
This method separates Alg. 1. and Alg. 2 such that the latter does not need any information of p. Both methods are interpreted as goodness-of-fit problems, for which the natural choice of ranking statistic is χ 2 . We will show first how to find the last round key k with Alg. 2.
6.4.1. Algorithm 2 with χ 2 Given empirical p.d.q k , the mark based on χ 2 -statistic is calculated from the data by where M = 2 m −1 is the number of degrees of freedom. Formula (21) can be interpreted as the 2 -distance between the empirical p.d. and the uniform p.d.. By Assumption 1, the right round key should produce data that are farthest away from the uniform p.d. and the last round key candidate k for which the mark S(k) is largest is chosen. Obviously, (21) simplifies to (ĉ k ) 2 , if m = 1.
The following theorem describes the relationship between the data complexity and the time complexity of the search phase. Theorem 6.1. Suppose the cipher satisfies Assumption 1 where q = θ and the p.d.'s p g , g ∈ F m 2 and θ are close to each other. Then the advantage of the χ 2 -method using marks (21) is given by where M = 2 m − 1, P S (> 0.5) is the probability of success, N is the amount of data used in the attack and C( p) and m (≥ 5) are the capacity and the dimension of the linear approximation (10), respectively.
Proof. According to [15], the random variable S(k 0 ) is distributed approximately as because of the symmetry property (11). Since m ≥ 5, we may approximate the distribution by a normal distribution with μ R = M + N C( p) and σ 2 R = 2(M + 2N C( p)). The parameters do not depend on g 0 or k 0 . For the wrong keys k = k 0 , we obtain by [15] that The mean and variance in (13) are When a < l is large, we use approximation (20) to get a ≈ b 2 by which μ a = √ 2a M + M. Further by approximation (1), we get φ(b) ≈ b2 −a . Then For small a, this estimate holds trivially. Then to proceed, we restrict to the case N C( p) < M/4. This is not essential restriction, since finally N C( p) will be close to a small constant multiple of √ M. Then we have and we use the upper bound as an estimate for σ 2 a + σ 2 R . Then we can solve data complexity from the formula (14) of the success probability to obtain the following estimate of the data complexity By solving for a, we get the estimated advantage as claimed.
Note that the normal approximation of the wrong key distribution is valid only when m > 5, that is, when the approximation of χ 2 -distribution by a normal distribution is valid. It is not possible to perform the theoretical calculations for small m as the χ 2distribution does not have a simple asymptotic form in that case and we cannot determine f W and F W in (13). Since for m = 1 our χ 2 -method reduces to the square of s(k) that was used by Selçuk, the theoretical distributions differ from our calculations and we get a slightly different formula for the advantage. Despite this difference, the methods are essentially equivalent for m = 1.
Keeping the capacity and advantage constant, we see by (23) that the data complexity increases exponentially as 2 m/2 as the dimension m of the linear approximation increases and is sufficiently large. Hence, in order to strengthen the attack by increasing m, the capacity should increase faster than 2 m/2 . This is a very strong condition and it suggests that in applications, only approximations with small m should be used with χ 2 -attack. The experimental results for different m presented in Sect. 6.7 as well as the theoretical curves depicted in Fig. 9 suggest that increasing m in the χ 2 -method does not necessarily mean improved performance for Alg. 2.
While (22) and (23) depend on the theoretical distribution p, the actual χ 2 -mark (21) is independent of p. Therefore, to realise the attack, we do not need to know p accurately. It suffices to find an approximation (10) that deviates as much as possible from the uniform distribution. On the other hand, if we use time and effort for computing an approximation of the theoretical p.d. p and if we may assume that the approximation is accurate, we would also like to exploit this knowledge for finding the right class of the key used in R first rounds using Alg. 1. As noted in Sect. 5, there are several ways to realising a multidimensional Alg. 1. Considering Alg. 1 as a χ 2 -based goodness-of-fit problem, we may combine Alg. 1 and Alg. 2 to be described next.

Combined Method and Discussion
The sums of squares of correlations used in the method of Biryukov et al. [4] are closely related to the sums of squares (15) and (21). Indeed, we could define a combined χ 2 -mark S by considering the sum of (15) and (21) and setting where s(k, g) is calculated from the empirical p.d.q k , k ∈ F l 2 by (15). The right key (k 0 , g 0 ) should give the smallest mark. If we approximate the denominators in (15) by 2 −m and scale by 2 m N , we obtain from S (k, g) the mark which is closely related to the one used in [4]: whereĉ k is the empirical correlation vector corresponding to the last round key candidate k and c g is the theoretical correlation vector corresponding to the key class candidate g. Indeed, if in (25) all correlation vectorsĉ k and c g contain correlations from all linear approximations then (25) becomes the same as 2 m B(k, g) as can be seen by Proposition 2.2 and Parseval's theorem. Initially, in the theoretical derivation of (25) only linearly and statistically independent approximations were included in the correlation vectors. However, in Sect. 3.4 of [4] it was proposed to take into account all linear approximations with strong correlations when forming the mark (25) in practice. In practical experiments by Collard et al. [10], this heuristic enhancement was demonstrated to improve the results for Serpent. In this paper, we have shown how to remove the assumption about independence of the linear approximations and offer the option of including all linear approximations that have sufficient contribution to the capacity.
Other possibilities for combining Alg. 1 and Alg. 2 based on χ 2 or its variants are also possible, with different weights on the terms of the sum in (24), for instance. However, the mathematically more straightforward way is to use the pure χ 2 -method defined by (21) and (15), as its statistical behaviour is well known. An even more efficient method can be developed based on LLR as will be shown next.

The LLR-Method
This method is also based on the same heuristic as the wrong key hypothesis: For k = k 0 , the distribution of the data should look uniform and for k 0 it should look like p g 0 , for some g 0 . Hence, for each k, the problem is to distinguish the uniform distribution from the discrete and known set p g , g ∈ F m 2 . Let us use the notation L(k, g) = LLR(q k , p g , θ). We propose to use the following mark Now k 0 should be the key for which this maximum over g is the largest, and ideally, then the maximum is achieved for g = g 0 . While the symmetry property (11) allows one to develop statistical theory without knowing g 0 , in practice one must search through F l 2 for k 0 and F m 2 for g 0 even if we are only interested in determining k 0 . Theorem 6.2 gives the trade-off between the time complexity of the search phase and the data complexity of the algorithm. Theorem 6.2. Suppose the cipher satisfies Assumption 1 where q = θ and the p.d.'s p g , g ∈ F m 2 and θ are close to each other. Suppose further that for all the wrong last round key candidates k = k 0 , the 2 m random variables L(k, g), g ∈ F m 2 are statistically independent. Then the advantage of the LLR-method for finding the last round key k 0 , assuming that it is paired with the right key class g 0 , is given by Here N is the amount of data used in the attack, P 12 (> 0.5) is the probability of success and C( p) and m are the capacity and the dimension of the linear approximation (10), respectively.
For fixed k = k 0 , the random variables L(k, g), g ∈ F m 2 are calculated using the same dataq k . Hence, they are actually statistically dependent. Similarly as in Sect. 5, we will assume that they are statistically independent to simplify calculations. The practical experiments concentrated on finding the right last round key k 0 such that we did not actually need the knowledge of g 0 or its data complexity. The following calculations give a theoretical model that can be used in describing how the LLR-method behaves especially when compared to other methods.
Proof. Using Proposition 2.1 and property (11), we can state Assumption 1 as follows: For the right key (k 0 , g 0 ), and for k = k 0 and any g ∈ F m 2 , Hence, μ R , σ 2 R , μ W and σ 2 W do not depend on g ∈ F m 2 . For fixed k = k 0 , the random variables L(k, g) for k = k 0 are identically normally distributed with mean μ W and variance σ 2 W . Assuming that for each k = k 0 , the random variables L(k, g)'s are independent, we obtain that the c.d.f. of their maximum is given by [14] F Let us fix the advantage a such that r = 2 l−a . The mean μ a of the r th wrong key statistic L (r ) can now be calculated from (13) to be and the variance is Let be the probability that we rank k 0 among the r highest ranking keys provided that we pair g 0 with k 0 . We can calculate P 12 using (14), (30) and (31) to obtain Hence, the amount of data is estimated to be sufficient. We approximate Then by replacing a by a + m in (20) we have a + m ≈ b 2 and can solve an estimate of advantage a as a function of N from (32) to get the claim.
be the probability that given k 0 , we choose g 0 , i.e. the probability of success of Alg. 1. with full advantage a = m. Let be the probability that we rank k 0 , paired with any g ∈ F m 2 , among the r highest ranking keys. Then P 2 = P 12 P 1 + Pr(L(k 0 ) > L (r ) | L(k 0 ) = LLR(k 0 , p g , θ), g = g 0 )(1 − P 1 ) ≥ P 12 P 1 .
Denote by N 1 and N 2 the data complexities needed to achieve success probabilities P 1 , and P 2 , respectively. The data complexity N 1 is given in (19). If we pair k 0 with g = g 0 , then L(k 0 ) ≥ L(k 0 , g 0 ) for a fixed empirical p.d.q k 0 , so that k 0 gets ranked higher than by using the correct g 0 . Hence, assuming that k 0 gets paired with g 0 only decreases P 2 so the theory predicts N 2 to be (slightly) larger than in reality. Then we can approximate N 2 from above by (32) and the corresponding advantage with success probability P 2 is approximately given by (27). The data complexity N 1 is an overestimate for the actual data complexity of Alg. 1 so in practice, N 2 > N 1 . Then if k 0 is ranked with advantage a and success probability P 2 > 0.5 among the 2 l−a highest ranking keys, it is also paired with the right key class g 0 . We have the following corollary: Corollary 6.1. Under the assumptions of Theorem 6.2, the data complexity of the LLRmethod for ranking the last round key k 0 among the r = 2 l−a highest ranking keys can be estimated as  and with this data complexity k 0 is paired with the right key class g 0 with a high success probability. On the other hand, given an amount N of data, the advantage of the LLRmethod is where P 2 (> 0.5) is the probability of success and C( p) and m are the capacity and the dimensions of the linear approximation (10), respectively.
With fixed N and capacity C( p), the advantage decreases linearly with m, whereas in (22), the logarithm of advantage decreases linearly with m. For fixed m and p, the advantage of the LLR-method seems to be larger than the advantage of the χ 2 -method. The experimental comparison of the methods is presented Sect. 6.7.

Algorithms and Complexities
Similarly as in Sect. 5.5, we divide Alg. 2 into online and offline phases. In the online phase, depicted in Fig. 6 we calculate the empirical p.d.'s for the round key candidates. The marks S(k) for the χ 2 -method and L(k) for the LLR-method are then assigned to the keys in the offline phase. The offline phases for χ 2 -method and LLR-method are depicted in Figs. 7 and 8, respectively. After the keys k are each given the mark, they can be ranked according to the mark. Givenq k 0 , the multidimensional Alg. 1 can be used for finding g 0 offline. The method used in Alg. 1 does not depend on the one used in Alg. 2, so, for example, the LLR ranking, see Fig. 3, can be used for finding g 0 .
For fixed advantage a, the data complexities of the online phase for the χ 2 and LLR are given in (23) and (33), respectively. Theoretical and practical results imply that the data complexities of different methods in Alg. 2 dominate the data complexity of Alg. 1 given in (19). Hence, the total data complexities are given by (23) and (33), even for deriving both k 0 and g 0 . The memory and time complexities for online phase are 2 l+m and 2 l m N , where N is the data complexity.
For the offline phase of LLR-based Alg. 2, the time and memory complexities are 2 l+m and 2 m max(2 l , 2 m ), respectively. The method obtained by a combination of χ 2based Alg. 2 and LLR-based Alg. 1 has the same complexities. Thus, we recommend using the LLR-method rather than χ 2 -method unless the accuracy of the p.d. p of the linear approximation (10) is a concern. If we use χ 2 for finding only the last round key k 0 , we have the same time complexity as for LLR but a reduced memory complexity 2 l , since we do not have to store the theoretical distributions p g . In some situations, it may even be advantageous to combine different methods. For example, we may first find, say, r best last round keys by χ 2 -based Alg. 2, such that the data complexity is given by (23), with advantage a = l − log r . Then we can proceed by applying Alg. 2 based on LLR-method in a reduced search space of size r < 2 l . Other similar variants are possible. Their usefulness depends on the cipher that is being analysed.

Experiments on Five-Round Serpent for Alg. 2
The purpose of the experiments was to test the accuracy of the derived statistical models and to demonstrate the better performance of the LLR-based method in practice. We take 5 rounds of Serpent, from the 4th to 8th round, m-dimensional linear approximation over four rounds, 4th to 7th, and searched for a 12-bit part of the round key used in the 8th round. Different m were used. Each experiment was performed for 16 different, randomly selected keys. Given a data sample of size N , we measured the experimental ranking of k 0 among all 2 12 candidates and computed the corresponding advantage. In the figures, the average advantage taken over 16 keys is depicted for each sample sizes N that are integer multiples of 2 21   We calculated the capacities for the approximation (10) over four-round Serpent for different m. The results in Sect. 5.6 showed that p g 's can be considered to be close to each other and θ.
The theoretical advantage of the χ 2 -method predicted in (22) has been plotted as a function of data complexity in Fig. 9. The figure shows that increasing m larger than 4, the attack is weakened. This suggests using m = 4 base approximations in the χ 2attack. Since we should have m at least 5 for the normal approximation of χ 2 M to hold, the theoretical calculations do not necessarily hold for small m. However, the experiments, presented in Fig. 10, are in accordance with the theory also for m = 1 and m = 4. The most efficient attack is obtained by using m = 4 equations. Increasing m (and hence the time and memory complexities of the attack) actually weakens the attack. The optimal choice of m depends on the cipher and the size of the correlations of available linear approximations.
The reason is the χ 2 -statistic itself: it only measures if the data follow a certain distribution, the uniform distribution in this case. The more the approximations we use, the larger the distributions become and the more the uncertainty we have about the "fitting" of the data. Small errors in experiments generate large errors in χ 2 as the fluctuations from the relative frequency 2 −m become more significant.
The theoretical advantage of the LLR-method (34) is plotted against the data complexity in Fig. 11 for different m. The empirical advantages for several different m are shown in Fig. 12 in the same experimental setting as was used for χ 2 . Unlike for χ 2 , we see that the method can be strengthened by increasing m, until the increase in the capacity C( p) becomes negligible compared to increase in m. For four-round Serpent, this happens when m ≈ 12. Experimental results presented in Figs. 10 an 12 confirm the theoretical prediction that the LLR-method is more powerful than the χ 2 -method.  Also the theoretical and empirical curves agree nicely. For example, the full advantage of 12 bits with m = 7 is achieved at log N = 27.5 for LLR, whereas for χ 2 -method log N = 28. Moreover, the LLR can be strengthened by increasing m. For m = 12, the empirical and theoretical data complexities are close to 2 26 .

Conclusions
We studied several ways of extending Matsui's Algorithms 1 and 2 to multiple dimensions. Using the advantage, we could compare the efficiency of the different methods in theory and in practice. The theory predicted that for both Alg. 1 and Alg. 2 the key ranking based on LLR is more efficient than goodness-of-fit tests using χ 2 (or G test). However, in practical experiments for Alg. 1, the methods seemed to perform equally well. It remains an open question how to avoid the assumption about statistical independence of the ranking values and make the theoretical prediction more accurate for Alg. 1.
For Alg. 2, both theoretical and practical results are in agreement. We showed that the χ 2 -based method is weaker than the LLR-based method. Hence, we recommend to use the LLR-method proposed in this paper rather than the χ 2 -method as long as the estimated p.d. of the multidimensional linear approximation can be assumed to be accurate for all keys. If the theoretical p.d. varies a lot with keys or otherwise cannot be determined accurately, then the χ 2 -method is preferred.

Later Developments
Let us conclude this revised version of the paper by giving a brief overview of the most important advances in multidimensional linear cryptanalysis during the decade that passed since the submission of this paper. At CT-RSA 2010, the authors presented an improved method for Alg. 1 called as the convolution method [19] and an attack on block cipher PRESENT [9]. The dependence of p.d. p of the key we observed in this paper has been considered, and more advanced statistical models have been developed. The variance due to random key was first integrated in the wrong key model of one-dimensional linear cryptanalysis [8], then to the model of the cipher [21] and finally for both [6] considering also explicitly the case of sampling without replacement. Since 2012, the hypothesis testing model has been adopted in many works on statistical cryptanalysis to determine distinguishing complexity [7] and data complexity for key recovery. Instead of relating the advantage a to ranking of the candidate keys, the quantity 2 −a is interpreted as the probability of the error of accepting a wrong candidate key. With the hypothesis testing model, the problem with statistically dependent order marks we had in this paper disappears, while the data complexity estimates remain essentially the same. Zero-correlation multidimensional linear cryptanalysis was invented in [7], where it was also shown to be linked with integral attacks. The corresponding mathematical link between general multidimensional linear cryptanalysis and truncated differential cryptanalysis was proved in [5]. It allows to transform differential-type attacks to lineartype attacks and vice versa. For example, the multidimensional linear attack presented in [9] gives also the best differential-type attack on this cipher. Finally, let us mention that the recently presented affine multidimensional linear cryptanalysis allows to remove all trivial linear approximations, e.g. where either the input mask or the output mask is equal to zero without introducing artificial independence assumptions [30]. The 2 m − 1 linear approximations over four rounds of SERPENT we used in this paper share the same output mask. Using the affine method, the trivial linear approximations (linear combinations of even number of base approximations) can be removed. The set of the remaining 2 m−1 linear approximations has the same capacity as the original set, but the variance of the related χ 2 -statistic is significantly reduced, thus potentially improving the power of the attack.