<figure><img src="image_1.png" /></figure>

HAL
open science

# The regular multivariate quadratic problem

Antoine Joux, Rocco Mora

► To cite this version:

Antoine Joux, Rocco Mora. The regular multivariate quadratic problem. Designs, Codes and Cryptography, 2025. hal-05449355

HAL Id: hal-05449355

https://inria.hal.science/hal-05449355v1

Submitted on 8 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# The regular multivariate quadratic problem

Antoine Joux, Rocco Mora

CISPA Helmholtz Center for Information Security - Germany

## Abstract

In this work, we introduce a novel variant of the multivariate quadratic problem, which is at the core of one of the most promising post-quantum alternatives: multivariate cryptography. In this variant, the solution of a given multivariate quadratic system must also be regular, i.e. if it is split into multiple blocks of consecutive entries with the same fixed length, then each block has only one nonzero entry. We prove the NP-completeness of this variant and show similarities and differences with other computational problems used in cryptography. Then we analyze its hardness by reviewing the most common solvers for polynomial systems over finite fields, derive asymptotic formulas for the corresponding complexities and compare the different approaches.

# 1 Introduction

**The MQ problem and its use in cryptography.** The problem of solving a multivariate quadratic system over a finite field takes the name of Multivariate Quadratic (MQ) problem and is known to be NP-complete. It also seems to be hard on average for a wide range of parameters, but despite its evident difficulty, numerous algorithms have been developed to solve it. The MQ problem is fundamental in cryptography, especially in the post-quantum realm. It serves as the foundation for designing multivariate public key cryptosystems, which form one of the most promising quantum-resistant alternatives. In particular, it has emerged over the years as being remarkably well-suited to the construction of digital signatures. It is even more pervasive in algebraic cryptanalysis, which has a crucial impact on many different families of primitives. Indeed, algebraic attacks aim at solving a polynomial system that models a cryptographic problem.

**A new variant of the MQ problem.** We present here a new variant of the MQ problem. Whilst the polynomial equations are sampled at random and do not present any particular structure, the solution is required to have a regular shape. Concretely, if the solution is a vector of $n = lw$ variables split into $w$ blocks each of length $l$, then it must have exactly one nonzero entry for each block. The NP-completeness of the multivariate quadratic problem is preserved when adding the regular condition to the solution.

The same regular constraint was originally introduced in the context of code-based cryptography in [AFS05]. Indeed, [AFS05] proposed a family of cryptographic hash functions based on the provably secure regular syndrome decoding (RSD) problem. More recently, the same problem received new attention thanks to its introduction in secure computation to design MPC-protocols [HOSSV22], pseudorandom correlation generators [BCGI18] and digital signatures from MPC-in-the-head [CCJ23, BCC$^{+}$24] or VOLE-in-the-head [CLY$^{+}$24]. These works suggest that also the RMQ problem could be used for similar constructions. On the other hand, there have been attempts at improving the cryptanalysis, either with combinatorial techniques [CCJ23, ES24] or with algebraic modelings [BØ23]. In continuity with its code-based predecessor, the new variant has been named the Regular Multivariate Quadratic (RMQ) problem.

**Our contribution.** After introducing the computational problem and overviewing its main features, we mainly focus on the study of the asymptotic complexity of the average RMQ instance in the setting of the uniqueness bound. We analyzed the most important families of techniques for solving boolean systems and adapted them to take advantage of the regularity of the solution. On one side, we explore the algebraic algorithms related to Gröbner basis computation [CKPS00, BFSS13] by modeling the regular constraint as additional equations to be satisfied, reviewing different ways of hybridization and estimating the solving degree, in a similar fashion to what has been done for the code-based counterpart in [BØ23]. On the other hand, we investigate how the research line of probabilistic polynomial methods initiated with [LPT$^{+}$17] applies and can be optimized to our context. We also show how an RMQ instance can be transformed into a standard MQ one with fewer variables and analyze the complexity of solving the new system. Finally, we provide comparisons among all the approaches examined. Despite the improvements presented, we conclude that the problem remains difficult in practice in the average case for a large range of parameters and therefore it is potentially suitable for cryptographic purposes.
# 2 Notation and preliminaries

## 2.1 Notation

**Sets and intervals.** The closed interval of real numbers between $a$ and $b$ is denoted by $[a, b]$, while the corresponding interval of integers by $\llbracket a, b \rrbracket$ and if $a = 1$ we simply write $\llbracket b \rrbracket$. The cardinality of a set $S$ is indicated with $|S|$.

**Vectors, matrices and polynomials over finite fields.** Given a prime power $q$, we denote by $\mathbb{F}_q$ the finite field of size $q$. The length-$n$ linear space over $\mathbb{F}_q$ is $\mathbb{F}_q^n$ and the vector space of matrices of size $n \times m$ over $\mathbb{F}_q$ is $\mathbb{F}_q^{n \times m}$. Vectors (respectively matrices) are written with lowercase (respectively uppercase) bold letters, while their entries are not bold, for instance $\mathbf{v} = (v_1, \dots, v_n) \in \mathbb{F}_q^n$ and $\mathbf{A} = (a_{i,j}) \in \mathbb{F}_q^{n \times m}$. We denote by $(\llbracket n \rrbracket)$ (respectively $(\downarrow_w^n)$) the sets of vectors in $\mathbb{F}_2^n$ of weight exactly $w$ (respectively at most $w$). It holds that $|\downarrow_w^n| = \sum_{i \in \llbracket 0, w \rrbracket} \binom{n}{i}$ and we denote such sum with $\binom{n}{\downarrow_w w}$.

Finally, the polynomial ring in a vector of $n$ variables $\mathbf{x} = (x_1, \dots, x_n)$ over $\mathbb{F}_q$ is written as $R = \mathbb{F}_q[\mathbf{x}]$ and its polynomials as $P(\mathbf{x})$. We use calligraphic upper-case letters for polynomial ideals and angle brackets to denote ideals generated by sequences or sets of polynomials, e.g. $\mathcal{I} = \langle S \rangle = \langle f_1, \dots, f_m \rangle$.

**Regular vectors.** Let $n = lw$, for two positive integers $w$ and $l$. An $l$-regular vector $\mathbf{v} \in \mathbb{F}_q^n$ is a vector such that

$$ \forall i \in \llbracket 0, w-1 \rrbracket, \quad |\{v_{il+j} | j \in \llbracket 1, l \rrbracket, v_{il+j} \neq 0\}| = 1, $$

i.e. such that in each of the $w$ windows of entries of length $l$, exactly one coordinate is non-zero. The Hamming weight of an $l$-regular vector is thus $w = n/l$. We denote by $\mathcal{R}_l(\mathbb{F}_2^n)$ the set of such vectors and more in general by $\mathcal{R}_l(S)$ the set of regular vectors within the set $S$. We also say that a vector $\mathbf{v} \in \mathbb{F}_q^n$ is $l$-at most regular if

$$ \forall i \in \llbracket 0, w-1 \rrbracket, \quad |\{v_{il+j} | j \in \llbracket 1, l \rrbracket, v_{il+j} \neq 0\}| \le 1, $$

i.e. such that in each of the $w$ windows of entries of length $l$, at most one coordinate is non-zero. The Hamming weight of an $l$-regular vector is thus $\le w = n/l$. The set of all such vectors is indicated with $\mathcal{AM}\mathcal{R}_l(\mathbb{F}_2^n)$ and more in general the set of regular vectors within the set $S$ is $\mathcal{R}_l(S)$. We clearly have $\mathcal{R}_l(\mathbb{F}_2^n) \subset \mathcal{AM}\mathcal{R}_l(\mathbb{F}_2^n)$ and $|\mathcal{R}_l(\mathbb{F}_2^n)| = l^{n/l}$ and $|\mathcal{AM}\mathcal{R}_l(\mathbb{F}_2^n)| = (l+1)^{n/l}$. From now on, $l$ will always denote the length of a window, so we simply say “(at most) regular”, simplify the notation with $\mathcal{R}$ and $\mathcal{AM}\mathcal{R}$, and assume that $l$ divides the vector length $n$.

For the sake of clarity, we also denote the coordinates of a vector of length $n = lw$, by using two indexes; one for the window (starting from 0) and one for the position within the window (starting from 1). In other words, given a vector $\mathbf{v} \in \mathbb{F}_q^{lw}$, we use the notation

$$ v_{i,j} \stackrel{\text{def}}{=} v_{il+j}, \quad \forall i \in \llbracket 0, w-1 \rrbracket, j \in \llbracket 1, l \rrbracket. $$

This notation is applied to vectors of variables as well. We also use the term “block” as a synonym for “window”.

**Complexity and asymptotic.** Given a function $f$, we use the standard big O notation $O(f(n))$ to provide the cost of an algorithm or one of its steps, meaning that the cost is upper bounded by $cf(n)$, for a positive real constant $c$ and for any $n \ge n_0$, where $n_0$ is a fixed integer. The metric we adopt uses the number of operations over the target finite field $\mathbb{F}_q$ as unit. For exponential costs, we also make use of the notation $O^*(f(n))$ to denote the asymptotic behavior, up to polynomial factors in $n$, i.e. the cost is upper bounded by $P(n)f(n)$ for a fixed polynomial $P$ and for any $n \ge n_0$, where $n_0$ is a fixed integer.

The logarithm $\log(\cdot)$ has to be intended in base 2, while $\ln(\cdot)$ is the natural logarithm. Let $H: [0, 1] \to [0, 1]$ be the binary entropy function defined as

$$ H(0) = H(1) = 0 $$

and

$$ H(p) \stackrel{\text{def}}{=} -p \log(p) - (1-p) \log(1-p) \quad \text{otherwise.} $$

The binary entropy function is increasing in $[0, 1/2]$, reaches its maximum $H(1/2) = 1$ and it is decreasing in $[1/2, 1]$, being null at its extremes. We also define the non-decreasing function

$$ H^*(p) \stackrel{\text{def}}{=} \begin{cases} H(p), & p \in [0, 1/2) \\ 1, & p \in [1/2, 1]. \end{cases} $$

We have that

$$ \binom{n}{w} = \mathcal{O}^*(2^{H(w/n)n}) \quad \text{and} \quad \binom{n}{\downarrow w} = \mathcal{O}^*(2^{H^*(w/n)n}). $$

Finally, the Greek letter $2 \le \omega \le 3$ denotes the linear algebra constant. When sparse linear algebra is possible, quadratic-time algorithms allow to replace $w$ with 2.
## 2.2 Algebraic cryptanalysis

Let $R = \mathbb{F}[x_1, \dots, x_n] = \mathbb{F}[\mathbf{x}]$ be the polynomial ring over the finite field $\mathbb{F}$ in $n$ variables. A polynomial $f \in R$ is said homogeneous if all its monomials have the same degree and affine otherwise. We denote by $R_d$ the vector space of all homogeneous polynomials of degree $d$, which is generated by all degree-$d$ monomials. We also consider ideals generated by a set (or sequences) of polynomials in $R$. If an ideal $\mathcal{I} \subset R$ admits a set of homogeneous generators we say that the ideal is homogeneous. In this case, we denote with $\mathcal{I}_d \stackrel{\text{def}}{=} \mathcal{I} \cap R_d$, i.e. the subspace of $A_d$ consisting of all homogeneous polynomials of degree $d$ in $\mathcal{I}$.

The Hilbert function and the Hilbert series describe the combinatorial structure of homogeneous ideals (or equivalently the corresponding quotient ring) and represent a powerful tool to study the complexity of solving algebraic systems.

**Definition 1 (Hilbert function/series).** Let $\mathcal{I} \subseteq R$ be a homogeneous ideal. The Hilbert function is defined as

$$ \mathrm{HF}_{R/\mathcal{I}}: \mathbb{N} \to \mathbb{N} \\ d \mapsto \dim_{\mathbb{F}} R_d/\mathcal{I}_d $$

and the Hilbert series as its generating series, i.e.

$$ \mathrm{HS}_{R/\mathcal{I}}(z) \stackrel{\text{def}}{=} \sum_{d=0}^{\infty} \mathrm{HF}_{R/\mathcal{I}}(d)z^d. $$

Given an ideal $\mathcal{I}$, its variety $\mathcal{V}(\mathcal{I})$ is the set of common zeros for all the elements in $\mathcal{I}$ whose coordinates lie in the algebraic closure of $\mathbb{F}$. However, for cryptographic purposes, one usually wants to find solutions over the field where the constraints are defined or at least over some finite extension. The constraints that force solutions to lie in the finite field of $q$ elements $\mathbb{F} = \mathbb{F}_q$ are the so-called field equations:

$$ x_i^q - x_i = 0 \quad \forall i \in [1, n]. $$

When an ideal $\mathcal{I}$ contains the field equations, then it is always zero-dimensional, i.e. the quotient $R/\mathcal{I}$ is a finite-dimensional vector space. In this case, the Hilbert series is simply a polynomial, that is there exists a positive integer $d$ such that $\mathcal{I}_d = R_d$. We call degree of regularity $d_{reg}$ the smallest integer $d$ that satisfies the equality.

In the most general case, the Hilbert series of an ideal is difficult to compute and the cost of this task is essentially the same as computing a Gröbner basis of the associated system. On the other hand, the Hilbert series can be predicted for systems of equations that satisfy a certain notion of regularity.

**Definition 2 (Regular system).** Let $m \le n$. The homogeneous system $S = \{f_1, \dots, f_m\} \subseteq R = \mathbb{F}[x_1, \dots, x_n]$ is regular if

$$ \forall i \in [1, m], \quad f_i \text{ is not a zero divisor in } R/\langle f_1, \dots, f_{i-1} \rangle. $$

**Remark 1.** The notion of regularity of an algebraic system should not be confused with that of regular constraint of a solution, with which it only shares the term.

**Proposition 1.** Let $S = \{f_1, \dots, f_m\} \subseteq R$ be a regular system with $\deg(f_i) = d_i$ for all $i \in [1, m]$. Then

$$ \mathrm{HS}_{R/\langle S \rangle}(z) = \frac{\prod_{i=1}^{m} (1 - z^{d_i})}{(1 - z)^n}. $$

The definition of regularity does not cover the overdetermined case, i.e. when $m > n$. The notion has first been extended to this setting in [Bar04].

**Definition 3 (Semi-regular sequence).** Let $S = \{f_1, \dots, f_m\} \subseteq R$ be a homogeneous sequence such that $\mathcal{I} = \langle S \rangle$ is zero-dimensional and has degree of regularity $d_{reg}$. The sequence $S$ is called semi-regular if $\mathcal{I} \neq R$ and if for $i \in [1, m]$, $g_i f_i = 0$ in $R/\langle f_1, \dots, f_{i-1} \rangle$ with $\deg(g_i f_i) < d_{reg}$ implies $g_i = 0$ in $R/\langle f_1, \dots, f_{i-1} \rangle$.

A special definition is given for boolean systems, by taking into account the Frobenius morphism. To this purpose, we define $R' = \mathbb{F}_2[\mathbf{x}]/\langle x_1^2, \dots, x_n^2 \rangle$.

**Definition 4 (Semi-regular sequence over $\mathbb{F}_2$).** Let $S = \{f_1, \dots, f_m\} \subseteq R'$ be a homogeneous sequence such that $\mathcal{I} = \langle S \rangle$ has degree of regularity $d_{reg}$. The sequence $S$ is called semi-regular over $\mathbb{F}_2$ if $\mathcal{I} \neq R'$ and if for $i \in [1, m]$, $g_i f_i = 0$ in $R'/\langle f_1, \dots, f_{i-1} \rangle$ with $\deg(g_i f_i) < d_{reg}$ implies $g_i = 0$ in $R'/\langle f_1, \dots, f_{i-1} \rangle$.

The Hilbert series is determined for semi-regular sequences as well:
**Proposition 2.** Let $S = \{f_1, \dots, f_m\}$ be a homogeneous semi-regular sequence with $\deg(f_i) = d_i$ for all $i \in [1, m]$. Then

$$HS_{R/\langle S \rangle}(z) = \left[ \frac{\prod_{i=1}^{m} (1 - z^{d_i})}{(1 - z)^n} \right]_+,$$

where $[\cdot]_+$ denotes the truncation of the series starting from the first non-positive coefficient.

Let $S = \{f_1, \dots, f_m\}$ be a homogeneous semi-regular sequence over $\mathbb{F}_2$ with $\deg(f_i) = d_i$ for all $i \in [1, m]$. Then

$$HS_{R'/\langle S \rangle}(z) = \left[ \frac{(1+z)^n}{\prod_{i=1}^{m} (1+z^{d_i})} \right]_+.$$

We remark that the fact that the Hilbert series is a polynomial (because of the truncation) is in accordance with the ideal being zero-dimensional.

## 2.3 Affine systems

In cryptography, we often deal with affine sequences. In this case, we say that a sequence is semi-regular if $S^{(h)} \stackrel{\text{def}}{=} \{f_1^{(h)}, \dots, f_m^{(h)}\}$ is, where $f^{(h)}$ denotes the part of highest degree of $f$. It is conjectured (see Fröberg conjecture [Frö85]) that most systems behave in this way. Informally speaking, if one picks a random polynomial system of $m \le n$ ($m > n$ respectively), then it is regular (semi-regular respectively) except with negligible probability. On the other hand, if a system has some non-random structure, then the results on the Hilbert series may not hold anymore. This is the case for the systems considered in this work, where the regular constraint on the solution can be translated in terms of quadratic monomial (and thus non-random) equations. As a workaround, we can split our system into two parts, in the same fashion as done in [BØ23]. One of them is considered to be generic, and the other, expressing the regular constraint, is analyzed separately. We thus make use of an assumption that captures the notion of semi-regularity in our context and allows to derive the Hilbert series for the final ideal.

To solve the affine system, we make use of the linearization technique of the XL algorithm [CKPS00], which consists in computing the Macaulay matrix of the system $S = \{f_1, \dots, f_m\}$ at some sufficiently high degree $d$ such that there are more independent equations than monomials. We recall the definition of the Macaulay matrix for both homogeneous and affine systems.

**Definition 5** (Macaulay Matrix [Mac94]). Let $S = \{f_1, \dots, f_m\} \subset \mathbb{F}[x]$ be an algebraic system such that $\deg(f_i) = d_i$ and let $d$ be a positive integer. If $S$ is homogeneous, the (homogeneous) Macaulay matrix $\mathrm{Mac}(S, d)$ of $S$ in degree $d$ is the matrix whose rows are each indexed by a polynomial $m_j f_i$, for any $f_i \in S$ and any monomial $m_j$ of degree $d - d_i$, and whose columns are indexed by all the monomials of degree $d$. If $S$ is an affine system, then the (affine) Macaulay matrix $\mathrm{Mac}(S, \le d)$ is the matrix whose rows are each indexed by a polynomial $m_j f_i$, for any $f_i \in S$ and any monomial $m_j$ of degree $\le d - d_i$, and whose columns are indexed by all the monomials of degree $\le d$.

Echelonizing the Macaulay matrix can be interpreted as taking polynomial combinations of elements in $S$ and it is not surprising that it plays an important role in Grobnër basis techniques as well, with which the XL algorithm is indeed strictly related. Moreover, we take advantage of the Wiedemann algorithm [Wie86, Cop94], which exploits the matrix sparseness. The cost of block Wiedemann XL algorithm is estimated in many works of this genre to be

$$3 \cdot r_{\mathrm{AVG}} \cdot \#\mathrm{cols}(\mathrm{Mac}(S, \le d))^2, \qquad (1)$$

where $r_{\mathrm{AVG}}$ is the average row weight of elements of the Macaulay matrix (which coincides with the average number of terms in polynomials of $S$). On the other hand, $\#\mathrm{cols}(\mathrm{Mac}(S, \le d))$ is the number of columns of the Macaulay matrix and it is upper bounded by $\sum_{i=0}^{d} \binom{n}{i}$ in the binary case (i.e. when reducing by the quadratic field equations), and by $\sum_{i=0}^{d} \binom{n+i-1}{i}$ for generic fields. From an asymptotic point of view, Equation (1) means that it is possible to replace the linear algebra constant $\omega$ with the value 2.

For affine systems, especially structured ones, the degree of regularity may not always capture the difficulty of finding a solution, exactly because it ignores the part of the system of lower degree. In this context, a more informative object is the witness degree, originally defined in the binary case in [BFSS13] but here generalized to whatever finite field, as done in [BØ23].

**Definition 6** (Witness degree (adapted from [BFSS13]))). Let $S = \{f_1, \dots, f_m\} \subseteq \mathbb{F}[x]$ be an (affine) sequence of polynomials and $\mathcal{I} = \langle S \rangle$ the ideal that it generates. Let $\mathcal{I}_{\le d}$ and $\mathcal{J}_{\le d}$ be the $\mathbb{F}$-vector spaces defined by

$$\mathcal{I}_{\le d} \stackrel{\text{def}}{=} \{f \in \mathcal{I} \mid \deg(f) \le d\}$$
and

$$ \mathcal{J}_{\le d} \stackrel{\text{def}}{=} \{f \in \mathcal{I} \mid f = \sum_{i=1}^{m} g_i f_i \wedge \forall i \in [1, m], \deg(f_i g_i) \le d\}. $$

The witness degree $d_{wit}$ of $S$ is the smallest integer $d_0$ such that $\mathcal{I}_{\le d} = \mathcal{I}_{\le d}$ and $\langle \text{LM}(f) \mid f \in \mathcal{I}_{\le d_0} \rangle = \text{LM}(\mathcal{I})$, where $\text{LM}(f)$ denotes the leading monomial of $f$ according to some graded order and $\text{LM}(\mathcal{I})$ is the set of leading monomials of all nonzero element in $\mathcal{I}$.

Taking the part of highest degree in each polynomial is one way to obtain a homogeneous system from an affine one. Another strategy consists of adding a variable $y$ and, given an affine polynomial $f$, defining the polynomial

$$ f^{(y)}(x_1, \dots, x_n, y) = y^{\deg(f)} f\left(\frac{x_1}{y}, \dots, \frac{x_n}{y}\right). $$

Then $f^{(y)}$ is a homogeneous polynomial of degree $\deg(f)$ in $n+1$ variables and we call $S^{(y)} = \{f_1^{(y)}, \dots, f_m^{(y)}\}$ the homogenization of $S$. Then the degree of regularity of $S^{(y)}$ (which is in general different from that of $S$) provides under some conditions an upper bound on the witness degree. Again, we provide the adaptation of [BØ23] to general finite fields of the original results of [BFSS13], even if our main interest is the field $\mathbb{F}_2$.

**Proposition 3.** Let $S = \{f_1, \dots, f_m, x_1^q - x_1, \dots, x_n^q - x_n\}$ be a sequence of polynomials that admits no solutions, and let $\mathcal{I}^{(y)} \stackrel{\text{def}}{=} \langle S^{(y)} \rangle$. Then $d_{wit}(S) \le d_{reg}(\mathcal{I}^{(y)})$.

The proposition above explains why we estimate the complexity of solving the RMQ problem by deriving the degree of regularity from the Hilbert series of the homogenized system $S^{(y)}$. The condition in Proposition 3 of having no solutions is in apparent contradiction with the problem. However, it is consistent with the hybrid strategy of the BooleanSolve algorithm [BFSS13] because, for the almost totality of partial guesses, the specialized systems do not admit any solution indeed. The same happens when applying a hybrid approach to an RMQ instance.

# 3 The RMQ problem

We introduce here a new computational problem, called regular multivariate quadratic problem, which is a variant of the classical multivariate problem, but where only regular solutions to the system are admitted. We provide both the decision and search versions of the problem

**Definition 7** (l-Regular Multivariate Quadratic (l-RMQ) problem). Let $l \ge 2, w \ge 1$ and $m \ge 1$ be three integers, $q$ a prime power and $n = wI$. Let $P_1, \dots, P_m$ be $m$ affine quadratic polynomial in $n$ variables $x = (x_1, \dots, x_n)$ over the finite field $\mathbb{F}_q$ given as input. Let

$$ \{P_i(x) = 0 \mid i \in [1, m]\} \qquad (2) $$

be the associated quadratic system. The l-RMQ problem is:

* (Decision version) Determine whether the system (2) has an l-regular solution over $\mathbb{F}_q$

* (Search version) Find an l-regular solution over $\mathbb{F}_q$ of (2), if there is any.

As for the vectors, we will remove the parameter $l$ and simply write “RMQ problem” whenever the window length is not relevant. In this article, we mainly focus on the case of boolean system, i.e. when the underlying finite field is $\mathbb{F}_2$ and we sometimes refer to the RMQ problem over $\mathbb{F}_2$ as simply the RMQ problem. Moreover, for what concerns the asymptotic analysis of the problem, we are interested in the setting where the length of the window $l$ is constant and hence the number of windows $w$ grows linearly with $n$, which is the regime where the problem is supposed to be most difficult.

## 3.1 On the shape of the quadratic equations for the RMQ problem

Let

$$ P(x) = \sum_{1 \le i < j \le n} a_{i,j} x_i x_j + \sum_{1 \le i \le n} a_i x_i + a $$

be a generic quadratic polynomial. When $P$ is evaluated on a regular (or even at most regular) vector $v$, the monomials $x_i x_j$'s, for any pair $(i, j)$ such that $(i-1) \div l = (j-1) \div l$, always vanish, regardless of the coefficients
$a_{i,j}$'s, because at least one between $v_i$ and $v_j$ is 0. It is then enough to consider the monomials other than those and redefine, by using the regular-friendly notation,

$$P(\boldsymbol{x}) = \sum_{\substack{0 \le i_1 < i_2 \le w-1 \\ 1 \le j_1, j_2 \le l}} a_{i_1, i_2}^{(j_1, j_2)} x_{i_1, j_1} x_{i_2, j_2} + \sum_{\substack{0 \le i \le w-1 \\ 1 \le j \le l}} a_i^{(j)} x_{i,j} + a. \quad (3)$$

The memory size of a single quadratic polynomial for the multivariate quadratic problem thus decreases from $\binom{n+1}{2} + 1$ to

$$\binom{n+1}{2} + 1 - \frac{n}{l} \binom{l}{2} = \frac{n(n-l+2)}{2} + 1.$$

## 3.2 On the generation of RMQ instances

It is possible to construct quadratic systems whose variety contains a given regular vector. Even more, we can construct such a system in a way that all the quadratic equations are indistinguishable from random ones with the shape discussed above. Indeed, let $\boldsymbol{v} \in \mathcal{R}(\mathbb{F}_2^n)$. Let us sample at random over $\mathbb{F}_2$ all the coefficients $a_{i_1, i_2}^{(j_1, j_2)}$'s and $a_i^{(j)}$'s of $P$ in (3). Then we define $a = \sum_{\substack{0 \le i_1 < i_2 \le w-1 \\ 1 \le j_1, j_2 \le l}} a_{i_1, i_2}^{(j_1, j_2)} v_{i_1, j_1} v_{i_2, j_2} + \sum_{\substack{0 \le i \le w-1 \\ 1 \le j \le l}} a_i^{(j)} v_{i,j}$, so that whatever is the value of $\sum_{1 \le i < j \le n} a_{i,j} v_i v_j + \sum_{1 \le i \le n} a_i v_i$, we obtain $P(\boldsymbol{v}) = 0$. One can easily verify that, if all the coefficients have been sampled at random then, for a fixed $\boldsymbol{v}$,

$$\mathrm{Pr}[a=0] = \mathrm{Pr} \left[ \sum_{1 \le i < j \le n} a_{i,j} v_i v_j + \sum_{1 \le i \le n} a_i v_i = 0 \right] = 1/2.$$

On the other hand, it is interesting to study for which parameters $n$ and $m$, a regular solution to a random quadratic system in $n$ variables and $m = \mu n$ equations is expected. Let us first recall this expectation in the standard setting. The overdetermined regime is when $m > n$. If $\mu > 1$, for $n \to \infty$, solutions are not attended, unless with exponentially small probability in $m-n$. On the opposite side, when the system is underdetermined, i.e. when $m < n$, an exponential number of solutions is expected when $n \to \infty$. The intermediate case $\mu = 1$ is when the expectation of the number of solutions is 1.

When instead we search for regular solutions, the threshold for the unique solution expected corresponds to a smaller ratio $\mu$ of equations/variables. More precisely, the probability that any vector is a solution of an equation is $1/2$, hence it has probability $1/2^m$ of being a solution of the system. However, the number of regular vectors over $\mathbb{F}_2$ is only $l^w$. The mentioned threshold is thus obtained by equating

$$\frac{l^w}{2^m} = 1,$$

which gives

$$\mu = \frac{m}{n} = \frac{\log(l^w)}{n} = \frac{w \log(l)}{n} = \frac{\log(l)}{l}.$$

Our analysis, especially the cryptanalysis through Gröbner basis computation, then mainly focuses on this case.

If instead we consider the same problem over a bigger field $\mathbb{F}_q$, then equating $((q-1)l)^w = q^m$ provides the determined case for

$$\mu = \frac{\log_q((q-1)l)}{l}.$$

## 3.3 NP-completeness of the RMQ problem

The decision problem of determining whether a boolean multivariate polynomial systems has solutions over $\mathbb{F}_2$ (respectively over its closure) is NP-complete (respectively NP-hard). This result has been proved in [FY79] by reducing the 3-SAT problem to it. The result can be easily adapted to polynomials of degree 2, i.e. the decision MQ problem is NP-complete as well. In the same work where it has been first introduced [AFS05], the regular syndrome decoding problem was also shown to be NP-complete. The polynomial-time reduction is an adaptation of the classic one from the Three Dimensional Matching problem [BMVT78].

It is then not surprising that the decision RMQ problem (and consequently the search counterpart) is NP-complete, too. We provide here an elementary reduction from the MQ problem itself in the boolean case.

**Proposition 4.** For any $l \ge 2$, the decision $l$-RMQ problem is NP-complete.
Proof. The decision l-RMQ problem is clearly in NP, as a candidate solution can be tested by 1) evaluating all the quadratic polynomials, and 2) checking the l-regularity. Both tasks can be done in polynomial time in the parameters.

On the other hand, we prove the completeness by reducing the MQ problem to the l-RMQ problem in polynomial time. Let $S = \{f_1, \dots, f_m\} \subseteq \mathbb{F}_2[x_1, \dots, x_n] = R$ be a multivariate quadratic system. Let $R' = \mathbb{F}_2[x_{1,1}, \dots, x_{1,l}, \dots, x_{n,1}, \dots, x_{n,l}]$ and consider the natural embedding $\phi: R \hookrightarrow R'$ that sends $x_i$ to $x_{i,1}$. Let $f'_j \stackrel{\text{def}}{=} \phi(f_j)$ for all $j \in [1, m]$ and $S' = \{f'_1, \dots, f'_m\} \subseteq R'$.

Now let $v = (v_1, \dots, v_n)$ be a solution for the system associated with $S$. Define

$$v' \stackrel{\text{def}}{=} (v_1, v_1 + 1, 0, \dots, 0, v_n, v_n + 1, 0, \dots, 0) \in \mathbb{F}_2^{nl}.$$

Then, for any $j \in [1, m]$, $f'_j(v') = f(v) = 0$, i.e. $v'$ is a solution for the system associated with $S'$. Moreover $v'$ is l-regular, since in each block

$$(v_i, v_i + 1, 0, \dots, 0)$$

exactly one entry (one between the first and the second) is nonzero. Vice versa, let $v' = (v_{1,1}, \dots, v_{1,l}, \dots, v_{n,1}, \dots, v_{n,l})$ be a regular solution for $S'$ and $v \stackrel{\text{def}}{=} (v_{1,1}, \dots, v_{n,1}) \in \mathbb{F}_2^n$. Then, by construction of $f'_j$, $f_j(v) = f'_j(v') = 0$, for any $j \in [1, m]$ which implies that $v$ is a solution for $S$. This shows that the decision MQ problem for $S$ and the decision RMQ problem for $S'$ have either both positive answer or both negative answer. Since $S'$ can be constructed in polynomial time from $S$, the reduction is polynomial-time. $\square$

**Remark 2.** The simple reduction provided introduces a linear blow-up in the window length $l$. This, however, is not indicative of how the hardness of the MQ problem diverges when adding the regular constraint and the blow-up can probably be improved. The next subsection shows that in the regime of unique solution, the exhaustive search changes only by a factor $l/\log(l)$ in the exponent of the complexity. The analysis of more advanced methods for the average-case problem will then appear in the following sections.

## 3.4 Exhaustive search

The most naive approach for solving an algebraic over a (small) finite field is to evaluate it over all elements in the vector space until a zero that is common to all equations is found. For standard multivariate quadratic systems the evaluation can be performed in an extremely efficient way using Gray codes [BCC+10] and the exponential cost is just given by the enumeration of all possible vectors: $O*(2^n/n_{\text{sol}})$, where $n_{\text{sol}}$ is the number of solutions.

We do not focus on the matter of adapting the efficient solution space enumeration to the regular setting, as the brute force approach only serves here as a benchmark for other algorithms. Admitting an efficient adaptation of the enumeration, the cost of such an approach in the regular framework amounts to $O*(l^w/n_{\text{sol}})$. In the case of a “randomly” sampled system with $m = \mu n$, $\mu = \frac{\log(l)}{l}$, equations, the complexity becomes

$$O*(2^{\mu n}) = O^*\left(2^{\frac{\log(l)}{l}n}\right).$$

More in general, the brute force approach over the field $\mathbb{F}_q$ has complexity

$$O^*\left(2^{\frac{\log((q-1)t)}{t}n}\right).$$

# 4 Algebraic Cryptanalysis

Besides brute force, the most common tool to solve a polynomial system and analyze the corresponding complexity is to run the XL algorithm/compute a Gröbner basis for the associated ideal. We first explain how to model the regular constraint with algebraic equations and then study relevant objects like the Hilbert series of the associated ideal and the degree of regularity to derive an estimation of the complexity of solving the system. Since we focus on the case of boolean systems, we also make use of and adapt in several ways the hybrid strategy to improve the algorithm.

## 4.1 The algebraic modeling

In addition to the initial system of equations obtained by equating to 0 the quadratic polynomials in $S_{\text{init}} = \{P_i(x) | i \in [1, m]\}$, we first refine the ideal by adding the $n$ polynomials corresponding to the field equations over $\mathbb{F}_2$

$$S_{\text{FE}} \stackrel{\text{def}}{=} \{x_{i,j}^2 - x_{i,j} | i \in [0, w-1], j \in [1, l]\}.$$

which force solutions found through the Gröbner basis computation to lie over $\mathbb{F}_2$ instead of over its closure.

The weight constraint of each window is modeled by two families of equations:
* The following $w(\binom{l}{2})$ homogeneous quadratic polynomials express the fact that each window has at most a nonzero entry

$$S_Q \stackrel{\text{def}}{=} \{x_{i,j_1} x_{i,j_2} \mid 0 \le i \le w-1, 1 \le j_1 < j_2 \le l\}.$$

We have seen that the system $S_{\text{init}}$ could be taken in a way that the quadratic monomials of the form $x_{i,j_1}x_{i,j_2}$ already miss from all the equations. However, the constraints from $S_Q$ are still necessary to remove monomials of higher degree divisible by such $x_{i,j_1}x_{i,j_2}$'s during the Gröbner basis computation.

* The following $w$ affine linear polynomials express the fact that the sum (over $\mathbb{F}_2$) of the entry in each window is 1 (i.e. the weight of each window is odd)

$$S_L \stackrel{\text{def}}{=} \left\{ \sum_{j=1}^{l} x_{i,j} + 1 \mid i \in [0, w-1] \right\}.$$

Note that each linear equation enables the elimination of one variable, so in practice we can set up a system with only $n-w = (l-1)w$ variables and recover the values of the $w$ eliminated variables by substitution once the system has been solved.

Since the only positive odd number smaller than 2 is 1, the two systems $S_Q$ and $S_L$ together model exactly the regular constraint of the solution. The equations from $S_Q$ alone instead would model the “at most regular” constraint. One can see that the equations in $S_Q$ and $S_L$ are exactly the same as those analyzed in [BØ23] for the regular syndrome decoding problem.

The system that models the RMQ problem over $\mathbb{F}_2$ is then the following.

**Modeling A (Over $\mathbb{F}_2$).** Given a system $S_{\text{init}} = \{P_i(\boldsymbol{x}) \mid i \in [1, m]\}$ of quadratic equations over $\mathbb{F}_2$, the algebraic system that models the RMQ problem over $\mathbb{F}_2$ is given by

* the system of $m+n+w(\binom{l}{2})$ quadratic and $w$ linear equations corresponding to the polynomials

$$S \stackrel{\text{def}}{=} S_{\text{init}} \cup S_{\text{FE}} \cup S_Q \cup S_L$$

* in the $n = lw$ variables $x_{i,j}$'s.

More in general, we define an algebraic modeling for big fields, where the field equations are not used and the nonzero entries are not necessarily equal to 1.

**Modeling B (Over a big field $\mathbb{F}_q$).** Given a system $S_{\text{init}} = \{P_i(\boldsymbol{x}) = 0 \mid i \in [1, m]\}$ of quadratic equations over $\mathbb{F}_q$, the algebraic system that models the RMQ problem over $\mathbb{F}_q$ is given by

* the system of $m+w(\binom{l}{2})$ quadratic equations corresponding to the polynomials

$$S \stackrel{\text{def}}{=} S_{\text{init}} \cup S_Q$$

* in the $n = lw$ variables $x_{i,j}$'s.

These are the modelings for the plain Gröbner basis computation, we will see afterward several choices for the hybrid approach. While the former can be analyzed for big fields, the latter are relevant for small fields only and we will then restrict the hybrid methods analysis to boolean functions.

## 4.2 The general proof framework for the Hilbert series derivation

For each plain or hybrid strategy described in this section, we provide the Hilbert series of the associated (possibly specialized) ideal. Instead of providing the proof in each case separately, we give here a general framework of the arguments that lead to these results and devote the space of each subsection to the computation of the asymptotic degree of regularity. The reasons behind this choice are twofold:

1. At a macro level, all these proofs follow the same pattern with just a slight adaptation of the initial assumption.

2. The full proof details have already been given in [BØ23]. In particular, the structured part of the system is exactly the same as for the regular syndrome decoding problem and we just recall the results holding for the structured system. The impact of the instance-specific quadratic equations can then be analyzed analogously to parity-check equations for the regular syndrome decoding problem according to a more general framework.

We then proceed by briefly recalling each step of this analysis.
### 4.2.1 Structured part

The structured (and fixed) part of the system is $S_{fix} = S_Q$ for big fields and $S_{fix} = S_Q \cup S_{FE} \cup S_L$ in the boolean case. The Hilbert series of the high-degree part of the structured subsystem alone can be derived as follows.

**Lemma 1** (from Lemma 6, [BØ23]). Let $R = \mathbb{F}_2[x]$. Let $S_L^{(h)} = \{\sum_{j=1}^l x_{i,j} | i \in [0, w-1]\}$ and $S_{FE}^{(h)} = \{x_{i,j}^2 | i \in [0, w-1], j \in [1, l]\}$. be the part of highest degree of $S_L$ and $S_{FE}$ respectively. Then

$$HS_{R/\langle S_Q \cup S_L^{(h)} \cup S_{FE}^{(h)} \rangle}(z) = (1 + (l-1)z)^w.$$

**Lemma 2** (Lemma 4, [BØ23]). Let $R = \mathbb{F}[x]$, with $\mathbb{F}$ a finite field. Then

$$HS_{R/\langle S_Q \rangle}(z) = \left(1 + l \frac{z}{1-z}\right)^w = \left(\frac{1 + (l-1)z}{1-z}\right)^w.$$

Despite not rewriting the proofs, we just recall that the series is obtained by first determining the Hilbert series “for one block” as $1 + (l-1)z$ or $1 + l \frac{z}{1-z}$ and then applying a standard result about how the Hilbert series (more generally the generating function) compose with Cartesian product [SF09]. Essentially, as the fixed part does not introduce any constraint between variables of different windows, the whole Hilbert series is obtained by multiplying the series of each block. This argument allows to determine the Hilbert series of a specialized system, where some of the variables have been fixed to 0. The most general specialization that we analyze is when any window potentially has a different number of entries set to 0. In particular, we fix to $0 \le l-l' \le 1$ variables (wlog $x_{i,1} = \dots = x_{i,l-l'} = 0$) for each window of a relative portion of $\gamma_{l'} \ge 0$ windows, for any $l' \in [1, l]$ in such a way that $\sum_{l'=1}^l \gamma_{l'} = 1$ and we call $S_{guess}$ the set of variables, which can be seen as homogeneous linear polynomials, that are fixed to 0. The same arguments used in [BØ23] and just recalled here readily imply, in the boolean case,

$$HS_{R/\langle S_Q \cup S_L^{(h)} \cup S_{FE}^{(h)} \cup S_{guess} \rangle}(z) = \prod_{l'=1}^{l} (1 + (l' - 1)z)^{\gamma_{l'} w}.$$

Adding the specialized equations to the structured part allows to unify all cases in a single assumption.

### 4.2.2 Assumptions mimicking semi-regularity

The computation of the Hilbert series for the entire system, i.e. including the non-structured part $S_{init}$, is based on the assumption that acts as that of semi-regularity, but with respect to the quotient ideal of the fixed part.

**Assumption 1.** Consider an instance $S_{init} = \{1, \dots, f_m\}$ of Modeling A or Modeling B and let $d_{reg}$ be the degree of regularity of $\mathcal{I} = \langle S_{init} \rangle$. Define the quotient ring $Q \stackrel{def}{=} R/\langle S_Q \cup S_L^{(h)} \cup S_{FE}^{(h)} \rangle$ in the non-hybrid boolean case, $Q \stackrel{def}{=} R/\langle S_Q \cup S_L^{(h)} \cup S_{FE}^{(h)} \cup S_{guess} \rangle$ in the hybrid boolean case or $Q \stackrel{def}{=} R/\langle S_Q \rangle$ in the (non-hybrid) big field case. We assume that for $i \le i \le m$, $g_i f_i^{(h)} = 0$ in $Q/\langle f_1^{(h)}, \dots, f_{i-1}^{(h)} \rangle$ with $\deg(g_i f_i^{(h)}) < d_{reg}$ implies $g_i = 0$ in $Q/\langle f_1^{(h)}, \dots, f_{i-1}^{(h)} \rangle$.

We let the reader appreciate the strong similarity of this assumption with Definitions 3 and 4. We recall that analogous hypotheses has been used in [BØ23] for the regular syndrome decoding problem, by just replacing the quadratic equations defining the instance with the parity-check equations. Moreover, the assumption, formalizing the fact that randomly generated quadratic polynomials usually behave nicely in the quotient ring with respect to the structural part, has been tested extensively for several parameters.

### 4.2.3 The random quadratic equations on top of the structured part

Retrieving the Hilbert series of $R/\langle S^{(h)} \rangle$ from the quotient $Q$ of Assumption 1 is simply a particular case of [Bar04, Sections 3.3.1, 3.3.2], cf. also [BØ23]. In particular, for each quadratic polynomial in $S_{init}$, the Hilbert series is multiplied by $\frac{1}{1+z^2}$ in the boolean case and by $1-z^2$ in the big field case. Moreover, the system is assumed to behave “semi-regularly” with respect to the extra variable $y$ used to homogenize it (cf. [BØ23, Section 4.3]), which provides a factor $\frac{1}{1-z}$. Finally, the series obtained has to be truncated at the first non-positive coefficient.

All the Hilbert series proposed in the next subsections can thus be determined in a straightforward way within this framework.
## 4.3 Plain Gröbner basis analysis

Let $x' \stackrel{\text{def}}{=} (x_1, \dots, x_n, y)$, $R \stackrel{\text{def}}{=} \mathbb{F}_2[x']$ and $\mathcal{I}$ be the polynomial ideal generated by $S^{(y)}$. The Hilbert series of $R/\mathcal{I}$ is

$$HS_{R/\mathcal{I}}(z) = \left[ \frac{(1 + (l-1)z)^w}{(1-z)(1+z^2)^m} \right]_+,$$

where the term $(1+z^2)^m$ accounts for the initial $m$ quadratic equations and the term $1-z$ for the homogenization variable $y$.

We make use of complex analysis techniques, namely the saddle-point method, to estimate asymptotically the degree of regularity of the system from its Hilbert series. We refer the reader to [Won01] for a broader general discussion of these techniques. These methods have already been widely used in the context of studying the Hilbert series. Two examples that are close to our analysis are [Bar04] for classic semi-regular systems and [BØ23] for the regular syndrome decoding modeling. Here we provide directly the formula for computing the asymptotic degree of regularity (and consequently the complexity) for the plain Gröbner basis approach over $\mathbb{F}_2$. All the necessary computation and explanations to get this value are moved to Appendix A.

The asymptotic relative degree of regularity is the smallest positive root $\bar{\delta}$ of the quartic polynomial

$$r_4\delta^4 + r_3\delta^3 + r_2\delta^2 + r_1\delta + r_0$$

in the variable $\delta$, whose coefficients are polynomial functions of $\mu$ and $l$:

$$\left\{
\begin{aligned}
r_4 &= ((l-1)^2 + 1)^2 \\
r_3 &= 2\mu((l-1)^2 + 1)((l-1)^2 + 3) - 4\frac{(l-1)^2((l-1)^2+1)}{l} \\
r_2 &= 4\mu^2(2(l-1)^2 + 3) - 2\mu\frac{(l-1)^2(3(l-1)^2-1)}{l^2} + 2\frac{(l-1)^2(3(l-1)^2+1)}{l^2} \\
r_1 &= 8\mu^3 + 20\mu^2\frac{(l-1)^2}{l} + 2\mu(l-1)^2\frac{3(l-1)^2-5}{l^2} - 4\frac{(l-1)^4}{l^3} \\
r_0 &= -\mu^2\frac{(l-1)^2}{l^2} - 2\mu\frac{(l-1)^4}{l^3} + \frac{(l-1)^4}{l^4}.
\end{aligned}
\right.$$

Hence, the final complexity becomes

$$O^*\left(\binom{n}{\bar{\delta}n}^\omega\right) = O^*\left(2^{\omega H(\bar{\delta})n}\right).$$

A closed-form expression for roots of quartic polynomials exists but, as the coefficients are complicated functions in $l$ (and possibly $\mu$), we do not provide the exact formula. Here and after, we instead proceed by fixing $l$ and approximate the zeros of the polynomials to find $\bar{\delta}$.

*Remark 3.* To be more precise, the cost of reducing the Macaulay matrix would be $\binom{n}{\bar{\delta}n}^\omega$. However, as it happens for the standard multivariate quadratic problem, $\bar{\delta}$ is small enough to have $O^*\left(\binom{n}{\bar{\delta}n}^\omega\right) = O^*\left(\binom{n}{\bar{\delta}n}^\omega\right)$. The same approximation holds a fortiori for the next hybrid approaches, as $\bar{\delta}$ generally decreases with respect to plain Gröbner basis.

*Remark 4.* Due to the regular constraint and the filed equations, monomials whose degree with restricted to any block is at least 2 are reduced. Hence the total degree of a polynomial is at most $w$, which implies an upper bound on the relative degree regularity $\bar{\delta} \le w/n = 1/l$ and consequently on the final complexity $O^*\left(\binom{n}{\bar{\delta}n}^\omega\right) = O^*\left(2^{\omega H(1/l)n}\right)$. It can be verified that for any $l$, the found degree of regularity correctly satisfies the upper bound, and it tends to it for $l \to \infty$.

### 4.3.1 Big fields

For completeness, we conclude the section with an analysis of Modeling B for a big field $\mathbb{F}_q$. If $\mathcal{I}$ is the ideal generated by the polynomials in $S^{(y)}$, then the Hilbert series is

$$HS_{R/\mathcal{I}}(z) = \left[ \frac{(1 + (l-1)z)^w(1-z^2)^m}{(1-z)^{w+1}} \right]_+.$$

As we did for the boolean case, we detail the computation needed to derive the asymptotic relative degree of regularity in Appendix B. It turns out that this is the smallest positive root $\bar{\delta}$ of the quartic polynomial

$$r_4\delta^4 + r_3\delta^3 + r_2\delta^2 + r_1\delta + r_0$$
in the variable $\delta$, whose coefficients are polynomial functions of $\mu$ and $l$:

$$ 
\begin{cases}
r_4 = 4l^2(2l+1)^2 \\
r_3 = -8\mu(l^3 - 4l^2 - 2l + 4) - 4(3l^2 - 8l + 8)(l-2) \\
r_2 = -16\mu^2(2l^2 - 4l - 1) + 8\mu(3l^3 - 14l^2 + 29l - 24) + 13l^2 - 48l + 48 \\
r_1 = -32\mu^3 - 16\mu^2(5l - 8) - 4\mu(6l^2 - 23 + 24) - 6(l-2) \\
r_0 = 4\mu^2 + 4\mu(2l-3) + 1.
\end{cases}
$$

Hence, the final complexity becomes

$$ \mathcal{O}^* \left( \left( \frac{n + \bar{\delta}n}{\bar{\delta}n} \right)^{\omega} \right) = \mathcal{O}^* (2^{\tau n}), $$

with $\tau = \omega(1 + \bar{\delta})H(\bar{\delta}/(1 + \bar{\delta}))$.

For a fixed $q$, the plain Gröbner basis computation is convenient with respect to brute force for small enough values of $l$ and whose maximum (slowly) increases with $q$. To give some examples,

* for $q \le 5$, brute force is always better;

* for $q \in \{7, 8\}$ plain Gröbner basis beats brute force iff $l = 2$;

* for $q \in \{31, 32\}$ plain Gröbner basis beats brute force iff $l \le 8$;

* for $q \in \{255, 256\}$ plain Gröbner basis beats brute force iff $l \le 31$.

## 4.4 Hybrid full-windows guess

A first idea of hybridization over $\mathbb{F}_2$ is to guess the nonzero entry of an amount of $\gamma w$ windows for some $\gamma \in [0, 1]$. The choice $\gamma = 0$ coincides with the plain Gröbner basis approach, while $\gamma = 1$ with exhaustive search. The probability that the guess is consistent with the unique solution is $1/l^{\gamma w}$. The factor that takes into account the cost of guessing is $l^{\gamma w}$.

For a fixed specialization and after partial evaluation, this method reduces the initial RMQ problem to another RMQ problem with $(1 - \gamma)w$ windows of length $l$ and hence with $(1 - \gamma)n$ variables. The corresponding Hilbert series becomes

$$ \mathrm{HS}_{R/\mathcal{I}}(z) = \left[ \frac{(1 + (l-1)z)^{(1-\gamma)w}}{(1-z)(1+z^2)^m} \right]_+ . $$

By repeating the computation from the plain Gröbner basis strategy shown in Appendix A, but replacing $w$ with $\gamma w$, we obtain the quartic equation $r_4\delta^4 + r_3\delta^3 + r_2\delta^2 + r_1\delta + r_0 = 0$ with

$$ 
\begin{cases}
r_4 = ((l-1)^2 + 1)^2 \\
r_3 = 2\mu((l-1)^2 + 1)((l-1)^2 + 3) - 4(1-\gamma)\frac{(l-1)^2((l-1)^2+1)}{l} \\
r_2 = 4\mu^2(2(l-1)^2 + 3) - 2\mu(1-\gamma)\frac{(l-1)^2(3(l-1)^2-1)}{l} + 2(1-\gamma)^2\frac{(l-1)^2(3(l-1)^2+1)}{l^2} \\
r_1 = 8\mu^3 + 20\mu^2(1-\gamma)\frac{(l-1)^2}{l} + 2\mu(1-\gamma)\frac{(l-1)^2(3(l-1)^2-5)}{l^2} - 4(1-\gamma)^3\frac{(l-1)^4}{l^3} \\
r_0 = -\mu^2(1-\gamma)^2\frac{(l-1)^2}{l^2} - 2\mu(1-\gamma)^3\frac{(l-1)^4}{l^3} + (1-\gamma)^4\frac{(l-1)^4}{l^4}.
\end{cases}
$$

If $\bar{\delta}$ is the smallest positive root of this quartic equation, then the final complexity of the hybrid full-windows method is

$$ \mathcal{O}^* \left( l^{\gamma w} \left( \frac{(1 - \gamma)n}{\bar{\delta}n} \right)^{\omega} \right) = \mathcal{O}^* (2^{\tau n}), $$

with $\tau = \gamma \frac{\log(l)}{l} + \omega(1 - \gamma)H(\frac{\bar{\delta}}{1-\gamma})$.

## 4.5 Hybrid partial-windows guess

Guessing the nonzero entry in a window can be equivalently seen as guessing the $l-1$ zero entries of the same window. Another possible way to hybridize over $\mathbb{F}_2$ is to guess only the same amount of zeros in each window, so that only $l'$ variables per block remain. This method contains the plain Gröbner basis one and the exhaustive search as subcases for $l' = l$ and $l' = 1$ respectively. The probability that the guess is consistent with the unique solution is $(l/l')^w$. Hence, the factor that takes into account the cost of guessing is $(l'/l)^w$. Differently from the real number $\gamma$ from the previous hybridization, the value $l'$ only assumes discrete values.
For a fixed specialization and after partial evaluation, this method reduces the initial RMQ problem to another RMQ problem with w windows of length l' and hence with $l'w = \frac{l'}{l}n$ variables. The corresponding Hilbert series becomes

$$HS_{R/I}(z) = \left[ \frac{(1 + (l' - 1)z)^w}{(1 - z)(1 + z^2)^m} \right]_+ .$$

By repeating the computation from the plain Gröbner basis strategy shown in Appendix A, but replacing l with l' (and paying attention to the fact that w remains 1/l and not 1/l'), we obtain the quartic equation $r_4\delta^4 + r_3\delta^3 + r_2\delta^2 + r_1\delta + r_0 = 0$ with

$$\left\{
\begin{aligned}
r_4 &= ((l' - 1)^2 + 1)^2 \\
r_3 &= 2\mu((l' - 1)^2 + 1)((l' - 1)^2 + 3) - 4(1 - \gamma)\frac{(l' - 1)^2((l' - 1)^2 + 1)}{l} \\
r_2 &= 4\mu^2(2(l' - 1)^2 + 3) - 2\mu(1 - \gamma)\frac{(l' - 1)^2(3(l' - 1)^2 - 1)}{l} + 2(1 - \gamma)^2\frac{(l' - 1)^2(3(l' - 1)^2 + 1)}{l^2} \\
r_1 &= 8\mu^3 + 20\mu^2(1 - \gamma)\frac{(l' - 1)^2}{l} + 2\mu(1 - \gamma)\frac{(l' - 1)^2(3(l' - 1)^2 - 5)}{l^2} - 4(1 - \gamma)^3\frac{(l' - 1)^4}{l^3} \\
r_0 &= -\mu^2(1 - \gamma)^2\frac{(l' - 1)^2}{l^2} - 2\mu(1 - \gamma)^3\frac{(l' - 1)^4}{l^3} + (1 - \gamma)^4\frac{(l' - 1)^4}{l^4}.
\end{aligned}
\right.$$

If $\bar{\delta}$ is the smallest positive root of this quartic equation, then the final complexity of the hybrid partial-windows method is

$$O^*\left(\left(\frac{l}{l'}\right)^w \binom{(l'/l)n}{\bar{\delta}n}\right)^\omega = O^*(2^{\tau n}),$$

$$\text{with } \tau = \frac{\log(l/l')}{l} + \omega \frac{l'}{l} H(\frac{\bar{\delta}l}{l}).$$

Remark 5. For all $4 \le l < 40$, our estimates show that the choice $l' = 2$ below provides the best trade-off. For $l \le \text{split} = 1000$, the optimal $l'$ never exceeds 3. More in general, the optimal $l'$ grows very slowly with respect to l, meaning the hybrid approach is quite unbalanced towards brute force rather than the plain Gröbner basis computation.

## 4.6 Hybrid different-windows guess

The most general hybrid approach is to allow different partial guessing of zeros for each window, i.e. guessing $l - l'$ variables for each window of a relative portion of $\gamma_{l'} \ge 0$ windows, for any $l' \in [1, l]$ in such a way that $\sum_{l'=1}^l \gamma_{l'} = 1$. The probability that the guess is consistent with the unique solution is $\prod_{l' \in [1, l]} (l/l')^{\gamma_{l'}w} = \prod_{l' \in [1, l-1]} (l/l')^{\gamma_{l'}w}$. Hence, the factor that takes into account the cost of guessing is $\prod_{l' \in [1, l-1]} (l'/l)^{\gamma_{l'}w}$. We remark that this hybridization can be seen as a natural generalization of all previous strategies:

* $\gamma_l = 1$ and $\gamma_{l'} = 0$ otherwise, corresponds to a plain Gröbner basis computation;

* $\gamma_1 = 1$ and $\gamma_{l'} = 0$ otherwise, corresponds to brute force;

* the minimum over the subset of tuples such that $\gamma_1 = \gamma$, $\gamma_l = 1 - \gamma$ and $\gamma_{l'} = 0$ otherwise, corresponds to hybrid full-window guess;

* the minimum over the subset of tuples such that $\gamma_{l'} = 1$ for some $\bar{l}'$ and $\gamma_{l'} = 0$ otherwise, corresponds to hybrid partial-window guess.

We also observe that, differently from the previous subcases, this additional layer of hybridization has not been considered in [BØ23] in the context of regular syndrome decoding. Moreover, a generic specialization of this method does not reduce to another RMQ problem as windows can have different lengths.

To help the reader navigate through this general hybrid approach, we start by providing the case studies for different values of l, postponing the computation necessary to obtain the tuple $(\gamma_1, \dots, \gamma_l)$ that minimizes the complexity.

### Full case study for $l = 2$

The regular multivariate quadratic problem with $l = 2$ directly reduces to an instance of the standard multivariate quadratic problem. Indeed the linear equation $x-1+x_2 = 1$ from $S_L$ (and corresponding to the first window) allows to eliminate the variable $x_2$ from the regular constraint $x_1x_2 = 0$ from $S_Q$, obtaining $x_1(x_1+1) = 0$. But the latter equation is exactly the field equation for $x_1$, so the regular constraint is redundant. The same happens for any other window. Moreover, the complexity is the same as that of the “hybrid partial-window guess” approach, as one can only choose to guess either zero or one variable in a window. Since we are left with half the variables, the coefficient at the exponent gets exactly halved with respect to [Bar04], leading to a complexity

$$O^*(2^{\frac{0.791}{2}n}) = O^*(2^{0.3955n}),$$

corresponding to the vector $(\gamma_1, \gamma_2) = (0.449, 0.551)$ and $\delta_{reg} = 0.0154$, for $\mu = \frac{\log(l)}{l} = \frac{1}{2}$.
## Full case study for $l = 3$

The case $l = 3$ is therefore the first non-trivial example. The best choice of $\gamma$ with the different-windows hybrid approach is $(\gamma_1, \gamma_2, \gamma_3) \approx (0.127, 0.873, 0)$, leading to $\delta_{reg} \approx 0.01622$ and a complexity of

$$O^*(2^{0.4179n}).$$

The fact that $\gamma_3 = 0$ means that at least 1 variable is thus guessed in each window and this results in finding the best trade-off for a random quadratic system in $m' = m$ unknowns and $n' = \frac{n}{3}$ variables. The classic hybrid approach applied for

$$\mu' = \frac{m'}{n'} = \frac{3m}{n} = 3\frac{\log(l)}{l} = \log(3)$$

is obtained for $\gamma \approx 0.127$, providing $\delta_{reg} \approx 0.0783$ and a complexity of $2^{0.6689n}$.

If we now multiply by the cost of guessing variables, we obtain

$$\left(\frac{l}{l-1}\right)^{n/3} \cdot 2^{0.6689n'}$$

that leads to

$$O^*\left(2^{\left(\frac{\log(3/2)}{3} + \frac{0.6689}{3}\right)n}\right) = O^*(2^{0.4179n})$$

i.e. we retrieve exactly the complexity found with our method, thus corroborating and verifying our result.

## Full case study for $l = 4, 5, 6$

For $l = 4, 5, 6$, the best choice of $\gamma$ with the different-windows hybrid approach is $\gamma_2 = 1$ and $\gamma_{l'} = 0$ otherwise, the RMQ problem is reduced to a random system of $m' = m$ equations in $n' = \frac{n}{l}$ variables. We thus have

$$\mu' = \frac{m'}{n'} = \log(l).$$

According to [BFSS13], whenever $0.55\mu' \ge 1$, the best algorithm is known to be a plain Gröbner basis computation (no hybrid). This is satisfied by any $l \ge 4$ and explains why, assuming $\gamma_{l'}$ for all $l' > 2$, then $(\gamma_1, \gamma_2) = (0, 1)$. We remark that, in these cases, the hybrid different-windows method boils down to the partial-windows one and it is very close to the brute force, as in each window only one variable less than the maximum is guessed.

## Large $l$

For large (constant) values of $l$, estimating the best parameters for the hybrid different-windows approach becomes very costly. As already said in Remark 5, if we restrict to the hybrid partial-window guess, the optimal choice seems to remain $\gamma_2 = 1$ and $\gamma_{l'} = 0$ otherwise. We expect the result to change very little or nothing when allowing for different amounts of guesses in each window.

## Hilbert series and general complexity

The Hilbert series corresponding to the specialized systems is

$$HS_{R/\mathcal{I}}(z) = \left[ \frac{\prod_{l'=1}^{l} (1 + (l' - 1)z)^{\gamma_{l'w}}}{(1-z)(1+z^2)^m} \right]_+ = \left[ \frac{\prod_{l'=2}^{l} (1 + (l' - 1)z)^{\gamma_{l'w}}}{(1-z)(1+z^2)^m} \right]_+ .$$

Consider the Cauchy integral

$$\mathcal{I}_d(n) \stackrel{\text{def}}{=} \frac{1}{2\pi i} \oint \frac{1}{z^{d+1}} \frac{\prod_{l'=2}^{l} (1 + (l' - 1)z)^{\gamma_{l'w}}}{(1-z)(1+z^2)^m} dz$$

and set $f(z)$ such that $e^{nf(z)} = \frac{1}{z^{d+1}} \frac{\prod_{l'=2}^{l} (1 + (l' - 1)z)^{\gamma_{l'w}}}{(1+z^2)^m}$ (as the term $(1-z)$ is asymptotically negligible), i.e.

$$f(z) \stackrel{\text{def}}{=} \left( \frac{1}{l} \sum_{l'=2}^{l} \gamma_{l'} \ln(1 + (l' - 1)z) \right) - \mu \ln(1 + z^2) - \delta \ln(z).$$

We compute the saddle points from the function

$$zf'(z) = \frac{1}{l} \left( \sum_{l'=2}^{l} \gamma_{l'} \frac{(l' - 1)z}{1 + (l' - 1)z} \right) - \frac{2\mu z^2}{1 + z^2} - \delta.$$
For each tuple $(\gamma_{l'})_{l'}$ and a fixed window-length $l$, we multiply $zf'(z)$ by $1+z^2$ and all the denominators $1+(l'-1)z$'s such that $\gamma_{l'} \neq 0$, in order to obtain a polynomial

$$g = \sum_{i=0}^{d_g} g_i z^i = \sum_{i=0}^{d_g} g_i(\delta) z^i$$

in the variable $z$ and parameter $\delta$ of degree $d_g = 1 + |\{l' | 2 \le l' \le l, \gamma'_{l'} \ne 0\}|$, from which we can obtain the saddle points. These are given by the roots of the discriminant of $g$, whose degree is $d_g - 2$. As long as at least two $\gamma_{l'}$ are nonzero, such a degree becomes $\ge 6$ and closed formula expressions for its root can not even be computed in theory. We recall that roots of the discriminant of $g$ are exactly the common roots of $g$ and its derivative $g'$, i.e. those of the resultant $\mathrm{Res}(g, g')$. Such resultant is a polynomial in $\delta$ of degree $d_g + (d_g - 1) = 2d_g - 1$. Hence, the asymptotic value of the degree of regularity can be equivalently taken as the smallest positive root of the resultant:

$$\bar{\delta} = \min\{\delta \in \mathbb{R} \mid (\mathrm{Res}(g, g'))(\delta) = 0, \ \delta > 0\}.$$

As the degree of the resultant potentially grows with $l$ and deriving a closed formula of its roots is not possible in general, we proceed by substituting all possible tuples $(\gamma_{l'})_{l'}$ according to some discrete split of the number of windows split. Thus, a priori, we need to compute the complexity for

$$\binom{\text{split} + l - 1}{l - 1}$$

tuple values, which quickly increases with $l$. We remark that this discrete splitting is also useful for concrete parameters where the number does not tend to infinity and therefore taking split = $w$ provides the best valid tuple $(\gamma_{l'})_{l'}$.

If $\bar{\delta}$ is the smallest positive root of the resultant, then the final complexity of the hybrid different-windows method is

$$O^*\left(\left(\prod_{l' \in [1, l-1]} \left(\frac{l}{l'}\right)^{\gamma_{l'} w}\right) \left(\frac{\sum_{l' \in [2, l]} \gamma_{l'} \frac{l'-1}{l} n}{\bar{\delta} n}\right)^{\omega}\right) = O^*(2^{\tau n}),$$

$$\text{with } \tau = \sum_{l' \in [1, l-1]} \gamma_{l'} \frac{\log(l/l')}{l} + \omega \sum_{l' \in [2, l]} \gamma_{l'} \frac{l'-1}{l} H\left(\frac{\bar{\delta} \cdot l}{\sum_{l' \in [2, l]} \gamma_{l'} (l'-1)}\right).$$

# 5 Cryptanalysis: the probabilistic polynomial method

We now adapt the polynomial method introduced in [LPT$^+$17] and then improved by [BKW19, Din21b, Din21a] to the regular setting.

We briefly and informally summarize the contributions brought by each article of this research line. Then, we provide a more detailed overview of their costs in standard case which serves as a comparison term for the subsequent improved adaptations to the regular setting. This family of algorithms, which are substantially different from algebraic approaches, relies on the strategy of representing a polynomial system through a single *identifying polynomial*. The mathematical aspect of this idea is credited to Razborov and Smolensky [Raz87, Smo87]. As we will recall, this approach is however impractical due to the very large degree attained by such polynomial. The probabilistic polynomial method consists in replacing the *identifying polynomial* with another one of much lower degree that represents the former with high probability.

This technique was previously developed for circuit design and its complexity analysis [Wil14], but used for the first time for polynomial system solving in [LPT$^+$17]. The algorithm of Lokshtanov et al. is the first one that achieves an exponential speedup over exhaustive search in the worst case setting (for instance without semiregular or other algebraic assumptions and without requiring the system to be over- or underdetermined). The algorithm reduces the search version of the problem of solving a system to a linear number of calls to its decision counterpart and makes use of the probabilistic construction to solve the latter faster.

The first improvement on this first milestone is due to Björklund et al., leading to a further exponential speedup. The core idea of [BKW19] is to isolate solutions in such a way that the decision version of the system solving problem can be in turn reduced to finding the parity of the number of solutions. We remark that in many cryptographic applications, algebraic systems are chosen in a way that only one solution exists. Our description thus immediately starts from Björklund et al.'s point of view. The other main contribution of [BKW19] is to observe that a certain evaluation steps can be accelerated by using a self-recursion to smaller parity-counting instances.

Building upon this version of the algorithm, Dinur remarked that the partial parities that need to be computed are strongly correlated with each other. Therefore, in [Din21b] the multiple parity-counting problem
is introduced, together with a procedure to solve it. By also making use of self-recursion, the algorithm form [Din21b] attains the best asymptotic complexity for solving binary polynomial systems among all known methods. Dinur later observed that the aforementioned reductions between problem add a large polynomial overhead that can be avoided completely. The same happens for the self-recursion, and [Din21a] hence proposes a variant of the algorithm that has a worse exponential complexity but is better for concrete parameters of cryptographic interest because it cut out a significant polynomial factor. Differently from the previous approaches, [Din21a] assumes that, for a given partition, an isolated solution exists with high probability. For cryptographic instances, this constraint is typically satisfied. We refer the reader to [BBSV22] for a survey about this line of research and for an implementation and practical comparison of the mentioned methods. In any case, as for the Gröbner basis approaches, our focus is mainly on the asymptotically relevant part and we provide the final complexity neglecting polynomial factors in n (i.e. using the $O^*$ notation).

## 5.1 Preliminaries

Given the initial system $S = \{P_1, \dots, P_m\}$, we define the polynomial

$$F(\boldsymbol{x}) = \prod_{i=1}^{m} (1 + P_i(\boldsymbol{x})).$$

We have that $\hat{\boldsymbol{x}}$ is a solution of the system, if and only if $F(\hat{\boldsymbol{x}}) = 1$. In order to simplify the discussion, we also assume that there is only one solution to the system. Indeed, the same assumption (which is common in cryptography) is made when adapting the method to the regular framework, as the value m is chosen in a way that the expected number of solutions is 1. The polynomial F is thus called identifying polynomial of S. However, the polynomial F has a too high degree and the cost of evaluations becomes forbidding. Therefore, F is replaced by a probabilistic polynomial $\tilde{F}$ constructed in the following way. Let $k < m$ be a positive integer. Pick uniformly at random a full rank matrix $\mathbf{A} = (a_{i,j}) \in \mathbb{F}_2^{k \times m}$ and define for any $i \in [1, k]$ the quadratic polynomial

$$R_i(\boldsymbol{x}) = \sum_{j=1}^{m} a_{i,j} P_j(\boldsymbol{x}).$$

Finally, we define the probabilistic polynomial

$$\tilde{F}(\boldsymbol{x}) = \prod_{i=1}^{k} (1 + R_i(\boldsymbol{x}))$$

which is the identifying polynomial of $\tilde{S} = \{R_1, \dots, R_k\}$ and whose degree is upper bounded by $2k$.

**Proposition 5.** For any $\hat{\boldsymbol{x}} \in \mathbb{F}_2^n$, if $F(\hat{\boldsymbol{x}}) = 1$, then $\tilde{F}(\hat{\boldsymbol{x}}) = 1$. If instead $F(\hat{\boldsymbol{x}}) = 0$, then $\mathrm{Pr}[\tilde{F}(\hat{\boldsymbol{x}}) = 0] \ge 1-2^{-k}$.

The approach of [BKW19] proceeds by reducing the problem of finding a solution to S to a parity-counting problem. The first step consists of reducing the search version to the decision version of the problem, by fixing one by one the variables and calling $\Theta(n)$ times the decision algorithm.

Then the decision version of the problem is reduced to the parity-counting problem. As the solution is unique, the output of the decision algorithm is positive if and only if the output of the parity-counting algorithm is 1. In the case where the solution is not unique, [BKW19] makes use of the Valiant-Vazirani affine hashing [VV85], which consists of adding random linear equations to the system in order to isolate a solution. In [BKW19], the parity of the number of solutions to S is computed in parts by partitioning the length-n vector of variables $\boldsymbol{x}$ into 2 vectors $\boldsymbol{y} = (y_1, \dots, y_{n-n_1})$ and $\boldsymbol{z} = (z_1, \dots, z_{n_1})$, for some parameter $\gamma$ such that $n_1 = \gamma n$, i.e.

$$\sum_{\hat{\boldsymbol{x}} \in \mathbb{F}_2^n} F(\hat{\boldsymbol{x}}) = \sum_{\hat{\boldsymbol{y}} \in \mathbb{F}_2^{n-n_1}} \sum_{\hat{\boldsymbol{z}} \in \mathbb{F}_2^{n_1}} F(\hat{\boldsymbol{y}}, \hat{\boldsymbol{z}}).$$

When replacing F with $\tilde{F}$, Proposition 5 and the union bound over all $\hat{\boldsymbol{z}} \in \mathbb{F}_2^{n_1}$ provide the following upper bound on the partial parity-sums for each $\hat{\boldsymbol{y}} \in \mathbb{F}_2^{n-n_1}$:

$$\mathrm{Pr} \left[ \sum_{\hat{\boldsymbol{z}} \in \mathbb{F}_2^{n_1}} \tilde{F}(\hat{\boldsymbol{y}}, \hat{\boldsymbol{z}}) = \sum_{\hat{\boldsymbol{z}} \in \mathbb{F}_2^{n_1}} F(\hat{\boldsymbol{y}}, \hat{\boldsymbol{z}}) \right] \ge 1 - 2^{n_1 - k}. \quad (4)$$

The integer k is thus chosen in an optimal way to obtain a probability non-negligibly larger than 1/2, that is $k = n_1 + 2$, i.e. $\lim_{n \to \infty} \frac{k}{n} = \gamma$.
A boolean function $F(y, z)$ can be written as

$$F(y, z) = z_1 \cdots z_{n_1} F_1(y) + F_2(y, z),$$

where at least a variable $z_i$ is missing in each monomial of $F_2(y, z)$. Then, because the field of coefficients is $\mathbb{F}_2$, the polynomial

$$G(y) \stackrel{\text{def}}{=} \sum_{\hat{z} \in \mathbb{F}_2^{n_1}} \tilde{F}(y, \hat{z}) \in \mathbb{F}_2[y]$$

is such that $G(y) = F_1(y)$. Therefore, the degree of $G(y)$ is at most

$$d_G \le 2k - n_1.$$

### 5.1.1 Non-recursive version

In order to interpolate $G(y)$ via the Möbius transform, it is then enough to evaluate it on $\binom{[n-n_1]}{\downarrow d_G} \times \mathbb{F}_2^{n_1}$, instead of on $\mathbb{F}_2^{n_1}$. Each evaluation corresponds to finding solutions to $\hat{S}$ via brute force, thus leading to a complexity of

$$O^*\left(\binom{[n-n_1]}{\downarrow d_G} 2^{n_1}\right).$$

The following step consists of evaluating $G(y)$ on all $y \in \mathbb{F}_2^{n-n_1}$ via Möbius transform to obtain the partial parity sums with probability of correctness at least $3/4$. Repeating $t = \Theta(n)$ times the procedure and applying a majority vote to choose the value of the partial parity leads to an error with exponentially small probability. This step trivially requires

$$O^*(2^{n-n_1})$$

operations. Finally the algorithms outputs the parity as $\sum_{\hat{y} \in \mathbb{F}_2^{n-n_1}} G(\hat{y})$.

Ignoring polynomial factors in $n$, the total complexity is

$$O^*\left(\binom{n-n_1}{\downarrow d_G} 2^{n_1} + 2^{n-n_1}\right).$$

Let us define the function $g: [0, 1] \to [0, 1]$

$$g(p) \stackrel{\text{def}}{=} (1-p)H^*\left(\frac{p}{1-p}\right). \qquad (5)$$

For a fixed $\gamma$, the asymptotic complexity can thus be written as

$$O^*(2^{\tau n}),$$

where

$$\tau \stackrel{\text{def}}{=} \max(g(\gamma) + \gamma, 1 - \gamma).$$

We want to minimize $\tau$ over $\gamma \in [0, 1]$. Moreover, we observe that, whenever $\gamma \ge 1/3$, we have $g(\gamma) = 1 - \gamma$, therefore $\tau \ge 1$ does not improve upon brute force. It turns out that the best choice to balance the two costs is $\gamma \approx 0.185 = 1 - 0.815$, which leads to a complexity of $O^*(2^{0.815n})$.

### 5.1.2 Recursion à la Björklund et al.

In [BKW19], a recursive self-reduction of the parity-counting problem is also used to improve the total complexity. Instead of evaluating over $\binom{[n-n_1]}{\downarrow d_G} \times \mathbb{F}_2^{n_1}$ for each element $y \in \binom{[n-n_1]}{\downarrow d_G}$, the recursion is applied on the smaller space $\mathbb{F}_2^{n_1}$. The recursion stops after a constant number of iterations $D$. In the last step, an exhaustive search is run. A parity-counting instance thus reduces to many smaller instances of the same problem, where at each step the number of variables decreases by a factor $\frac{n_1}{n} = \gamma$. If we fix $\gamma$ and denote the complexity coefficient at depth $i$ of the recursion as $\tau_i$, the recursive formula thus gives

$$\tau_i = \gamma^i \max(g(\gamma) + \tau_{i+1}, 1 - \gamma) \quad \forall i \in [0, D-1] \qquad (6)$$

and

$$\tau_D = \gamma^D.$$

We want to determine the value of $\gamma$ that optimizes $\tau_0$. By taking enough large $D$, [BKW19] shows that

$$\tau_0 \le 1 - \gamma \quad \text{whenever} \quad H^*\left(\frac{\gamma}{1-\gamma}\right) < 1 - \gamma,$$
i.e.

$$ \tau_0 \le 1 - \gamma \quad \text{whenever} \quad g(\gamma) < (1 - \gamma)^2. $$

In order to optimize the final complexity, the best choice consists of finding the value for which the increasing function $H^* (\frac{\gamma}{1-\gamma})$ equals the decreasing function $1 - \gamma$. In particular, we can select $\gamma \approx 0.197 = 1 - 0.803$, leading to a complexity of $O^*(2^{0.803n})$.

### 5.1.3 Recursion à la Dinur

In [BKW19], the many parity-counting instances are solved independently. However, this approach is suboptimal as such instances are related. Dinur introduces in [Din21b] the multiple parity-counting problem which solves all parity-counting instances at once. When recursion is applied, the multiple parity-counting instance is reduced to only a few smaller instances of the same problem. Moreover, the multiple counting algorithm is used to enumerate all the solutions of a polynomial system.

The multiple parity-counting instance must output the vector

$$ (G(\mathbf{y}))_{\mathbf{y} \in (\begin{smallmatrix} n-n_1 \\ \downarrow d_G \end{smallmatrix})}. $$

In addition to the parameter $\gamma = n/n_1$, a parameter $\eta$ is introduced in order to further split the space of the last $n_1$ variables into a right part of $n_2 = (\gamma - \eta)n$ variables and a left part of $n_1 - n_2 = \eta n$ variables. In order to compute the multiple parity-counting instance, Dinur's algorithm applies a self-reduction to only one multiple parity-counting instance that outputs the vector whose coordinates are

$$ G(\mathbf{y}, \mathbf{u}) = \sum_{\hat{\mathbf{v}} \in \mathbb{F}_2^{n_2}} F(\mathbf{y}, \mathbf{u}, \hat{\mathbf{v}}), $$

with $(\mathbf{y}, \mathbf{u}) \in \mathbb{F}_2^{n_1-n_2} \times \mathbb{F}_2^{n_1-n_2}$ such that the total weight is upper bounded by a certain quantity that depends on the rank $k'$ of the randomized subsystem that needs to be taken. Using the same argument derived from (4) as before, we can choose $k' = n_2 + 2$, i.e. $\lim_{n \to \infty} \frac{k'}{n} = \gamma - \eta$, so that $\deg G(\mathbf{y}, \mathbf{u}) \le 2k' - n_2 = n_2 + 4 \sim n_2$. The following step is the interpolation of $G(\mathbf{y}, \mathbf{u})$ with the Möbius transform from its low-weight evaluations. The corresponding cost is

$$ \mathcal{O}^* \left( \left( \begin{smallmatrix} n - n_2 \\ \downarrow \deg G(\mathbf{y}, \mathbf{u}) \end{smallmatrix} \right) \right) = \mathcal{O}^* \left( \left( \begin{smallmatrix} n - n_2 \\ \downarrow n_2 \end{smallmatrix} \right) \right). \quad (7) $$

Finally $G(\mathbf{y}, \mathbf{u})$ is evaluated on all the elements of

$$ \left( \begin{smallmatrix} [n - n_1] \\ \downarrow \deg G(\mathbf{y}) \end{smallmatrix} \right) \times \mathbb{F}_2^{n_1 - n_2}, \quad (8) $$

so that the initial multi-parity counting output

$$ (G(\mathbf{y}))_{\mathbf{y} \in (\begin{smallmatrix} n-n_1 \\ \downarrow d_G \end{smallmatrix})} $$

can be reconstructed. The algorithm calls itself until depth $D$. If we denote with $n_1^{(i)}$ and $n_2^{(i)}$ the values of $n_1$ and $n_2$ at step $i \ge 0$ respectively, we obtain $n_1^{(0)} = n_1 = \gamma n$ and $n_2^{(i)} = n_1^{(i+1)} \sim n_1^{(i)} - \eta n \sim n_1^{(0)} - i\eta n \sim (\gamma - i\eta)n$, for any $i \in [0, D-1]$, where $D$ is the integer such that $n_1^{(D)} \le \eta n$. Then we have $n_2^{(D)} = 0$ and the algorithm runs brute force instead.

We want to analyze the cost of each step of the recursion in order to find the time complexity of the whole algorithm. At each step $1 \le i < D$, the cost (7) coming from the interpolation can be written as

$$ \mathcal{O}^* \left( \left( \begin{smallmatrix} n - n_2^{(i)} \\ \downarrow n_2^{(i)} \end{smallmatrix} \right) \right) = \mathcal{O}^* (2^{g(\gamma - i\eta)n}), $$

while the evaluation step over (8) takes time

$$ \mathcal{O}^* \left( \left( \begin{smallmatrix} n - n_1^{(i)} \\ \downarrow n_1^{(i)} \end{smallmatrix} \right) 2^{n_1^{(i)} - n_2^{(i)}} \right) = \mathcal{O}^* (2^{(\eta + g(\gamma - (i-1)\eta))n}). $$

We then have to add the cost of brute force at depth $D$, which is

$$ \mathcal{O}^* \left( \left( \begin{smallmatrix} n - n_1^{(D)} \\ \downarrow n_1^{(D)} \end{smallmatrix} \right) 2^{n_1^{(D)}} \right) \le \mathcal{O}^* (2^{(\eta + g(\gamma - (D-1)\eta))n}). $$
and the cost of the evaluation of $G(\mathbf{y})$ for all elements $\mathbf{y} \in \mathbb{F}_2^{n_1}$, i.e.

$$ \mathcal{O}^*(2^{(1-\gamma)n}). $$

We recall now how to choose the parameters that provide the total complexity claimed in [Din21b]. We fix $\eta > 0$ arbitrarily small, define $\tau = \max_{p \in [0, 1/3]} g(p) = \log \psi \approx 0.6943$, where $\psi \approx 1.618$ is the golden ratio, and choose $\gamma = 1 - \tau$. In this way, each recursive step has an upper bound on the time complexity given by

$$ \mathcal{O}^*(2^{(\tau+\epsilon)n}). $$

Moreover, since $n_2^{(D)} = 0$, we derive that $D\eta \le \gamma \le (D+1)\eta$ and hence $D \sim \frac{\gamma}{\eta} = O(1)$ can be taken constant so that the final complexity remains

$$ \mathcal{O}^*(2^{0.6943n}). $$

## 5.2 Adapting to the regular framework

The search to decision problem reduction can be easily transferred to the regular setting. Assuming a solution exists, for the first window all the $l$ possibilities are tested. Assuming the uniqueness of the solution, only one test gives a positive answer. Then the algorithm proceeds recursively with the next window. Thus, the number of calls needed is upper bounded by $wl = n$.

We define the polynomials $R_i$'s and $F$ in the same way. Proposition 5 still holds. However, since we want to find the unique regular solution, we need to consider the parity of solutions over regular elements, i.e.

$$ \sum_{\hat{x} \in \mathcal{R}(\mathbb{F}_2^n)} F(\hat{x}) = \sum_{\hat{y} \in \mathcal{R}(\mathbb{F}_2^{n-1})} \sum_{\hat{z} \in \mathcal{R}(\mathbb{F}_2^{n_1})} F(\hat{y}, \hat{z}) \text{ and we define} $$

$$ G(\mathbf{y}) = \sum_{\hat{z} \in \mathcal{R}(\mathbb{F}_2^{n_1})} F(\mathbf{y}, \hat{z}). $$

We want $\hat{y} \in \mathbb{F}_2^{n_1}$ to be an isolated solution, and we have

$$ \mathrm{Pr} \left[ \sum_{\hat{z} \in \mathcal{R}(\mathbb{F}_2^{n_1})} \tilde{F}(\hat{y}, \hat{z}) = \sum_{\hat{z} \in \mathcal{R}(\mathbb{F}_2^{n_1})} F(\hat{y}, \hat{z}) \right] \ge 1 - \frac{\lceil \gamma w \rceil}{2^k}. $$

Therefore, we can pick $k = \lceil \gamma w \log l \rceil + 2$ so that the probability above remains above 1/2 by a constant factor. Asymptotically, this choice corresponds to

$$ k \sim \gamma w \log l. $$

Let $\delta_v$ be the set of indices where $v$ is not zero, i.e.

$$ \delta_v \stackrel{\text{def}}{=} \{i \in \llbracket 1, n \rrbracket \mid v_i = 1\}. $$

We can write the polynomial $F$ as

$$ F(\mathbf{y}, \mathbf{z}) = \sum_{v \in \mathcal{R}(\mathbb{F}_2^{n_1})} \left( \prod_{i \in \delta_v} z_i \right) F_v(\mathbf{y}) + F'(\mathbf{y}, \mathbf{z}), $$

where each monomial appearing in $F'(\mathbf{y}, \mathbf{z})$ misses completely at least one window of variables in $\mathbf{z}$. Hence, if $l$ is even, then

$$ G(\mathbf{y}) = \sum_{z \in \mathcal{R}(\mathbb{F}_2^{n_1})} \left[ \sum_{v \in \mathcal{R}(\mathbb{F}_2^{n_1})} \left( \prod_{i \in \delta_v} z_i \right) F_v(\mathbf{y}) + F'(\mathbf{y}, \mathbf{z}) \right] = \sum_{v \in \mathcal{R}(\mathbb{F}_2^{n_1})} F_v(\mathbf{y}), $$

and its degree is $d_G \le 2k - \gamma w \sim \gamma w(2 \log(l) - 1)$. If instead $l$ is odd, then $F'$ does not vanish when summing over because its monomials may appear an odd number of times. To circumvent this issue, we can guess a 0 position in each block and reduce the regular multivariate quadratic problem to a new instance of the new problem with $w$ windows of length $l-1$ (which is even). More in general, we can choose to initially fix an even number of zeros in each window when $l$ is even and an odd number of zeros in each window when $l$ is odd. In both cases, we are left with an even number $l' < l$ of variables for each window. From now on, we assume that $l$ is even, and take into account the cost of guessing the zeros in the final complexity.

In the standard setting, the efficiency of the interpolation and evaluation steps is guaranteed by the fact the Möbius transform allows to reconstruct the algebraic normal form (ANF) of a polynomial of degree $d$ from evaluations of vectors of weight at most $d$ only. In the regular setting, we can show that evaluations of regular vectors can be reconstructed from the knowledge of low-weight at most regular vectors. More precisely, we have the following results that is exploited throughout the algorithms to take advantage of the regular structure and improve on the complexity.
**Lemma 3.** Let $\mathbf{y}_1 \in \mathbb{F}_2^n$ and $\mathbf{y}_2 \in \mathcal{AMR}(\mathbb{F}_2^n)$, such that $\text{supp}(\mathbf{y}_1) \subseteq \text{supp}(\mathbf{y}_2)$. Then $\mathbf{y}_1 \in \mathcal{AMR}(\mathbb{F}_2^n)$.

*Proof.* In each window, $\mathbf{y}_2$ has at most one coordinate equal to 1. Since $\text{supp}(\mathbf{y}_1) \subseteq \text{supp}(\mathbf{y}_2)$, the same holds for $\mathbf{y}_1$. Therefore $\mathbf{y}_1$ is also at most regular. $\square$

**Proposition 6.** Let $F: \mathbb{F}_2^n \to \mathbb{F}_2$ be a boolean function of degree $d$. Let $\mathbf{x} \in \mathcal{R}(\mathbb{F}_2^n)$. Then the evaluations $F(\mathbf{x})$ can be computed from the evaluations of $F$ over all $\mathbf{y} \in \mathcal{AMR}(\binom{[n]}{d})$.

*Proof.* Let $S \subseteq \llbracket 1, n \rrbracket$ and $\mathbf{x}_S$ be the monomial that is the product of all variables whose index belongs to $S$. We can thus write

$$F(\mathbf{x}) \stackrel{\text{def}}{=} \sum_{|S| \le d} a_S \mathbf{x}^S.$$

We first prove that from the evaluations of all at most regular vectors of weight at most $d$, we obtain the value of each $a_S$ such that $\mathbf{x}_S$ is at most regular. We prove this first result by induction on the cardinality of $S$.

**Base case:** if $|S| = 0$, then $S = \emptyset$ and $a_0 = F(\mathbf{0})$ is known.

**Inductive step:** Assume we know $a_{\hat{S}}$ for all $\hat{S}$ such that $\mathbf{x}_{\hat{S}} \in \mathcal{AMR}(\binom{[n]}{d})$. Let $S$ be such that $\mathbf{x}_S \in \mathcal{AMR}(\binom{[n]}{w+1})$. We have

$$F(\mathbf{x}_S) = \sum_{\hat{S} \subseteq S} a_{\hat{S}} = a_S + \sum_{\hat{S} \subset S} a_{\hat{S}}.$$

For any $\hat{S} \subset S$, by the previous lemma $\mathbf{x}_{\hat{S}} \in \mathcal{AMR}(\mathbb{F}_2^n)$, moreover $\hat{S}$ has weight at most $w+1-1=w$, hence $\mathbf{x}_{\hat{S}} \in \mathcal{AMR}(\binom{[n]}{d})$ and $a_{\hat{S}}$ is known. Therefore $a_S$ can also be computed.

Now let $\mathbf{x}_S \in \mathcal{R}(\mathbb{F}_2^n)$. We have

$$F(\mathbf{x}_S) = \sum_{\hat{S} \subseteq S} a_{\hat{S}} = \sum_{\hat{S} \subset S, |\hat{S}| \le d} a_{\hat{S}},$$

where the last equality follows from the fact that $\text{deg} F = d$ and thus $a_{\hat{S}} = 0$ whenever $|\hat{S}| > d$. Again by the previous lemma, all subsets $\hat{S} \subseteq S$ are such that $\mathbf{x}_{\hat{S}}$ is at most regular. Hence $F(\mathbf{x}_S)$ can be reconstructed because all the relevant $a_{\hat{S}}$'s are known. $\square$

If $2k - \gamma w \ge (1 - \gamma)w$, we then need to compute $G(\mathbf{y})$ for all $\mathbf{y} \in \mathcal{AMR}(\mathbb{F}_2^{n-n_1})$, so that $G(\mathbf{y})$ can be interpolated. In this case, the parity counting computation already requires a number of operations given by $(l+1)^{(1-\gamma)w}l^{\gamma w} > l^w$ and is therefore worse than a brute force approach. We can thus restrict the (asymptotic) analysis to the case where

$$\gamma w(2 \log(l) - 1) < (1 - \gamma)w \\ \iff \gamma < \frac{1}{2 \log(l)}.$$

In such a case, the regular vectors $\mathbf{y}$'s for which $G$ needs to be calculated are only those with total weight at most $d_G$, whose number is

$$\sum_{i=0}^{2k-\gamma w} \binom{(1-\gamma)w}{i} l^i,$$

Then $G(\mathbf{y})$ is interpolated from these values and finally evaluated on all the $l^{(1-\gamma)w}$ regular vectors. The final complexity for the non-recursive version of the algorithm thus becomes

$$\mathcal{O}^*\left(\left(\sum_{i=0}^{2k-\gamma w} \binom{(1-\gamma)w}{i} l^i\right) l^{\gamma w} + l^{(1-\gamma)w}\right). \qquad (9)$$

We now have to determine the relevant function $g: [0, 1] \times Z_+ \to [0, 1]$ (that depends on $\gamma$ and $l$) for the regular setting. It is defined as the function such that

$$\mathcal{O}^*(2^{g(p,l)n}) = \mathcal{O}^*\left(\sum_{i=0}^{pw(2\log(l)-1)} \binom{(1-p)w}{i} l^i\right).$$

Since the sum on the right-hand side has a linear number of terms in $n$, the asymptotic is given by

$$\max_{i \in \llbracket 0, pw(2 \log(l)-1) \rrbracket} \binom{(1-p)w}{i} l^i.$$
The argmax can be determined as the $i$ for which two consecutive terms are as close as possible. By imposing

$$1 \approx \frac{\binom{(1-p)w}{i+1} l^{i+1}}{\binom{(1-p)w}{i} l^i} = \frac{l((1-p)w - i)}{i+1},$$

we obtain

$$i \sim \frac{l(1-p)w}{l+1} = \frac{(1-p)}{l+1}n.$$

Therefore, whenever $pw(2\log(l) - 1) \le \frac{l}{l+1}(1-p)w$, the maximum term is the last one, i.e.

$$\binom{(1-p)w}{pw(2\log(l) - 1)} l^{pw(2\log(l) - 1)}.$$

Otherwise, we obtain

$$2^{g(p,l)n} \sim \binom{(1-p)w}{\frac{(1-p)}{l+1}n} l^{\frac{(1-p)}{l+1}n} = (l+1)^{(1-p)w}.$$

This leads to the definition

$$g(p, l) \stackrel{\text{def}}{=} \begin{cases} \frac{1}{l} \left( (1-p)H\left(\frac{p(2\log(l)-1)}{1-p}\right) + p(2\log(l)-1)\log(l) \right), & \text{if } p(2\log(l)-1) \le \frac{l}{l+1}(1-p) \\ (1-p)^{\frac{\log(l+1)}{l}}, & \text{otherwise.} \end{cases}$$

For fixed $l$ and $\gamma$, the asymptotic exponent coefficient of our algorithm becomes

$$\max \left( g(\gamma, l) + \gamma \frac{\log(l)}{l}, (1-\gamma) \frac{\log(l)}{l} \right).$$

The probability of guessing correctly $l - l'$ zeros in each window amounts to $(l/l')^w$ and the partial evaluations on the guessed zeros for each system equation can be done in time $O^*(1)$. Therefore, the initial problem reduces in polynomial time to a new instance of regular multivariate quadratic with $w$ windows of length $l'$, so that the total number of variables is $n' = l'w = \frac{l'}{l}n$. Finally, when taking into account the initial guess of zeros in each window, the asymptotic complexity turns out to be

$$O^*(2^{\tau n})$$

where

$$\tau \stackrel{\text{def}}{=} \min \left\{ \frac{\log(l/l')}{l} + \frac{l'}{l} \max \left( g(\gamma, l') + \gamma \frac{\log(l')}{l'}, (1-\gamma) \frac{\log(l')}{l'} \right) \mid l' \text{ even}, 0 \le l' \le l, 0 \le \gamma \le \frac{1}{2\log(l')} \right\}.$$

### 5.2.1 Recursion à la Björklund et al.

We adapt now the recursion argument from [BKW19] in the context of regular solutions. For a fixed $y$ of low weight, instead of partially evaluating $F(y, z)$ over all $z \in \mathcal{R}(\mathbb{F}_2^{n_1})$ to obtain $G(y)$, many parity counting instances with a number of variables that decreases by a factor $\gamma$ are solved recursively. The term $l^{\gamma w}$ in (9) is thus replaced with the cost of the recursive call. Hence, fixed $\gamma$ and called $\tau_i$ the cost at depth $i$, the recursive formula (6) translates in the regular setting to

$$\tau_i = \gamma^i \max \left( g(\gamma, l) + \tau_{i+1}, (1-\gamma) \frac{\log(l)}{l} \right) \quad \forall i \in [0, D-1] \quad (10)$$

and

$$\tau_D = \gamma^D \frac{\log(l)}{l}.$$

In order to find the optimal value for $\gamma$, in a way analogous to [BKW19], we prove the following result.

**Proposition 7.** Let $0 \le \gamma \le \frac{1}{2\log(l)-1}$ be fixed and $g(\gamma, l) \le (1-\gamma)^2 \frac{\log(l)}{l}$. Then

$$\tau_0 \le (1-\gamma) \frac{\log(l)}{l}.$$
Proof. We proceed by induction. For sufficiently big D, we have

$$
\begin{align*}
\tau_{D-1} &= \gamma^{D-1} \max \left( g(\gamma, l) + \gamma_D \frac{\log(l)}{l}, (1-\gamma) \frac{\log(l)}{l} \right) \\
&= \gamma^{D-1} \max \left( (1-\gamma)^2 \frac{\log(l)}{l} + \gamma_D \frac{\log(l)}{l}, (1-\gamma) \frac{\log(l)}{l} \right) \\
&= \gamma^{D-1} \frac{\log(l)}{l} \max(1-2\gamma+\gamma^2+\gamma^D, 1-\gamma) \\
&= \gamma^{D-1} (1-\gamma) \frac{\log(l)}{l}.
\end{align*}
$$

Let us now assume that $\tau_{i+1} \le \gamma^{i+1}(1-\gamma) \frac{\log(l)}{l}$. The inductive step is to prove that

$$
\tau_i \le \gamma^i (1 - \gamma) \frac{\log(l)}{l},
$$

while the case $i+1 = D-1$ shown above represents the base case. From (10),

$$
\begin{align*}
\tau_i &= \gamma^i \max \left( \gamma^{i+1} (1-\gamma) \frac{\log(l)}{l} + \tau_{i+1}, (1-\gamma) \frac{\log(l)}{l} \right) \\
&\le \gamma^i \max \left( (1-\gamma)^2 \frac{\log(l)}{l} + \gamma^{i+1} (1-\gamma) \frac{\log(l)}{l}, (1-\gamma) \frac{\log(l)}{l} \right) \\
&= \gamma^i (1-\gamma) \frac{\log(l)}{l} \max((1-\gamma) + \gamma^{i+1}, 1) \\
&= \gamma^i (1-\gamma) \frac{\log(l)}{l}.
\end{align*}
$$

Hence, it follows by induction that $\tau_0 \le (1 - \gamma) \frac{\log(l)}{l}$. $\square$

Depending on the initial amounts $l - l'$ of fixed zeros per block, the final complexity can thus be upper bounded by

$$
\mathcal{O}^*(2^{\tau n})
$$

where

$$
\tau \stackrel{\text{def}}{=} \min \left\{ \frac{\log(l/l')}{l} + (1 - \gamma_{l'}) \frac{\log(l')}{l} \mid l' \text{even}, 0 \le l' \le l \right\}.
$$

and

$$
\gamma_{l'} \stackrel{\text{def}}{=} \max \left\{ \gamma \mid g(\gamma, l') \le (1-\gamma)^2 \frac{\log(l')}{l'} \right\}.
$$

We remark that from our experiments, the optimal value of $l'$ is always $l$ or $l-1$ depending on the parity of $l$, i.e. at most one 0 per window is guessed.

5.2.2 Recursion à la Dinur

Analogously to [Din21b], we define a multiple parity-counting problem that outputs

$$
(G(\mathbf{y}))_{\mathbf{y} \in \mathcal{AM}\mathcal{R}(\mathbb{F}_2^{(1-\gamma)n})}.
$$

Then, we further split the block $z$ of length $n_1 = \gamma n$ into two blocks $\mathbf{u}$ and $\mathbf{v}$ of length $n_1 - n_2 = \eta n$ and $n_2 = (\gamma - \eta)n$ respectively. In this way, a self-reduction to a new multiple parity-counting instance is applied and such instance outputs the vectors with coordinates

$$
G(\mathbf{y}, \mathbf{u}) = \sum_{\hat{\mathbf{v}} \in \mathcal{R}(\mathbb{F}_2^{n_2})} F(\mathbf{y}, \mathbf{u}, \hat{\mathbf{v}}),
$$

with $(\mathbf{y}, \mathbf{u}) \in \mathcal{AM}\mathcal{R}(\mathbb{F}_2^{n-n_1} \times \mathbb{F}_2^{n_1-n_2})$ of weight upper bounded by $\deg G(\mathbf{y}, \mathbf{u}) \le 2k' - \frac{2\log(l)-1}{l}n_2 \sim (\gamma - \eta)\frac{2\log(l)-1}{l}n$. The interpolating step through the Möbius transform has asymptotic cost

$$
\mathcal{O}^*(2^{g(\frac{n_2}{n}, l)n}). \qquad (11)
$$

Finally $G(\mathbf{y}, \mathbf{u})$ is evaluated on all the elements of

$$
\mathcal{AM}\mathcal{R}\left(\left(\begin{array}{c}
n-n_{1} \\
\downarrow \operatorname{deg} G(\mathbf{y})\end{array}\right)\right) \times \mathcal{R}(\mathbb{F}_{2}^{n_{1}-n_{2}}),
\qquad(12)
$$
so that the initial multi-parity counting output

$$ (G(\mathbf{y}))_{\mathbf{y} \in \binom{[n-n_1]}{\downarrow d_G}} $$

can be reconstructed.

The same reasoning exhibited before allows to evaluate the cost of each recursive call.

With the notation already introduced, at each step $1 \le i < D$, the cost of interpolation can be written as

$$ |\mathcal{AMR}\left(\left(\begin{array}{c}[n-n_2^{(i)}] \\ \downarrow n_2^{(i)}/l\end{array}\right)\right)| = \mathcal{O}^*(2^{g(\gamma-i\eta,l)n}), $$

while the evaluation step takes time

$$ |\mathcal{AMR}\left(\left(\begin{array}{c}[n-n_1^{(i)}] \\ \downarrow n_1^{(i)}/l\end{array}\right)\right)| \cdot |\mathcal{R}(\mathbb{F}_2^{n_1^{(i)}-n_2^{(i)}})| = \mathcal{O}^*(2^{(\eta \frac{\log(l)}{l} + g(\gamma-(i-1)\eta,l))n}). $$

We then have to add the cost of brute force at depth $D$, which is

$$ |\mathcal{AMR}\left(\left(\begin{array}{c}[n-n_1^{(i)}] \\ \downarrow n_1^{(i)}/l\end{array}\right)\right)| \cdot |\mathcal{R}(\mathbb{F}_2^{n_1^{(i)}-n_2^{(i)}})| \le \mathcal{O}^*(2^{(\eta \frac{\log(l)}{l} + g(\gamma-(D-1)\eta,l))n}). $$

and the cost of the evaluation of $G(\mathbf{y})$ for all elements $\mathbf{y} \in \mathcal{R}(\mathbb{F}_2^{n-n_1})$, i.e.

$$ \mathcal{O}^*(2^{(1-\gamma)\frac{\log(l)}{l}n}). $$

By taking into account the initial guess of zeros, the complexity exponent coefficient then becomes

$$ \mathcal{O}^*(2^{(\tau+\eta\frac{\log(l')}{l})n}). $$

where

$$ \tau \stackrel{\text{def}}{=} \min \left\{ \max \left\{ \frac{\log(l/l')}{l} + \frac{l'}{l} g(\gamma, l') \mid \gamma \in [0, \hat{\gamma}_{l'}) \right\} \mid l' \text{even}, 0 \le l' \le l \right\}. $$

and

$$ \hat{\gamma}_{l'} \stackrel{\text{def}}{=} \min \left\{ \gamma \mid g(\gamma, l') \ge (1-\gamma) \frac{\log(l')}{l'} \right\}. $$

Again, the term $\eta$ can be chosen arbitrarily small (and we thus neglect it in the comparison with other approaches) while keeping constant the number of iterations $D$. Moreover, as well as for the recursion à la Björklund et al., the optimal value of $l'$ is either $l$ or $l-1$, depending on the parity of $l$.

# 6 An alternative interpretation of boolean regular multivariate quadratic systems

Despite each entry of the solution vector can assume two values, it is evident that the information content/entropy carried out by the corresponding variable is lower than that for a random system with no constraints on the solution. It is then natural to wonder whether the solution can be written with less bits of information while preserving the possibility of expressing it as the solution of an algebraic system. We present here an alternative way to describe an instance of the regular multivariate quadratic problem that uses less variables, indeed. Let us consider a *binary* multivariate quadratic system whose solution of interest is regular and with the usual parameters. For the sake of simplicity, we restrict to the case where the window length is a power of 2, i.e. $l = 2^s$ for some positive integer $s$. The same description could be carried out to the more general case with a slightly suboptimal algebraic system.

We observe that the set of variables of $\mathbb{F}_2[x_1, \dots, x_l]$ can be mapped injectively into the set of monic polynomials of degree $s = \log(l)$ in $\mathbb{F}_2[x'_1, \dots, x'_s]$. A canonical example of such map is

$$ g: x_j \mapsto \prod_{a=1}^{s} (x_a' + \text{bin}_a(j-1)), $$

where $\text{bin}_a(b)$ denotes the $a$-th least significant bit of the binary representation of $b$. Then, $\tilde{g}$ induces a bijection between the set of binary vectors of length $l$ and weight 1 and the vector space $\mathbb{F}_2^s$: $v$ is sent to the unique vector $v' \in \mathbb{F}_2^s$ such that for all $i \in [1, l]$,

$$ \text{ev}_{x_i}(v) = v_i = \text{ev}_{g(x_i)}(v'), $$
where $ev_{P(x)}(v)$ is the evaluation $P(v)$.

The map $g$ extends via multiplication to monomials (including those formed by variables from different blocks). We split the $x'$ vector into $w$ blocks of length $s$ and denote its entries with two indexes, so that we define

$$ \bar{g}: \prod_{i=0}^{w-1} \prod_{j=1}^{l} x_{i,j}^{\alpha_{i,j}} \mapsto \prod_{i=0}^{w-1} \prod_{j=1}^{l} \prod_{a=1}^{s} (x'_{i,a} + \text{bin}_a(j-1))^{\alpha_{i,j}} $$

We are interested in evaluation of regular binary vectors, for our purposes it is thus enough to restrict the map $\bar{g}$ to polynomials of degree at most 1 in each block of variables:

$$ \bar{g}: \prod_{i=0}^{w-1} x_{i,j_i} \mapsto \prod_{i=0}^{w-1} \prod_{a=1}^{s} (x'_{i,a} + \text{bin}_a(j_i - 1)). \quad (13) $$

If we now consider an instance of the multivariate quadratic problem, the map $\bar{g}$ sends, by linearity, any multivariate quadratic equation into an equation of degree $2s$ in $sw$ variables. We thus reduced the multivariate quadratic problem to a multivariate degree-(2s) problem in less variables. Notice that even if the instance of the RMQ problem is picked randomly, the new system obtained still preserves some structure. Indeed the only terms that can appear are those coming from the expansion of

$$ \prod_{a=1}^{s} (x'_{i_1,a} + b_{i_1,a})(x'_{i_2,a} + b_{i_2,a}), \quad (14) $$

for some $b_{i_1,a}, b_{i_2,a} \in \mathbb{F}_2$ and $i_1 \neq i_2$. We will discuss how this structure affects the Hilbert series of the associated ideal and consequently the resolution of the system via Gröbner basis techniques later. On the other hand, the threshold ratio of equations over variables that corresponds to the unicity of the solution becomes

$$ \mu' = \frac{m'}{n'} = \frac{m}{sw} = 1. $$

This value resembles that of random systems and reflects the fact that no conditions are put on the shape of the solution anymore. We remove the prime symbol from now on.

## 6.1 Algebraic cryptanalysis of the alternative modeling

We want to study the Hilbert series of the homogenization of the system obtained from the construction above. We call $\mathcal{I}$ the ideal generated by it. First of all, we characterize the admissible terms in the Macaulay matrix at a given degree $d$. After homogenization, the only terms that can appear in a degree-(2s) equation can be derived by homogenizing (14):

$$ \prod_{a=1}^{s} (x_{i_1,a} + b_{i_1,a}y)(x_{i_2,a} + b_{i_2,a}y) $$

and then expanding it, i.e.

$$ \left\{ \left( \prod_{a_1 \in A_1} x_{i_1, a_1} \cdot \prod_{a_2 \in A_2} x_{i_2, a_2} \right) y^{2s - |A_1| - |A_2|} \mid i_1 \neq i_2, A_1, A_2 \subseteq [1, l] \right\}. $$

The characterization of admissible terms in the Macaulay matrix at a given degree $d \ge 2s$ boils down to listing all the terms that are contained in the ideal generated by the monomials above. It is straightforward to see that a monomial

$$ \prod_{i=0}^{w-1} \prod_{j=1}^{l} x_{i,j}^{\alpha_{i,j}} \cdot y^{\alpha_h} $$

can appear if and only if

$$ \exists i_1 \neq i_2 \text{ s.t. } \left( \sum_{j=1}^{l} \alpha_{i_1, j} + \alpha_{i_2, j} \right) + \alpha_h \geq 2s. $$

Let $\mathcal{M}_{NA}(d)$ be the set of non-admissible monomials at degree $d$,

$$ \mathcal{M}_{NA}(d) \stackrel{\text{def}}{=} \left\{ \prod_{i=0}^{w-1} \prod_{j=1}^{l} x_{i,j}^{\alpha_{i,j}} \cdot y^{\alpha_h} \mid \sum_{i,j} \alpha_{i,j} + \alpha_h = d \wedge \forall i_1 \neq i_2, \left( \sum_{j=1}^{l} \alpha_{i_1,j} + \alpha_{i_2,j} \right) + \alpha_h < 2s \right\}. $$

We are mainly interested in determining the maximum degree at which such set are not empty, because this implies a corresponding non-zero Hilbert function.
**Proposition 8.** For any $d > (s-1)w + 1$, $\mathcal{M}_{NA}(d) = \emptyset$.

*Proof.* Let us consider the degree-$d$ monomial $\prod_{i=0}^{w-1} \prod_{j=1}^{l} x_{i,j}^{\alpha_{i,j}} \cdot y^{\alpha_h}$, let $d_i = \sum_{j+1}^{l} \alpha_{i,j}$ be its degree restricted to the $i$-th window and $D_1 \ge D_2$ be the first and second maximums of $\{d_i\}$. The condition of non-admissibility requires $D_1 + D_2 + d_h \le 2s-1$, which implies $D_2 \le \frac{D_1+D_2}{2} \le \frac{D_1+D_2}{2} = d_h \le s-1$. Therefore

$$d = d_h + \sum_i d_i \le D_1 + D_2 + (w-2)D_2 + d_h \le 2s-1 + (w-2)(s-1) = (s-1)w + 1.$$

Hence, for any $d > (s-1)w + 1$, $\mathcal{M}_{NA}(d)$ does not contain any monomial. $\square$

It immediately follows from the definition of the Hilbert series that

$$\mathrm{HF}_{R/\mathcal{I}}(d) \ge |\mathcal{M}_{NA}(d)|.$$

Experimentally, we found out that, beside non-admissible monomials, the homogenized system $S^{(y)}$ behaves as a semi-regular one. More precisely,

**Heuristic 1.** Let $S_{\mathrm{init}} = \{f_1, \dots, f_m\} \subseteq \mathbb{R} = \mathbb{F}_2[x_1, \dots, x_n]$ be a randomly picked system, with $n = 2^s w$. Let $S'_{\mathrm{init}} = \{\bar{g}(f_1), \dots, \bar{g}(f_m)\} \subseteq R' = \mathbb{F}_2[x'_1, \dots, x'_{sw}]$ be set of polynomials obtained by applying the map $\bar{g}$ of Equation (13) to $S_{\mathrm{init}}$ and $S'_{\mathrm{FE}}$ the set of polynomials corresponding to field equations for $x'$'s entries. Define the homogeneous ideal $\mathcal{I} = (\langle S'_{\mathrm{init}} \rangle^{(y)} \cup \langle S'_{\mathrm{FE}} \rangle^{(y)}\rangle$. Then with high-probability

$$\mathrm{HF}_{R'/\mathcal{I}}(d) = \max \left( [z^d] \left[ \frac{(1+z)^{sw}}{(1-z)(1+z^{2s})^{sw}} \right]_+, |\mathcal{M}_{NA}(d)| \right),$$

where $[z^d]h(z)$ is the coefficient of the term of degree $d$ of the series $h(z)$.

The existence of non-admissible monomials potentially increases the degree of regularity of the homogenized system. In practice, it is reasonable to consider that this does not reflect on a higher solving degree. Indeed, when solving the affine system with XL algorithm or through a Gröbner basis computation. On the opposite, we experimentally notice that for small concrete parameters, the degree falls may even let the solving degree decrease. Table 1 provides a comparison of the experimental solving degree and times for some small parameters between random instances of degree-2s polynomial systems and systems originated from boolean regular multivariate quadratic instances with the procedure explained above. We first remark that for random degree-2s systems, the solving degree always coincides with the degree of regularity of the system before homogenizing, i.e. with the degree of regularity of the part of highest degree, which is consistent with the possibility of exploiting degree falls during the Gröbner basis computation. Moreover, we can see that the smaller solving degree is also generally reflected in a smaller time complexity. An exception is given by the case $s=1$, i.e. $l=2$, for which the results coincide. This is not surprising, for $l=2$ the solution has the same number of 0s and 1s, which means that the RMQ problem could be rewritten as the standard MQ problem with half the number of variables. A smaller complexity with our alternative approach would thus imply a better Gröbner basis algorithm for the multivariate quadratic problem.

In any case, we expect the degree of regularity to grow asymptotically in a very similar manner to the one for a random system of equations of degree 2s with $m$ equations and $sw$ unknowns, and we use it to benchmark the cost of this alternative approach. In other words, we expect our analysis to give a pretty tight upper bound on the effective complexity of this method.

### 6.1.1 Hybrid approaches

As well as for the standard modeling, the Gröbner basis computation can be sped up by specializing some variables, at the cost of repeating the process until a good guess is found. The analogy goes further, as there exist a correspondence between specializing a variable in the alternative modeling and fixing several variables in the standard modeling according to the methods explained before. Indeed, let us restrict to the analysis of one window. From the definition of the bijective map $g: \mathbb{F}_2[x_1, \dots, x_l] \to \mathbb{F}_2[x'_1, \dots, x'_s]$, we have the following equivalences:

$$\begin{align*}
x'_a &= 0 &\iff \forall j \text{ s.t. } \mathrm{bin}_a(j-1) &= 0, x_j = 0, \\
x'_a &= 1 &\iff \forall j \text{ s.t. } \mathrm{bin}_a(j-1) &= 1, x_j = 0,
\end{align*}$$

i.e. specializing one variable in the alternative modeling coincides with guessing $2^{s-1} = l/2$ zeros in a same window for the standard modeling, as the $a$-th binary digits of an element in $[0, 2^s - 1]$ is equally distributed between the values 0 and 1. Moreover, these $s$ distributions for $a \in [1, s]$ are mutually independent, hence specializing $s-s'$ unknowns from a same block corresponds to choosing $\sum_{i=1}^{s'} 2^{s-i}$ zero positions. In particular,
<table><caption>Table 1: Comparison of concrete complexities of the plain Gröbner basis approach for the standard quadratic modeling, the alternative degree-(2s) modeling and a random degree-(2s) system with *lw* unknowns. For a comparison with the alternative modeling we only pick regular instances with *w* blocks of length *l* = 2s. The ratio *μ* = 1.2 is slightly above the uniqueness bound to guarantee with high probability the uniqueness of the (forced) regular solution. The random system instead is forced to have a random solution. The "d<sub>solv</sub>" column stands for the maximum degree reached by the MAGMA built-in Gröbner basis function. The time is computed in seconds and the column "Rows" denotes the largest number of rows of a (Macaulay) matrix computed during the Gröbner basis computation. The calculations have been run on one core of an AMD EPYC 9374F processor with the MAGMA version V2.28-7.</caption><thead><tr><th rowspan="2">(s, w)</th><th colspan="3">Quadratic modeling</th><th colspan="3">Alternative modeling</th><th colspan="3">Random degree-(2s)</th></tr><tr><th>d<sub>solv</sub></th><th>Time</th><th>Matrix</th><th>d<sub>solv</sub></th><th>Time</th><th>Matrix</th><th>d<sub>solv</sub></th><th>Time</th><th>Matrix</th></tr></thead><tbody><tr><td>(1, 20)</td><td>5</td><td>0.92</td><td>1.1 · 10<sup>5</sup></td><td>5</td><td>0.33</td><td>2.3 · 10<sup>4</sup></td><td>5</td><td>0.34</td><td>2.4 · 10<sup>4</sup></td></tr><tr><td>(1, 25)</td><td>5</td><td>270</td><td>3.8 · 10<sup>5</sup></td><td>5</td><td>9.1</td><td>6.0 · 10<sup>4</sup></td><td>5</td><td>9.1</td><td>5.8 · 10<sup>4</sup></td></tr><tr><td>(1, 30)</td><td>6</td><td>7.1 · 10<sup>3</sup></td><td>3.1 · 10<sup>6</sup></td><td>6</td><td>410</td><td>6.3 · 10<sup>5</sup></td><td>6</td><td>420</td><td>6.3 · 10<sup>5</sup></td></tr><tr><td>(2, 5)</td><td>4</td><td>0.01</td><td>2.0 · 10<sup>3</sup></td><td>5</td><td>0</td><td>868</td><td>6</td><td>0.01</td><td>640</td></tr><tr><td>(2, 6)</td><td>4</td><td>0.01</td><td>7.4 · 10<sup>3</sup></td><td>6</td><td>0.06</td><td>3.2 · 10<sup>3</sup></td><td>7</td><td>0.1</td><td>3.3 · 10<sup>3</sup></td></tr><tr><td>(2, 7)</td><td>5</td><td>0.25</td><td>3.1 · 10<sup>4</sup></td><td>6</td><td>0.35</td><td>6.1 · 10<sup>3</sup></td><td>8</td><td>1.23</td><td>1.010<sup>4</sup></td></tr><tr><td>(2, 8)</td><td>5</td><td>1.0</td><td>1.0 · 10<sup>5</sup></td><td>6(7)</td><td>2.1</td><td>3.7 · 10<sup>4</sup></td><td>8</td><td>22</td><td>7.1 · 10<sup>4</sup></td></tr><tr><td>(2, 9)</td><td>5</td><td>4.8</td><td>8.9 · 10<sup>4</sup></td><td>7</td><td>40</td><td>1.3 · 10<sup>5</sup></td><td>9</td><td>310</td><td>1.5 · 10<sup>5</sup></td></tr><tr><td>(2, 10)</td><td>5</td><td>63</td><td>4.3 · 10<sup>5</sup></td><td>7</td><td>360</td><td>1.2 · 10<sup>5</sup></td><td>9</td><td>1.0 · 10<sup>4</sup></td><td>8.6 · 10<sup>5</sup></td></tr><tr><td>(2, 11)</td><td>6</td><td>530</td><td>1.0 · 10<sup>6</sup></td><td>8</td><td>5.8 · 10<sup>3</sup></td><td>8.6 · 10<sup>5</sup></td><td></td><td></td><td></td></tr><tr><td>(2, 12)</td><td>6</td><td>2.8 · 10<sup>3</sup></td><td>1.7 · 10<sup>6</sup></td><td>8</td><td>5.0 · 10<sup>4</sup></td><td>3.6 · 10<sup>5</sup></td><td></td><td></td><td></td></tr><tr><td>(2, 13)</td><td>6</td><td>1.7 · 10<sup>4</sup></td><td>2.6 · 10<sup>6</sup></td><td>≥ 8</td><td></td><td>≥ 4.0 · 10<sup>5</sup></td><td></td><td></td><td></td></tr><tr><td>(3, 4)</td><td>4</td><td>0.1</td><td>1.9 · 10<sup>4</sup></td><td>8</td><td>0.21</td><td>5.4 · 10<sup>3</sup></td><td>9</td><td>0.28</td><td>5.1 · 10<sup>3</sup></td></tr><tr><td>(3, 5)</td><td>5</td><td>15</td><td>2.4 · 10<sup>5</sup></td><td>9</td><td>15</td><td>4.8 · 10<sup>4</sup></td><td>10</td><td>20</td><td>3.0 · 10<sup>4</sup></td></tr><tr><td>(3, 6)</td><td>5</td><td>170</td><td>6.7 · 10<sup>5</sup></td><td>9</td><td>710</td><td>2.4 · 10<sup>5</sup></td><td>11</td><td>4.3 · 10<sup>3</sup></td><td>2.2 · 10<sup>5</sup></td></tr><tr><td>(3, 7)</td><td>6</td><td>3.1 · 10<sup>4</sup></td><td>7.4 · 10<sup>6</sup></td><td>≥ 9</td><td></td><td>≥ 12 · 10<sup>5</sup></td><td>≥ 12</td><td></td><td></td></tr><tr><td>(4, 2)</td><td>3</td><td>0.01</td><td>6.7 · 10<sup>3</sup></td><td>8</td><td>0</td><td>300</td><td>8</td><td>0.01</td><td>170</td></tr><tr><td>(4, 3)</td><td>4</td><td>1.2</td><td>9.5 · 10<sup>4</sup></td><td>9</td><td>0.27</td><td>3.3 · 10<sup>3</sup></td><td>10</td><td>0.3</td><td>4.4 · 10<sup>3</sup></td></tr><tr><td>(4, 4)</td><td>4</td><td>190</td><td>3.7 · 10<sup>5</sup></td><td>10</td><td>86</td><td>5.5 · 10<sup>4</sup></td><td>12</td><td>180</td><td>4.6 · 10<sup>4</sup></td></tr><tr><td>(5, 2)</td><td>3</td><td>0.31</td><td>7.4 · 10<sup>4</sup></td><td>10</td><td>0.02</td><td>1.1 · 10<sup>3</sup></td><td>10</td><td>0.02</td><td>980</td></tr><tr><td>(5, 3)</td><td>4</td><td>82</td><td>1.5 · 10<sup>6</sup></td><td>11</td><td>32</td><td>4.9 · 10<sup>4</sup></td><td>13</td><td>36</td><td>2.6 · 10<sup>4</sup></td></tr></tbody></table>
guessing all s variables reduces to guess all the $2^s - 1$ zeros, i.e. the nonzero position. The guessing strategies are therefore very similar, however the alternative modeling admits a less refined choice for partial-window guess, since only a number of variables that is a power of 2 can left unguessed in each window.

The hybrid approach corresponding to full-window guess boils down, for the alternative modeling, to determine the optimal value $0 \le \gamma \le 1$ such that specializing all $\gamma sw$ variables from $\gamma w$ blocks minimizes the complexity of solving the system.

Therefore we study the Hilbert series

$$HS_{R/\mathcal{I}}(z) = \left[ \frac{(1+z)^{(1-\gamma)sw}}{(1-z)(1+z^{2s})^m} \right]_+ .$$

We can thus write the Hilbert function at degree $d$ as

$$HF_{R/\mathcal{I}}(d) = I_d(n)$$

for the Cauchy integrals

$$I_d(n) \stackrel{\text{def}}{=} \frac{1}{2\pi i} \oint \frac{1}{z^{d+1}} \frac{(1+z)^{(1-\gamma)sw}}{(1-z)(1+z^{2s})^m} dz.$$

We set $f(z)$ such that $e^{nf(z)} = \frac{1}{z^{d+1}} \frac{(1+z)^{(1-\gamma)sw}}{(1-z)(1+z^{2s})^m}$, i.e.

$$f(z) \stackrel{\text{def}}{=} \frac{(1-\gamma)sw}{n} \ln(1+z) - \mu \ln(1+z^{2s}) - \frac{1}{n} \ln(1-z) - \frac{d+1}{n} \ln(z),$$

with $\mu = m/n$. We compute the saddle points as

$$zf'(z) = \frac{(1-\gamma)sw}{n} \frac{z}{1+z} - 2s\mu \frac{z^{2s}}{1+z^{2s}} + \frac{1}{n} \frac{z}{1-z} - \frac{d+1}{n}.$$

We observe that the term $\frac{1}{n} \frac{z}{1-z}$ is negligible when $n \to \infty$. Hence, after multiplying by all denominators and defining the degree-(2 log(l) + 1) polynomial

$$\begin{aligned} g(z) &= \frac{(1-\gamma)\log(l)w}{n}z(1+z^{2\log(l)}) - 2\log(l)\mu z^{2\log(l)}(1+z) - \frac{d+1}{n}(1+z)(1+z^{2\log(l)}) \\ &= \left( \frac{(1-\gamma)\log(l)w}{n} - 2\log(l)\mu - \frac{d+1}{n} \right) z^{2\log(l)+1} + \left( -2\log(l)\mu - \frac{d+1}{n} \right) z^{2\log(l)} \\ &\quad + \left( \frac{(1-\gamma)\log(l)w}{n} - \frac{d+1}{n} \right) z - \frac{d+1}{n}, \end{aligned}$$

the asymptotic relative degree of regularity $\delta \stackrel{\text{def}}{=} \frac{d+1}{n}$ is given by the smallest positive root of the resultant Res(g, g'), whose degree is $2\deg(g) - 1 = 4\log(l) + 1$:

$$\bar{\delta} = \min\{\delta \in \mathbb{R} \mid (\text{Res}(g, g'))(\delta) = 0, \delta > 0\}.$$

The final asymptotic complexity thus becomes

$$O^* \left( 2^{\gamma \frac{\log(l)}{l} n} \cdot \left( \frac{\frac{(1-\gamma)\log(l)}{l}}{\bar{\delta} n} n \right)^{\omega} \right) = O^* (2^{\tau n}),$$

$$\text{with } \tau = \gamma \frac{\log(l)}{l} + \omega(1-\gamma) \frac{\log(l)}{l} H \left( \frac{\bar{\delta} \cdot l}{(1-\gamma) \log(l)} \right).$$

For the partial-window strategy, we notice that specializing variables decreases the degree of the terms involving that variables blocks. Since the same amount $s - s'$ of variables is guessed in each block, the degree decreases from $2s$ to $2s'$. The relevant Hilbert series becomes instead

$$HS_{R/\mathcal{I}}(z) = \left[ \frac{(1+z)^{s'w}}{(1-z)(1+z^{2s'})^m} \right]_+ .$$

We can thus write the Hilbert function at degree $d$ as

$$HF_{R/\mathcal{I}}(d) = I_d(n)$$

for the Cauchy integrals

$$I_d(n) \stackrel{\text{def}}{=} \frac{1}{2\pi i} \oint \frac{1}{z^{d+1}} \frac{(1+z)^{s'w}}{(1-z)(1+z^{2s})^m} dz.$$
We set $f(z)$ such that $e^{nf(z)} = \frac{1}{z^{d+1}} \frac{(1+z)^{s'w}}{(1-z)(1+z^{2s'})^m}$, i.e.

$$f(z) \stackrel{\text{def}}{=} \frac{s'w}{n} \ln(1+z) - \mu \ln(1+z^{2s'}) - \frac{1}{n} \ln(1-z) - \frac{d+1}{n} \ln(z),$$

with $\mu = m/n$. We compute the saddle points as

$$zf'(z) = \frac{s'w}{n} \frac{z}{1+z} - 2s'\mu \frac{z^{2s'}}{1+z^{2s'}} + \frac{1}{n} \frac{z}{1-z} - \frac{d+1}{n}.$$

We observe that the term $\frac{z}{1+z}$ is negligible when $n \to \infty$. Let $l' = s^{s'}$. After multiplying by all denominators and defining the degree-(2 log(l') + 1) polynomial

$$\begin{aligned} g(z) &= \frac{\log(l')w}{n}z(1+z^{2\log(l')}) - 2\log(l')\mu z^{2\log(l')}(1+z) - \frac{d+1}{n}(1+z)(1+z^{2\log(l')}) \\ &= \left(\frac{\log(l')w}{n} - 2\log(l')\mu - \frac{d+1}{n}\right)z^{2\log(l') + 1} + \left(-2\log(l')\mu - \frac{d+1}{n}\right)z^{2\log(l')} \\ &\quad + \left(\frac{\log(l')w}{n} - \frac{d+1}{n}\right)z - \frac{d+1}{n}, \end{aligned}$$

the asymptotic relative degree of regularity $\delta \stackrel{\text{def}}{=} \frac{d+1}{n}$ is given by the smallest positive root of the resultant Res(g, g'), whose degree is $2\deg(g) - 1 = 4\log(l') + 1$:

$$\bar{\delta} = \min\{\delta \in \mathbb{R} \mid (\text{Res}(g, g'))(\delta) = 0, \ \delta > 0\}.$$

The final asymptotic complexity thus becomes

$$\mathcal{O}^* \left( 2^{\frac{\log(l/l')}{l}} n \cdot \left( \frac{\log(l')}{\bar{\delta}n} \right)^{\omega} \right) = \mathcal{O}^* (2^{\tau n}),$$

$$\text{with } \tau = \frac{\log(l')}{l} + \omega \frac{\log(l')}{l} H \left( \frac{\bar{\delta} \cdot l}{\log(l')} \right).$$

## 6.2 Polynomial method on the alternative modeling

In our previous analysis of the probabilistic polynomial method, we focused on the case of quadratic systems. The general idea of the algorithm, however, works for equations of higher degree as well and each of the already mentioned works [LPT+17, BKW19, Din21b, Din21a] provides the estimated complexity in function of the degree. The aim of this subsection is to provide the complexity of the polynomial method to the system of $m$ equations of degree $2s$ in $sw$ variables turning out from the application of the map $\bar{g}$. Since the solutions are no longer subject to the regular constraints, the cost estimate is a mere application of the best result among all the probabilistic method variants with respect to polynomials of degree $2s$. In particular the recursive version of [Din21b] holds the asymptotic complexity record of $\mathcal{O}^*(2^{(1-1/(2d))n})$ for any degree $d > 2$. We recall that [Din21b] was also the best for degree 2. When applied to our system, we thus obtain a complexity of

$$\mathcal{O}^*(2^{(1-1/(4s))sw}) = \mathcal{O}^*\left(2^{\left(\frac{\log(l)}{l} - \frac{1}{4l}\right)n}\right).$$

## 7 Comparison

We provide here below figures with the asymptotic complexities for the plain and hybrid Gröbner basis approaches assuming $\omega = 2$ and $\mu = \log(l)/l$, as well as for the polynomial methods and we provide the exhaustive search as a benchmark. We show both the absolute complexity and the relative complexity with respect to brute force, i.e. where we reparametrize in such a way that brute force has unitary cost. In the latter case, we note that plain Gröbner basis does not tend to 1. The following explanation provides an intuition of this behavior. The extremal case is when there is only one window and the total weight is 1. Indeed, the absolute degree of regularity tends to 2 and thus the (nonconstant) complexity coefficient is given by $\omega H(2/n)$, while for brute force it is clearly $\log(n)$. We thus obtain

$$\lim_{n \to \infty} \frac{\omega H(2/n)}{\log(n)} = \lim_{n \to \infty} \frac{\omega(2 \log(n) - 2 - (n-2) \log(n) + (n-2) \log(n-2))}{\log(n)} = 2\omega = 4.$$

Note that for $l = 2$ the different window guess coincides with the full window guess, because either 0 or 1 entries per window are guessed. For the same reason, the hybrid partial-window guess is nothing but the best between
<figure><img src="image_3.png" /><figcaption>Figure 1: Absolute complexity, comparison among algebraic and probabilistic methods</figcaption></figure>

a plain Gröbner basis and the exhaustive search, with the former being better for $l = 2$. However, for $l \ge 3$, the hybrid partial-window guess becomes better than the full-window guess approach. Overall, the best asymptotic approach is the hybrid approach with a different number of guesses per window, as it is the one that makes use of more parameters to optimize the cost among the Gröbner basis methods. The drawback is that even finding the best parameters may become non-trivial as $l$ increases. On the other hand, we observe that already for $l \ge 4$, the method seems to coincide (or almost coincide) with the sub-case of hybrid partial-window guess, and we can thus consider the latter to quite reliably estimate the asymptotic complexity for a fixed $l$.

We remark that the probabilistic methods can still be relevant in regimes other than the one analyzed here for $\mu = \frac{\log(l)}{l}$. Indeed, differently from a Gröbner basis computation, whose complexity highly varies with the number of equations, the formers are mostly independent from them. Therefore, if the system is enough underdetermined, then the probabilistic methods may still outperform asymptotically Gröbner basis approaches in finding a regular solution.

Regarding the complexity for the solvers of the alternative modeling (Figures 3 and 4, we see that a plain Gröbner basis beats that of original modeling for big enough windows (the smallest power of 2 for which this is realized is $l = 2^5$). The best value of $\gamma$ for the hybrid approach quickly tends to 1 with $l$ growing, making this method coincide with brute force and the two 2 curves almost overlap already starting from $l = 8$. This is due to the growth of the polynomial degrees, the algebraic part becomes less competitive as the latter increase. For $l = \{4, 8, 16, 32\}$, the two modelings with hybrid partial-window specialization coincide. This follows from Remark 5 and the fact that the RMQ problem with $l' = 2$ reduces to the standard MQ problem with half the number of variables. We also note that, when $s' = 1$, the (specialized) alternative modeling has just a random structure, because of degree $2s' = 2$, i.e. it can be treated as an overdetermined random quadratic system.

Due to all what have been said so far, we conclude that a pretty simple and elegant, yet accurate, estimation of the complexity of an RMQ instance at the uniqueness bound is obtained by guessing $l-2$ in each block and thus reducing the problem to many MQ instances in $w$ unknowns and $\log(l)w$ quadratic equations, i.e. with $\mu = \log(l)$. For $0.55\mu \ge 1$, which corresponds to $l \ge 4$, we have already said that the best trade-off consists of a plain Gröbner basis computation. The simplicity of this analysis stems from the existence of a closed form asymptotic expression for the relative degree of regularity $\delta$. We recall indeed that [Bar04, Corollary 4.4.3] determines the asymptotic of $\delta$ for a semi-regular quadratic system over $\mathbb{F}_2$ with ratio $\mu$ as

$$ \delta \to \bar{\delta} \stackrel{\text{def}}{=} -\mu + \frac{1}{2} + \frac{1}{2} \sqrt{2\mu^2 - 10\mu - 1 + 2(\mu + 2)\sqrt{\mu(\mu + 2)}}. \quad (15) $$

Hence an asymptotic upper bound on the complexity of solving a boolean RMQ instance with $w$ blocks of length $l$ is

$$ \mathcal{O}^* \left( \left( \left( \frac{l}{l'} \right)^w \cdot \left( \frac{w}{\bar{\delta}w} \right)^\omega \right) = \mathcal{O}^* \left( 2^{\left( \frac{\log(1/l') + \omega H(\delta)}{l} \right) n} \right), $$
<figure><img src="image_4.png" /><figcaption>Figure 2: Complexity relative to brute force, comparison among algebraic and probabilistic methods</figcaption></figure>

where $\bar{\delta}$ is defined as in (15) with $\mu = \log(l)$.

Finally, our estimates show that the partial/different-window hybrid technique remains asymptotically the most efficient for higher values of $l$, even compared to other methods. Indeed, for $l = 2$ the recursion à la Dinur in the probabilistic method for the alternative modeling beats any other approach, which is consistent with the probabilistic method being asymptotically better than Gröbner basis techniques for MQ instances at the uniqueness bound. However, the former, when applied to either the standard or the alternative modeling, becomes more costly already for $l \ge 4$.

# 8 Conclusions

We presented a new NP-complete variant of the multivariate quadratic problem where, instead of considering structured polynomials, the additional property of having weight-1 blocks is required from the solution. We reviewed all the most competitive techniques for solving polynomial systems over finite fields. By doing this, we devised tailored approaches that take into account the regular constraint. Despite the non-negligible improvements over the standard multivariate quadratic problem, our analysis shows that, even in the average case, the problem remains hard. Indeed, the complexity is still exponential and, for block lengths of cryptographic interest, has a not-so-small constant at the exponent.

Algebraic cryptanalysis resulted in being, asymptotically, the most efficient method. This is surprising to some extent because it is not what happens in the usual MQ setting, where the best polynomial methods beat the hybrid Gröbner basis computation for sufficiently large parameters, at least in theory. Moreover, as remarked in the corresponding section, the best hybrid strategy is very close in essence to an exhaustive search rather than a plain Gröbner basis computation, meaning that the best trade-off is achieved by specializing almost all variables in each block. This fact is an additional piece of evidence that the problem is difficult in practice.

Even more so for concrete instances, we expect hybrid algebraic cryptanalysis to perform better than polynomial methods, as the latter have not found much use in practice for cryptographic parameters, not even for the MQ problem. In this regard, an interesting open problem would be to study the efficiency of the Crossbred algorithm with respect to RMQ and possibly adapt it to this scenario, as this is the only work that defeats fast exhaustive search for practical parameters.

The RMQ problem naturally inherits properties from MQ and RSD problems, and both revealed to be well suited for cryptographic constructions, for instance for designing digital signatures. It is then tantalizing to conjecture that the RMQ problem could have many cryptographic applications as well, and we leave for future work devising competitive protocols based on the corresponding hardness assumption.
<figure><img src="image_5.png" /><figcaption>Figure 3: Absolute complexity, comparison with the methods for the alternative system description</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Complexity relative to brute force, comparison with the methods for the alternative system description</figcaption></figure>
References

[AFS05] Daniel Augot, Matthieu Finiasz, and Nicolas Sendrier. A family of fast syndrome based cryptographic hash functions. In *Progress in Cryptology–Mycrypt 2005: First International Conference on Cryptology in Malaysia*, Kuala Lumpur, Malaysia, September 28-30, 2005. *Proceedings 1*, pages 64–83. Springer, 2005.

[Bar04] Magali Bardet. *Étude des systèmes algébriques surdéterminés. Applications aux codes correcteurs et à la cryptographie*. PhD thesis, Université Pierre et Marie Curie-Paris VI, 2004.

[BBSV22] Stefano Barbero, Emanuele Bellini, Carlo Sanna, and Javier Verbel. Practical complexities of probabilistic algorithms for solving boolean polynomial systems. *Discrete Applied Mathematics*, 309:13–31, 2022.

[BCC$^{+}$10] Charles Bouillaguet, Hsieh-Chung Chen, Chen-Mou Cheng, Tung Chou, Ruben Niederhagen, Adi Shamir, and Bo-Yin Yang. Fast exhaustive search for polynomial systems in $f^2$. In *International Workshop on Cryptographic Hardware and Embedded Systems*, pages 203–218. Springer, 2010.

[BCC$^{+}$24] Dung Bui, Eliana Carozza, Geoffroy Couteau, Dahmun Goudarzi, and Antoine Joux. Faster signatures from mpc-in-the-head. In *ASIACRYPT 2024-International Conference on the Theory and Application of Cryptology and Information Security*, 2024.

[BCGI18] Elette Boyle, Geoffroy Couteau, Niv Gilboa, and Yuval Ishai. Compressing vector ole. In *Proceedings of the 2018 ACM SIGSAC Conference on Computer and Communications Security*, pages 896–912, 2018.

[BFSS13] Magali Bardet, Jean-Charles Faugère, Bruno Salvy, and Pierre-Jean Spaenlehauer. On the complexity of solving quadratic boolean systems. *Journal of Complexity*, 29(1):53–75, 2013.

[BKW19] Andreas Björklund, Petteri Kaski, and Ryan Williams. Solving systems of polynomial equations over gf (2) by a parity-counting self-reduction. In *International Colloquium on Automata, Languages and Programming*, pages 1–13. Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2019.

[BMVT78] Elwyn Berlekamp, Robert McEliece, and Henk Van Tilborg. On the inherent intractability of certain coding problems (corresp.). *IEEE Transactions on Information theory*, 24(3):384–386, 1978.

[BØ23] Pierre Briaud and Morten Øygarden. A new algebraic approach to the regular syndrome decoding problem and implications for pcg constructions. In *Annual International Conference on the Theory and Applications of Cryptographic Techniques*, pages 391–422. Springer, 2023.

[CCJ23] Eliana Carozza, Geoffroy Couteau, and Antoine Joux. Short signatures from regular syndrome decoding in the head. In *Annual International Conference on the Theory and Applications of Cryptographic Techniques*, pages 532–563. Springer, 2023.

[CKPS00] Nicolas Courtois, Alexander Klimov, Jacques Patarin, and Adi Shamir. Efficient algorithms for solving overdefined systems of multivariate polynomial equations. In *International Conference on the Theory and Applications of Cryptographic Techniques*, pages 392–407. Springer, 2000.

[CLY$^{+}$24] Hongrui Cui, Hanlin Liu, Di Yan, Kang Yang, Yu Yu, and Kaiyi Zhang. Resolved: Shorter signatures from regular syndrome decoding and vole-in-the-head. In *IACR International Conference on Public-Key Cryptography*, pages 229–258. Springer, 2024.

[Cop94] Don Coppersmith. Solving homogeneous linear equations over gf (2) via block wiedemann algorithm. *Mathematics of Computation*, 62(205):333–350, 1994.

[Din21a] Itai Dinur. Cryptanalytic applications of the polynomial method for solving multivariate equation systems over gf (2). In *Annual International Conference on the Theory and Applications of Cryptographic Techniques*, pages 374–403. Springer, 2021.

[Din21b] Itai Dinur. Improved algorithms for solving polynomial systems over gf (2) by multiple parity-counting. In *Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA)*, pages 2550–2564. SIAM, 2021.

[ES24] Andre Esser and Paolo Santini. Not just regular decoding: asymptotics and improvements of regular syndrome decoding attacks. In *Annual International Cryptology Conference*, pages 183–217. Springer, 2024.
[Frö85] Ralf Fröberg. An inequality for hilbert series of graded algebras. *Mathematica Scandinavica*, 56(2):117–144, 1985.

[FY79] Aviezri S Fraenkel and Yaacov Yesha. Complexity of problems in games, graphs and algebraic equations. *Discrete Applied Mathematics*, 1(1-2):15–30, 1979.

[HOSSV22] Carmit Hazay, Emmanuela Orsini, Peter Scholl, and Eduardo Soria-Vazquez. Tinykeys: A new approach to efficient multi-party computation. *Journal of Cryptology*, 35(2):13, 2022.

[LPT$^+$17] Daniel Lokshtanov, Ramamohan Paturi, Suguru Tamaki, Ryan Williams, and Huacheng Yu. Beat- ing brute force for systems of polynomial equations over finite fields. In *Proceedings of the Twenty- Eighth Annual ACM-SIAM Symposium on Discrete Algorithms*, pages 2190–2202. SIAM, 2017.

[Mac94] Francis Sowerby Macaulay. *The algebraic theory of modular systems*. Cambridge University Press, 1994.

[Raz87] Alexander A Razborov. Lower bounds on the size of bounded depth circuits over a complete basis with logical addition. *Mat. Zametki*, 41(4):598–607, 1987.

[SF09] Robert Sedgewick and P Flajolet. *Analytic combinatorics*. Cambridge University Press Cambridge, 2009.

[Smo87] Roman Smolensky. Algebraic methods in the theory of lower bounds for boolean circuit complexity. In *Proceedings of the nineteenth annual ACM symposium on Theory of computing*, pages 77–82, 1987.

[VV85] Leslie G Valiant and Vijay V Vazirani. Np is as easy as detecting unique solutions. In *Proceedings of the seventeenth annual ACM symposium on Theory of computing*, pages 458–463, 1985.

[Wie86] Douglas Wiedemann. Solving sparse linear equations over finite fields. *IEEE transactions on information theory*, 32(1):54–62, 1986.

[Will14] Richard Ryan Williams. The polynomial method in circuit complexity applied to algorithm design (invited talk). In *34th International Conference on Foundation of Software Technology and Theoretical Computer Science (FSTTCS 2014)*. Schloss-Dagstuhl-Leibniz Zentrum für Informatik, 2014.

[Won01] Roderick Wong. *Asymptotic approximations of integrals*. SIAM, 2001.

A Plain Gröbner basis analysis over $\mathbb{F}_2$

Let us define the Cauchy integrals

$$I_d(n) \stackrel{\text{def}}{=} \frac{1}{2\pi i} \oint \frac{1}{z^{d+1}} \frac{(1 + (l-1)z)^w}{(1-z)(1+z^2)^m} dz$$

for all non-negative integers $d$. Each coefficient of the Hilbert series can be written as the corresponding integral

$$\mathrm{HF}_{R/\mathcal{I}}(d) = I_d(n).$$

We set $f(z)$ such that $e^{nf(z)} = \frac{1}{z^{d+1}} \frac{(1+(l-1)z)^w}{(1-z)(1+z^2)^m}$, i.e.

$$f(z) \stackrel{\text{def}}{=} \frac{w}{n} \ln(1 + (l - 1)z) - \frac{m}{n} \ln(1 + z^2) - \frac{1}{n} \ln(1 - z) - \frac{d+1}{n} \ln(z).$$

We compute the saddle points as

$$\begin{aligned} z f'(z) &= \frac{(l-1)w}{n} \frac{z}{1+(l-1)z} - \frac{2m}{n} \frac{z^2}{1+z^2} + \frac{1}{n} \frac{z}{1-z} - \frac{d+1}{n} \\ &= \frac{(l-1)}{l} \frac{z}{1+(l-1)z} - 2\mu \frac{z^2}{1+z^2} + \frac{1}{n} \frac{z}{1-z} - \frac{d+1}{n}. \end{aligned}$$

We observe that the term $\frac{1}{n} \frac{z}{1-z}$ is negligible when $n \to \infty$. Hence, the computation leads to the condition

$$\frac{(l-1)}{l} z(1+z^2) - 2\mu z^2(1+(l-1)z) - \frac{d+1}{n}(1+z^2)(1+(l-1)z) = 0$$
and we want to find the value $\delta \stackrel{\text{def}}{=} \frac{d+1}{n}$ which provides the asymptotic of the relative degree of regularity. The cubic equation in $z$ above can be rewritten in the canonical form

$$P(z) = p_3 z^3 + p_2 z^2 + p_1 z + p_0 = 0$$

where

$$\left\{
\begin{aligned}
p_3 &= (l-1)(2\mu + \delta - \frac{1}{l}) \\
p_2 &= 2\mu + \delta \\
p_1 &= (l-1)(\delta - \frac{1}{l}) \\
p_0 &= \delta.
\end{aligned}
\right.$$

We want to find coalescent saddle points, i.e. a root of multiplicity 2 for $P$. For cubic polynomials, this can be found by forcing the discriminant to vanish. The discriminant is defined as

$$\Delta(P) = 18p_3p_2p_1p_0 - 4p_2^3p_0 + p_2^2p_1^2 - 4p_3p_1^3 - 27p_3^2p_0^2.$$

Equating the discriminant to 0 (and dividing by -4) provides a quartic equation

$$r_4\delta^4 + r_3\delta^3 + r_2\delta^2 + r_1\delta + r_0 = 0 \qquad (16)$$

in the variable $\delta$, whose coefficients are polynomial functions of $\mu$ and $l$:

$$\left\{
\begin{aligned}
r_4 &= ((l-1)^2 + 1)^2 \\
r_3 &= 2\mu((l-1)^2 + 1)((l-1)^2 + 3) - 4\frac{(l-1)^2((l-1)^2+1)}{l} \\
r_2 &= 4\mu^2(2(l-1)^2 + 3) - 2\mu\frac{(l-1)^2(3(l-1)^2-1)}{l} + 2\frac{(l-1)^2(3(l-1)^2+1)}{l^2} \\
r_1 &= 8\mu^3 + 20\mu^2\frac{(l-1)^2}{l} + 2\mu(l-1)^2\frac{3(l-1)^2-5}{l^2} - 4\frac{(l-1)^4}{l^3} \\
r_0 &= -\mu^2\frac{(l-1)^2}{l^2} - 2\mu\frac{(l-1)^4}{l^3} + \frac{(l-1)^4}{l^4}.
\end{aligned}
\right.$$

The smallest positive root $\bar{\delta}$ of (16) gives the good asymptotic for the relative degree of regularity.

## B Plain Gröbner basis approach over $\mathbb{F}_q$

We can write the Hilbert function at degree $d$ as

$$\mathrm{HF}_{R/\mathcal{I}}(d) = I_d(n)$$

for the Cauchy integrals

$$I_d(n) \stackrel{\text{def}}{=} \frac{1}{2\pi i} \oint \frac{1}{z^{d+1}} \frac{(1 + (l-1)z)^w (1-z^2)^m}{(1-z)^{w+1}} dz.$$

We set $f(z)$ such that $e^{nf(z)} = \frac{1}{z^{d+1}} \frac{(1+(l-1)z)^w (1-z^2)^m}{(1-z)^{w+1}}$, i.e.

$$f(z) \stackrel{\text{def}}{=} \frac{w}{n} \ln(1 + (l-1)z) - \frac{m}{n} \ln(1-z^2) - \frac{w+1}{n} \ln(1-z) - \frac{d+1}{n} \ln(z).$$

We compute the saddle points as

$$\begin{aligned}
z f'(z) &= \frac{(l-1)w}{n} \frac{z}{1+(l-1)z} - \frac{2m}{n} \frac{z^2}{1-z^2} + \frac{w+1}{n} \frac{z}{1-z} - \frac{d+1}{n} \\
&= \frac{(l-1)}{l} \frac{z}{1+(l-1)z} - 2\mu \frac{z^2}{1-z^2} + \frac{1}{l} \frac{z}{1-z} + \frac{1}{n} \frac{z}{1-z} - \frac{d+1}{n}.
\end{aligned}$$

We observe that the term $\frac{1}{n} \frac{z}{1-z}$ is negligible when $n \to \infty$. Hence, the computation leads to the condition

$$\frac{(l-1)}{l} z(1-z^2) - 2\mu z^2(1+(l-1)z) + \frac{1}{l} z(1+z)(1+(l-1)z) - \frac{d+1}{n}(1-z^2)(1+(l-1)z) = 0$$

and we want to find the value $\delta \stackrel{\text{def}}{=} \frac{d+1}{n}$ which provides the asymptotic of the relative degree of regularity. The cubic equation in $z$ above can be rewritten in the canonical form

$$P(z) = p_3 z^3 + p_2 z^2 + p_1 z + p_0 = 0$$
where

$$
\begin{cases}
p_3 = (l-1)(2\mu - \delta) \\
p_2 = 2\mu - \delta - 1 \\
p_1 = \delta(l-1) - 1 \\
p_0 = \delta.
\end{cases}
$$

We want to find coalescent saddle points, i.e. a root of multiplicity 2 for P. For cubic polynomials, this can be found by forcing the discriminant to vanish. The discriminant is defined as

$$
\Delta(P) = 18p_3p_2p_1p_0 - 4p_2^3p_0 + p_2^2p_1^2 - 4p_3p_1^3 - 27p_3^2p_0^2.
$$

Equating the discriminant to 0 provides a quartic equation

$$
r_4 \delta^4 + r_3 \delta^3 + r_2 \delta^2 + r_1 \delta + r_0 = 0 \quad (17)
$$

in the variable $\delta$, whose coefficients are polynomial functions of $\mu$ and $l$:

$$
\begin{cases}
r_4 = 4l^2(2l+1)^2 \\
r_3 = -8\mu(l^3 - 4l^2 - 2l + 4) - 4(3l^2 - 8l + 8)(l-2) \\
r_2 = -16\mu^2(2l^2 - 4l - 1) + 8\mu(3l^3 - 14l^2 + 29l - 24) + 13l^2 - 48l + 48 \\
r_1 = -32\mu^3 - 16\mu^2(5l-8) - 4\mu(6l^2 - 23l + 24) - 6(l-2) \\
r_0 = 4\mu^2 + 4\mu(2l-3) + 1.
\end{cases}
$$

The smallest positive root $\bar{\delta}$ of (17) gives the good asymptotic for the relative degree of regularity.