Algebraic Attacks against Grendel: An Arithmetization-Oriented Primitive with the Legendre Symbol

Abstract

The rise of modern cryptographic protocols such as Zero-Knowledge proofs and secure Multi-party Computation has led to an increased demand for a new class of symmetric primitives. Unlike traditional platforms such as servers, microcontrollers, and desktop computers, these primitives are designed to be implemented in arithmetical circuits. In terms of security evaluation, arithmetization-oriented primitives are more complex compared to traditional symmetric cryptographic primitives. The arithmetization-oriented permutation Grendel employs the Legendre Symbol to increase the growth of algebraic degrees in its nonlinear layer. To analyze the security of Grendel thoroughly, it is crucial to investigate its resilience against algebraic attacks. This paper presents a preimage attack on the sponge hash function instantiated with the complete rounds of the Grendel permutation, employing algebraic methods. A technique is introduced that enables the elimination of two complete rounds of substitution permutation networks (SPN) in the sponge hash function without significant additional cost. This method can be combined with univariate root-finding techniques and Gröbner basis attacks to break the number of rounds claimed by the designers. By employing this strategy, our attack achieves a gain of two additional rounds compared to the previous state-of-the-art attack. With no compromise to its security margin, this approach deepens our understanding of the design and analysis of such cryptographic primitives.

1. Introduction

Arithmetization-oriented primitives have recently been widely employed in advanced cryptographic protocols, including Fully Homomorphic Encryption (FHE) protocols, Multiparty Computation (MPC) protocols, and Zero-Knowledge (ZK) proofs. These advanced cryptographic protocols employ arithmetic to convert normal calculations into a sequence of finite field operations, such as addition and multiplication over a large finite field ${\mathbb{F}}_p$, where p is a big prime integer greater than or equal to $2^{63}$. To characterize these finite field operations, arithmetization-oriented primitives are created. The design criterion for arithmetical primitives is to lessen the complexity of multiplication in cryptographic algorithms, as the primary resource consumption in advanced cryptographic protocols comes from the multiplication operation. Using the low-degree round function is an easy route to accomplishing this objective.

Various arithmetization-oriented primitives have been developed, including MiMC [1], GMiMC [2], HadesMiMC/Poseidon [3,4], Masta [5], Pasta [6], Ciminion [7], Chaghri [8], and Neptune [9]. These primitives directly use a low-degree round function as power maps $x\mapsto x^d$. More complex ones such as Rescue use the low-degree power map $x\mapsto x^3$ and its inverse $x\mapsto x^{1/3}$ as round functions. A new arithmetization-oriented primitive, Grendel [10], is designed for zero-knowledge proof systems. Grendel uses the Legendre symbol to enhance the round functions in combination with the SHARK-like construction. Here, ${\chi }_p(·):{\mathbb{F}}_p\mapsto \{-1,0,1\}$ is defined as ${\chi }_p\left(x\right):=x^{\frac{p-1}{2}}(modp)$ for the Legendre symbol. The application of the Legendre symbol in cryptography dates back to 1988. Ivan Damgård [11] proposed a new problem of predicting consecutive Legendre (Jacobi) symbols modulo a prime. This problem can be utilized to construct a cryptographically strong pseudorandom bit generator. This concept is closely tied to the distribution of quadratic residues and nonresidues modulo a prime number [12]. In 1997, Mauduit and Sárközy [13] introduced a range of metrics for quantifying the pseudo-randomness of binary sequences. Pseudo-random number generators and pseudo-random bit sequence generators, also known as pseudo-random sequence generators, have widespread applications in numerous scientific, technological, and industrial fields. They are utilized for process modeling, industrial problem solving, and cybersecurity purposes, serving as essential tools in these domains [14]. Tóth [15] and Gyarmati et al. [16] introduced new measures of pseudorandomness (avalanche effect and cross-correlation), and have asserted that those values in the Legendre symbol sequence (known as the Legendre symbol PRF, $x\mapsto {\chi }_p^k\left(x\right):={\chi }_p(x+k)$, where k is the private key) are high. In [17], Khovratovich developed a birthday-bound attack for the security analysis of the Legendre symbol PRF. This attack was subsequently enhanced by Beullens et al. [18] and Kaluderovic et al. [19]. According to more recent research by Seres et al., key-recovery attacks against the Legendre symbol PRF may be converted into the solution of a certain set of multivariate quadratic equation systems over a prime field [20].

In symmetric cryptographic schemes, the incorporation of the Legendre symbol into a round function requires the resulting nonlinear layer to be invertible. In [21], the authors proposed the construction of an invertible function using the Legendre symbol as follows: $x\mapsto x·({\chi }_p\left(x\right)+\alpha )$. The resulting function is invertible when ${\chi }_p({\alpha }^2-1)=1$. By combining the Legendre symbol with the power map, the map $x\mapsto x^d·{\chi }_p\left(x\right)$ is obtained, which is invertible when $gcd(d+(p-1)/2,p-1)=1$. In [22], Grassi et al. conducted a further analysis on the generalization of $x\mapsto x·({\chi }_p\left(x\right)+\alpha )$ to $x\mapsto x^d·({\chi }_p\left(x\right)+\alpha )$, building upon the foundations of [10,21]. They proposed new invertible functions that combine the Legendre symbol, and analyzed their statistical and algebraic properties.

When operating on large finite fields, arithmetization-oriented ciphers are less susceptible to statistical attacks such as differential [23] and linear [24] attacks. However, they are more vulnerable to algebraic attacks. For example, the cipher Jarvis [25] was found to be vulnerable to Gröbner basis attacks. The high-order differential attack is an effective method, as demonstrated in [26] for the high-order differential attack on GMiMC and in [27], where Eichlseder et al. first applied the high-order differential attack on MiMC. Subsequently, Bouvier et al. [28] and Cui et al. [29] analyzed the upper bounds on the algebraic degrees of MiMC, reevaluating its security margin against high-order differentials using different approaches. In [30], Liu et al. proposed an innovative technique called the coefficient grouping technique, which reduces the evaluation of algebraic degrees to a well-structured optimization problem. They applied this technique to launch a high-order differential attack on Chaghri, a fully homomorphic encryption scheme. Exploring the application of algebraic methods further to analyze arithmetization-oriented ciphers remains an interesting avenue for future investigation.

• Related Works. In [21], the authors proposed the construction of an invertible function using the Legendre symbol, although it was not utilized for cryptographic design. Recently, Grassi et al. [31] re-proposed the use of Legendre symbols for secure multi-party computation (MPC) applications. In the original security analysis of Grendel proposed by the designers in [10], the utilization of S-boxes based on the Legendre symbol was highlighted as a notable advantage. This choice allows a higher algebraic degree to be achieved within a relatively small number of rounds, providing resilience against high-order differential [32,33] and interpolation attacks [34]. Consequently, the focus of the analysis shifted towards the Gröbner basis attack, which presents two distinct approaches for constructing equation systems:
  • The first approach involves the attacker guessing all the Legendre symbols used in the scheme. Subsequently, they can solve the resulting system of equations and verify the correctness of the guessed symbols based on the obtained solution. The complexity of this attack increases by approximately a factor of 2 for each correctly guessed symbol, considering a probability of accurately guessing of around $1/2$.
  • On the other hand, the second approach avoids guessing the Legendre symbols and instead relies on the introduction of auxiliary variables to facilitate the establishment of the equation system. For more detailed information, please refer to [10].

Additionally, after guessing all Legendre symbols, the S-boxes in Grendel exhibit a low degree. As a result, it is not necessary to introduce intermediate variables in each round to mitigate the degree of growth. Instead, the attacker can directly solve a higher-degree system of equations. This alternative approach has already been used to attack the full hash function of Grendel [22]. In Section 3.1, we provide a detailed description of this attack. Despite proposing a full-round preimage attack on the hash function of Grendel, Ref. [22] only employed basic algebraic attacks. They treated the hash preimage as an unknown variable x, constructed a single-variable equation, and then solved it.

• Our Contribution. In this paper, we further analyze the hash function of Grendel on the basis of [10,22].
  • We introduce the Constrained Input/Constrained Output (CICO) problem [4] and exploit its solution to obtain preimages of the hash function of Grendel. In this way, we extend the previously proposed technique in [35] and improve the preimage attack by bypassing two additional rounds of the SPN structure. By introducing the CICO problem, our attack is capable of attacking two additional rounds compared to the attack presented in [22], as shown in

    Table 1. For a security level of $s=128$ and a modulus $p=2^{256}$, the number of rounds that can be attacked using the univariate root-finding method in different instances of the hash function Grendel.

    Instance $(\mathit{d},\mathit{n})$	Attacked Rounds in [10]	Attacked Rounds in [22]	Our Result
    (2, 3)	28	25	27
    (2, 4)	21	20	22
    (2, 8)	11	12	14
    (2, 12)	7	8	10
    (3, 3)	22	22	24
    (3, 4)	16	18	20
    (3, 8)	8	11	13
    (3, 12)	6	8	10
    (5, 3)	16	19	21
    (5, 4)	12	16	18
    (5, 8)	6	10	12
    (5, 12)	4	7	9

    .
  • Additionally, we employ the CICO problem to formulate a system of multivariate equations for the hash function of Grendel. Through an analysis of intermediate variable introduction and the core intricacies of Gröbner basis attacks, we enhance the understanding of constructing equation systems and executing Gröbner basis attacks.

• Organization. Here, we provide a brief overview of the organization of this paper. We present preliminaries in the following section. Section 3 covers the Algebraic Cryptanalysis of the Grendel hash function. In Section 3.3, we discuss our preimage attack on the Grendel hash function, utilizing a construction of a univariate equation system. Section 3.4 delves into the investigation of the Grendel hash function’s security through Gröbner basis attacks. Finally, we conclude the paper in Section 4.

2. Preliminaries

2.1. Notation

In the following, let p be a prime number and let ${\mathbb{F}}_p$ be a finite field with p elements. Let ${\mathbb{F}}_p^n$ denote a vector space with n elements and with each element in ${\mathbb{F}}_p$. The notation ${\mathbf{0}}^u$ represents a vector of length u in ${\mathbb{F}}_p^u$, where all components are zero. Considering $\mathit{X}\in {\mathbb{F}}_p^n$, we denote by $X_i$ its i-th component for each $i\in \{0,1,\dots ,n-1\}$, i.e., $\mathit{X}=(X_0,X_1,\dots ,X_{n-1})$.

Let ${\mathbb{F}}_p^n$ be a vector space with standard basis $\{e_0,e_1,\dots ,e_{n-1}\}$. A vector subspace $V_u$ of ${\mathbb{F}}_p^n$ can be represented as the span of a subset of a standard basis $\{e_0,e_1,\dots e_{u-1}\}$, where $0<u<n$. For ease of reference, please review the abbreviated

Table 2. List of abbreviations.

Abbreviation	Explanation
${\mathbb{F}}_p$	Finite field with p elements
${\mathbb{F}}_p^n$	Vector space with n elements over ${\mathbb{F}}_p$
${\mathbf{0}}^u$	Zero vector in ${\mathbb{F}}_p^u$ of length u
$\mathit{X}=(X_0,X_1,\dots ,X_{n-1})$	Vector in ${\mathbb{F}}_p^n$
$X_i$	i-th component of vector $\mathit{X}$
$V_u$	Vector subspace of ${\mathbb{F}}_p^n$ spanned by $\{e_0,e_1,\dots ,e_{u-1}\}$

.

Definition 1

(The Legendre Symbol). The Legendre symbol ${\chi }_p(·)$ is a function ${\chi }_p:{\mathbb{F}}_p\mapsto \{-1,0,1\}$, defined as

$$\begin{array}{c}\hfill {\chi }_p\left(x\right)=\left\{\begin{array}{cc}1\hfill & ifxisanonzeroquadraticresiduemodulop,\hfill \\ -1\hfill & ifxisaquadraticnon\text{-}residuemodulop,\hfill \\ 0\hfill & x=0.\hfill \end{array}\right.\end{array}$$

Proposition 1

([36]). For two prime integers $p,q\ge 3$, the Legendre symbol has the following properties:

• If $x\equiv y(modp)$, then ${\chi }_p\left(x\right)={\chi }_p\left(y\right)$.
• ${\chi }_p(x·y)={\chi }_p\left(x\right)·{\chi }_p\left(y\right)$.
• ${\chi }_p\left(q\right)·{\chi }_q\left(p\right)={(-1)}^{\frac{p-1}{2}·\frac{q-1}{2}}$.

2.2. CICO Problem

In the cryptanalysis of traditional symmetric schemes, the goal is to recover the key (or some subkeys) with complexity lower than $2^k$. However, the security of arithmetization-oriented hash functions such as the Grendel hash function relies on the computational infeasibility of solving the CICO problem.

Definition 2

(CICO Problem [4]). Let $F:{\mathbb{F}}_p^n\to {\mathbb{F}}_p^n$ be a permutation, and let $0<u<n$ be an integer. For given $(a_0,\dots a_{n-u-1}),(b_0,\dots ,b_{n-u-1})\in {\mathbb{F}}_p^{n-u}$, the CICO problem aims to find $(X_0,\dots ,X_{u-1}),(Y_0,\dots ,Y_{u-1})\in {\mathbb{F}}_p^u$ such that

$$\begin{array}{c}\hfill F(X_0,\dots ,X_{u-1},a_0,\dots ,a_{n-u-1})=(Y_0,\dots ,Y_{u-1},b_0,\dots ,b_{n-u-1}).\end{array}$$

A simpler version of the CICO problem can be defined as follows: when $n=2$ and $u=1$, the goal is to find $(X,Y)\in {\mathbb{F}}_p^2$ such that $F(X,0)=(Y,0)$. It can be observed that both the input and output of the permutation belong to the same vector subspace $V_u$. The CICO problem is highly relevant to the security of hash functions. Therefore, if the adversary has the ability to solve the problem with a complexity of less than $p^{n-u}$ permutation calls, it is possible to find a preimage or collision of the hash function under the sponge structure. The CICO problem can usually be modelled as a system of equations and solved algebraically.

2.3. Solve the Systems of Algebraic Equations

Our attack is based on modelling cryptographic primitives as a system of polynomial equations. In this section, we present the methods and complexities of solving some univariate and multivariate equations, which are then used to attack the hash function Grendel.

We assume that the cryptographic primitive is represented as a well-defined system, specifically, a system of m polynomial equations consisting of n variables $\mathit{X}=(X_0,\dots ,X_{n-1})\in {\mathbb{F}}_p^n$,

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{F}_0(X_0,\dots X_{n-1})=0\hfill \\ F_1(X_0,\dots X_{n-1})=0\hfill \\ ⋮\hfill \\ F_{m-1}(X_0,\dots X_{n-1})=0.\hfill \end{array}\right.\end{array}$$

Then, our purpose is to find the ordinary solution of the equation in the hope of obtaining the round key of the encryption schemes or the preimage of the hash function.

2.3.1. Solve a System of Univariate Equations

A univariate equation has only one variable and an equation of $F\left(x\right)=0$. Solving the given system is equivalent to finding the roots of the univariate polynomial $F\in {\mathbb{F}}_p\left[x\right]$ with degree $\mathcal{D}$ of F. Because all operations are performed on the finite field ${\mathbb{F}}_p$, the computational complexity is measured in terms of field operations.

• Compute $G=x^p-x(modF)$.
  The computation of $x^p(modF)$ requires $\mathcal{O}\left(\mathcal{D}·log\left(p\right)·log\left(\mathcal{D}\right)·\log \left(log\left(\mathcal{D}\right)\right)\right)$ field operations with a double-and-add algorithm.
• Compute $H=gcd(F,G)$.
  H has the same roots as F in ${\mathbb{F}}_p$, as $H=gcd(F,x^p-x)$; however, its degree is likely much lower. This step [37] requires $\mathcal{O}\left(\mathcal{D}·log{\left(\mathcal{D}\right)}^2\right)$ field operations.
• Factor H.
  In general, the polynomial H has only a few roots in ${\mathbb{F}}_p$. Thus, this step is negligible in complexity.

• This root-finding approach using GCD computations is provided in [37], and the final complexity of the algorithm is estimated by

  $$\begin{array}{c}\hfill \mathcal{O}\left(\mathtt{M}\left(\mathcal{D}\right)log\left(\mathcal{D}\right)log(\mathcal{D}·p)\right),\end{array}$$

  (1)

  where $\mathtt{M}\left(\mathcal{D}\right):=63.43·\mathcal{D}log\left(\mathcal{D}\right)log(log(\mathcal{D}\left)\right)+\mathcal{O}(\mathcal{D}log(\mathcal{D}\left)\right)$ is the complexity of multiplying two polynomials with a degree of at most $\mathcal{D}$ over ${\mathbb{F}}_p$.

2.3.2. Solve a System of Multivariate Equations

The Gröbner basis attack is a method of recovering a secret from a system of polynomial equations. The first step is to convert the primitive into a system of multivariate equations. Then, a Gröbner basis is computed for the ideal generated by these polynomials. Finally, the Gröbner basis is utilized to compute the target variables in the given system. This attack method involves the following three phases:

• When launching a Gröbner basis attack, the first step is to construct a set of polynomial equations describing the primitive. After that, a Gröbner basis for the ideal generated by these equations is computed, usually concerning the degrevlex ordering for better efficiency. The algorithm used for the computation of the Gröbner basis could be Buchberger’s algorithm [38], F4 [39], or F5 [40].
• After computing the Gröbner basis for the given system of polynomial equations, the next step is to perform a change of term order to facilitate the computation of the elimination ideals and the elimination of variables. This is typically achieved by going from the degrevlex term order to the lex one using an algorithm such as FGLM [41]. It is worth noting that in many applications, including those in cryptography, the systems of algebraic equations results in zero-dimensional ideals, meaning that they have only finite solutions.
• The final step of a Gröbner basis attack is to solve the univariate equation for the last variable using a polynomial factoring algorithm. This allows the specific value of the last variable to be obtained; this can then be substituted into the remaining equations to obtain the full solution of the system. This step can use the algorithm mentioned above to find the univariate equation system. When the polynomial has been factored, its roots can be easily found, and correspond to the possible values of the last variable. By substituting each root into the remaining equations, it is possible to obtain all possible solutions to the system of equations.

• Cost of Gröbner Basis Computation. For a system of m polynomial equations and n variables, we have

  $$F_0(X_0,\dots ,X_{n-1})=F_1(X_0,...,X_{n-1})=\dots =F_{m-1}(X_0,...,X_{n-1})=0,$$

  where $F_i\in {\mathbb{F}}_p[X_0,...,X_{n-1}]$, $0\le i\le m$. The complexity of computing a Gröbner basis in degrevlex term order [42] is

  $$\begin{array}{c}\hfill \mathcal{O}\left({\left(\genfrac{}{}{0pt}{}{n+D_{reg}}{D_{reg}}\right)}^{\omega }\right).\end{array}$$

  (2)
  In [43], another bound for the complexity of computing the Gröbner basis was provided:

  $$\begin{array}{c}\hfill \mathcal{O}\left(n{D}_{reg}·{\left(\genfrac{}{}{0pt}{}{n+D_{reg}-1}{D_{reg}}\right)}^{\omega }\right),\end{array}$$

  (3)

  where $2\le \omega <2.3727$ is the linear algebra constant representing the complexity of matrix multiplication and $D_{reg}$ is the degree of regularity. By further comparing these two complexities and computing their ratio, we can observe that

  $$\begin{array}{c}\hfill \frac{{\left(\genfrac{}{}{0pt}{}{n+D_{reg}}{D_{reg}}\right)}^{\omega }}{n{D}_{reg}·{\left(\genfrac{}{}{0pt}{}{n+D_{reg}-1}{D_{reg}}\right)}^{\omega }}=\frac{(n+D_{reg}^{\omega })}{n^{\omega +1}·D_{reg}}.\end{array}$$
  When n is small and $D_{reg}$ is large, the complexity calculation of Formula (3) provides a tighter bound. However, the authors of [22] found that when n is small, the complexity of computing the Gröbner basis is asymptotically smaller than the complexity of the FGLM algorithm. Therefore, for small values of n the complexity of the Gröbner basis attack depends on the complexity of the FGLM algorithm. On the other hand, Formula (2) becomes more restrictive when n has a comparatively larger value.
• Cost of FGLM algorithm. When employing the FGLM [41] algorithm, the complexity of converting the degrevlex order to the lex order is

  $$\begin{array}{c}\hfill \mathcal{O}(n·{\mathcal{D}}_{\mathcal{I}}^3),\end{array}$$

  where n represents the number of variables and ${\mathcal{D}}_{\mathcal{I}}$ denotes the degree of the zero-dimensional ideal.

Assuming that the polynomial system is a regular system, we consider n polynomials with the same degree $d_i=\delta ,i\in [1,n]$ and n variables. In this case, $D_{reg}$ can be estimated as $1+{\sum }_{i=1}^n{d}_i-1$. If the polynomial system is not regular, then its $D_{reg}$ is less than $1+{\sum }_{i=1}^n{d}_i-1$, a bound known as Macaulay’s bound.

2.4. Description of the Grendel Hash Function

The hash function Grendel is composed of the Grendel permutation combined with a sponge structure. Let $p\ge 3$ be a prime number and let $n\ge 2$ be an integer. The Grendel permutation $\mathcal{P}:{\mathbb{F}}_p^n\to {\mathbb{F}}_p^n$ is obtained by iteratively applying the Grendel round function $\mathcal{F}:{\mathbb{F}}_p^n\to {\mathbb{F}}_p^n$ for R rounds. The state size is n. Each round employs a distinct round constant. Each round function consists of three parts: a nonlinear layer, a linear layer, and adding round constants, respectively, denoted as $\mathcal{NL}$, $\mathcal{L}$, and $\mathcal{AC}$.

• The Nonlinear Layer: let $\mathit{X}=(X_0,\dots ,X_{n-1})\in {\mathbb{F}}_p^n$; then, $\mathcal{NL}:{\mathbb{F}}_p^n\to {\mathbb{F}}_p^n$ consists of independent n identical S-boxes $\mathcal{NL}\left(X\right)=(S\left(X_0\right),S\left(X_1\right),\dots ,S\left(X_{n-1}\right))$, where

  $$\begin{array}{c}\hfill S\left(X\right)=X^d·{\chi }_p\left(X\right),\end{array}$$

  with ${\chi }_p$ being the Legendre symbol. Here, $d\ge 2$ is an integer that satisfies $gcd(\frac{2d+p-1}{2},p-1)=1$.
• The Linear Layer: $\mathcal{L}:{\mathbb{F}}_p^n\to {\mathbb{F}}_p^n$ is an $n\times n$ MDS matrix $\mathcal{M}\in {\mathbb{F}}_p^{n\times n}$.
• The Adding Round Constants Step: this involves the utilization of round constants $c_j^i\in {\mathbb{F}}_p$, where $0\le i\le R-1$ and $0\le j\le n-1$.

• The Grendel round function $\mathcal{F}$ consists of three parts, and can be described as $\mathcal{F}(·)=\mathcal{AC}\circ \mathcal{L}\circ \mathcal{NL}(·)$. The Grendel permutation $\mathcal{P}$ is iterated over R rounds by $\mathcal{F}$, which can be expressed as $\mathcal{P}(·)=\underset{R}{\underbrace{\mathcal{F}\circ \dots \circ \mathcal{F}(·)}}$. The pseudocode describing the Grendel permutation is shown in Algorithm 1.

Algorithm 1: The Grendel Permutation $\mathcal{P}$

Input: $\mathit{X}=(X_0,X_1,\dots ,X_{n-1})\in {\mathbb{F}}_p^n$;    Output: $\mathit{Y}=(Y_0,Y_1,\dots ,Y_{n-1})\in {\mathbb{F}}_p^n$.   1: for $r=0$ to $R-1$ do   2:       for $i=0$ to $n-1$ do   3:             $X_i\leftarrow X^d·{\chi }_p\left(X_i\right)$;   4:       end for   5:       $\mathit{X}\leftarrow \mathcal{M}·\mathit{X}$;   6:       for $i=0$ to $n-1$ do   7:             $X_i\leftarrow X_i+c_i^r$;   8:       end for   9: end for 10: $\mathit{Y}\leftarrow \mathit{X}$; 11: return $\mathit{Y}$;

The sponge construction (Figure 1) [44,45] is a cryptographic framework that utilizes an internal cryptographic permutation or function. It provides versatility in achieving different objectives, including encryption, authentication, and hashing. The construction is based on the concept of a sponge, which consists of an internal permutation that operates on a fixed-size state. By appropriately configuring the sponge, it can be adapted for various cryptographic applications, providing security and flexibility. In this paper, we make slight modifications to the original approach in order to operate on elements of ${\mathbb{F}}_p$ instead of ${\mathbb{F}}_2$. Both the input and the output may be of arbitrary size. The state size is $n=r+c$, where r denotes the rate and c denotes capacity. To process a message m which consists of elements from the field ${\mathbb{F}}_p$, we utilize the following operations:

• Padding: if the length of the message is already a multiple of r, no padding is necessary; however, if the length is not a multiple of r, we first append $1\in {\mathbb{F}}_p$ to the message, then pad the message with 0 until its length becomes a multiple of r.
• Absorption: the message is divided into blocks of size r. Each block is added to the first r blocks of the state using the addition operation. Afterwards, the entire state is processed by applying the permutation function $\mathcal{P}$. Repeat the above operation until all the messages are absorbed.
• Squeezing: in each iteration of the squeezing phase, a block of length r is squeezed out, then the permutation function $\mathcal{P}$ is applied to the entire state and the squeezed block is extracted. This process is repeated until the squeezing phase is completed.

Security. According to the proof presented in [45], when the inner permutation bears resemblance to a random permutation, the sponge construction is indistinguishable from a random oracle up to approximately $p^{c/2}$ queries. Equivalently, in order to provide s bits of security, we need $p^{c/2}\ge 2^s$ and $p^r\ge 2^s$, i.e., $c\ge \lceil 2s·{log}_p\left(2\right)\rceil$. For such a hash function $\mathcal{H}:{\mathbb{F}}_p^{\ast }\to {\mathbb{F}}_p^{\infty }$, it is hard to find

• collision resistance: $x,x^{\prime }\ne x$ such that $\mathcal{H}\left(x\right)=\mathcal{H}\left(x^{\prime }\right)$
• preimage resistance: x, given y such that $\mathcal{H}\left(x\right)=y$
• second-preimage resistance: $x^{\prime }$, given $x\ne x^{\prime }$ such that $\mathcal{H}\left(x^{\prime }\right)=\mathcal{H}\left(x\right)$.

We assume an output of at least $\lceil 2s/{log}_2\left(p\right)\rceil$ elements to prevent birthday bound attacks. Furthermore, we require $c\ge \lceil 2s/{log}_2\left(p\right)\rceil$ for an s-bit security level.

3. Algebraic Cryptanalysis of Grendel Hash Function

In this section, we simply review the preimage attack proposed by [22] on a sponge hash function instantiated with the Grendel permutation presented in Section 3.1. We then introduce the CICO problem and provide a further analysis of the security of the Grendel hash function.

3.1. Preimage Attack on Hash Function Grendel in [22]

Let s denote the security level and let p represent the prime that defines the field. The authors of [22] focused on the case where $p\ge 2^s$, allowing for $r\ge 1$. Here, r defines the rate of the sponge hash function. In this scenario, the hash function can output a single element from ${\mathbb{F}}_p$, which aligns with common practice, as p is typically chosen to be large.

For a hash digest $h\in {\mathbb{F}}_p$, the objective is to find a preimage. For cases where $r\ge 2$, the authors of [22] started by fixing $r-1$ input elements. In contrast to the analysis in [10], they avoided introducing intermediate variables. Instead, they [22] employed polynomials of degree $d^R$ with fixed Legendre symbols, where R denotes the number of attacked rounds. Essentially, the attack on R-round construction involves the following steps:

• Iterating over all possible sets of Legendre symbols. The probability of a Legendre symbol being $\pm 1$ is approximately $\frac{1}{2}$, while the probability of it being 0 is $\frac{1}{p}$. Consequently, the probability that l Legendre symbols are different from zero can be calculated as ${(1-\frac{1}{p})}^l$. For a large number of rounds, if p is approximately $2^{32}$, this probability exceeds $99.99\%$. In their attack, $l=nR-(n-1)=n(R-1)+1$. In the first round, it is possible to compute $n-1$ Legendre symbols deterministically because there is no linear layer before the initial application of the S-boxes.
• Solving the resulting univariate equation to identify a preimage. The authors focused on the case in which the number of hash output elements is 1. By fixing all Legendre symbols, there is only a single unknown (the input variable) and a single equation of degree at most $d^R$ in the end. The equation system consists of only one univariate equation, and can be solved by applying a root-finding algorithm to this equation.
• Verifying whether the obtained solution is a valid preimage. With the roots discovered, the authors proceeded to verify the validity of the obtained solution. They did this by comparing the computed Legendre symbols to the fixed ones for the given instance. If a inconsistency was found between the computed symbol using their solution and the fixed symbol, they promptly terminated the verification process, indicating an invalid trial. Considering that we only need to compute the first Legendre symbol in each instance with a probability of $50\%$, the first two symbols with a probability of $25\%$, etc., we can expect to compute an average of 3 Legendre symbols for each trial before encountering an inconsistency.

3.2. Techniques to Skip SPN Rounds

In this section, we introduce a trick proposed by [35] which can help us skip two rounds without additional consumption when analyzing the permutation based on the SPN structure using the CICO problem.

Let permutation $\mathcal{P}:{\mathbb{F}}_p^n\to {\mathbb{F}}_p^n$ be s-secure against the CICO problem. We can split it into two permutations, ${\mathcal{F}}_0$ and ${\mathcal{F}}_1$, i.e., $\mathcal{P}={\mathcal{F}}_1\circ {\mathcal{F}}_0(·)$. Here, $V_u$ is a vector subspace spanned by $\{e_0,\dots ,e_{u-1}\}$. We denote the input state and output state of $\mathcal{P}$ as

$$\begin{array}{cc}\hfill & \hfill \mathit{X}=(X_0,X_1,\dots ,X_{u-1},A_0,\dots ,A_{n-u-1})\in V_u,\hfill \\ \hfill & \hfill \mathit{Z}=(Z_0,Z_1,\dots ,Z_{u-1},C_0,\dots ,C_{n-u-1})\in V_u,\hfill \end{array}$$

respectively. Here, $(A_0,\dots ,A_{n-u-1})\in {\mathbb{F}}_p^{n-u}$ and $(C_0,\dots ,C_{n-u-1})\in {\mathbb{F}}_p^{n-u}$ are fixed constants. According to the definition of the CICO problem, if we can find a pair of inputs $\mathit{X}$ and $\mathit{Z}$ such that $\mathcal{P}\left(\mathit{X}\right)=\mathit{Z}$ with a complexity less than $2^s$, we can conclude that the security margin of the permutation is insufficient.

We denote $\mathit{Y}=(Y_0,Y_1,\dots ,Y_{n-1})\in {\mathbb{F}}_p^n$ as the intermediate variable after ${\mathcal{F}}_0$. If $\mathit{Y}$ can be found to belong to the vector subspace $V_u$, a polynomial system with $n-u$ outputs can be constructed through ${\mathcal{F}}_1$, meaning that we can find its root. Ultimately, the value of $\mathit{X}$ can be obtained based on the value of $\mathit{Y}$, which is sufficient to solve the CICO problem. Then, to solve the CICO problem of permutation $\mathcal{P}$, only the ${\mathcal{F}}_1$ part needs to be dealt with, not the whole permutation $\mathcal{P}$.

To provide a detailed description of this technique, we assume that the permutation $\mathcal{P}$ corresponds to the Grendel permutation. ${\mathcal{F}}_0$ consists of two nonlinear layers, one linear layer, and one round key addition in the Grendel round function, while ${\mathcal{F}}_0$ can be expressed as ${\mathcal{F}}_0(·)=\mathcal{NL}\circ \mathcal{AC}\circ \mathcal{L}\circ \mathcal{NL}(·)$ and ${\mathcal{F}}_1$ can be regarded as an $R-2$ round Grendel round function with a linear layer and a round key addition; moreover, S represents the S-box while $S^{-1}$ denotes its inverse.

The S-box needs to satisfy the following property:

$$\begin{array}{c}\hfill S(A·X)=S\left(A\right)·S\left(X\right),\end{array}$$

(4)

where $A,X\in {\mathbb{F}}_p$.

Let the linear layer MDS matrix $\mathcal{M}$ satisfy

$$\begin{array}{c}\hfill {\mathcal{M}}^{-1}=\left[\begin{array}{ccccc}{m}_{0,0}& m_{0,1}& m_{0,2}& \dots & m_{0,n-1}\\ m_{1,0}& m_{1,1}& m_{1,2}& \dots & m_{1,n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ m_{n-1,0}& m_{n-1,1}& m_{n-1,2}& \dots & m_{n-1,n-1}\end{array}\right].\end{array}$$

The round constant is denoted as $c_j^i$ ($0\le i\le R-1,0\le j\le n-1$). Next, we show how to construct univariate equations with the CICO problem. For the following discussion, we set $u=n-1$ at all times.

When $n=3$, we set $u=n-1=2$; then, $V_u$ is a vector subspace spanned by $\{e_0,e_1\}$. Let the input states of ${\mathcal{F}}_0$ be $\mathit{X}=(X_1,X_2,A_0)\in V_u$, where $A_0$ is a fixed constant, and let the states after ${\mathcal{F}}_0$ be $\mathit{Y}=(Y_0,Y_1,Y_2)\in {\mathbb{F}}_p^3$. When passing through the first nonlinear layer of ${\mathcal{F}}_0$, we have

$$\begin{array}{cc}\hfill S\left(A_0\right)& =m_{2,0}(S^{-1}\left(Y_0\right)-c_0^0)+m_{2,1}(S^{-1}\left(Y_1\right)-c_1^0)+m_{2,2}(S^{-1}\left(Y_2\right)-c_2^0)\hfill \\ \hfill & =m_{2,0}{S}^{-1}\left(Y_0\right)+m_{2,1}{S}^{-1}\left(Y_1\right)+m_{2,2}{S}^{-1}\left(Y_2\right)-(m_{2,0}{c}_0^0+m_{2,1}{c}_1^0+m_{2,2}{c}_2^0).\hfill \end{array}$$

(5)

We fix $Y_2$ to a constant value $B_0=S\left(m_{2,2}^{-1}(m_{2,0}{c}_0^0+m_{2,1}{c}_1^0+m_{2,2}{c}_2^0+S\left(A_0\right))\right)$. Then, we can simplify Equation (5) as follows:

$$\begin{array}{cc}\hfill & m_{2,0}{S}^{-1}\left(Y_0\right)+m_{2,1}{S}^{-1}\left(Y_1\right)=0\hfill \\ ⟺\hfill & m_{2,0}{S}^{-1}\left(Y_0\right)=-m_{2,1}{S}^{-1}\left(Y_1\right)\hfill \\ ⟺\hfill & S\left(m_{2,0}{S}^{-1}\left(Y_0\right)\right)=S(-m_{2,1}{S}^{-1}\left(Y_1\right))\hfill \\ ⟺\hfill & S\left(m_{2,0}\right)Y_0=S\left(m_{2,1}\right)Y_1.\hfill \end{array}$$

(6)

The S-box must satisfy Formula (4) for the above Equation (6) to be established successfully. We find that $A_0$ and $B_0$ are fixed and that $Y_1$ can be represented by $Y_0$ as $Y_1=\frac{S\left(m_{2,0}\right)}{S\left(m_{2,1}\right)}{Y}_0$. Then, we have

$$\begin{array}{cc}\hfill & \mathit{X}=(X_0,X_1,A_0)\in V_2,\hfill \\ \hfill & \mathit{Y}=Y_0(1,\frac{S\left(m_{2,0}\right)}{S\left(m_{2,1}\right)},0)+(0,0,B_0)\in V_2.\hfill \end{array}$$

When $n=4$, we set $u=n-1=3$; then, $V_u$ is a vector subspace spanned by $\{e_0,e_1,e_2\}$. As seen in Figure 2, we denote the input and output states of ${\mathcal{F}}_0$ as $\mathit{X}=(X_0,X_1,X_2,A_0)\in V_u$ and $\mathit{Y}=(Y_0,Y_1,Y_2,Y_3)\in {\mathbb{F}}_p^4$, respectively, where $A_0$ are fixed constants. When passing through the first nonlinear layer of ${\mathcal{F}}_0$, we have

$$\begin{array}{cc}S\left(A_0\right)\hfill & =m_{3,0}(S^{-1}\left(Y_0\right)-c_0^0)+m_{3,1}(S^{-1}\left(Y_1\right)-c_1^0)+m_{3,2}(S^{-1}\left(Y_2\right)-c_2^0)+m_{3,3}(S^{-1}\left(Y_3\right)-c_3^0)\hfill \\ \hfill & =m_{3,0}{S}^{-1}\left(Y_0\right)+m_{3,1}{S}^{-1}\left(Y_1\right)+m_{3,2}{S}^{-1}\left(Y_2\right)+m_{3,3}{S}^{-1}\left(Y_3\right)\hfill \\ & -(m_{3,0}{c}_0^0+m_{3,1}{c}_1^0+m_{3,2}{c}_2^0+m_{3,3}{c}_3^0)\hfill \\ \hfill & =\sum _{i=0}^3{m}_{3,i}{S}^{-1}\left(Y_i\right)-\sum _{i=0}^3{m}_{3,i}{c}_i^0.\hfill \end{array}$$

(7)

We fix $Y_3$ to a constant denoted as $B_0$, while $Y_3$ satisfies

$$\begin{array}{c}\hfill m_{3,3}{S}^{-1}\left(Y_3\right)=\sum _{i=0}^3{m}_{3,i}{c}_i^0+S\left(A_0\right).\end{array}$$

Then, we can obtain

$$\begin{array}{cc}\hfill & \hfill m_{3,0}{S}^{-1}\left(Y_0\right)+m_{3,1}{S}^{-1}\left(Y_1\right)+m_{3,2}{S}^{-1}\left(Y_2\right)=0.\hfill \end{array}$$

(8)

In order to simplify the equation, we set $(Y_1,Y_2)=(S\left(Q_1\right)Y_0,S\left(Q_2\right)Y_0)$ and bring $(Y_1,Y_2)$ into Equation (8); then, we have

$$\begin{array}{c}\hfill S^{-1}\left(Y_0\right)(m_{3,0}+m_{3,1}{Q}_1+m_{3,2}{Q}_2)=0.\end{array}$$

Therefore, if $(Q_1,Q_2)$ and $Y_3$ satisfy

$$\begin{array}{c}\left\{\begin{array}{cc}\hfill & m_{3,0}+m_{3,1}{Q}_1+m_{3,2}{Q}_2=0\hfill \\ \hfill & Y_3=S\left(m_{3,3}^{-1}(\sum _{i=0}^3{m}_{3,i}{c}_i^0+S\left(A_0\right))\right),\hfill \end{array}\right.\hfill \end{array}$$

we have

$$\begin{array}{cc}\hfill & \mathit{X}=(X_0,X_1,X_2,A_0)\in V_3,\hfill \\ \hfill & \mathit{Y}=Y_0(1,Q_1,Q_2,0)+(0,0,0,B_0)\in V_3.\hfill \end{array}$$

When $n\ge 4$, we set $u=n-1$ in general, while $V_{n-1}$ is a vector subspace spanned by $\{e_0,e_1,\dots ,e_{n-2}\}$. Similarly, the input and output states of ${\mathcal{F}}_0$ are in the form of $\mathit{X}=(X_0,X_1,\dots ,X_{n-2},A_0)$ and $\mathit{Y}=(Y_0,Y_1,\dots ,Y_{n-1})\in {\mathbb{F}}_p^n$, respectively. Let $A_0\in {\mathbb{F}}_p$ be a fixed constant. When passing through the first nonlinear layer of ${\mathcal{F}}_0$, we have

$$\begin{array}{c}\hfill S\left(A_0\right)=\sum _{i=0}^{n-1}{m}_{n-1,i}(S^{-1}\left(Y_i\right)-c_i^0)=\sum _{i=0}^{n-1}{m}_{n-1,i}{S}^{-1}\left(Y_i\right)-\sum _{i=0}^{n-1}{m}_{n-1,i}{c}_i^0.\end{array}$$

(9)

We can fix $Y_{n-1}$ to a constant denoted as $B_0$; then, $Y_{n-1}$ fulfills

$$\begin{array}{c}\hfill m_{n-1,n-1}{S}^{-1}\left(Y_{n-1}\right)=\sum _{i=0}^{n-1}{m}_{n-1,i}{c}_i^0+S\left(A_{n-1}\right).\end{array}$$

(10)

Just as for $n=3$ and $n=4$, we set $(Y_1,Y_2,\dots ,Y_{n-2})=(S\left(Q_1\right)Y_0,S\left(Q_2\right)Y_0,\dots ,S\left(Q_{n-2}\right)Y_0)$. By bringing $(Y_1,Y_2,\dots ,Y_{n-1})$ and the constant $Y_{n-1}$ back into the Equation (9), we can obtain

$$\begin{array}{c}\hfill S^{-1}\left(Y_0\right)(m_{n-1,0}+\sum _{i=1}^{n-1}{m}_{n-1,i}{Q}_i)=0.\end{array}$$

Therefore, if $(Q_1,Q_2,\dots ,Q_{n-1})$ and $Y_{n-1}$ satisfy

$$\begin{array}{c}\left\{\begin{array}{cc}\hfill & m_{n-1,0}+\sum _{i=1}^{n-1}{m}_{n-1,i}{Q}_i=0\hfill \\ \hfill & Y_{n-1}=S\left(m_{n-1,n-1}^{-1}(\sum _{i=0}^{n-1}{m}_{n-1,i}{c}_i^0+S\left(A_0\right))\right),\hfill \end{array}\right.\hfill \end{array}$$

we have

$$\begin{array}{cc}\hfill & \mathit{X}=(X_0,X_1,\dots ,X_{n-2},A_0)\in V_{n-1},\hfill \\ \hfill & \mathit{Y}=Y_0(1,Q_1,Q_2,\dots ,Q_{n-2},0)+({\mathbf{0}}^{n-1},B_0)\in V_{n-1}.\hfill \end{array}$$

Let $\mathit{Y}$ be the input to ${\mathcal{F}}_1$, where $Y_0$ is the only unknown variable. We define the output of ${\mathcal{F}}_1$ as $\mathit{Z}=(Z_0,Z_1,\dots ,Z_{n-2},C_0)\in V_{n-1}$, with $C_0\in {\mathbb{F}}_p$ being a fixed constant. Considering the final position of the output from ${\mathcal{F}}_1$, a univariate equation is constructed with $Y_0$ as its variable, taking the form of

$$\begin{array}{c}\hfill F\left(Y_0\right)=C_0.\end{array}$$

(11)

With a valid $Y_0$, we can invariably infer an input $\mathit{X}$ specifically tailored for the R-round permutation $\mathcal{P}$ that projects onto the vector subspace $V_{n-1}$.

3.3. Application to Grendel Hash Function

In this section, we build upon the full-round preimage attack on the Grendel hash function presented in [22] by employing the trick described in the previous section to decrease the degree and complexity of the polynomial system. Consider a security level denoted by s and a prime p that defines the field. We limit ourselves here to the case in which $p\ge 2^s$. The following are the details of our attack:

• We first divide the Grendel permutation into two parts, ${\mathcal{F}}_0$ and ${\mathcal{F}}_1$, as before. The Grendel permutation consists of R rounds. Consider the Grendel hash function with the parameters $n=r+c$. The Grendel S-box, denoted as $S\left(x\right):x\mapsto x^d·{\chi }_p\left(x\right)$, satisfies Formula (4), which can be easily proven using Proposition 1. Similarly, we set $u=n-1$ and let $V_u$ be a vector subspace. The Grendel permutation takes an input $\mathit{X}=(X_0,\dots ,X_{r-1},{\mathbf{0}}^c)$, where $X_0,\dots ,X_{r-1}$ represent the input messages, and produces an output $\mathit{Z}=(Z_0,\dots ,Z_{r-1},{\mathbf{0}}^c)$. The initial value IV of Grendel is set to all zeros, and the last c elements of the output $\mathit{Z}$ are also zeros. Consequently, when $\mathit{X}\in V_u$ passes through ${\mathcal{F}}_0$ it results in $\mathit{Y}=(Y_0,Y_1,\dots ,Y_{n-1})\in V_u$. As stated in the previous section, we have $Y_{n-1}$, and $(Q_1,\dots ,Q_{n-1})$ satisfy

  $$\begin{array}{c}\left\{\begin{array}{c}{m}_{n-1,0}+\sum _{i=1}^{n-1}{m}_{n-1,i}{Q}_i=0\hfill \\ Y_{n-1}={\left(m_{n-1,n-1}^{-1}\left(\sum _{i=0}^{n-1}{m}_{n-1,i}{c}_i^0\right)\right)}^3·{\chi }_p\left(m_{n-1,n-1}^{-1}\left(\sum _{i=0}^{n-1}{m}_{n-1,i}{c}_i^0\right)\right).\hfill \end{array}\right.\hfill \end{array}$$
  Thus, for ${\mathcal{F}}_1$ there is only one unknown input variable $Y_0$. The subsequent processing can be carried out in a similar manner as described in [22].
• According to [22], it can be observed that when $p\ge 2^{32}$, the probability of the Legendre symbol being $\pm 1$ is greater than $99.99\%$. Therefore, we only consider guessing $\pm 1$. Based on the previous step, ${\mathcal{F}}_1$ has an input $\mathit{Y}$ with only one unknown variable $Y_0$. The ${\mathcal{F}}_1$ has $R-2$ rounds; we must guess the number of Legendre symbols, which is provided by $l=n(R-3)+1=nR-3n+1$. Because the Legendre symbol of $Y_0$ only needs to be guessed in the first round of ${\mathcal{F}}_1$, while the other values are constant, the Legendre symbol is known. Consequently, there are at most $2^l=2^{nR-3n+1}$ distinct sets of Legendre symbols to guess until the correct set of Legendre symbols is found.
• After fixing the Legendre symbols, we can construct a polynomial with $Y_0$ as an unknown variable. The polynomial equation, as defined in Formula (11), has a degree of $\mathcal{D}=d^{R-2}$. To determine the specific value of $Y_0$ we can employ the root-finding algorithm in Section 2.3.1. The complexity $T_1$ of the root-finding algorithm is

  $$\begin{array}{cc}\hfill & T_1=\mathcal{O}(\mathtt{M}\left(d^{R-2}\right)log\left(d^{R-2}\right)log(d^{R-2}·p)),\hfill \\ \hfill & \mathtt{M}\left(d^{R-2}\right)=63.43·d^{R-2}log\left(d^{R-2}\right)log(log\left(d^{R-2}\right))+\mathcal{O}(d^{R-2}log\left(d^{R-2}\right)).\hfill \end{array}$$
• Upon obtaining the value of $Y_0$, we need to verify its validity. This requires checking the correctness of each guessed Legendre symbol. According to [22], for each set of guessed Legendre symbols we only need to verify three of them to exclude an invalid set. The complexity $T_2$ of computing a Legendre symbol [46] is evaluated as $\mathcal{O}(\sigma {(log\sigma )}^{2}log(log\left(\sigma \right)))$ for $\sigma =log\left(p\right)$; therefore, the complexity of this step is $3·T_2$.

• Upon obtaining a valid $Y_0$ according to the CICO definition, we can always deduce $\mathit{X}$ such that they are mapped to the vector subspace $V_u$ through the Grendel permutation. The overall computational complexity of this attack, denoted as T, is represented by

  $$\begin{array}{c}\hfill T=(T_1+3·T_2)·2^{nR-3n+1}.\end{array}$$

  Therefore, this particular instance is vulnerable to attack if $T\le 2^s$. As summarized in Table 1, for a security level of $s=128$ it can be observed that under different parameter settings we are able to perform two additional attacking rounds compared to the previous work [22]. However, our advances do not exceed the newly established security margin in [22].

In our investigation of the Grendel hash function, we have ascertained that capitalizing on the CICO problem to devise a univariate equation is feasible solely under the conditions $u=r=n-1$. This premise holds because in this specific scenario it enables the generation of an intermediate variable intimately associated with the vector subspace $V_u$, ensuring that the hash output remains confined to this subspace.

3.4. The Gröbner Basis Attacks for the Grendel Hash Function

In this section, we employ the CICO problem to construct a multivariate equation system for the Grendel hash function instantiation. Similarly, we consider the message absorption size to be r, resulting in r hash digests being squeezed out after a Grendel permutation. To further analyze the complexity, we utilize the Gröbner basis attack method described in Section 2.3.2 and incorporate insights gained from our experimental observations.

Building upon our previous assumption of guessing the Legendre symbols, we delve deeper into the analysis by considering the introduction of intermediate variables to reduce the degree of the polynomials. Based on the presence of intermediate variables, we categorize our attacks into two scenarios: one without the introduction of intermediate variables, and another in which intermediate variables are introduced in each round.

Considering the Grendel hash function with input messages $(X_1,X_2,\dots ,X_{r-1})$ of size r and an IV set to all zeros, the Grendel permutation takes an input $\mathit{X}=(X_1,X_2,\dots ,X_{r-1},{\mathbf{0}}^c)$, where ${\mathbf{0}}^c$ denotes a vector of zeros with length c. The resulting output $\mathit{Z}=(Z_0,\dots ,Z_{n-1})$ is subject to the CICO problem, where the input $\mathit{X}$ belongs to the vector subspace $V_r$ spanned by $\{e_0,\dots ,e_{r-1}\}$. To satisfy this condition, the last c positions of $\mathit{Z}$, denoted as $C_0,\dots ,C_{c-1}$, are fixed constants. In the following attacks, we always set $r=c=\frac{2}{n}$. Consequently, we can construct a system of multivariate equations.

• Without Intermediate Variables. Let $\mathit{X}$ and $\mathit{Z}$ be the input and output of the permutation. No additional intermediate variables, such as $\mathit{Y}$, are introduced. We can build an equation system with c variables and c equations:

  $$\begin{array}{c}\left\{\begin{array}{cc}\hfill & F_0(X_0,X_1,\dots ,X_{c-1})=C_0\hfill \\ \hfill & F_1(X_0,X_1,\dots ,X_{c-1})=C_1\hfill \\ \hfill & ⋮\hfill \\ \hfill & F_{c-1}(X_0,X_1,\dots ,X_{c-1})=C_{c-1}.\hfill \end{array}\right.\hfill \end{array}$$

  We obtain a system of equations with c equations and c variables, where each equation has a degree of ${\mathcal{D}}_i=d^R,0\le i\le c-1$. It is evident that the degrees of the equations is much larger than the number of variables. Therefore, Formula (3) is used to calculate the complexity of the Gröbner basis algorithm. Specifically, when setting $d=2$, we can compare the computational complexities of the Gröbner basis and FGLM algorithms. For a system of c equations in which each equation has a degree of $2^R$, we can compute the upper bound on the regularity degree ${\mathcal{D}}_{reg}$ of the equation system and the zero-dimensional ideal ${\mathcal{D}}_{\mathcal{I}}$ as follows:

  $$\begin{array}{ccc}\hfill & \hfill {\mathcal{D}}_{reg}\le 1+\sum _{i=0}^{c-1}({\mathcal{D}}_i-1)=(2^R-1)·c+1,\hfill & \hfill {\mathcal{D}}_{\mathcal{I}}\le \prod _{i=0}^{c-1}{\mathcal{D}}_i=2^{Rc}.\end{array}$$

  Computing the Gröbner basis with respect to the grevlex term order using Formula (3) exhibits asymptotic complexity:

  $$\begin{array}{c}\hfill T_G=n{D}_{reg}·{\left(\genfrac{}{}{0pt}{}{n+D_{reg}-1}{D_{reg}}\right)}^{\omega }\le c·((2^R-1)·c+1)·\left(\genfrac{}{}{0pt}{}{2^R·c}{2^R·c-c+1}\right).\end{array}$$

  Then, applying a fast variant of the FGLM algorithm to perform the change of term order exhibits asymptotic complexity:

  $$\begin{array}{c}\hfill T_F={\left({\mathcal{D}}_{\mathcal{I}}\right)}^{\omega }=2^{Rc\omega }.\end{array}$$

  By evaluating $T_G$ and $T_F$ for $c=2$, we find that

  $$\begin{array}{cc}\hfill & T_G=(2^{R+2}-2)\times \left(\genfrac{}{}{0pt}{}{2^{R+1}}{2^{R+1}-1}\right)={\left(2^{R+1}\right)}^{\omega }\times (2^{R+2}-2)\le 2^{Rw+w+r+2},\hfill \\ \hfill & T_F=2^{2Rw}.\hfill \end{array}$$

  Based on this, it is clear that $T_F$ is larger than $T_G$, implying that the FGLM algorithm becomes the bottleneck in the computation of Gröbner bases. As shown in

  Table 3. For a security margin of $s=128$ and a modulus $p=2^{256}$, by setting c = r = 2/n we can evaluate the number of rounds that can be susceptible to Gröbner basis attacks with varying computational complexities.

  Instance $(\mathit{d},\mathit{n})$	Attacked Rounds with ${\mathit{T}}_{\mathit{F}}$	Attacked Rounds with ${\mathit{T}}_{\mathit{G}}$
  (2, 4)	16	21
  (2, 8)	8	10
  (2, 12)	5	7
  (3, 4)	12	17
  (3, 8)	6	8
  (3, 12)	4	5
  (5, 4)	9	14
  (5, 8)	4	7
  (5, 12)	3	4

  for $p=2^{256}$ and $s=128$, by setting $c=r=2/n$, we can evaluate the number of rounds that can be susceptible to Gröbner basis attacks with varying computational complexities. Our practical tests regarding the actual degrees reached in the computation are given in Figure 3.

• Intermediate Variables. Using the input and output of the Grendel permutation directly may be infeasible due to the high degree and dense nature of the polynomials involved. To overcome this challenge, one possible strategy is to introduce intermediate variables. This approach reduces the degrees in the equation system, thereby reducing the number of monomials, although it introduces additional variables. In each round of Grendel permutation, we introduce new variables to prevent an increase in degrees (Figure 4). Let $\mathit{X}$ and ${\mathit{Y}}^0=(Y_0^0,\dots ,Y_{n-1}^0)\in {\mathbb{F}}_p^n$ represent the input and output of the nonlinear layer in the first round, respectively. The relationship between $\mathit{X}$ and ${\mathit{Y}}^0$ can be expressed through r equations of degree d and c equations of degree 1. Specifically, $Y_0^0=X_0^d$, in accordance with the definition of the S-box (excluding the consideration of the Legendre symbol, as it is determined based on the conjecture); consequently, we add n variables in each round. Then, we simply use the output values $C_0,\dots ,C_{c-1}$ to construct the system of equations except for the last one, meaning that we have $r+Rn$ variables and the same number of equations. Among these equations, there are $Rn$ equations with a degree of d and r equations with a degree of 1.

When n is large, indicating the presence of a greater number of intermediate variables, then $D_{reg}$ becomes relatively small. Therefore, we utilize Formula (2) to evaluate the complexity of the Gröbner basis attack.

In summary, the complexity of the Gröbner basis attack on the Grendel hash function can be divided into three parts. The first part is the complexity of guessing the Legendre symbols, denoted as $T_{guess}$. The second part is the complexity of the Gröbner basis attack, denoted as $T_{GB}$ (with specific calculations selected from the various scenarios mentioned earlier). The third part is the complexity of verifying the Legendre symbols, denoted as $T_{verify}$. Lastly,

$$\begin{array}{cc}\hfill & T_{guess}=2^{(R-1)n+c},\hfill \\ \hfill & T_{verify}=3·\mathcal{O}(\sigma {(log\sigma )}^{2}log(log\left(\sigma \right)))\mathrm{for}\sigma =log\left(p\right).\hfill \end{array}$$

The overall complexity of the Gröbner basis attack for the Grendel hash function can be evaluated by

$$\begin{array}{c}\hfill (T_{GB}+T_{verify})·T_{guess}.\end{array}$$

Our practical tests regarding the actual degrees reached in the computation are given in Figure 5.

4. Conclusions

In this paper, we propose a preimage attack on the sponge hash function implemented with full rounds of the Grendel permutation, utilizing algebraic approaches. By introducing the CICO problem, we investigate the construction of univariate and multivariate equation systems for the Grendel hash function and employ different algorithms to solve these equations, resulting in new analytical findings. This provides additional insights into the factors that designers should consider when developing arithmetization-oriented cryptographic primitives in response to the CICO problem. Moreover, our research highlights the influence that the selection of distinct algebraic methods for equation construction can have on the security analysis of cryptographic primitives.

• Further Discussion. It is worthwhile to investigate the potential of merging various algebraic methods for the analysis of arithmetization-oriented cryptographic primitives in the future. Similar to the algebraic techniques used in this study, we introduce the CICO problem and employ specific strategies to bypass one round of the cryptographic algorithm, enabling further equation construction or utilization of alternative techniques to build equation systems. This study has focused solely on the analysis of permutation in SPN structures; however, in the future it may be possible to extend this method to Feistel structures.