A new hybrid method combining search and direct based construction ideas to generate all 4 × 4 involutory maximum distance separable (MDS) matrices over binary field extensions

This article presents a new hybrid method (combining search based methods and direct construction methods) to generate all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $4 \times 4$\end{document}4×4 involutory maximum distance separable (MDS) matrices over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathbf{F}_{2^m}$\end{document}F2m. The proposed method reduces the search space complexity at the level of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$\sqrt n $$\end{document}n, where n represents the number of all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $4 \times 4$\end{document}4×4 invertible matrices over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathbf{F}_{2^m}$\end{document}F2m to be searched for. Hence, this enables us to generate all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $4 \times 4$\end{document}4×4 involutory MDS matrices over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathbf{F}_{2^3}$\end{document}F23 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathbf{F}_{2^4}$\end{document}F24. After applying global optimization technique that supports higher Exclusive-OR (XOR) gates (e.g., XOR3, XOR4) to the generated matrices, to the best of our knowledge, we generate the lightest involutory/non-involutory MDS matrices known over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathbf{F}_{2^3}$\end{document}F23, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathbf{F}_{2^4}$\end{document}F24 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathbf{F}_{2^8}$\end{document}F28 in terms of XOR count. In this context, we present new \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $4 \times 4$\end{document}4×4 involutory MDS matrices over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathbf{F}_{2^3}$\end{document}F23, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathbf{F}_{2^4}$\end{document}F24 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathbf{F}_{2^8}$\end{document}F28, which can be implemented by 13 XOR operations with depth 5, 25 XOR operations with depth 5 and 42 XOR operations with depth 4, respectively. Finally, we denote a new property of Hadamard matrix, i.e., (involutory and MDS) Hadamard matrix form is, in fact, a representative matrix form that can be used to generate a small subset of all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $2^k\times 2^k$\end{document}2k×2k involutory MDS matrices, where k > 1. For k = 1, Hadamard matrix form can be used to generate all involutory MDS matrices.


INTRODUCTION
Two important properties used in the design of block ciphers and defined by Shannon (1949) are confusion and diffusion.These properties are respectively satisfied by substitution boxes (or shortly S-boxes) and linear transformations in a round function of a block cipher.Maximum distance separable (MDS) matrices are derived from MDS codes and provide the maximum diffusion.They are used as the core component of diffusion layers/linear transformations in the design of cryptographic primitives like block ciphers and hash functions.MDS matrices have the maximum branch number, which is an important cryptographic criterion used for defining diffusion rate.Branch number also helps to measure security against some well-known attacks like differential (Biham & Shamir, 1991) and linear cryptanalysis (Matsui, 1994).Block ciphers (or a cryptographic primitive) generally use MDS matrices over F 2 3 , F 2 4 and F 2 8 in their diffusion layers according to the design strategies and considering implementation issues of a block cipher.Also, using involutory MDS matrices in the design of a diffusion layer of a block cipher has the advantage of reusing the same circuit in the decryption process and helps to implement a block cipher at close encryption and decryption costs.In this context, we present our experimental results by considering the finite fields F 2 3 , F 2 4 and F 2 8 .Note that our construction technique can easily be used for designing new involutory MDS matrices over finite field F 2 m especially for m 8.
Generally, the construction techniques of MDS matrices can be categorized into three groups: direct construction methods, search based methods, and hybrid methods (combining search based methods and direct construction methods).The direct construction methods include the methods such as Cauchy matrices (Youssef, Mister & Tavares, 1997;Cui, Jin & Kong, 2015), Vandermonde matrices (Sajadieh et al., 2012b) and Companion matrices (Gupta & Ray, 2013).It should also be noted that Cauchy and Vandermonde matrices are generally not efficient for low-cost implementations (Gupta et al., 2019).Recently, unlike the other direct construction methods, a new direct construction method (or a new matrix form) to generate all 3 Â 3 involutory MDS matrices has been given in Güzel et al. (2019).On the other hand, search based methods to find MDS matrices are based on using hybrid structures (Sim et al., 2015), recursive structures (Sajadieh et al., 2012a(Sajadieh et al., , 2012b)), and searching some special matrix forms like circulant and Hadamard matrix forms, which have some advantages in the implementation phase and have a higher probability of finding MDS matrix when compared to a randomized square matrix.Moreover, the other special matrix forms used to find (lightweight) MDS matrices are as follows: circulant-like (Gupta & Ray, 2015), Toeplitz, Toeplitz-like, and Hankel matrices (Gupta et al., 2019).The search based construction methods are useful for finding MDS matrices with small orders.But, finding MDS matrices with higher orders is exactly an NP-complete problem (computation cost for checking a matrix to be MDS is still too expensive) (Sim et al., 2015).To handle this problem, the Generalized Hadamard (shortly GHadamard) matrix form, a hybrid construction method, was proposed in Pehlivanoğlu et al. (2018).Overall, in Sakallı et al. (2020), the authors proposed a complementary method for the current construction methods in the literature, which generates isomorphic k Â k MDS matrices (new MDS matrices from the implementation point of view) from any existing k Â k MDS matrix (due to its ground field structure).All these methods can be evaluated within the local optimization category that focuses on the coefficients of a given matrix.In recent years, global optimization techniques have been proposed to construct smaller diffusion layer circuits for involutory/non-involutory MDS matrices.This challenging problem is also known as the problem of finding the shortest linear straight-line program (SLP) and optimizing circuits globally.In this context, two main global optimization techniques can be given as cancellation-free programs (e.g., Paar's algorithms (Paar, 1997) and heuristic techniques (e.g., Boyar-Peralta (BP) algorithm (Boyar & Peralta, 2010;Boyar, Matthews & Peralta, 2012) RNBP, A1, A2 (Tan & Peyrin, 2019).All these heuristics support 2-input XOR gates (XOR2) only, but in Baksi et al. (2021) the authors by inspiring the idea given in Banik, Funabiki & Isobe (2019) proposed a new version of the original BP heuristic, called BDKCI, that supports higher input XOR gates such as 3-input XOR (XOR3) and 4-input XOR (XOR4) gates.Utilizing higher input XOR gates can result in a lower cost in specific ASIC libraries, so in this article, we use BDKCI heuristic for finding the efficient implementation of a given matrix.Moreover, we consider not only gate count (GC) and circuit depth but also gate equivalent (GE) metric for lightweight implementations.BDKCI heuristic directly calculates the total GEs for each ASIC library (STM 90 nm (ASIC1), STM 65nm (ASIC2), TSMC 65 nm (ASIC3), and STM 130 nm (ASIC4)), more technical details can be found in Baksi et al. (2021).
To the best of our knowledge, there is no known method in the literature to represent/ generate all involutory MDS matrices over F 2 m .In this study, we present a new hybrid method to generate all 4 Â 4 involutory MDS matrices over F 2 m .We consider the problem of finding lightweight involutory/non-involutory MDS matrices with low implementation costs.In this context, we generate all 4 Â 4 involutory MDS matrices over F 2 3 and F 2 4 and evaluate these matrices up to a threshold value with respect to naive XOR count (d-XOR) (Khoo et al., 2014;Jean et al., 2017), and then we look for the lightest hardware circuits known for 4 Â 4 involutory MDS matrices in terms of XOR counts and circuit depths by using BDKCI algorithm.Note that when obtaining the lightest 4 Â 4 non-involutory MDS matrix over F 2 4 , we take the benefit of the elements of the lightest 4 Â 4 involutory MDS matrix over F 2 4 generated.
In this article, we propose a new hybrid method to generate all 4 Â 4 involutory MDS matrices over F 2 m .The proposed hybrid method consists of two parts: the first part includes searching for all 4 Â 4 representative involutory MDS matrices and the second part includes generating all 4 Â 4 involutory MDS matrices using all 4 Â 4 representative involutory MDS matrices found by search in the first part of the proposed method.Hence, we present the number of all 4 Â 4 involutory MDS matrices over F 2 3 and F 2 4 .In addition, we give the lightest known involutory/non-involutory MDS matrices over F 2 3 , F 2 4 and F 2 8 obtained by using the global optimization algorithm BDKCI, which is an improved version of BP algorithm.

Motivation and our contribution
The demand for low-cost security design targeting resource-constrained devices has triggered the exploration of lightweight diffusion layers (diffusion layers with low implementation costs).In particular, the circuit area is a crucial criterion for lightweight cryptography in terms of hardware.This means that reducing the cost of hardware implementation is to minimize the number of expensive logical gates (especially XOR operations) and the depth of the circuit (for low latency), which is the number of gates on the longest circuit path.However, it is not an easy task to construct lightweight MDS matrices, especially involutory ones.Involutory MDS matrices use the same circuit and have the same implementation costs in the encryption and decryption phases.Therefore, the designers have considered the problem of building optimal or close to optimal linear circuits for involutory/non-involutory MDS matrices.
Previous studies mainly focused on the problem of constructing a lightweight MDS matrix from two perspectives: building a locally optimized implementation of an MDS matrix consisting of the coefficients that are easy to evaluate or finding the globally optimized implementation of a given MDS matrix.However, while searching lightweight MDS matrices, most studies use some heuristic searching methods that cannot find all MDS matrices especially involutory ones in a specific finite field (except for the matrix form given in Güzel et al. (2019) used for generating all 3 Â 3 involutory MDS matrices).In this article, for the first time in the literature, we focus on a hybrid method generating all 4 Â 4 involutory MDS matrices over F 2 m with the lightest known implementation costs because many block ciphers use 4 Â 4 and 8 Â 8 (in the form of 2 k Â 2 k ) MDS matrices.For example, the Advanced Encryption Standard (AES) (Daemen & Rijmen, 2002) uses a 4 Â 4 MDS matrix as the main part of its diffusion layer.Moreover, (lightweight) block ciphers to be designed in the future may likely use 4 Â 4 involutory/non-involutory MDS matrices over F 2 m with nice hardware implementation costs.In this article, we combine local optimization (new hybrid construction method) and global optimization (BDKCI algorithm) techniques.The main contributions of this article can be given as follows: A new hybrid construction method to generate all 4 Â 4 involutory MDS matrices over F 2 m is proposed.The proposed method reduces the search space complexity (the search space defining the total number of all invertible 4 Â 4 matrices over F 2 m ) at the level of ffiffiffi n p , where n represents the number of all 4 Â 4 invertible matrices to be searched for.It is shown that the number of all 4 Â 4 involutory MDS matrices over F 2 m can be calculated by the formula #RIM Â ð2 m À 1Þ 3 , where #RIM represents the number of all 4 Â 4 representative involutory MDS matrices.Hence, there are, respectively, 16,464 (¼ 48 Â ð2 3 À 1Þ 3 ) and 242,514,000 (¼ 71; 856 Â ð2 4 À 1Þ 3 ) 4 Â 4 involutory and MDS matrices over F 2 3 and F 2 4 .The new lightest involutory/non-involutory MDS matrices over F 2 3 , F 2 4 and F 2 8 achieved by the proposed hybrid construction method are presented: -A new 4 Â 4 involutory MDS matrix over F 2 4 which can be implemented by only 25 XOR operations with depth 5 is presented, whereas the previously known lightest one (Sarkar & Syed, 2016) requires 26 XOR operations to be implemented with the same depth using XOR2 and XOR3 gates (please see the matrix M 4 ).-A new 4 Â 4 involutory MDS matrix over F 2 3 is presented, which can be implemented by only 13 XOR operations with depth 5 using XOR3 and XOR4 gates (please see the matrix M 3 ).-A new 4 Â 4 non-involutory MDS matrix over F 2 4 which can be implemented by only 19 XOR operations with depth 3 is presented, whereas previously known lightest one (Sarkar & Syed, 2016) costs 20 XOR operations the same depth using XOR2, XOR3 and XOR4 gates (please see the matrix M 5 ).-A new 4 Â 4 involutory MDS matrix over F 2 8 can be implemented by only 52 XOR operations with depth 5 using XOR2 and XOR3 gates.Moreover, the same matrix requires only 42 XOR operations with depth 4 using XOR2, XOR3, and XOR4 gates (please see the matrix M 6 ).While compare to the previous best-known implementations of linear layers of some block ciphers, it improves the results in GCs and GEs.
These findings clearly establish that the circuit implementations of matrices M 3 , M 4 , M 5 and M 6 over the specified fields are superior to those given in previous studies.This aspect also underscores one of the key novelties in our article.
A new property of Hadamard matrix is denoted, i.e., (involutory and MDS) Hadamard matrix form is, in fact, a representative matrix form that can be used to generate a small amount of all 2 k Â 2 k involutory MDS matrices, where k .1.For k ¼ 1, Hadamard matrix form can be used to generate all involutory MDS matrices.For 4 Â 4 Hadamard matrices, we present the matrix form R 1 (which has the generic properties of Hadamard matrix, i.e., XOR sum of the elements in any row or column of a Hadamard matrix is equal to 1 and XOR sum of the elements in the main diagonal is equal to 0) given in "Proposed Method" section for finding all 4 Â 4 representative involutory MDS matrices by search.The matrix form R 1 can also be adaptable to 2 k Â 2 k Hadamard matrices for k .2. Hence, one can generate new representative involutory MDS matrices over any finite field by using the matrix form R 1 or adapted versions of R 1 for 2 k Â 2 k Hadamard matrices for k .2, which may help us find new involutory MDS matrices (also, with the help of b i parameters given in Theorem 1) with better implementation properties.All the optimization results of matrices are available at https://github.com/mkurtpehlivanoglu/Hybrid_Method.

Organization
This article is organized as follows: We give some notations, properties of MDS matrices, and two metrics used for identifying lightweightness of an MDS matrix in the "Preliminaries" section.In the "Proposed Method" section, we propose a new hybrid method to generate all 4 Â 4 matrices over F 2 m .Experimental results on the number of all 4 Â 4 involutory MDS matrices and some examples for 4 Â 4 involutory MDS matrices over F 2 3 , F 2 4 and F 2 8 with the lowest XOR counts are presented in "Experimental Results" section.Finally, we conclude the article in the "Conclusion" section.

PRELIMINARIES
In this section, we describe the mathematical background needed throughout the article.In this context, some notations and definitions are presented.
The finite field F 2 m consisting of 2 m elements is defined by an irreducible polynomial pðxÞ of degree m over F 2 and is denoted by F 2 ½x=ðpðxÞÞ.The elements of the finite field F 2 m can be represented as polynomials over F 2 ½x=pðxÞ, i.e., P mÀ1 i¼0 a i a i , where a i 2 F 2 and a is a root (and a primitive element) of F 2 m .For simplicity, we denote F 2 m defined by irreducible polynomial pðxÞ as F 2 m =pðxÞ (the finite field with 2 m elements) and use hexadecimal notation to represent the elements of F 2 m and the irreducible polynomial pðxÞ used for defining F 2 m .As an example, the four-bit string 1110 which can also be represented by 0xe in hexadecimal notation corresponds to the polynomial a 3 þ a 2 þ a in F 2 4 .In the same manner, 0x13 used when denoting F 2 4 =0x13 stands for the irreducible polynomial where d, n, and k represent the length and the number of rows of the generating matrix of the code C, respectively).Moreover, generator matrices that produce MDS codes are called MDS matrices.MDS matrices derived from MDS codes provide the maximum diffusion in a block cipher and have the maximum differential and linear branch number (k þ 1 for k Â k for MDS matrices), which are two important cryptographic criteria for linear transformations.The followings are the properties of an MDS matrix: Let M be a k Â k square matrix.M is an MDS matrix, if and only if every square submatrix of M is non-singular.If M is a k Â k MDS matrix, the transpose matrix of M (M T ) is also an MDS matrix.
Let c 2 F 2 m be a non-zero constant and let M be a k Â k MDS matrix, then the multiplication of a row (or column) of M by c does not affect the MDS property of M.
If a square matrix A is its own inverse (i.e., A = A À1 ), then the matrix A is called an involutory matrix.Equivalently, the matrix A is an involution if and only if A 2 = I, where I is the identity matrix.Definition 1.A k Â k Finite Field Hadamard matrix (simply Hadamard matrix) H over F 2 m with k ¼ 2 t for t .0 can be expressed as follows: where sub-matrices A 0 and A 1 are also 2 tÀ1 Â 2 tÀ1 Hadamard matrices.
When denoting a k Â k Hadamard matrix, we use the notation had(a 0 , a 1 ,…, a kÀ1 ), where a i 2 F 2 m for 0 i k À 1.In this respect, a 4 Â 4 Hadamard matrix H can be given as follows: (2) Some important properties of a k Â k Hadamard matrix H over F 2 m are as follows (Pehlivanoğlu et al., 2018): given above for a k Â k Hadamard matrix and preserves involutory and MDS properties of a given k Â k involutory and MDS Hadamard matrix.The idea in preserving MDS property of GHadamard matrix form is based on the last property for defining an MDS matrix.A k Â k GHadamard matrix GH is generated by using the combination of nonzero k À 1 more b i parameters and their inverses with a k Â k Hadamard matrix H over F 2 m .In this context, a 4 Â 4 GHadamard matrix GH can be denoted as follows: We estimate the hardware cost of an (involutory) MDS matrix (i.e., linear transformation) with the number of XOR operations required in hardware implementation which can be described as x i x a i È x b i with a i ; b i , i, where x a i and x b i are inputs and some subset of x i s are outputs.In the literature, there are two important approximations of the implementation cost in terms of XOR count: direct XOR (d-XOR) count (Khoo et al., 2014) and sequential XOR (s-XOR) count (Beierle, Kranz & Leander, 2016).While d-XOR count is defined as the Hamming weight (the number of 1 bits) of the corresponding mk Â mk binary matrix (transformed from k Â k matrix over F 2 m ) minus mk, s-XOR count is defined as the minimum number of XOR operations needed to implement the mk Â mk binary matrix with in-place operations and without extra intermediate computations.Although s-XOR count seems like a better approximation, it causes a high computational cost for optimizing full MDS matrices (Duval & Leurent, 2018).
Local and global optimization techniques are used to find optimized implementations of MDS matrices in terms of the required number of XOR operations.The main difference between these two techniques is that the global optimization technique focuses on the optimization of a linear Boolean function of a whole matrix circuit while the local optimization technique focuses on the evaluation of diffusion matrix coefficients.Local optimization techniques do not guarantee finding efficient circuit implementation of an MDS matrix and therefore, more recently, global optimization techniques have been addressed to find well optimized circuits.The idea is based on the reduction of XOR count which is extracted by the naive implementation of an MDS matrix that contains a lot of repeated calculations.
In this article, while performing global optimization of an (involutory) MDS matrix, the circuit constructed by only XOR gates is handled as a linear Boolean function consisting of n input signals fx 0 ; x 1 ; . . .; x n g and m output signals fy 0 ; y 1 ; . . .; y m g (also called target signals).We also use BDKCI algorithm to build a globally optimized implementation of a 4 Â 4 (involutory) MDS matrix generated by using the proposed new hybrid method.The aim of BDKCI algorithm is to find efficient circuits with intermediate variables t i s (calculated once) for other computation sequences.These circuits can reuse intermediate variables t i s that lead to reducing the number of gates required.More details about BDKCI heuristic can be found in Baksi et al. (2021).

PROPOSED METHOD
In this section, we present a new method to generate all 4 Â 4 involutory MDS matrices over F 2 m .The proposed method is a hybrid construction method and is based on the combination of search based and direct based construction methods.The main idea of the proposed method is first to generate all 4 Â 4 representative involutory MDS matrices by search, and then to obtain all 4 Â 4 involutory MDS matrices by applying three more nonzero parameters and their inverses to these representative matrices.When generating representative involutory MDS matrices, we take the benefit of generic properties of a Hadamard matrix satisfying the involutory property, i.e., XOR sum of the elements in any row or column of a Hadamard matrix is equal to 1 and XOR sum of the elements in the main diagonal is equal to 0. Note that these properties also force XOR sum of the elements in the anti-diagonal (counter diagonal) to be equal to 0. Then, we can easily define the matrix form R 1 in order to search for all representative involutory MDS matrices over F 2 m as follows: The matrix form R 1 above for finding representatives involutory MDS matrix is defined by eight parameters ðr 11 ; r 12 ; r 13 ; r 21 ; r 22 ; r 31 ; r 32 ; r 33 Þ over F 2 m À f0g.Note that when defining the matrix form R 1 , the finite field element 0 is not considered because of the fact that all entries of an MDS matrix should be non-zero by the first property used for defining an MDS matrix.Hence, the search space for finding representative involutory MDS matrices over F 2 m is approximately obtained as ð2 m À 1Þ 8 (by omitting non-invertible or singular matrices).For example, for finding all 4 Â 4 representative involutory MDS matrices over F 2 4 , the search space is approximately ð2 4 À 1Þ 8 % 2 31:25 .Before proceeding with the details on the matrix form R 1 , we introduce the matrix form A 2Â2 used for generating all 2 Â 2 involutory MDS matrices over F 2 m in Remark 2 and Lemma 3 from which the main idea of the article comes.
Remark 2. Consider 2 k Â 2 k involutory MDS matrices, where k is a positive integer.For k ¼ 1, we show in Lemma 3 that one can directly generate all 2 Â 2 involutory MDS matrices by the following matrix form: Proof.Let A ¼ r 11 r 12 r 21 r 22 !be a 2 Â 2 involutory matrix with the restriction r 11 6 ¼ 0. Let c ij denote the elements of A 2 for i; j 2 f1; 2g, i.e., c ij ¼ P 2 k¼1 r ik r kj .Since A 2 ¼ I, where I is the 2 Â 2 identity matrix, we obtain the following equations (by considering if i ¼ j then c ij ¼ 1 and if i 6 ¼ j then c ij ¼ 0): After adding the Eqs.( 4) and ( 7) given above, we obtain the equation r 2 11 ¼ r 2 22 .The equation r 2 11 ¼ r 2 22 can be rewritten as ðr 11 þ r 22 Þ 2 ¼ 0, and therefore we obtain r 11 ¼ r 22 because the operations are performed in the finite field F 2 m .On the other hand, we have 2 from the Eq. ( 4).Then, r 12 and r 21 depending on r 11 and a new parameter b 1 are respectively obtained as r 12 ¼ ð1 þ r 11 Þb 1 and r 21 ¼ ð1 þ r 11 Þb À1 1 with the restrictions b 1 2 F 2 m À f0g and r 11 2 F 2 m À f0; 1g so that the matrix form can also satisfy the MDS property (i.e., the determinant of the matrix form A 2Â2 is not equal to 0).Remark 4 Lemma 3 shows that the parameter b 1 and its inverse (b À1 1 ) are hidden parameters keeping involutory and MDS property for the matrix form A 2Â2 , and also that involutory (MDS) matrices are formed as distinct classes in the search space.Similarly, generating all 4 Â 4 involutory MDS matrices becomes an easier problem i.e., the problem is first to focus on searching for and finding all representative involutory MDS matrices by using the matrix form R 1 , and then generate all 4 Â 4 involutory MDS matrices by using representative involutory MDS matrices obtained and the hidden parameters (b 1 , b 2 and b 3 with their inverses b À1 1 , b À1 2 and b À1 3 ) given in Theorem 8.In Lemma 6, we give the mathematical background needed to define matrix form R 1 used for finding representative involutory MDS matrices, which constitutes the search side of the proposed method and we give the characteristics of 4 Â 4 representative involutory MDS matrices with Definition 5.Then, in Theorem 8, we present the mathematical background needed for the direct construction side of the proposed method, which is based on three non-zero parameters preserving involutory and MDS property of a 4 Â 4 representative involutory MDS matrix given.
By Lemma 6, representative involutory MDS matrices are found by searching through all possible candidates using the matrix form R 1 .This allows us to use 8 parameters (r 11 , r 12 , r 13 , r 21 , r 22 , r 31 , r 32 and r 33 ) (instead of 16 parameters needed for defining all 4 Â 4 matrices over F 2 m ) in the search phase of the matrix R 1 and thus enables us to reduce the search space complexity from ð2 m À 1Þ 16 to ð2 m À 1Þ 8 , which is approximately at the level of ffiffiffi n p , where n represents the number of all invertible 4 Â 4 matrices over F 2 m .In this context, we present the matrix form R 1 again below: Definition 5. A representative involutory MDS matrix (or shortly RIM) is a 4 Â 4 involutory matrix R and satisfies the following conditions: XOR sum of all elements of its main diagonal is equal to 0, XOR sum of any rows (and columns) of R is equal to 1 and the MDS property given in "Preliminaries" section.Lemma 6.A representative involutory matrix satisfies the first two conditions given in Definition 5 and these two conditions used to define the matrix form R 1 guarantee to find all 4 Â 4 representative involutory matrices.Proof.Let R ¼ ½r ij be a 4 Â 4 involutory and representative matrix such that all r ij 6 ¼ 0 and let c ij ¼ P 4 k¼1 r ik r kj denote the elements of R 2 for i; j ¼ 1; 2; 3; 4. Then R 2 ¼ I, where I is the 4 Â 4 identity matrix.If i ¼ j, then c ij ¼ 1.Otherwise, c ij ¼ 0. Hence, the following equations are satisfied: By adding the Eqs.( 8)-( 11) side by side, the following equation is obtained: Because we study in the finite field F 2 m the Eq. ( 24) can be rewritten as follows: Hence, we verify the first statement of the lemma.To show the proof of the second statement, we assume that the sum of elements of any row of R is equal to 1.We prefer to study with rows instead of columns of R without breaking the generality.Then, the equations corresponding to our assumption are obtained as follows: By adding the Eqs.( 12)-( 14) side by side, we get the following equality: From the assumption Eqs. ( 26)-( 28), it is clear that: By replacing the expressions r 12 þ r 13 þ r 14 , r 22 þ r 23 þ r 24 , r 32 þ r 33 þ r 34 and r 42 þ r 43 þ r 44 in the Eq. ( 30) with the expressions r 11 þ 1, r 21 þ 1, r 31 þ 1 and r 41 þ 1, respectively, the Eq. ( 35) and then the Eq. ( 36) are obtained: We can rewrite the Eq. ( 36) as follows: Then, by using the assumption (26), we can obtain the Eq. ( 8) as follows: In a similar manner, the Eqs.( 9)-( 11) can be obtained by using the equations from ( 15) to ( 23) and the assumptions ( 26)-( 29).
Similar assumptions can be given when proving the second part of the lemma, related to the sum of the elements of any column, as follows: One can easily show that the Eqs.( 8)-( 11) can again be obtained by using the equations from ( 12) to ( 23) and the assumption Eqs. ( 39)-( 42).For instance, the Eq. ( 8) can be obtained by applying similar operations shown in proving the first part of the lemma and by using the assumption equations from (39) to (42) and the Eqs.( 15), ( 18) and ( 21).Hence, the second part of the lemma is satisfied.Remark 7. Lemma 6 does not guarantee that any matrix satisfying assumptions is exactly a representative involutory matrix.Theorem 8. Let RIM ¼ ½r ij be a 4 Â 4 representative involutory MDS matrix, then the matrix form A ¼ ½a ij obtained by the matrix RIM and some parameters (b i parameters) is also involutory MDS in the following form: where b i s for i 2 f1; 2; 3g are the elements of F 2 m À f0g.
Proof.Since the representative MDS matrix RIM is involutory, it satisfies all equations from (8) to ( 23) given in Lemma 6.If we square the matrix A ¼ ½a ij , we obtain the following expressions corresponding to c ij 's, where c ij denote the elements of A 2 for i; j 2 f1; 2; 3; 4g: needed to define a 4 Â 4 finite field matrix, which makes the search space of a conventional search ð2 m À 1Þ 16 .Thus, our hybrid method reduces the search space complexity from approximately ð2 m À 1Þ 16 to ð2 m À 1Þ 8 , which is approximately at the level of ffiffiffi n p , where n represents the number of all invertible 4 Â 4 matrices over F 2 m .By Theorem 8, all 4 Â 4 involutory MDS matrices are generated directly by using representative involutory MDS matrices found by search in the previous stage and b i parameters.
In Examples 10 and 11, we obtain 4 Â 4 involutory MDS matrices over F 2 4 by using the proposed method, which is also given in the literature.Example 10.Let F 2 4 be generated by the primitive element a which is a root of the primitive polynomial x 4 þ x þ 1 (0x13).Consider the 4 Â 4 h-circulant involutory MDS matrix M 1 recently given in Cauchois & Loidreau (2019).
1 a 14 a 7 a 14 a 2 1 a 13 a 11 a 13 a 4 1 1 a 7 a 11 a 8 2 6 6 4 3 7 7 5 over F 2 4 =0x13.In fact, the involutory MDS matrix M 1 belongs to a class of which representative involutory MDS matrix is as follows: RIM 1 ¼ a a 7 a 5 a 11 a 7 a 2 a 14 a 10 a 5 a 14 a 4 a 13 a 11 a 10 a 13 a 8 2 6 6 4 3 7 7 5 which is symmetric and is also 4 Â 4 h-circulant involutory MDS matrix.The matrix M 1 can easily be obtained by applying the parameters b 1 ¼ a 8 , b 2 ¼ a 9 and b 3 ¼ a 11 (and their inverses) to representative involutory MDS matrix RIM 1 .Example 11.Let F 2 4 be generated by the primitive element a which is a root of the primitive polynomial x 4 þ x þ 1 (0x13).Consider the 4 Â 4 involutory MDS matrix M 2 given in Pehlivanoğlu et al. (2018).
a a 2 a 5 a 7 a 10 a a 5 a 9 a 2 a 10 1 2 6 6 4 3 7 7 5 over F 2 4 =0x13.In fact, the involutory MDS matrix M 2 belongs to a class of which representative involutory MDS matrix is as follows: RIM 2 ¼ 1 a 6 a 2 a 3 a 9 a a 13 a 2 a 5 a 14 a a 6 a 6 a 5 a 9 1 The involutory MDS matrix M 2 can easily be obtained by applying the parameters b 1 ¼ a 9 , b 2 ¼ a 13 and b 3 ¼ a 12 (and their inverses) to representative involutory MDS matrix RIM 2 .

EXPERIMENTAL RESULTS
In this section, we generate all 4 Â 4 involutory MDS matrices over F 2 3 and F 2 4 by finding all representative involutory MDS matrices over these finite fields.In order to generate all representative involutory MDS matrices, we use the the matrix form R 1 given in the "Proposed Method" section.
From the experimental results given above, it is shown that (involutory and MDS) Hadamard matrix form, which is also representative (involutory and MDS) matrix form, can approximately generate 2.1% of all involutory MDS matrices over F 2 4 .It is likely that this percentage will reduce for involutory MDS matrices over larger finite fields.Since Lemma 6 and Theorem 8 given in "Proposed Method" section can be updated for 2 k Â 2 k involutory matrices, it is clear that (involutory and MDS) Hadamard matrix form can be used to generate a small subset of all 2 k Â 2 k involutory MDS matrices, where k .1, especially over larger finite fields.Remark 12.As stated above, there are 71,856 4 Â 4 representative involutory MDS matrices over F 2 4 , 1,512 of which can also be obtained by 4 Â 4 Hadamard matrix (by search).Then, 70,344 (¼ 71; 856 À 1; 512) 4 Â 4 representative involutory MDS matrices cannot be generated by 4 Â 4 Hadamard matrix form.That means, by using the method given in Sakallı et al. (2020), one can map these representative involutory MDS matrices in order to obtain directly isomorphic counterparts over F 2 8 .Note that there are 4 isomorphisms from the finite field F 2 4 (defined by any irreducible polynomial) to the finite field F 2 8 (defined by any irreducible polynomial).Hence, by using the parameters b 1 , b 2 , b 3 2 F 2 m À f0g and their inverses, one can directly generate 4,665,600,972,000 (= ð255Þ 3 Á 70; 344 Á 4 % 2 42:08 Þ 4 Â 4 involutory MDS matrices over F 2 8 (defined by any irreducible polynomial), which cannot also be generated by GHadamard matrix.
In Tables 1 and 2, we present the number of occurrences of 1s in all 4 Â 4 involutory MDS matrices over F 2 3 and F 2 4 , respectively.Remark 13.In Junod & Vaudenay (2005b), the maximum number of occurrences of 1s in 4 Â 4 MDS matrices is shown to be 9.In this article, we modify this result by showing that the maximum number of occurrences of 1s for 4 Â 4 involutory and MDS matrices is also 9 (see "Appendix" section).
In this article, we generate the lightest 4 Â 4 involutory MDS matrices over F 2 3 , F 2 4 and F 2 8 in terms of XOR count.To do so, we first generate involutory MDS matrices with low naive XOR counts.In this context, we consider all 4 Â 4 involutory MDS matrices over F 2 3 and consider 4 Â 4 involutory MDS matrices over F 2 4 with up to a threshold naive XOR count.For 4 Â 4 involutory MDS matrices over F 2 8 , we consider anti-diagonal symmetric matrices with up to a threshold naive XOR count and choose among them the ones whose results are up to a maximum of 110 XOR operations after optimizing with PAAR 1. Finally, 4 Â 4 involutory MDS matrices suitable for the given criteria are optimized with BDKCI algorithm.
In Example 14, we present a new 4 Â 4 involutory MDS matrix over F 2 3 generated by using the proposed method, which can be implemented by only 13 XOR operations with depth 6. Example 14.Let F 2 3 be generated by the primitive element a which is a root of the primitive polynomial x 3 þ x 2 þ 1 (0xd).Consider the 4 Â 4 involutory MDS matrix M 3 as given below: over F 2 3 =0xd.In fact, the involutory MDS matrix M 3 belongs to a class of which representative involutory MDS matrix is as follows: RIM 3 ¼ 1 a 6 a 5 a 3 a 6 a a a 4 a 5 a a 2 a 2 a 3 a 4 a 2 a 4 2 6 6 4 3 7 7 5 The involutory MDS matrix M 3 with d-XOR count 56 (¼ 20 þ 4 Á 3 Á 3) can easily be obtained by applying the parameters b 1 ¼ a, b 2 ¼ a 2 and b 3 ¼ a 4 (and their inverses) to representative involutory MDS matrix RIM 3 .After applying BDKCI heuristic to the matrix M 3 , we find the circuits for the matrix M 3 with 13 XORs and depth 5 using XOR3 and XOR4 gates (please see Table 5 for details).
In Example 15, we present a new 4 Â 4 involutory MDS matrix over F 2 4 by using the proposed method, which is better than the matrices given in the literature.It can be implemented by only 25 XOR operations with depth 5, whereas the previously known lightest one (Sarkar & Syed, 2016) requires 26 XOR operations with the same depth using XOR2 and XOR3 gates.Example 15.Let F 2 4 be generated by the primitive element a which is a root of the primitive polynomial x 4 þ x þ 1 (0x13).Consider 4 Â 4 involutory MDS matrix M 4 as given below: M 4 ¼ 1 a 4 a 14 a 1 a 14 a a 14 a 8 a a 14 a 4 a 2 a 8 1 1 2 6 6 4 3 7 7 5 over F 2 4 =0x13.In fact, the involutory MDS matrix M 4 belongs to a class of which representative involutory MDS matrix is as follows: a 9 a 7 a 4 a 14 1 a 9 a 13 a 2 a 14 1 a 11 a 13 a 4 1 2 6 6 4 3 7 7 5 The involutory MDS matrix M 4 with d-XOR count 79 ð¼ 31 þ 4 Â 3 Â 4Þ can easily be obtained by applying the parameters b 1 ¼ a 4 , b 2 ¼ a 5 and b 3 ¼ a 9 (and their inverses) to representative involutory MDS matrix RIM 4 .After applying BDKCI heuristic to the matrix M 4 , we find the circuits for the matrix M 4 with 25 XORs and depth 5 using XOR2 and XOR3 gates (please see Table 6 for details).
In Example 16, we present a new 4 Â 4 non-involutory MDS matrix over F 2 4 .We have found this matrix by searching through all anti-diagonal symmetric matrices with the same elements of the 4 Â 4 involutory MDS matrix presented in Example 15.This noninvolutory MDS matrix can be implemented by only 19 XOR operations with depth 3. Example 16.Let F 2 4 be generated by the primitive element a which is a root of the primitive polynomial x 4 þ x þ 1 (0x13).Consider the 4 Â 4 anti-diagonal symmetric non-involutory MDS matrix M 5 as given below:     As shown in Table 3, our proposed method leads to better results than the best-known XOR GCs and GEs given in the literature.In the case of involutory 4 Â 4 MDS matrices over F 2 4 with depth 5, the matrix M 4 presented in Example 15 offers 1 XOR operations improvement from the best-known result given in Sarkar & Syed (2016).On the other hand, the lightest 4 Â 4 involutory MDS matrix M 3 over F 2 3 with the lowest known XOR count 13 (the lowest XOR count known so far) with depth 5 is presented in Table 3.In non-involutory MDS matrices section in Table 3, we compare the matrix M 5 presented in Example 16 with the best known XOR count results in the literature.The matrix M 5 as depth 3 is the lightest XOR count result with 19 XOR operations offering 1 XOR operations improvement.
Operation No. Operation No.

Operation
No. Operation Example 17.Let F 2 8 be generated by the primitive element a which is a root of the primitive polynomial The involutory MDS matrix M 6 with d-XOR count 162 (¼ 66 þ 4 Á 3 Á 8) can easily be obtained by applying the parameters b 1 ¼ a 146 , b 2 ¼ a 127 and b 3 ¼ a 18 (and their inverses)  to representative involutory MDS matrix RIM 6 .After applying BDKCI heuristic to the matrix M 6 , we find the circuit for the matrix M 6 with 52 XOR operations and depth 5 using XOR2 and XOR3 gates.Moreover, the same matrix M 6 requires only 42 XOR operations with depth 4 using XOR2, XOR3 and XOR4 gates (please see Table 7 for details).
In Table 4, we compare the matrix M 6 with the circuit implementations of ciphers given in the literature.When comparing the matrix M 6 with the previous best-known implementations of linear layers of some block ciphers, the involutory MDS matrix M 6 improves the results in GCs and GEs with the minimum circuit depths.
Note that in Tables 3 and 4 we do not only consider the construction of circuits with minimum GCs (i.e., XOR gate counts) and GEs but also take into account the minimum circuit depth for low latency criterion.Since we stopped the algorithm after four hours of runtime for each matrix, lighter circuits can be found because the proposed new hybrid method allows us to generate all 4 Â 4 involutory MDS matrices over F 2 m .For further improvements, more runtime may be required.

CONCLUSION AND FUTURE WORKS
In this article, we proposed a new hybrid method to generate all 4 Â 4 involutory MDS matrices over F 2 m .The proposed method reduces the search space complexity at the level of ffiffiffi n p by searching and finding representative involutory MDS matrices, where n represents the number of all invertible 4 Â 4 matrices over F 2 m .In this respect, we were able to generate all 4 Â 4 involutory MDS matrices over F 2 3 and F 2 4 .For the finite field F 2 8 , the search space to generate all representative involutory MDS matrices is approximately 2 64 , which is still too high to be searched for.Nevertheless, one can easily generate 4 Â 4 involutory and MDS matrices over F 2 8 by focusing on representative involutory MDS matrices.In the future, if the number of parameters in R 1 presented in the "Proposed Method" section is reduced by eliminating some parameters and a new search form is obtained, then the search space for finding representative involutory MDS matrices for larger finite fields can be reduced.However, this new form is highly possible to be a more complex structure than R 1 .As a result, we believe that our proposed method is not only useful for generating all involutory MDS matrices over F 2 m , but it is also useful for other methods in the literature like the methods used for constructing lightweight involutory MDS matrices over the general linear groups GL (m,F 2 ).

DECLARATIONS
All authors have read and agreed to the submitted version of the manuscript.Funding is not applicable.The authors declare no conflict of interest.Data sharing is not applicable to this article as no data sets were generated or analyzed during the current study.

APPENDIX
Let F 2 4 be generated by the primitive element a which is a root of the primitive polynomial x 4 þ x þ 1 (0x13).Consider 4 Â 4 representative involutory MDS matrix as given below: The number of 1s in the matrix M 7 is 9, which is the maximum number of occurrences of 1s in 4 Â 4 MDS matrices.
Tuncay et al. (2023), PeerJ Comput.Sci., DOI 10.7717/peerj-cs.157714/25 Tuncay et al. (2023), PeerJ Comput.Sci., DOI 10.7717/peerj-cs.157721/25RIM 7 ¼ 1 a 4 a 11 a 13 a 11 a 12 a 11 a 11 a 3 a 8 a 12 a 4 a 5 a 3 a 11 =0x13.Then, one can generate 4 Â 4 involutory and MDS matrix M 7 by applying the parameters b 1 ¼ a 11 , b 2 ¼ a 4 and b 3 ¼ 1 (and their inverses) to RIM 7 as follows: Pehlivanoğlu et al. (2018) m is equal to 1, then H is involutory matrix, i.e., H 2 ¼ I and H ¼ H À1 , where I is identity matrix and H À1 is the inverse of H. On the other hand, GHadamard matrix form presented inPehlivanoğlu et al. (2018)satisfies the last property À Tuncay et al. (2023), PeerJ Comput.Sci., DOI 10.7717/peerj-cs.15776/25 H i;j ¼ a iÈj , where a i parameters are the entries of the first row of a k Â k Hadamard matrix H. H is a bi-symmetric matrix, namely, H ¼ H T and HJ ¼ JH where J is a k Â k exchange matrix ðJ i;kÀiþ1 ¼ 1 and other elements of J are 0), H 2 ¼ c 2 Â I, where c ¼ È kÀ1 i¼0 a i and I is the k Â k identity matrix.If XOR sum of the first row elements of a k Â k f0g.Note that the following involutory MDS Hadamard matrix form RIM 2Â2 : is also representative involutory MDS matrix form for 2 Â 2 matrices over F 2 m with the restriction r 11 2 F 2 m À f0; 1g.Hence, the representative involutory MDS matrix form RIM 2Â2 can be used to generate all involutory MDS matrices over F 2 m .
!be a 2 Â 2 matrix over F 2 m .If the matrix A is involutory and MDS, then all 2 Â 2 involutory MDS matrices can be generated by the following matrix form:

Table 1
The number of occurrences of 1s in all 4 Â 4 involutory MDS matrices over F 2 3.

Table 3
Comparison of 4 Â 4 involutory and non-involutory MDS matrices in view of XOR counts with different depths obtained by using BDKCI heuristic.

Table 4
Comparison of best known implementation cost of few 32 Â 32 matrices (linear layers of block ciphers) in ASIC libraries.
Bold values indicate the best results.