On the construction of ultra-light MDS matrices
Yu Tian Xiutao Feng Guangrong Li


Article Title]On the construction of ultra-light MDS matrices[This research was supported by National Key Research and Development Project under Grant No. 2018YFA0704705 and CAS Project for Young Scientists in Basic Research (Grant No. YSBR-035).]










1]Yu Tian

[2]Xiutao Fengfengxt@amss.ac.cn


1]Guangrong Li

[1]College of Artificial Intelligence, Nanning Vocational and Technical University  Naning,   China

*[2]Key Laboratory of Mathematics Mechanization, Academy of Mathematics and Systems Science, Chinese Academy of Sciences,
 Beijing,   China







In recent years, the Substitution-Permutation Network has emerged as a crucial structure for constructing symmetric key ciphers. Composed primarily of linear matrices and nonlinear S-boxes, it offers a robust foundation for cryptographic security. Among the various metrics used to assess the cryptographic properties of linear matrices, the branch number stands out as a particularly important index. Matrices with an optimal branch number are referred to as MDS matrices and are highly prized in the field of cryptography. In this paper we delve into the construction of lightweight MDS matrices. We commence implementation trees of MDS matrices, which is a vital tool for understanding and manipulating their implementations, and then present an algorithm that efficiently enumerates all the lightest MDS matrices based on the word representation. As results, we obtain a series of ultra-lightweight 4× 4 MDS matrices, remarkably, 4-bit input MDS matrices with 35 XOR operations and 8-bit input ones with 67 XOR operations . These matrices represent the most comprehensive lightweight MDS matrices available to date. Furthermore, we craft some involution 4× 4 MDS matrices with a mere 68 XOR gates.To our best knowledge, they are the best up to date. In the realm of higher-order MDS matrices, we have successfully constructed 5× 5 and 6× 6 matrices with 114 and 148 XOR gates respectively. These findings outperform the current state-of-the-art.



















[
    [
    September 5, 2024
=====================





§ INTRODUCTION


The Substitution-Permutation Network (SPN, for brevity) has emerged as a pivotal component in the construction of symmetric ciphers, particularly in recent years. Due to the prevalence of the block cipher AES <cit.>, this structure is extensively studied and typically comprises three primary elements: a key schedule, a small-sized (usually 4- or 8-bits) nonlinear function known as the S-box, and a larger (often 32-, 64- or 128-bits) linear function. The S-box serves to mix the bits within a 4- or 8-bit word, while the linear function mixes words together.

According to the wide trail design strategy <cit.>, the resilience of SPN ciphers against classical attacks, particularly differential and linear attacks, can be assessed by examining their constituent components. Specifically, for a nonlinear S-box, a high nonlinearity and low differential uniformity are desirable; S-boxes achieving optimal differential uniformity are termed Almost Perfect Nonlinear (APN). As for the linear function, a high branch number is expected to establish strong diffusion between input and output words; matrices corresponding to linear functions with an optimal branch number are known as MDS matrices (related to Maximum Distance Separable codes). Besides their use in SPN ciphers like AES, MDS matrices also find application in Feistel ciphers (e.g., Camellia <cit.>, SMS4 <cit.>), hash functions (e.g., Whirlpool <cit.>, Grostl <cit.>), and stream ciphers (e.g., ZUC <cit.>).

With the increasing miniaturization of electronic devices handling sensitive data, there is a growing need for novel cryptographic primitives that offer low implementation costs. In contrast to standardized, robust primitives like AES, lighter cryptographic solutions are often preferred in resource-constrained environments. To reduce the cost of the SPN structure, one approach is to optimize its main components: the S-box and the linear function. Numerous studies have explored lightweight alternatives for these components, as evidenced by works on S-boxes <cit.> and linear functions <cit.>. These findings have facilitated the design of various new cipher proposals aimed at balancing cost and security in constrained settings, such as Present <cit.>, KATAN <cit.>, LED <cit.>, LBlock <cit.>, Prince <cit.>, Skinny <cit.>, among others. These examples underscore the advantages of low-cost cryptographic primitives.

MDS matrices, with their optimal branch number and excellent diffusion properties, are well-suited for linear diffusion layers and are widely used in cryptographic designs. However, it's important to note that MDS matrices can be relatively dense, resulting in higher hardware implementation costs(larger circuit areas). Therefore, to incorporate MDS matrices as linear diffusion layers in lightweight ciphers, it becomes necessary to optimize existing MDS matrices or construct new lightweight MDS matrices that strike a balance between security and efficiency.



 §.§ Related works


There has been a significant amount of research dedicated to constructing lightweight MDS matrices or near-MDS matrices.

One common approach for constructing MDS matrices involves selecting a matrix with a specific structure, such as Hadamard, Cauchy, Vandermonde, Toeplitz, circulant, or other matrices with special properties. Then, specific values are assigned to the variable parameters of the matrix to ensure that all its submatrices are nonsingular, resulting in an MDS matrix. This method reduces the search space by focusing on one or a few specific types of matrices. Additionally, for matrices with special structures, there are often many repeated values in the determinants of their submatrices, which reduces the computational complexity when verifying the nonsingularity of the submatrices. Remarkably, for Cauchy matrices, their structure guarantees that they are MDS matrices regardless of the parameter values. In summary, searching for MDS matrices within matrices with special structures is computationally efficient and often successful. There is a considerable amount of literature on this topic, including works based on Cauchy matrices <cit.>, <cit.>, <cit.>, <cit.>, <cit.>,  Vandermonde matrices  <cit.>, <cit.>, <cit.>, <cit.>, Hadamard matrices <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>,  Toeplitz matrices <cit.>, <cit.>,  <cit.>, <cit.>, <cit.>, and circulant matrices <cit.>, <cit.>, <cit.>,  <cit.>, <cit.>.

However, this approach of limiting the search to specific matrix structures has some drawbacks. The implementation cost of the resulting matrices is often difficult to control. Early methods focused on optimizing the implementation of matrix elements, known as local optimization. For example, in a circulant matrix, where the entire matrix is determined by its first row consisting of k elements (which are n× n invertible matrices over F_2), local optimization involves optimizing the implementation of these k matrices while preserving the XOR operations between them. This approach cannot avoid the overhead of (k-1)kn XOR operations for the entire matrix. Relevant works include <cit.>, <cit.>, <cit.>. Specifically, in <cit.>, Kölsch provided the specific structure of all matrices that can be implemented using only 2 XORs for element (treated as a matrix) multiplication with a vector in a field of characteristic 2. Local optimization methods have a smaller computational footprint because they only focus on optimizing matrix elements, but their effectiveness is limited as they do not consider the overall matrix structure.

To improve optimization efficiency, a natural approach is to treat a k× k matrix with  n× n invertible matrices as elements as a single kn× kn matrix over F_2 and optimize the entire matrix, known as global optimization <cit.>. It involves in the so-called SLP (Straight Line Program) problem. In <cit.>, Boyar et al. reduced the point set covering problem to a subproblem of the SLP problem (SLP problem with Hamming weight limited to 3 for each row) and proved that the SLP problem is NP-hard. This indicates that finding a universally applicable and efficient algorithm for it is challenging. Currently, heuristic algorithms are commonly used to address this issue, including the Paar algorithm <cit.>, the BP algorithm <cit.> and its variants <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>.

Recently, Yang et al. introduced a new search method for involutory MDS matrices and obtained a 4× 4 involutory MDS matrix implementable with only 35 XORs (4-bit input) and 70 XORs (8-bit input) <cit.>. Here it should be pointed out that their 35-XOR result has being optimal in the 4-bit input. Pehlivano?lu et al. further considered 3-input and higher XOR gates in their variant of the BP algorithm and achieved MDS matrices with optimized circuit area <cit.>.

An alternative approach to constructing lightweight MDS matrices is to limit the implementation cost of the matrix, then derive the matrix from its implementation, and finally check if it satisfies the MDS condition, i.e., all submatrices are nonsingular. This method was first introduced by Duval and Leurent in <cit.> and extended by Zeng et al <cit.>. Roughly speaking,  it considers matrices of the form M = A_1A_2⋯ A_t with an upper bound t, where A_1,A_2,⋯,A_t have low implementation costs. By means of careful selection for A_i, some 4× 4 MDS matrix with 35 XOR gates (4-bit input) and 67 XOR gates (8-bit input) are constructed, which have been proven to be optimal under a word-based structure <cit.>. In <cit.>, Sajadieh and Mousavi further applied it to search higher-order MDS matrices with a General Feistel Structure and obtained some  6× 6 and 8× 8 MDS matrices with relatively low costs.




 §.§ Our contributions


In this study, our main objective is to construct the lightest weight MDS matrices. We first investigate the structures of MDS matrices with low implementation costs and give all possible structures with optimal costs. Based on these structures, we then generate a significant number of new ultra-light MDS matrices. 

More speaking, we introduce the concept of an "implementation tree". Based on the implementation tree, we develop a new method that starts with several unit vectors and iteratively generates rows for the MDS matrix. As each new row is generated, we check if it satisfies the MDS condition. This verification step is crucial as it eliminates a large number of unproductive branches during the search process. Through this refined approach, we are able to identify all MDS matrix structures with optimal costs (based on word metrics). By leveraging these structures, we successfully find all 60 nonequivalent 4× 4 MDS matrices that can be implemented with just 67 XOR operations on the ring F_2[x]/(x^8 + x^2 + 1) and 35 XOR operations on the ring F_2[x]/(x^4 + x + 1). Notably, our results encompass the 10 matrices in <cit.> and the 52 matrices in <cit.>. Here it should be pointed that our method is also capable of producing a significant quantity of 4× 4 MDS matrices with 36-41 XOR gates for 4-bit input and 68-80 XOR gates for 8-bit input by modestly relaxing the implementation cost constraints. Due to the limit of the paper, we only list the counts of those with the lowest cost for various depths in Table <ref>. Furthermore, we identify 4× 4 involutory MDS matrices on the same ring that require only 68 XOR operations for implementation. To the best of our knowledge, it is the most optimal result currently available.

To demonstrate the effectiveness of our algorithm, we conduct a comparative analysis with the methodology described in <cit.>. The results of this comparison are presented in Table <ref>.



0.7
[htbp]
	
	
	The number of MDS matrices of different depths
	
		 
		
			
			Ring     Depth     Cost     <cit.>     <cit.>     Our results  

			
			F_2[α]     -      67     2     52     60

			    5     67     2     4     4

			    4     69     1     4     7

			    3     77     1     2     3 

			F_2[β]     -      35     2     52     60

			    5     35     2     4     4

			    4     37     1     4     7

			    3     41     1     2     3 

			
		
			
  * α is the companion matrix of x^8+x^2+1, β is the companion matrix of x^4+x+1.
		
	







[htbp]
	Comparison of the efficiency of our algorithm with previous results
	
		2*Ring     2*Cost     2*Depth     2cPrevious Algorithm     2cOur Algorithm     2*Ref. 

		4-56-7
		            Memory     Time     Memory     Time     

		F_2[α]     67     5     30.9G     19.5h     4M     <10ms     <cit.>

		F_2[α]     68     5     24.3G     2.3h     4M     <10ms     <cit.>

		F_2[α]     69     4     274G     30.2h     4M     <10ms     <cit.>

		
	
		
  *     α is the companion matrix of x^8+x^2+1.
	




As for higher-order MDS matrices, we successfully constructs 5× 5 and 6× 6 matrices for 8-bit inputs with 114 and 148 XOR gates respectively, which are the best results up to date, see Table <ref>.




[H]
	
	
	
	The XOR gates required to implement MDS matrices of different orders
	
		       
		
		
			
			Size of matrix     Previous results     Our results     Ref. 

			
			5     129      114     <cit.> 

			6     156     148     <cit.> 

			
	






 §.§ Organization of the paper


In Section 2, we introduce several essential conceptions that serve as the foundation for our subsequent discussions. In Section 3, we propose a novel algorithm aimed at efficiently searching for ultra-light 4× 4 MDS (Maximum Distance Separable) matrices and involutory MDS matrices. In Section 4, we extend our focus to the construction of higher-order MDS matrices, especially for 5× 5 and 6× 6 MDS matrices. Finally, in Section 5 we summarize our key findings and conclusions.




§ PRELIMINARIES


Since linear functions defined over vector spaces can be equivalently expressed as matrix-vector multiplications, we restrict our attention to the corresponding matrix representations in the following discussion. Let n and k be positive integers, and let F_2 denote the binary field consisting of the elements 0 and 1.

In the context of SPN ciphers, linear functions often operate on the concatenated outputs of multiple S-boxes. Here, we specifically consider a k× k matrix, whose entries are n× n invertible matrices over F_2. In this setting, k represents the number of S-boxes, while n represents the size of a single S-box. Alternatively, this structure can be interpreted as a larger nk× nk matrix defined over F_2. For simplicity and consistency, we further constrain these n× n matrices to belong to a commutative ring R_n that is a subset of GL(n,2) - the general linear group of n× n invertible matrices over the binary field F_2. This restriction applies throughout the remainder of this paper. As a concrete example, one might consider the case where R_n=F_2[α] and α is the companion matrix of a polynomial of degree n. We denote the set of all k× k invertible matrices with entries in R_n by M_k(R_n).



 §.§ MDS matrix



The branch number is a significant cryptographic parameter used to measure the diffusion properties of linear functions. Given that any arbitrary linear function L can be represented in matrix form, i.e., there exists a matrix M such that L(x) = Mx, the concept of the branch number can be extended to any matrix. In the following, we introduce the notion of the branch number for matrices.


	Let M ∈ M_k(R_n) be a given matrix. Its differential and linear branch numbers are defined respectively as:	
	
    B_d(M) = min_x ≠ 0{wt(x) + wt(Mx)} and  B_l(M) = min_x ≠ 0{wt(x) + wt(M^Tx)},

	where wt(x) denotes the weight of x, i.e., the number of nonzero components in x, and M^T represents the transpose of M.


It is evident that for any matrix M in M_k(R_n), its differential and linear branch numbers are at most k+1. If either of these branch numbers equals k+1, the matrix is designated as an MDS matrix. Notably, if the differential branch number of a matrix achieves optimality, its linear branch number does so as well, and vice versa.

For a given matrix M, the following theorem provides a necessary and sufficient condition for M to be an MDS matrix.


	Let R_n be a commutative ring and M be a k × k matrix over R_n. Then M is an MDS matrix if and only if all of its minors are invertible.


Let 𝐞_i denote the i-th unit vector, that is, 𝐞_i = (0,…,0,i1,0,…,0) with 1 in the i-th position. A permutation matrix is an invertible matrix whose rows consist of unit vectors.


	Let M be an MDS matrix and P, Q be two permutation matrices. Then M' = PMQ is also an MDS matrix.


For any given MDS matrix M, an MDS matrix M' is said to be equivalent to M if there exist two permutation matrices P and Q such that M' = PMQ, denoted by M' ∼ M. Since the inverse of a permutation matrix is also a permutation matrix, this equivalence relation is symmetric, meaning if M' ∼ M, then M ∼ M'. Let [M] represent the set of all MDS matrices equivalent to M. It is straightforward to observe that [M] comprises (k!)^2 MDS matrices.



 §.§ Cost and depth of binary matrices


The differential and linear branch numbers constitute two pivotal cryptographic indicators of a matrix. When a matrix realized as a circuit, its depth and cost emerge as two additional essential parameters. Subsequently, we study the circuit realization of a general binary matrix, emphasizing its depth and cost implications.


	Consider any nonzero binary m× n matrix M. An implementation I of M refers to a circuit constructed solely using XOR gates that can compute y=Mx for any input vector x, where x = (x_1,x_2,…,x_n)∈ F_2^n and y = (y_1,y_2,…,y_m)∈ F_2^m.


As an illustrative example, let us consider a 4×4 binary matrix

    M=[ 1 0 0 0; 1 1 0 0; 1 1 1 0; 1 1 1 1 ].

Its implementation I can be expressed as:

    y_1    = x_1, 
    
    		y_2    = x_1 ⊕ x_2, 
    
    		y_3    = (x_1 ⊕ x_2) ⊕ x_3, 
    
    		y_4    = (x_1 ⊕ x_2) ⊕ (x_3 ⊕ x_4).


Figure <ref> depicts the structure of I, where intermediate terms x_5 = x_1 ⊕ x_2 and x_7 = x_3 ⊕ x_4 are introduced for clarity. Since I comprises solely of XOR gates, it gives rise to a binary tree B_i rooted at y_i for each y_i (1 ≤ i ≤ m), known as the spanning tree of y_i. In B_i, nodes other than the root y_i are denoted as descendants of y_i.








	For a given binary matrix M with an implementation I, the cost of I is defined as the total count of XOR gates employed in I.



	Given a binary m× n matrix M and its implementation I, the depth of I is determined as the maximum height attained among all the trees B_i rooted at y_i (1 ≤ i ≤ m).


Considering the aforementioned 4×4 binary matrix M and its implementation I, we observe that the depth and cost of I are 2 and 4, respectively. Notably, a binary matrix M may admit multiple distinct implementations, each potentially exhibiting different costs and depths. For instance, Figure <ref> presents an alternative implementation I' of the same 4×4 matrix M, with a depth and cost of 3. In our work, we prioritize implementations with the lowest cost for any given matrix M.







Let M and M' be matrices related by M' = PMQ, where P and Q are permutation matrices, and let I be an implementation of M. An implementation I' of M' is said to be a corresponding implementation of I if it is obtained solely by rearranging the input and output variables in accordance with P and Q. Thus we have the following conclusion.


	Given equivalent matrices M and M', with I being an implementation of M and I' being the corresponding implementation of I for M', it follows that I and I' have identical cost and depth.


By Theorem <ref>, it becomes evident that all matrices belonging to the same equivalence class share a common cost and depth. Therefore, in our search for ultra-light MDS matrices, we can confine our attention to the distinct equivalence classes.



 §.§ Metric


In this section, we discuss the method to compute the cost of a k × k matrix M over R_n:

    [ y_1; y_2;   ⋮; y_k ]
    	=
    	[ a_11 a_12    ⋯ a_1k; a_21 a_22    ⋯ a_2k;    ⋮    ⋮    ⋱    ⋮; a_k1 a_k2    ⋯ a_kk ][ x_1; x_2;   ⋮; x_k ],

where x_i, y_i ∈ F_2^n, a_ij∈ R_n, 1 ≤ i, j ≤ k. Here, we only consider block implementations of M, treating each x_i and y_i as individual blocks.

For a given matrix M and its block implementation I, two fundamental operations exist:

	
  * Word-wise XOR: X  X_1 ⊕ X_2,
	
  * Scalar multiplication over R_n: X  aX,


where X, X_1, X_2 are n-dimensional vectors over F_2 and a ∈ R_n. Word-wise XOR requires n XOR gates; the cost of scalar multiplication aX over R_n depends on the specific value of a.

Any block implementation I of matrix M can be expressed as a linear straight-line program, as exemplified below.


	
    [ y_1; y_2; y_3; y_4 ]
    		=
    		[ 1 α 1 0; 1 0 1 1; 0 α 1 1; 1 1 1 0 ][ x_1; x_2; x_3; x_4 ],

	where x_i, y_i ∈ F_2^8 (1 ≤ i ≤ 4), and α is the companion matrix of x^8 + x^2 + 1:
	
	
		
    α =
    			[ 0 0 0 0 0 0 0 1; 1 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 1; 0 0 1 0 0 0 0 0; 0 0 0 1 0 0 0 0; 0 0 0 0 1 0 0 0; 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 1 0;   ].

	
	
One implementation of M is given by:

    y_1    = (x_1 ⊕ x_3) ⊕α x_2, 
    
    		y_2    = (x_1 ⊕ x_3) ⊕ x_4, 
    
    		y_3    = (x_3 ⊕ x_4) ⊕α x_2, 
    
    		y_4    = (x_1 ⊕ x_3) ⊕ x_2.

This can be rearranged into a linear straight-line program I as follows:

    x_5    = x_1 ⊕ x_3, 
    
    	x_6    = x_3 ⊕ x_4, 
    
    	x_7    = x_5 ⊕α x_2 = y_1, 
    
    	x_8    = x_5 ⊕ x_4 = y_2, 
    
    	x_9    = x_6 ⊕α x_2 = y_3, 
    
    	x_10   = x_5 ⊕ x_2 = y_4.


Note that each step of I has the form:

    X = a · X_1 ⊕ b · X_2,

where X_1, X_2 are either intermediate results x_i(i≥ 5) calculated previously or input variables x_i (1 ≤ i ≤ 4).

Two distinct types of operations are observed in this form: word-wise XOR, denoted as word-xor(expressed as XOR) to distinguish it from bit-wise xor(expressed as xor), and scalar multiplication, which refers to matrix-vector multiplication over R_n.

To compute the cost of I, we proceed as follows:

First, the linear straight-line program comprises 6 steps, yielding an XOR cost of 6 × 8 = 48.

Second, considering scalar multiplication, the non-trivial multiplication in I is α x_2, which appears twice. Since the latter occurrence of α x_2 can reuse the result of the first, only one calculation is needed when assessing cost. Therefore, the cost of scalar multiplication is 1.

In conclusion, the total cost of I is 48 + 1 = 49 xor operations.


From the above example, it is evident that when determining the cost of an implementation I of a k × k matrix M over R_n, we follow a two-step process. Firstly, we ascertain the number of steps s in I as a linear straight-line program, dictating the number of XORs necessary in I. Secondly, we identify all unique and non-trivial scalar multiplications in I and tally the total number of xors t required for these scalar multiplications. Ultimately, the cost of I is given by ns + t.



§ THE CONSTRUCTION OF LIGHTWEIGHT 4×4 MDS MATRICES


In this section, we mainly focus on the construction of the 4× 4 MDS matrix. Because of its small size, we can discuss it in more details. First, we will introduce the concept of implementation tree.



 §.§ Implementation tree


In Example <ref>, we observe that under the word-based structure, a matrix implementation can be conceptualized as a linear straight-line program. Each such program uniquely determines a matrix. The lightweight MDS matrix can be constructed by searching for the linear straight-line program with the lowest cost capable of implementing an MDS matrix. Next, we investigate the properties of these linear straight-line programs.

Consider a linear straight-line program realizing a matrix:

    t_1 = a_1 t_1,1⊕ b_1 t_1,2, t_2 = a_2 t_2,1⊕ b_2 t_2,2, …, t_m = a_m t_m,1⊕ b_m t_m,2,

where t_i,j (1 ≤ i ≤ m, 1 ≤ j ≤ 2) are either initial variables x_p or intermediate results t_q (1 ≤ q < i). We designate t_q (1 ≤ q ≤ m) as a term of the program, and the number of terms m as its length.

A partial order relation can be established among the terms t_1, …, t_m: if t_p = a_pt_q ⊕ b_pt_p,2 or t_p = a_pt_p,1⊕ b_pt_q appears in the program, then t_q ≺ t_p. Furthermore, if t_q ≺ t_p and t_r ≺ t_q, then t_r ≺ t_p. For instance, if t_5 = at_3 + bt_4 and t_4 = ct_1 + dt_2, then t_i ≺ t_5 for 1 ≤ i ≤ 4. The relation t_q ≺ t_p signifies that t_q is an essential component of t_p.

Given the partial order relationship, suppose the following linear straight-line program can implement some k × k matrix on R_n, with a_pt_p,1⊕ b_pt_p,2 = t_p where a_p ≠ 0 and b_p ≠ 0:


    a_1t_1,1⊕ b_1t_1,2   = t_1, 
       ⋮
    
    	a_i_1t_i_1,1⊕ b_i_1t_i_1,2   = t_i_1 = y_1, 
    
    	a_i_1+1t_i_1+1,1⊕ b_i_1+1t_i_1+1,2   = t_i_1+1, 
       ⋮
    
    	a_i_2t_i_2,1⊕ b_i_2t_i_2,2   = t_i_2 = y_2, 
       ⋮
    
    	a_i_kt_i_k,1⊕ b_i_kt_i_k,2   = t_i_k = y_k.


This program exhibits the following properties:


	
  * If t_p is not an output (i.e., t_p ≠ y_i for all 1 ≤ i ≤ k), then t_p ≺ y_i for some 1 ≤ i ≤ k.
	
  * If 1 ≤ q < p ≤ i_k, then either t_q ≺ t_p or t_q and t_p are incomparable. If t_q and t_p are incomparable, their order in the program can be swapped.


Let us elaborate on these properties:


	
  * The partial order relation t_b ≺ t_a indicates that t_b is necessary for computing t_a, that is to say, t_a cannot be computed before t_b. Conversely, if t_b ⊀t_a, whether t_b has been generated is irrelevant to the computation of t_a. For Property 1, if there exists a non-output t_p such that t_p ⊀y_i for all 1 ≤ i ≤ k, this implies that t_p is not required for computing any outputs y_i (1 ≤ i ≤ k). Therefore, removing t_p from the matrix implementation would not affect its functionality, and we can reduce the program's length and simplify the implementation. Thus we always assume that no such a t_p exists in a matrix implementation.
	
	
  * Consider t_q and t_p (1 ≤ q < p ≤ i_k). If t_p ≺ t_q, it means that t_p must be used to compute t_q. However, since t_p is generated after t_q, this leads to a contradiction. Therefore, their relationship can only be t_q ≺ t_p or they are incomparable. If t_p and t_p+1 are incomparable, swapping their positions in the program has no impact on the overall implementation. This is because there are no terms between t_p and t_p+1, so the swap does not affect the generation of these two terms or any subsequent terms. However, it should be noted that if t_p and t_p+1 are incomparable and p - q ≥ 2, such a swap may not be feasible. For example, if there exists q < s < p such that t_s ≺ t_p, swapping t_q and t_p would place t_s after t_p, preventing the generation of t_p and leading to a contradiction.


For brevity, we represent the above linear straight-line program as [t_1, …, t_i_1 = y_1, …, t_i_k = y_k]. Applying properties 1 and 2, we arrive at the following theorem:


	For any given k × k MDS matrix M on R_n, an implementation I = [t_1, …, t_i_1 = y_1, …, t_i_k = y_k] of M can be rearranged to I' = [t_1', …, t_i_1'' = y_1, …, t_i_k' = y_k] by swapping the relative positions of t_p (1 ≤ p ≤ i_k) such that I' satisfies:
	
    t_1'    ≺ y_1, …, t_i_1'-1' ≺ y_1, 
    
    		t_i_1'+1'    ≺ y_2, …, t_i_2'-1' ≺ y_2, 
       ⋮
    
    		t_i_k-1'+1'    ≺ y_k, …, t_i_k-1' ≺ y_k.

  

	Consider the implementation I of M, which we divide into k segments as follows:
	
    t_1,   ⋯,t_i_1-1,t_i_1=y_1,
    
    		t_i_1+1,   ⋯,t_i_2-1,t_i_2=y_2,
       ⋮
    
    		t_i_k-1+1   ⋯,t_i_k-1,t_i_k=y_k.

	
	Utilizing Property 2, we commence the adjustment process segment by segment, starting from the first.
	
	In the first segment, we initiate the downward search from t_i_1-1. Suppose the first term encountered that is incomparable with y_1 is t_j_1. Then, under the partial order relation ≺, t_j_1 is incomparable with t_j_1+1,⋯,t_i_1=y_1. This arises because t_j_1 is the inaugural term discovered in the descending search that cannot be compared with y_1. If there existed some p such that t_j_1≺ t_p (j_1 < p < i_1) held true, then since t_p≺ y_1, we would have t_j_1≺ y_1?a contradiction. Hence, t_j_1 is indeed incomparable with t_j_1+1,⋯,t_i_1=y_1. By Property 2, we can successively exchange t_j_1 upward, akin to bubble sorting, until it is positioned behind y_1. This process is repeated until the remaining terms t_1',⋯,t_i_1'-1' satisfy t_s≺ y_1 for all 1≤ s < i_1'. Thus, the first segment is appropriately adjusted.
	
	Segments 2 to k-1 are adjusted analogously to the first segment. After the preceding adjustments, any term t_p in segment k necessarily satisfies t_p⊀y_s for all 1≤ s≤ k-1. By virtue of Property 1, t_p≺ y_k ensures that the last segment naturally conforms to the requirements subsequent to the adjustments of the preceding k-1 segments.
	
	After the aforementioned adjustments, I' fulfills the theorem's conditions, thereby the conclusion follows.


Below we designate I' as a normal implementation. Notably, our adjustments solely alter the order of terms within the linear straight-line program. Therefore they are identical when they are realized as circuits.

To illustrate Theorem <ref>, we employ the matrix and implementation from Example <ref>. In this instance, the implementation I takes the following form:

    x_1⊕ x_3    =x_5,
    
    	x_3⊕ x_4    =x_6,
    
    	x_5⊕α x_2    =x_7=y_1,
    
    	x_5⊕ x_4    =x_8=y_2,
    
    	x_6⊕α x_2   =x_9=y_3,
    
    	x_5⊕ x_2    =x_10=y_4.

The partial order relationships among the terms are as follows:

    x_5≺ x_7,x_5≺ x_8,x_5≺ x_10,x_6≺ x_9.

Adjusting I according to Theorem <ref>: Since x_6 is incomparable with x_7=y_1, we exchange x_6 and x_7 and finish the first segment's adjustment. Next, as x_6 is incomparable with x_8, we swap x_6 and x_8 in the second segment's adjustment. For the third segment, the relationship x_6≺ x_9=y_3 already exists, we do nothing. Similarly, the fourth segment requires no modification.

After renumbering the terms, we obtain I' as follows:

    x_1⊕ x_3    =x_5,
    
    	x_5⊕α x_2    =x_6=y_1,
    
    	x_5⊕ x_4    =x_7=y_2,
    
    	x_3⊕ x_4    =x_8,
    
    	x_8⊕α x_2   =x_9=y_3,
    
    	x_5⊕ x_2    =x_10=y_4.


The distinction between I and I' lies solely in the ordering of terms, as both ultimately express the same set of relationships:

    y_1=(x_1⊕ x_3)⊕α x_2,
    
    		y_2=(x_1⊕ x_3)⊕ x_4,
    
    		y_3=(x_3⊕ x_4)⊕α x_2,
    
    		y_4=(x_1⊕ x_3)⊕ x_2.


Based on the preceding discussion, when an MDS matrix is implemented, it suffices to consider a normal implementation akin to the form I' in Theorem <ref>. We call the vector (i_1', i_2' - i_1', …, i_k' - i_k-1') as a characteristic of I'. Once the order among y_i is fixed, its characteristic remains constant.

To further illustrate the concept, let us consider the MixColumn matrix of the AES algorithm as an example. We introduce the notion of an implementation tree through this example.


	The MixColumn matrix of AES is given by:
	
	
    M_AES =
    	[ [   α α+1   1   1;   1   α α+1   1;   1   1   α α+1; α+1   1   1   α ]],

	where α denotes the companion matrix of x^8 + x^4 + x^3 + x + 1.
	
	An implementation I of M_AES can be expressed as follows:
	
	[H]
		130mm! 
			
				x_5 = 1 · x_1 ⊕ 1 · x_2,     x_6 = 1 · x_2 ⊕α· x_5,     x_7 = 1 · x_3 ⊕ 1 · x_4,     x_8 = 1 · x_6 ⊕ 1 · x_7=y_1, 

				x_9 = 1 · x_2 ⊕ 1 · x_3,     x_10 = 1 · x_3 ⊕α· x_9,     x_11 = 1 · x_1 ⊕ 1 · x_4,     x_12 = 1 · x_10⊕ 1 · x_11=y_2, 

				x_13 = 1 · x_4 ⊕α· x_7,     x_14 = 1 · x_5 ⊕ 1 · x_13=y_3,         

				 x_15 = 1 · x_1 ⊕α· x_11,     x_16 = 1 · x_9 ⊕ 1 · x_15=y_4.         

		
	
	
	The cost associated with I is 12 × 8 + 4 × 3 = 108. Notably, the cost of XOR operations is 96 and comprises the majority of the total implementation cost. If we abstract away the specific values of scalar multiplications, the implementation I can be generalized as follows:
	
	[H]
		130mm!
			
				x_5 = a_1 · x_1 ⊕ a_2 · x_2,     x_6 = a_3 · x_2 ⊕ a_4 · x_5,     x_7 = a_5 · x_3 ⊕ a_6 · x_4,     x_8 = a_7 · x_6 ⊕ a_8 · x_7=y_1, 

				x_9 = a_9 · x_2 ⊕ a_10· x_3,     x_10 = a_11· x_3 ⊕ a_12· x_9,     x_11 = a_13· x_1 ⊕ a_14· x_4,     x_12 = a_15· x_10⊕ a_16· x_11=y_2, 

				x_13 = a_17· x_4 ⊕ a_18· x_7,     x_14 = a_19· x_5 ⊕ a_20· x_13=y_3,         

				x_15 = a_21· x_1 ⊕ a_22· x_11,     x_16 = a_23· x_9 ⊕ a_24· x_15=y_4.         

		
	
	
	This abstraction highlights the role of XOR operations in the implementation. It can be visually represented using a binary tree, as shown in Figure <ref>.
	
	


From Figure <ref>, we can observe the following:

	
  * The tree is regular, meaning each node is either a leaf node or has two non-empty subtrees.
	
  * Each node in the tree is generated from two nodes with smaller serial numbers, defining the order of implementation.



	Consider a given k × k MDS matrix M on R_n. Let I = [t_1, …, t_i_1 = y_1, …, t_i_k = y_k] be a normal implementation of M. We stipulate that x_1 = t_0, x_2 = t_-1, …, x_k = t_-(k-1). For any term t_p = a_p t_m ⊕ b_p t_n in I, it induces a triple T_p = (m, n, p). If t_s = y_i, we mark its corresponding triple as T_s → y_i. By arranging such triples in the order of their appearance in the original implementation I, we obtain a sequence T = [T_1, …, T_i_1→ y_1, …, T_i_k→ y_k].
	
	We refer to T as the implementation tree of M. Each triple T_p (1 ≤ p ≤ i_k) is called a node of T. The vector (i_1, i_2 - i_1, …, i_k - i_k-1) is designated as the type of T, and i_k is termed the capacity of T. Among all k × k MDS matrices on R_n, the implementation tree with the minimum capacity is deemed the simplest tree.
	


In Definition <ref>, the implementation tree is represented as triples in the form T_p = (m,n,p). Alternatively, it can be expressed using XOR operations with parameters. Let x_1=T_0, x_2=T_-1, …, x_k=T_-(k-1) and the implementation tree be denoted as T=[T_1, …, y_1=T_i_1, …, y_k=T_i_k], where T_p = a_pT_m ⊕ b_pT_n corresponds to T_p=(m,n,p). Here, a_p and b_p are formal parameters without specific values. Notably, the partial order relationship between terms in I can be seamlessly translated to the nodes in T.

We define the implementation tree as follows: matrix → matrix implementation → implementation tree of matrix. Initially, we have the matrix, from which we derive the word-based implementation of the matrix, and finally extract the concept of the implementation tree. When constructing a matrix, we follow the reverse process: first, we construct an implementation tree, then determine the parameter values for each scalar multiplication to achieve a concrete implementation, and finally derive the matrix from this implementation.

Before presenting the specific algorithm, we estimate the capacity of the MDS matrix implementation tree.


	For a given k × k MDS matrix M on R_n, let I=[t_1, …, t_i_1=y_1, …, t_i_k=y_k] be an implementation of M. When fully expanding t_i_p=y_p (1 ≤ p ≤ k) to the input variable x_q (1 ≤ q ≤ k), each x_q must appear.



	Each t_i_p=y_p (1 ≤ p ≤ k) corresponds to row p of M. Since M is an MDS matrix, the elements in row p are non-zero. Therefore, each x_q (1 ≤ q ≤ k) must appear in the expansion of t_i_p=y_p.


The following theorem provides a lower bound for the capacity of MDS matrix implementations.


	For a given k × k (k ≥ 3) MDS matrix M on R_n, let T=[T_1, …, T_i_1→ y_1, …, T_i_k→ y_k] be an implementation tree of M. Then, i_k ≥ 2k-1.



	Assume an implementation I=[t_1, …, t_i_1=y_1, …, t_i_k=y_k] of M corresponds to T. By Lemma <ref>, we have i_1 ≥ k-1. Since each output variable t_i_p=y_p (2 ≤ p ≤ k) must be generated, it follows that i_p - i_p-1≥ 1 (2 ≤ p ≤ k). Thus, i_k ≥ 2k-2. Equality holds only if i_1 = k-1 and i_p - i_p-1 = 1 (2 ≤ p ≤ k). We show that this situation does not arise, implying i_k ≥ 2k-1.
	
	If i_1 = k-1 and i_p - i_p-1 = 1 (2 ≤ p ≤ k), by Lemma <ref>, each x_q appears exactly once in the expansion of t_i_1=y_1. Consider y_2 = t_i_2 = a_i_2t_i_2,1⊕ b_i_2t_i_2,2 with i_2,1 < i_2,2. Then, t_i_2,2≺ y_1 or t_i_2,2 = y_1 since y_2 is generated immediately after y_1. We analyze two cases:
	
	
		
  * Consider the case where t_i_2,1 and t_i_2,2 are not comparable under the partial order relation ≺. It follows that t_i_2,2≺ y_1. Suppose that y_1 = t_i_1 = a_i_1t_i_2,1⊕ b_i_1t_i_2,2. If this is not the case, then y_1 = t_i_1 = a_i_1t_i_1,1⊕ b_i_1t_i_1,2, which implies that either t_i_2,1≺ t_i_1,1 or t_i_2,2≺ t_i_1,1 or t_i_2,1≺ t_i_1,2 or t_i_2,2≺ t_i_1,2. This would lead to the omission of some x_q when expanding y_2 in terms of the input variables, contradicting Lemma <ref>. Assuming y_1 = t_i_1 = a_i_1t_i_2,1⊕ b_i_1t_i_2,2, since k ≥ 3, at least one of t_i_2,1 or t_i_2,2 must contain more than two input variables. Let x_m and x_n be two such variables. Considering the 2 × 2 matrix M_2 of M corresponding to x_m, x_n; y_1, y_2, we find that the determinant of M_2 is 0, contradicting the fact that M is an MDS matrix. Therefore, the first case does not hold.
		
		
  * Now consider the case where t_i_2,1≺ t_i_2,2. In this scenario, it must be that t_i_2,2 = y_1. Otherwise, there would be some x_q missing when expanding y_2 in terms of the input variables, contradicting Lemma <ref>. Assuming t_i_2,2 = y_1, if t_i_2,1 contains more than two input variables, let x_m and x_n be two such variables. If not, choose x_m and x_n from the missing terms of t_i_2,1. Considering the 2 × 2 submatrix M_2' of M corresponding to x_m, x_n; y_1, y_2, we find that the determinant of M_2' is 0, contradicting the fact that M is an MDS matrix. Therefore, the second case does not hold.
	
	
	In both cases, we reach a contradiction. Therefore, i_k ≥ 2k-1.




 §.§ An algorithm for searching the simplest tree


In the preceding section, we introduced the concept of the implementation tree. We can leverage this structure to construct the k × k MDS matrix over R_n. As exemplified in Example <ref>, the cost of the word-xor operation constitutes a significant portion of the overall cost in implementing an MDS matrix. Therefore, our approach prioritizes searching for implementation trees with minimal capacity. Subsequently, we assign parameters for scalar multiplications to derive the specific implementation, and ultimately yield a low-cost MDS matrix.

Our algorithm is divided into two main parts:

	
  * Firstly we treat the parameters a, b of each node Y = a · Y_1 ⊕ b · Y_2 in an implementation tree as undetermined coefficients. We perform a traversal search to identify all the simplest trees which have the potential to form MDS matrices.
	
	
  * Secondly we assign values to the undetermined parameters within the simplest trees. This instantiation process will result in a specific MDS matrix.


By Theorem <ref>, for a given k × k (k ≥ 3) MDS matrix M over R_n, any implementation T = [T_1, ⋯, y_1 = T_i_1, ⋯, y_k = T_i_k] of M satisfies the conditions i_k ≥ 2k - 1, i_1 ≥ k - 1, and i_p - i_p-1≥ 1 (2 ≤ p ≤ k). Therefore, we commence our search with i_k = 2k - 1 and look for the simplest trees with this capacity. If there does not exist such a simplest tree, we increment i_k and repeat the search.

Next, we outline the procedure for finding a possible implementation tree with a capacity of i_k as below:

	
  * We solve the indefinite equation s_1 + s_2 + ⋯ + s_k = i_k under the constraints s_1 ≥ k - 1 and s_i ≥ 1 (2 ≤ i ≤ k) to obtain the set of all possible types, denoted by I_i_k, where each type represents the number of non-leaf nodes in its corresponding implementation tree.
	
	
  * For each type (i_1, i_2 - i_1, ⋯, i_k - i_k-1) ∈ I_i_k, we initiate the search with input variables x_1, x_2, ⋯, x_k. Leveraging the property that each node T_p in T = [T_1, ⋯, y_1 = T_i_1, ⋯, y_k = T_i_k] is the sum of preceding nodes, we perform a traversal search and sequentially generate k binary trees (corresponding to rows of the matrix) of the specified type. Concurrently, we verify if the matrix associated with T is an MDS matrix. For simplicity, we set x_1 = T_0, x_2 = T_-1, ⋯, x_k = T_-(k-1).


To generate the first binary tree, rooted at y_1:

	
  * We recursively generate the tree starting from T_1. Each T_1 is expressed as T_1 = a_1T_p_1⊕ b_1T_p_2 with -(k-1) ≤ p_1 < p_2 ≤ 0. This process yields C_k^2 distinct T_1 options. Each T_1 can then serve as a child node in subsequent trees.
	
	
  * We proceed to generate T_2 using a similar approach: T_2 = a_2T_p_1⊕ b_2T_p_2 with -(k-1) ≤ p_1 < p_2 ≤ 1. For each T_1, we obtain C_k+1^2 different T_2 options, resulting in a total of C_k^2 · C_k+1^2 unique (T_1, T_2) sequences.
	
	
  * We repeat this process until we generate y_1 = T_i_1. At this point, we have generated C_k^2 · C_k+1^2 ⋯ C_k+i_1-1^2 sequences of the form (T_1, T_2, ⋯, y_1 = T_i_1).


At this stage, each generated sequence ensures that each term T_p is expressed as the sum of its predecessors. However, we must further verify the following conditions to determine if these sequences can constitute the first tree in the implementation tree:


	
  * T_1 ≺ y_1, T_2 ≺ y_1, ⋯, T_i_1-1≺ y_1. This condition ensures the normality of the implementation and can be checked relatively easily.
	
  * Whether the 1 × k matrix corresponding to y_1 satisfies the MDS condition (i.e., all submatrix determinants are non-zero). Since the scalar multiplications in the implementation tree lack specific values, the elements of the 1 × k matrix are polynomials in a_1, b_1, ⋯, a_i_1, b_i_1. In practice, we verify if the determinants of all submatrices of the 1 × k matrix are not divisible by 2 by means of symbolic calculation. If a determinant is divisible by 2, it implies that regardless of the specific values assigned to the parameters a_1, b_1, ⋯, a_i_1, b_i_1, they cannot form an MDS matrix. Such sequences are subsequently discarded.


Let S_1 denote the set of all sequences that pass the above test. Each sequence (T_-(k-1),⋯,T_i_1) in S_1 serves as the foundation for generating the second tree in the implementation tree hierarchy. 
Following a similar approach to the generation of the first tree, we obtain | S_1|· C_k+i_1^2⋯ C_k+i_1+i_2-1^2 candidate sequences. Subsequently, we evaluate whether they meet two conditions:

	
  * T_i_1+1≺ y_2, T_i_1+2≺ y_2, ⋯, T_i_1+i_2-1≺ y_2
	
  * Whether the 2× k matrix associated with y_1, y_2 satisfies the MDS condition, i.e., all submatrix determinants are nonzero.


The set of sequences that fulfill the above two conditions is denoted as S_2. By extension, all subsequent trees in the implementation tree hierarchy can be generated analogously. This iterative process terminates in two scenarios:

	
  * At row p, if S_p = ∅, it indicates the absence of an implementation tree corresponding to the type (i_1, i_2-i_1, ⋯, i_k-i_k-1) (or, more precisely, the absence of an implementation tree for types whose first p terms equal (i_1, i_2-i_1, ⋯, i_p-i_p-1)).
	
  * If S_k ≠∅, then we obtain a successful construction of an implementation tree for the MDS matrix of the given type. Since the search is exhaustive, we can obtain all implementation trees corresponding to the type (i_1, i_2-i_1, ⋯, i_k-i_k-1).


Finally, by iterating over all possible values of (i_1, i_2-i_1, ⋯, i_k-i_k-1) in I_i_k, we either obtain all implementation trees with capacity i_k or determine that no such trees exist. At this point, we increment i_k by 1 and repeat all the above procedures.





 §.§ Feasible structures


We executed Algorithm <ref> on matrices of sizes 3× 3, 4× 4, and 5× 5 defined over the ring R_n and obtained the following results.



Of particular interest is the 4× 4 MDS matrix over R_n. Guided by Theorem <ref>, our traversal search began with an implementation tree of capacity 7. Utilizing Algorithm <ref>, we swiftly discovered, within milliseconds, that no tree of this capacity existed. Subsequently, we tested the trees with a capacity of 8 and successfully located some viable types. A comprehensive list of all the simplest trees with a capacity of 8 is provided in Appendix A.



 §.§ The 4× 4 MDS matrices with the lowest cost on R_8


In this section we establish the ring R_8 as F_2[α], where α denotes the companion matrix associated with the polynomial α^8 + α^2 + 1. Explicitly, α is represented by the following matrix:

    α =
    	[ 0 0 0 0 0 0 0 1; 1 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 1; 0 0 1 0 0 0 0 0; 0 0 0 1 0 0 0 0; 0 0 0 0 1 0 0 0; 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 1 0;   ].


The choice of this ring is motivated by the fact that it contains a unique zero divisor, α^4 + α + 1, and the multiplication by α has a minimal cost of just one xor operation.

To construct MDS matrices with the lowest possible cost, we must assign values to the simplest trees described in section 3.3. During the assignment process, we restrict the values of all parameters to the set {1, α, α^2, α^3, α^-1, α^-2, α^-3}, whose elements have the cost of no more than 2 xor operations. We go through all cases with low costs and obtain 60 distinct 4 × 4 MDS matrices over F_2[α] with 67 xor operations. 30 of them are listed in the following tables, and another 30 MDS matrices can be got by replacing α with α^-1.








		

	
	
	
	
	
	


	The two optimal 4× 4 matrices, as cited in <cit.>, are encompassed within matrices labeled 21 and 22; the ten matrices documented in <cit.> are integrated into matrices numbered 17, 19, 20, 21, and 22; the 52 matrices detailed in <cit.> are represented in the matrices numbered 1-8 and 13-30. The rest eight matrices labeled 9-12 represent newly reported findings.


Furthermore, if we replace α with the companion matrix of the polynomial x^4 + x + 1, we can naturally obtain 60 new MDS matrices with a word length of 4 that require only 35 XOR operations for their implementation.

It is important to note that substituting α with any invertible 8 × 8 matrix over the finite field F_2, whose minimal polynomial is α^8 + α^2 + 1, might yield some new results. However, as long as the substituted matrix can be implemented using a single xor operation, the resulting matrices will retain their MDS property and can still be implemented using 67 xor operations. It should be emphasized that the matrix obtained through such a substitution is not equivalent to the original matrix. Furthermore, while there exist numerous substitution methods that satisfy these conditions, we refrain from enumerating them here due to their extensive nature.

By slightly relaxing the requirements on the implementation cost, our method can yield a large number of 4× 4 MDS matrices that can be implemented with 36-41 XOR gates (for 4-bit input) and 68-80 XOR gates (for 8-bit input). If we limit the circuit depth, we can obtain MDS matrices with optimal costs for different depths. These matrices are included in Appendix B.



 §.§ Construction of 4× 4 Involutory MDS Matrices


In this section, we endeavor to construct lightweight involutory MDS matrices by assigning specific parameters to the eight simplest trees obtained in Section 3.3. As a result, we have identified six involutory MDS matrices that can be efficiently implemented using only 68 XOR gates. To the best of our knowledge, this represents the current optimal result in terms of implementation cost.

Since row permutations and row multiplications affect the involutory properties of a matrix, we must consider these operations during our search. Specifically, when searching for an involutory matrix, we must account for the original 8× 2 = 16 parameters, as well as additional 4 parameters to represent row multiplications and 4! = 24 possible row orders.


	It is important to note that the order of the matrix columns does not need to be considered. If P and Q are permutation matrices such that (PAQ)^2 = I, then right-multiplying both sides of the equation by Q^-1 and left-multiplying by Q yields QPAQPA = I. This means that if there exists an involutory matrix in the equivalence class of A, it can be obtained simply by swapping the rows of A.


We define the equivalence class of an involutory MDS matrix as follows: If A is an involutory MDS matrix and A' = P^-1AP, where P is a permutation matrix, then A and A' belong to the same equivalence class. It is evident that the implementation of A' differs from that of A only in the order of the input and output variables, and both matrices are either involutory or non-involutory.

Our search is restricted to the ring F[x]/(x^8 + x^2 + 1), and the parameters are powers of α, where α is the companion matrix of x^8 + x^2 + 1. We calculate the cost of α^n as |n| and consider the sum of these costs as the total cost t. Scalar multiplications can be reused, potentially reducing the practical cost below t. During the search, we select s (where s ≤ 6) out of 20 parameters for assignment, fixing the unselected parameters at 1. We limit the total cost of the parameters to t ≤ 8.

Under these conditions, we find that only the simplest trees 3 and 4 can produce involutory matrices. These involutory matrices satisfy s ≥ 4 and t ≥ 5. Focusing on the lightest involutory MDS matrices with t = 5, we identify 6 matrices from the simplest tree 3 and 12 matrices from the simplest tree 4. We provide these 18 matrices satisfying t = 5 below.

Simplest tree 3:

    y_1=(a_5(a_3(a_0x_1⊕ a_1x_2)⊕ a_2x_3)⊕ a_4x_4)· a_16,
    
    		y_2=(a_10x_3⊕ a_11(a_9(a_7(a_3(a_0x_1⊕ a_1x_2)⊕ a_2x_3)⊕ a_6x_1)⊕ a_8y_1))· a_17,
    
    		y_3=(a_12(a_0x_1⊕ a_1x_2)⊕ a_13(a_9(a_7(a_3(a_0x_1⊕ a_1x_2)⊕ a_2x_3)⊕ a_6x_1)⊕ a_8y_1))· a_18,
    
    		y_4=(a_15y_3⊕ a_14(a_7(a_3(a_0x_1⊕ a_1x_2)⊕ a_2x_3)⊕ a_6x_1))· a_19.


The six MDS matrices, associated with the simplest tree 3 and adhering to the condition t=5, are designated as follows:

	
  * a_6=α,a_7=α^2,a_12=α^-1,a_13=α.
	
  * a_3=α,a_9=α^2,a_12=α^-1,a_14=α.
	
  * a_3=α,a_6=α,a_7=α,a_9=α,a_12=α^-1.
	
  * a_7=α,a_9=α,a_12=α^-1,a_13=α,a_14=α.
	
  * a_0=α^-1,a_6=α,a_7=α,a_12=α^-1,a_15=α^-1.
	
  * a_0=α^-1,a_9=α,a_12=α^-1,a_14=α,a_15=α^-1.

The aforementioned six matrices exhibit a specific order among their rows, which is as follows:

    P=[ [ 0 0 1 0; 0 0 0 1; 1 0 0 0; 0 1 0 0 ] ].


The resulting matrix M is expected to have the order (y_3, y_4, y_1, y_2). Among the six involutory MDS matrices with t=5 mentioned previously, the fourth and sixth assignments demonstrate scalar multiplication reuse, which are listed below. Since the terms multiplied by these scalars are identical, the multiplication cost is calculated only once. Therefore, the total cost of the entire matrix operation is determined to be 8 × 8 + 4 = 68.

    4:[ [ α^3+α^2+α+α^-1     α^3+α+α^-1          α^3+α              α;       α^3+α^-1 α^3+α^2+α+α^-1      α^3+α^2+α              α;              1              1              1              1;        α^2+α+1          α^2+1            α^2              1 ] ],



    6:[ [    α+1+α^-1+α^-2         α+1+α^-1              α+1                1; α+α^-1+α^-2+α^-3    α+1+α^-1+α^-2         α+1+α^-1             α^-1;             α^-1                1                1                1;         α+1+α^-1              α+1                α                1 ] ].


Simplest tree 4:

    y_1=(a_4(a_0x_1⊕ a_1x_2)⊕ a_5(a_2x_3⊕ a_3x_4))· a_16, 
    
    		y_2=(a_11(a_9(a_6x_1⊕ a_7(a_2x_3⊕ a_3x_4))⊕ a_8x_3)⊕ a_10(a_0x_1⊕ a_1x_2))· a_17, 
    
    		y_3=(a_12y_1⊕ a_13(a_9(a_6x_1⊕ a_7(a_2x_3⊕ a_3x_4))⊕ a_8x_3))· a_18, 
    
    		y_4=(a_15y_2⊕ a_14(a_6x_1⊕ a_7(a_2x_3⊕ a_3x_4)))· a_19.


Here we present the twelve assignments of MDS matrices corresponding to the simplest tree 4 satisfying the condition (t=5). These assignments are enumerated as follows:

	
  * a_4=α,a_5=α^-1,a_8=α,a_9=α^2.
	
  * a_6=α,a_7=α^2,a_10=α^-1,a_11=α.
	
  * a_5=α^-1,a_7=α,a_11=α^2,a_13=α.
	
  * a_4=α,a_9=α^2,a_10=α^-1,a_14=α.
	
  * a_4=α,a_6=α,a_7=α,a_9=α,a_10=α^-1.
	
  * a_5=α^-1,a_7=α,a_8=α,a_9=α,a_11=α.
	
  * a_2=α^-1,a_5=α^-1,a_8=α,a_9=α,a_12=α^-1.
	
  * a_4=α,a_5=α^-1,a_9=α,a_11=α,a_13=α.
	
  * a_2=α^-1,a_5=α^-1,a_11=α,a_12=α^-1,a_13=α.
	
  * a_7=α,a_9=α,a_10=α^-1,a_11=α,a_14=α.
	
  * a_0=α^-1,a_6=α,a_7=α,a_10=α^-1,a_15=α^-1.
	
  * a_0=α^-1,a_9=α,a_10=α^-1,a_14=α,a_15=α^-1.

The order between the rows for the above assignments is given by the permutation matrix:

    P=[ [ 0 0 1 0; 0 0 0 1; 1 0 0 0; 0 1 0 0 ] ],

which implies that the resulting matrix M should have the order (y_2, y_4, y_1, y_3).

Among the twelve involutory MDS matrices presented above, the 8th, 9th, 10th, and 12th assignments exhibit scalar multiplication reuse and have the cost of 68:


    [ [      α^2+1          1      α^2+α        α^2;        α^2          1    α^2+α+1      α^2+1;          α          α       α^-1       α^-1;      α^2+α          α α^2+α+α^-1   α^2+α^-1 ]]    [ [      α+1        1      α+1        α;        α        1 α+1+α^-1      α+1;        1        1     α^-2     α^-1;   α+α^-1     α^-1 α+1+α^-3   α+α^-2 ] ], 
    
    		(8)                                         (9)



    [ [   α^2+α^-1       α^-1      α^3+α        α^3; α^2+α+α^-1       α^-1  α^3+α^2+α    α^3+α^2;          1          1          1          1;        α+1          1        α^2      α^2+1 ] ]    [ [   α+α^-2     α^-1      α+1        α; α+1+α^-3     α^-2 α+1+α^-1      α+1;     α^-1        1        1        1;   α+α^-1        1        α      α+1 ] ]. 
    
    		(10)                                         (12)




§ CONSTRUCTION OF HIGHER ORDER MDS MATRICES


Due to the high efficiency of our algorithm, we are able to search for higher-order MDS matrices. Below, we present the implementation tree structures for MDS matrices of order 5 and 6, along with the actual matrices. Specifically, we demonstrate a 5th-order MDS matrix with a word length of 8 that can be implemented using 114 xor operations, and a 6th-order MDS matrix with a word length of 8 that can be implemented using 148 xor operations.



 §.§ 5th-order MDS Matrix


Using Algorithm <ref>, we can obtain the implementation tree of a 5th-order MDS matrix. Figure <ref> illustrates such an implementation tree for a 5× 5 MDS matrix.



When dealing with matrices of higher orders, we arrive the case that the number of parameters to be assigned increases significantly, so we necessitate a filtering process. Consider the implementation tree of the 5× 5 MDS matrix depicted in Figure <ref> as an example. This tree is comprised of 12 XOR operations, and totally 24 parameters (a_1, a_2, …, a_24) need to be assigned, where each a_i represents an n× n matrix over F_2 and n denotes the word length:

    x_6 = a_1· x_1 ⊕ a_2· x_2, x_7 = a_3· x_3 ⊕ a_4· x_4, x_8 = a_5· x_6 ⊕ a_6· x_7, x_9 = a_7· x_5 ⊕ a_8· x_8 = y_1,


    x_10 = a_9· x_1 ⊕ a_10· x_3, x_11 = a_11· x_5 ⊕ a_12· x_8, x_12 = a_13· x_2 ⊕ a_14· x_11,


    x_13 = a_15· x_7 ⊕ a_16· x_12 = y_2 ,


    x_14 =a_17· x_4 ⊕ a_18· x_11, x_15 = a_19· x_6 ⊕ a_20· x_14 = y_3,


    x_16 = a_21· x_8 ⊕ a_22· x_12 = y_4,


    x_17 = a_23· x_8 ⊕ a_24· x_14 = y_5,


To minimize the implementation cost of the MDS matrix, we aim to maximize the number of identity matrices among the assigned parameters (a_1, a_2, …, a_24). This simplification process involves searching for the minimum number of non-identity matrices required to form an MDS matrix implementation tree. Initially, we assume that at least s non-identity matrices are needed. We start with s = 1 and iterate through all C_24^s possible combinations of positions for these non-identity matrices. If a valid MDS matrix can be formed with this assignment, we have a simplified implementation tree. Otherwise, we increment s and repeat the process until a valid tree is found.

In practice, it is found that an implementation tree for a 5× 5 MDS matrix can be constructed using only 5 non-identity matrices among the 24 positions. Such a tree is given by the following equations:

    x_6    = x_1 ⊕ b_1 · x_2,     x_7    = x_3 ⊕ x_4, 
    
    	x_8    = b_2 · x_6 ⊕ b_3 · x_7,     x_9    = x_5 ⊕ x_8 = y_1, 
    
    	x_10   = x_1 ⊕ x_3,     x_11   = x_5 ⊕ x_8, 
    
    	x_12   = x_2 ⊕ b_4 · x_11,     x_13   = x_7 ⊕ x_12 = y_2, 
    
    	x_14   = x_4 ⊕ b_5 · x_11,     x_15   = x_6 ⊕ x_14 = y_3, 
    
    	x_16   = x_8 ⊕ x_12 = y_4,     x_17   = x_8 ⊕ x_14 = y_5.


The algorithm for simplifying the MDS matrix implementation tree structure is outlined in Algorithm <ref>.



From the preceding example, it is evident that Algorithm <ref> significantly reduces the number of parameters that require assignment. For these reduced parameters, we can directly assign values. Consider the implementation tree of the 5× 5 MDS matrix presented in the previous section (Figure <ref>) as a case study. Assuming a word length of 8, each parameter takes an 8× 8 matrix over F_2.

Initially, we select the matrix α as our basis, whose minimum polynomial is x^8+x^6+x^5+x^3+1:

	
    α =
    		[ 0 1 0 0 0 0 0 1; 1 0 0 0 0 0 0 0; 0 1 0 0 1 0 0 0; 0 0 1 0 0 0 0 0; 0 0 0 1 0 0 0 0; 0 0 0 0 1 0 0 0; 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 1 0;   ].


This polynomial is primitive of degree 8 and requires only 2 xor operations for implementation. To minimize the cost of matrix implementation, we opt for parameter values from lower powers of α, such as α^-4,α^-3,α^-2,α^-1,α^1,α^2,α^3,α^4. For a given 5× 5 MDS matrix, we need to go through 8^5 = 2^15 possible assignments.



Running Algorithms <ref> and <ref>, we get a concrete implementation of the 5× 5 MDS matrix on GL(2,8):

    x_6    = x_1 ⊕α· x_2,     x_7    = x_3 ⊕ x_4, 
    
    	x_8    = α^-2· x_6 ⊕α^2 · x_7,     x_9    = x_5 ⊕ x_8 = y_1, 
    
    	x_10   = x_1 ⊕ x_3,     x_11   = x_5 ⊕ x_8, 
    
    	x_12   = x_2 ⊕α^-3· x_11,     x_13   = x_7 ⊕ x_12 = y_2, 
    
    	x_14   = x_4 ⊕α· x_11,     x_15   = x_6 ⊕ x_14 = y_3, 
    
    	x_16   = x_8 ⊕ x_12 = y_4,     x_17   = x_8 ⊕ x_14 = y_5,

where α is as defined earlier, with a minimal polynomial of x^8+x^6+x^5+x^3+1. The overall implementation cost is 12 × 8 + 9 × 2 = 114 xor gates.



 §.§ 6th-order MDS Matrix


Similarly to the 5th-order matrix, we do the same for a 6th-order matrix. An implementation tree of 6th-order MDS matrix is depicted below:

Its corresponding matrix is:

    x_7 =  x_1 ⊕  x_2, x_8 =  x_3 ⊕α· x_7, x_9 =  x_4 ⊕α^3· x_8, x_10 =  x_5 ⊕α^-1· x_9, x_11 =  x_6 ⊕  x_10 = y_1,


    x_12 =  x_6 ⊕  x_7, x_13 =  x_9 ⊕  x_12, x_14 =  x_2 ⊕α^-1· x_10, x_15 =  α^2· x_13⊕  x_14,


    x_16 =  x_4 ⊕  x_15 = y_2,


    x_17 = x_3 ⊕ x_15, x_18 =  x_5 ⊕α^2· x_17, x_19 =  x_14⊕  x_18 = y_3,


    x_20 =  x_12⊕  x_17 = y_4,


    x_21 =  x_8 ⊕  x_18 = y_5,


    x_22 =  x_13⊕  x_21 = y_6,

where

	
    α=[ [ 0 1 0 0 0 0 0 1; 1 0 0 0 0 0 0 0; 0 1 0 0 1 0 0 0; 0 0 1 0 0 0 0 0; 0 0 0 1 0 0 0 0; 0 0 0 0 1 0 0 0; 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 1 0;   ]]


and its minimal polynomial is α is x^8+x^6+x^5+x^3+1. The cost of the entire implementation is 16× 8 + 10× 2 = 148 xor gates.



§ CONLUSION


In this paper we present a traversal algorithm tailored for the discovery of lightweight 4× 4 MDS matrices. As results, we successfully generate all implementation trees for 4× 4 MDS matrices using only 8 word-XOR operations. Based on this systematic approach, we derive a series of MDS matrices that exhibit the lowest computational cost. Additionally, we obtain the lowest-cost 4× 4 involutory MDS matrix currently known over the finite field GL(2,8). Looking beyond the scope of 4× 4 matrices, we extend our method to higher-order MDS matrices. Specifically, we employ our method to obtain the lowest-cost 5× 5 and 6× 6 MDS matrices as well.




	
	

§ THE SIMPLEST TREES WITH CAPACITY 8

	
	Tree 1-5 are have the type of (3, 3, 1, 1) and tree 6-8 have the type (4, 2, 1, 1).
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	

§ THE 4× 4 MDS MATRICES FOR DIFFERENT DEPTHS

	
	
	
	
	In Table <ref>, α represents the companion matrix of either x^8+x^2+1 or. The 4× 4 MDS matrices presented in Table <ref> for depth 4 can be implemented using 69 XOR operations (for 8-bit input) or 37 XOR operations (for 4-bit input), representing the most cost-effective solution at this depth.
	
	As the depth decreases to 3, the minimum cost increases to 77 XOR operations (for 8-bit input) or 41 XOR operations (for 4-bit input). Owing to the strict limitation on depth, the available parameter options become constrained. For depth 3, there are three matrices that achieve the lowest cost. These matrices are enumerated below.
	
    1:[ [      1      1   α^-1 α^-1+α;      1    1+α      α      α;      α   α^-1 1+α^-1      1;    1+α      1      1    1+α ] ]

	
	
    x_5 =  x_1 ⊕  x_2, x_6 = α· x_4 ⊕  x_5, x_7 =  x_3 ⊕  x_4,  x_8 =  x_6 ⊕α^-1· x_7 = y_1.

	
    x_9 =    x_2 ⊕  x_3,  x_10 = x_6 ⊕α· x_9 = y_2.

	
    x_11 =  α· x_1 ⊕  x_7, x_12 = α^-1· x_9 ⊕  x_11 = y_3.

	
    x_13 =   x_6 ⊕  x_11 = y_4.

	
	
    2:[ [      1      1      α α^-1+α;      1 1+α^-1   α^-1   α^-1;   α^-1      α    1+α      1; 1+α^-1      1      1 1+α^-1 ] ]

	
	
    x_5 =  x_1 ⊕  x_2, x_6 = α^-1· x_4 ⊕  x_5, x_7 =  x_3 ⊕  x_4,  x_8 =  x_6 ⊕α· x_7 = y_1

	
    x_9 =    x_2 ⊕  x_3,  x_10 = x_6 ⊕α^-1· x_9 = y_2

	
    x_11 =  α^-1· x_1 ⊕  x_7, x_12 = α· x_9 ⊕  x_11 = y_3

	
    x_13 =   x_6 ⊕  x_11 = y_4

	
	
    3:[ [   α^-1      α      1   α^-1; α^-1+1    1+α      1      1;      1      1 α^-1+1 α^-2+1;    1+α      1    1+α      1 ] ]

	
	
    x_5 = α^-1· x_1 ⊕α· x_2, x_6 =  x_3 ⊕α^-1· x_4, x_7 =  x_5 ⊕  x_6=y_1

	
    x_8 =    x_1 ⊕  x_3,  x_9 = x_2 ⊕α^-1· x_4, x_10 =  x_8 ⊕  x_9, x_11 =  x_5 ⊕  x_10 = y_2

	
    x_12 = α^-1· x_6 ⊕  x_10 = y_3

	
    x_13 =  α· x_8 ⊕  x_11 = y_4




