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

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 ${\mathbf{F}}_{2^3}$, ${\mathbf{F}}_{2^4}$ and ${\mathbf{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 ${\mathbf{F}}_{2^3}$, ${\mathbf{F}}_{2^4}$ and ${\mathbf{F}}_{2^8}$. Note that our construction technique can easily be used for designing new involutory MDS matrices over finite field ${\mathbf{F}}_{2^m}$ especially for $m\le 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\times 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, 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\times k$ MDS matrices (new MDS matrices from the implementation point of view) from any existing $k\times 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 ${\mathbf{F}}_{2^m}$. In this study, we present a new hybrid method to generate all $4\times 4$ involutory MDS matrices over ${\mathbf{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\times 4$ involutory MDS matrices over ${\mathbf{F}}_{2^3}$ and ${\mathbf{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\times 4$ involutory MDS matrices in terms of XOR counts and circuit depths by using BDKCI algorithm. Note that when obtaining the lightest $4\times 4$ non-involutory MDS matrix over ${\mathbf{F}}_{2^4}$, we take the benefit of the elements of the lightest $4\times 4$ involutory MDS matrix over ${\mathbf{F}}_{2^4}$ generated.

In this article, we propose a new hybrid method to generate all $4\times 4$ involutory MDS matrices over ${\mathbf{F}}_{2^m}$. The proposed hybrid method consists of two parts: the first part includes searching for all $4\times 4$ representative involutory MDS matrices and the second part includes generating all $4\times 4$ involutory MDS matrices using all $4\times 4$ representative involutory MDS matrices found by search in the first part of the proposed method. Hence, we present the number of all $4\times 4$ involutory MDS matrices over ${\mathbf{F}}_{2^3}$ and ${\mathbf{F}}_{2^4}$. In addition, we give the lightest known involutory/non-involutory MDS matrices over ${\mathbf{F}}_{2^3}$, ${\mathbf{F}}_{2^4}$ and ${\mathbf{F}}_{2^8}$ obtained by using the global optimization algorithm BDKCI, which is an improved version of BP algorithm.

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 ${\mathbf{F}}_{2^m}$ consisting of $2^m$ elements is defined by an irreducible polynomial $p(x)$ of degree $m$ over ${\mathbf{F}}_2$ and is denoted by ${\mathbf{F}}_2[x]/(p(x))$. The elements of the finite field ${\mathbf{F}}_{2^m}$ can be represented as polynomials over ${\mathbf{F}}_2[x]/p(x)$, i.e., ${\sum }_{i=0}^{m-1}{a}_i{\alpha }^i$, where $a_i$ $\in$ ${\mathbf{F}}_2$ and $\alpha$ is a root (and a primitive element) of ${\mathbf{F}}_{2^m}$. For simplicity, we denote ${\mathbf{F}}_{2^m}$ defined by irreducible polynomial $p(x)$ as ${\mathbf{F}}_{2^m}/p(x)$ (the finite field with $2^m$ elements) and use hexadecimal notation to represent the elements of ${\mathbf{F}}_{2^m}$ and the irreducible polynomial $p(x)$ used for defining ${\mathbf{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 ${\alpha }^3+{\alpha }^2+\alpha$ in ${\mathbf{F}}_{2^4}$. In the same manner, 0x13 used when denoting ${\mathbf{F}}_{2^4}/\mathsf{0}\mathsf{x}\mathsf{13}$ stands for the irreducible polynomial $p(x)=x^4+x+1$.

If an $[n,k,d]$ code C reaches the Singleton bound $d=n-k+1$, then C is called an MDS code (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\times k$ for MDS matrices), which are two important cryptographic criteria for linear transformations. The followings are the properties of an MDS matrix:

• Let $\mathit{M}$ be a $k\times k$ square matrix. $\mathit{M}$ is an MDS matrix, if and only if every square sub-matrix of $\mathit{M}$ is non-singular.

• If $\mathit{M}$ is a $k\times k$ MDS matrix, the transpose matrix of $\mathit{M}$ ( ${\mathit{M}}^T$) is also an MDS matrix.

• Let $c\in$ ${\mathbf{F}}_{2^m}$ be a non-zero constant and let $\mathit{M}$ be a $k\times k$ MDS matrix, then the multiplication of a row (or column) of $\mathit{M}$ by $c$ does not affect the MDS property of $\mathit{M}$.

If a square matrix $\mathit{A}$ is its own inverse (i.e., $\mathit{A}$ = ${\mathit{A}}^{-1}$), then the matrix $\mathit{A}$ is called an involutory matrix. Equivalently, the matrix $\mathit{A}$ is an involution if and only if ${\mathit{A}}^2$ = $\mathit{I}$, where $\mathit{I}$ is the identity matrix.

When denoting a $k\times k$ Hadamard matrix, we use the notation had( $a_0$, $a_1$,…, $a_{k-1}$), where $a_i$ $\in$ ${\mathbf{F}}_{2^m}$ for $0\le i\le k-1$. In this respect, a $4\times 4$ Hadamard matrix $\mathit{H}$ can be given as follows:

(2) $$\mathit{H}=had(a_0,a_1,a_2,a_3)=\left[\begin{array}{cccc}{a}_0& a_1& a_2& a_3\\ a_1& a_0& a_3& a_2\\ a_2& a_3& a_0& a_1\\ a_3& a_2& a_1& a_0\end{array}\right]$$

Some important properties of a $k\times k$ Hadamard matrix $\mathit{H}$ over ${\mathbf{F}}_{2^m}$ are as follows (Pehlivanoğlu et al., 2018):

• ${\mathit{H}}_{i,j}=a_{i\oplus j}$, where $a_i$ parameters are the entries of the first row of a $k\times k$ Hadamard matrix $\mathit{H}$.

• $\mathit{H}$ is a bi-symmetric matrix, namely, $\mathit{H}={\mathit{H}}^T$ and $\mathit{H}\mathit{J}=\mathit{J}\mathit{H}$ where $\mathit{J}$ is a $k\times k$ exchange matrix $({\mathit{J}}_{i,k-i+1}=1$ and other elements of $\mathit{J}$ are 0),

• ${\mathit{H}}^2=c^2\times \mathit{I}$, where $c={\oplus }_{i=0}^{k-1}{a}_i$ and $\mathit{I}$ is the $k\times k$ identity matrix.

If XOR sum of the first row elements of a $k\times k$ Hadamard matrix $\mathit{H}$ over ${\mathbf{F}}_{2^m}$ is equal to 1, then $\mathit{H}$ is involutory matrix, i.e., ${\mathit{H}}^2=\mathit{I}$ and $\mathit{H}={\mathit{H}}^{-1}$, where $\mathit{I}$ is identity matrix and ${\mathit{H}}^{-1}$ is the inverse of $\mathit{H}$. On the other hand, GHadamard matrix form presented in Pehlivanoğlu et al. (2018) satisfies the last property $(i.e.,{\mathit{H}}^2=c^2\times \mathit{I}$, where $c={\oplus }_{i=0}^{k-1}{a}_i)$ given above for a $k\times k$ Hadamard matrix and preserves involutory and MDS properties of a given $k\times 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\times k$ GHadamard matrix $\mathit{G}\mathit{H}$ is generated by using the combination of non-zero $k-1$ more $b_i$ parameters and their inverses with a $k\times k$ Hadamard matrix $\mathit{H}$ over ${\mathbf{F}}_{2^m}$. In this context, a $4\times 4$ GHadamard matrix $\mathit{G}\mathit{H}$ can be denoted as follows:

(3) $$\mathit{G}\mathit{H}=Ghad(a_0,a_1;b_1,a_2;b_2,a_3;b_3)=\left[\begin{array}{cccc}{a}_0& a_1{b}_1& a_2{b}_2& a_3{b}_3\\ a_1{b}_1^{-1}& a_0& a_3{b}_1^{-1}{b}_2& a_2{b}_1^{-1}{b}_3\\ a_2{b}_2^{-1}& a_3{b}_2^{-1}{b}_1& a_0& a_1{b}_2^{-1}{b}_3\\ a_3{b}_3^{-1}& a_2{b}_3^{-1}{b}_1& a_1{b}_3^{-1}{b}_2& a_0\end{array}\right]$$

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\leftarrow x_{a_i}\oplus 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\times mk$ binary matrix (transformed from $k\times k$ matrix over ${\mathbf{F}}_{2^m}$) minus $mk$, s-XOR count is defined as the minimum number of XOR operations needed to implement the $mk\times 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 $\{x_0,x_1,\dots ,x_n\}$ and $m$ output signals $\{y_0,y_1,\dots ,y_m\}$ (also called target signals). We also use BDKCI algorithm to build a globally optimized implementation of a $4\times 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\times 4$ involutory MDS matrices over ${\mathbf{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\times 4$ representative involutory MDS matrices by search, and then to obtain all $4\times 4$ involutory MDS matrices by applying three more non-zero 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 ${\mathbf{F}}_{2^m}$ as follows:

$$R_1=\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& r_{11}+r_{12}+r_{13}+1\\ r_{21}& r_{22}& r_{12}+r_{13}+r_{21}+r_{31}+r_{32}& r_{12}+r_{13}+r_{22}+r_{31}+r_{32}+1\\ r_{31}& r_{32}& r_{33}& r_{31}+r_{32}+r_{33}+1\\ r_{11}+r_{21}+r_{31}+1& r_{12}+r_{22}+r_{32}+1& r_{12}+r_{21}+r_{31}+r_{32}+r_{33}+1& r_{11}+r_{22}+r_{33}\end{array}\right].$$

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 ${\mathbf{F}}_{2^m}-\{0\}$. 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 ${\mathbf{F}}_{2^m}$ is approximately obtained as $(2^m-1{)}^8$ (by omitting non-invertible or singular matrices). For example, for finding all $4\times 4$ representative involutory MDS matrices over ${\mathbf{F}}_{2^4}$, the search space is approximately $(2^4-1{)}^8\approx 2^{31.25}$. Before proceeding with the details on the matrix form $R_1$, we introduce the matrix form $A_{2\times 2}$ used for generating all $2\times 2$ involutory MDS matrices over ${\mathbf{F}}_{2^m}$ in Remark 2 and Lemma 3 from which the main idea of the article comes.

Proof. Let $A=\left[\begin{array}{cc}{r}_{11}& r_{12}\\ r_{21}& r_{22}\end{array}\right]$ be a $2\times 2$ involutory matrix with the restriction $r_{11}\ne 0$. Let $c_{ij}$ denote the elements of $A^2$ for $i,j\in \{1,2\}$, i.e., $c_{ij}={\sum }_{k=1}^2{r}_{ik}{r}_{kj}$. Since $A^2=I$, where I is the $2\times 2$ identity matrix, we obtain the following equations (by considering if $i=j$ then $c_{ij}=1$ and if $i\ne j$ then $c_{ij}=0$):

(4) $$r_{11}^2+r_{12}{r}_{21}=1$$

(5) $$r_{11}{r}_{12}+r_{12}{r}_{22}=0$$

(6) $$r_{21}{r}_{11}+r_{22}{r}_{21}=0$$

(7) $$r_{21}{r}_{12}+r_{22}^2=1$$

After adding the Eqs. (4) and (7) given above, we obtain the equation $r_{11}^2=r_{22}^2$. The equation $r_{11}^2=r_{22}^2$ 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 ${\mathbf{F}}_{2^m}$. On the other hand, we have $r_{12}{r}_{21}=r_{11}^2+1=(r_{11}+1{)}^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\in {\mathbf{F}}_{2^m}-\{0\}$ and $r_{11}\in {\mathbf{F}}_{2^m}-\{0,1\}$ so that the matrix form can also satisfy the MDS property (i.e., the determinant of the matrix form $A_{2\times 2}$ is not equal to 0).

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\times 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\times 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\times 4$ matrices over ${\mathbf{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 $\sqrt{n}$, where $n$ represents the number of all invertible $4\times 4$ matrices over ${\mathbf{F}}_{2^m}$. In this context, we present the matrix form $R_1$ again below:

$$R_1=\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& r_{11}+r_{12}+r_{13}+1\\ r_{21}& r_{22}& r_{12}+r_{13}+r_{21}+r_{31}+r_{32}& r_{12}+r_{13}+r_{22}+r_{31}+r_{32}+1\\ r_{31}& r_{32}& r_{33}& r_{31}+r_{32}+r_{33}+1\\ r_{11}+r_{21}+r_{31}+1& r_{12}+r_{22}+r_{32}+1& r_{12}+r_{21}+r_{31}+r_{32}+r_{33}+1& r_{11}+r_{22}+r_{33}\end{array}\right].$$

Proof. Let $R=[r_{ij}]$ be a $4\times 4$ involutory and representative matrix such that all $r_{ij}\ne 0$ and let $c_{ij}={\sum }_{k=1}^4{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\times 4$ identity matrix. If $i=j$, then $c_{ij}=1$. Otherwise, $c_{ij}=0$. Hence, the following equations are satisfied:

(8) $${r_{11}}^2+r_{12}{r}_{21}+r_{13}{r}_{31}+r_{14}{r}_{41}=1$$

(9) $$r_{21}{r}_{12}+{r_{22}}^2+r_{23}{r}_{32}+r_{24}{r}_{42}=1$$

(10) $$r_{31}{r}_{13}+r_{32}{r}_{23}+{r_{33}}^2+r_{34}{r}_{43}=1$$

(11) $$r_{41}{r}_{14}+r_{42}{r}_{24}+r_{43}{r}_{34}+{r_{44}}^2=1$$

(12) $$r_{11}{r}_{12}+r_{12}{r}_{22}+r_{13}{r}_{32}+r_{14}{r}_{42}=0$$

(13) $$r_{11}{r}_{13}+r_{12}{r}_{23}+r_{13}{r}_{33}+r_{14}{r}_{43}=0$$

(14) $$r_{11}{r}_{14}+r_{12}{r}_{24}+r_{13}{r}_{34}+r_{14}{r}_{44}=0$$

(15) $$r_{21}{r}_{11}+r_{22}{r}_{21}+r_{23}{r}_{31}+r_{24}{r}_{41}=0$$

(16) $$r_{21}{r}_{13}+r_{22}{r}_{23}+r_{23}{r}_{33}+r_{24}{r}_{43}=0$$

(17) $$r_{21}{r}_{14}+r_{22}{r}_{24}+r_{23}{r}_{34}+r_{24}{r}_{44}=0$$

(18) $$r_{31}{r}_{11}+r_{32}{r}_{21}+r_{33}{r}_{31}+r_{34}{r}_{41}=0$$

(19) $$r_{31}{r}_{12}+r_{32}{r}_{22}+r_{33}{r}_{32}+r_{34}{r}_{42}=0$$

(20) $$r_{31}{r}_{14}+r_{32}{r}_{24}+r_{33}{r}_{34}+r_{34}{r}_{44}=0$$

(21) $$r_{41}{r}_{11}+r_{42}{r}_{21}+r_{43}{r}_{31}+r_{44}{r}_{41}=0$$

(22) $$r_{41}{r}_{12}+r_{42}{r}_{22}+r_{43}{r}_{32}+r_{44}{r}_{42}=0$$

(23) $$r_{41}{r}_{13}+r_{42}{r}_{23}+r_{43}{r}_{33}+r_{44}{r}_{43}=0$$

By adding the Eqs. (8)–(11) side by side, the following equation is obtained:

(24) $${r_{11}}^2+{r_{22}}^2+{r_{33}}^2+{r_{44}}^2=0$$

Because we study in the finite field ${\mathbf{F}}_{2^m}$ the Eq. (24) can be rewritten as follows:

(25) $$r_{11}+r_{22}+r_{33}+r_{44}=0$$

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:

(26) $$r_{11}+r_{12}+r_{13}+r_{14}=1$$

(27) $$r_{21}+r_{22}+r_{23}+r_{24}=1$$

(28) $$r_{31}+r_{32}+r_{33}+r_{34}=1$$

(29) $$r_{41}+r_{42}+r_{43}+r_{44}=1$$

By adding the Eqs. (12)–(14) side by side, we get the following equality:

(30) $$r_{11}\left(r_{12}+r_{13}+r_{14}\right)+r_{12}\left(r_{22}+r_{23}+r_{24}\right)+r_{13}\left(r_{32}+r_{33}+r_{34}\right)+r_{14}\left(r_{42}+r_{43}+r_{44}\right)=0$$

From the assumption Eqs. (26)–(28), it is clear that:

(31) $$r_{11}+1=r_{12}+r_{13}+r_{14}$$

(32) $$r_{21}+1=r_{22}+r_{23}+r_{24}$$

(33) $$r_{31}+1=r_{32}+r_{33}+r_{34}$$

(34) $$r_{41}+1=r_{42}+r_{43}+r_{44}$$

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:

(35) $$r_{11}(r_{11}+1)+r_{12}(r_{21}+1)+r_{13}(r_{31}+1)+r_{14}(r_{41}+1)=0$$

(36) $${r_{11}}^2+r_{11}+r_{12}{r}_{21}+r_{12}+r_{13}{r}_{31}+r_{13}+r_{14}{r}_{41}+r_{14}=0$$

We can rewrite the Eq. (36) as follows:

(37) $${r_{11}}^2+r_{12}{r}_{21}+r_{12}+r_{13}{r}_{31}+r_{14}{r}_{41}=r_{11}+r_{12}+r_{13}+r_{14}$$

Then, by using the assumption (26), we can obtain the Eq. (8) as follows:

(38) $${r_{11}}^2+r_{12}{r}_{21}+r_{12}+r_{13}{r}_{31}+r_{14}{r}_{41}=1$$

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:

(39) $$r_{11}+r_{21}+r_{31}+r_{41}=1$$

(40) $$r_{12}+r_{22}+r_{32}+r_{42}=1$$

(41) $$r_{13}+r_{23}+r_{33}+r_{43}=1$$

(42) $$r_{14}+r_{24}+r_{34}+r_{44}=1$$

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.

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\in \{1,2,3,4\}$:

(43) $$c_{11}={r_{11}}^2+r_{12}{r}_{21}+r_{13}{r}_{31}+r_{14}{r}_{41}$$

(44) $$c_{22}=r_{21}{r}_{12}+{r_{22}}^2+r_{23}{r}_{32}+r_{24}{r}_{42}$$

(45) $$c_{33}=r_{31}{r}_{13}+r_{32}{r}_{23}+{r_{33}}^2+r_{34}{r}_{43}$$

(46) $$c_{44}=r_{41}{r}_{14}+r_{42}{r}_{24}+r_{43}{r}_{34}+{r_{44}}^2$$

(47) $$c_{12}=b_1(r_{11}{r}_{12}+r_{12}{r}_{22}+r_{13}{r}_{32}+r_{14}{r}_{42})$$

(48) $$c_{13}=b_2(r_{11}{r}_{13}+r_{12}{r}_{23}+r_{13}{r}_{33}+r_{14}{r}_{43})$$

(49) $$c_{14}=b_3(r_{11}{r}_{14}+r_{12}{r}_{24}+r_{13}{r}_{34}+r_{14}{r}_{44})$$

(50) $$c_{21}=b_1^{-1}(r_{21}{r}_{11}+r_{22}{r}_{21}+r_{23}{r}_{31}+r_{24}{r}_{41})$$

(51) $$c_{23}=b_1^{-1}{b}_2(r_{21}{r}_{13}+r_{22}{r}_{23}+r_{23}{r}_{33}+r_{24}{r}_{43})$$

(52) $$c_{24}=b_1^{-1}{b}_3(r_{21}{r}_{14}+r_{22}{r}_{24}+r_{23}{r}_{34}+r_{24}{r}_{44})$$

(53) $$c_{31}=b_2^{-1}(r_{31}{r}_{11}+r_{32}{r}_{21}+r_{33}{r}_{31}+r_{34}{r}_{41})$$

(54) $$c_{32}=b_2^{-1}{b}_1(r_{31}{r}_{12}+r_{32}{r}_{22}+r_{33}{r}_{32}+r_{34}{r}_{42})$$

(55) $$c_{34}=b_2^{-1}{b}_3(r_{31}{r}_{14}+r_{32}{r}_{24}+r_{33}{r}_{34}+r_{34}{r}_{44})$$

(56) $$c_{41}=b_3^{-1}(r_{41}{r}_{11}+r_{42}{r}_{21}+r_{43}{r}_{31}+r_{44}{r}_{41})$$

(57) $$c_{42}=b_3^{-1}{b}_1(r_{41}{r}_{12}+r_{42}{r}_{22}+r_{43}{r}_{32}+r_{44}{r}_{42})$$

(58) $$c_{43}=b_3^{-1}{b}_2(r_{41}{r}_{13}+r_{42}{r}_{23}+r_{43}{r}_{33}+r_{44}{r}_{43})$$

The expressions from (43) to (46) are equal to 1 and the expressions from (47) to (58) are equal to 0 because the expressions from (43) to (46) are, respectively, the same with the left side of the equations from (8) to (11) given in Lemma 2, the expressions from (47) to (58) contain the results of multiplying 0 by non-zero finite field element(s) because the same expressions from (12) to (23) in Lemma 6, which are equal to 0, appear respectively in the expressions from (47) to (58). Hence, the involutory property of the matrix RIM is preserved. On the other hand, $b_i$ parameters in the matrix form A can be obtained by multiplying the constant elements $c_1$, $c_2$, $c_3$ and $c_4$ $\in$ ${\mathbf{F}}_{2^m}-\{0\}$ with the first column, the second column, the third column and the fourth column, respectively, and by multiplying the inverses of these elements $c_1^{-1}$, $c_2^{-1}$, $c_3^{-1}$ and $c_4^{-1}$ with the first row, the second row, the third row and the fourth row respectively. As given in “Preliminaries” section, the MDS property is invariant under the multiplication of a row/column with a non-zero constant. Hence, the MDS property of the matrix RIM is also preserved.

The proposed method can be divided into two stages as follows:

• By Lemma 6, representative involutory MDS matrices are found by searching through all possible candidates using the matrix form $R_1$, which constitutes the search side of the proposed method. This allows us to use 8 parameters ( $r_{11}$, $r_{12}$, $r_{13}$, $r_{21}$, $r_{22}$, $r_{31}$, $r_{32}$ and $r_{33}$) in the search phase of the matrix $R_1$. Hence, the search space for finding representative involutory MDS matrices over ${\mathbf{F}}_{2^m}$ is $(2^m-1{)}^8$ by excluding the finite field element 0 in all parameters (the finite field element 0 is ignored in the search space calculation because all entries of an MDS matrix should be non-zero). However, when searching for finding all $4\times 4$ involutory MDS matrices in a conventional manner, 16 parameters are needed to define a $4\times 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 $\sqrt{n}$, where $n$ represents the number of all invertible $4\times 4$ matrices over ${\mathbf{F}}_{2^m}$.

• By Theorem 8, all $4\times 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\times 4$ involutory MDS matrices over ${\mathbf{F}}_{2^4}$ by using the proposed method, which is also given in the literature.