# Contributions to polynomial system solving:
## Recurrences and Gröbner bases

Jérémy Berthomieu

► To cite this version:

Jérémy Berthomieu. Contributions to polynomial system solving: Recurrences and Gröbner bases. Symbolic Computation [cs.SC]. Sorbonne Université, 2023. tel-04289532

HAL Id: tel-04289532

https://hal.sorbonne-universite.fr/tel-04289532

Submitted on 16 Nov 2023

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Sorbonne Université

Mémoire d'Habilitation à Diriger des Recherches

Mention Informatique

# Contributions to polynomial system solving:
Recurrences and Gröbner bases

Par : Jérémy Berthomieu

## RAPPORTEURS :

M. BERNARD MOURRAIN,

M. CORDIAN RIENER,

M. GILLES VILLARD,

Directeur de recherche,

Professeur,

Directeur de recherche,

INRIA & Université Côte d'Azur

Universitetet i Tromsø – Norges arktiske universitet

CNRS & École Normale Supérieure de Lyon

## EXAMINATEURS :

M. ALIN BOSTAN,

M. MANUEL KAUERS,

MME FATEMEH MOHAMMADI,

M. MOHAB SAFEY EL DIN,

Directeur de recherche,

Professeur,

Professeure,

Professeur,

INRIA

Johannes Kepler Universität Linz

Katholieke Universiteit Leuven, Présidente

Sorbonne Université

Date de soutenance : 21 septembre 2023
This work would not have been possible without the support of many people.

I would like first to thank Mohab SAFEY EL DIN, for giving me so much advice which helped me a lot and of course to be in the jury of this habilitation thesis. I also want to thank all my past and current colleagues in the POLSYS group with whom it is a pleasure to collaborate and work: Christian EDER, Jean-Charles FAUGÈRE, Vincent NEIGER, Ludovic PERRET and Guénaël RENAULT.

Many thanks to Bernard MOURRAIN, Cordian RIENER and Gilles VILLARD for reviewing this manuscript!

I would also like to thank Alin BOSTAN for accepting to be in this jury and his advice and support together with the other members of this jury: Manuel KAUERS and Fatemeh MOHAMMADI.

I am lucky to be part of a very nice work environment. Thank you to all the other members of the POLSYS group in these last few years: Lorenzo BALDI, Olive CHAKRABORTY, Andrew FERGUSON, Jorge GARCÍA FONTÁN, Sriram GOPALAKRISHNAN, Dimitri LESNOFF, Rafael MOHR, Hadrien NOTARANTONIO, Pierre PÉBEREAU, Rémi PRÉBET, Georgy SCHOLTEN, Thi Xuan Vu, Trung-Hieu Vu. And thank you to all the former and current members of the PEQUAN team, in particular Stef GRAILLAT for all his support and the fun talking about our children, and also Dominique BÉRÉZIAT, Pierre FORTIN, Fabienne JÉZÉQUEL, Thibault HILAIRE Christoph LAUTER, Théo MARY, Valérie MÉNISSIER-MORAIN and Marc MEZZAROBBA.

Thank you to all the administrative staff of LIP6 for their help.

A little bit further, thank you to the members of MATHEXP and in particular Frédéric CHYZAK and Pierre LAIREZ.

More generally, thank you the whole Computer Algebra community, especially the members of the *GT Calcul Formel*.

This works has been partly supported by joint ANR-FWF ANR-19-CE48-0015 ECARP and ANR-22-CE91-0007 EAGLES projects, the ANR grants ANR-18-CE33-0011 SESAME and ANR-19-CE40-0018 DE RERUM NATURA projects, European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement N. 813211 (POEMA), grant DIM-RFSI 2021-02-C21/1131 of the Paris Île-de-France Region, grant FA8665-20-1-7029 of the EOARD-AFOSR and PGMO grant CAMISADo.

On a personal level, I would like to thank my friends, my family in law, Amélie, Benjamin and Nathalie, and my family, Arnaud, Géraldine, Juliette, Joséphine, Chloé, Thibault, Antoine, Baptiste and my parents for their constant support.

Finally, my deepest thank and love for everything to Camille and our two beautiful sons, Tristan and Paul.
<table>
  <tbody>
    <tr>
      <td>1</td>
      <td>Introduction</td>
      <td>1</td>
    </tr>
    <tr>
      <td>1.1</td>
      <td>Gröbner bases: a fundamental object</td>
      <td>1</td>
    </tr>
    <tr>
      <td>1.2</td>
      <td>Solving polynomial systems</td>
      <td>2</td>
    </tr>
    <tr>
      <td>1.3</td>
      <td>Sequences</td>
      <td>2</td>
    </tr>
    <tr>
      <td>1.4</td>
      <td>Contributions</td>
      <td>4</td>
    </tr>
    <tr>
      <td>1.4.1</td>
      <td>Problem 1: Complexity of SPARSE-FGLM for generic determinantal systems</td>
      <td>4</td>
    </tr>
    <tr>
      <td>1.4.2</td>
      <td>Problem 2: Complexity of Gröbner bases change of order</td>
      <td>5</td>
    </tr>
    <tr>
      <td>1.4.3</td>
      <td>Problem 3: Computation of Gröbner bases of saturated ideals</td>
      <td>5</td>
    </tr>
    <tr>
      <td>1.4.4</td>
      <td>Problem 4: Polynomial arithmetic for guessing C-relations</td>
      <td>6</td>
    </tr>
    <tr>
      <td>1.4.5</td>
      <td>Problem 5: Minimization of number of queries</td>
      <td>6</td>
    </tr>
    <tr>
      <td>1.4.6</td>
      <td>Problem 6: Structured guessing</td>
      <td>7</td>
    </tr>
    <tr>
      <td>1.4.7</td>
      <td>MSOLVE</td>
      <td>7</td>
    </tr>
    <tr>
      <td>1.5</td>
      <td>Organization</td>
      <td>7</td>
    </tr>
    <tr>
      <td>2</td>
      <td>MSOLVE</td>
      <td>9</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>State-of-the-art algorithms</td>
      <td>9</td>
    </tr>
    <tr>
      <td>2.1.1</td>
      <td>The F<sub>4</sub> algorithm</td>
      <td>9</td>
    </tr>
    <tr>
      <td>2.1.2</td>
      <td>The SPARSE-FGLM algorithm</td>
      <td>11</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Experimental results</td>
      <td>14</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Computing Gröbner bases</td>
      <td>17</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Gröbner bases change of order algorithms</td>
      <td>17</td>
    </tr>
    <tr>
      <td>3.1.1</td>
      <td>Complexity of SPARSE-FGLM</td>
      <td>17</td>
    </tr>
    <tr>
      <td>3.1.2</td>
      <td>A polynomial-matrix algorithm</td>
      <td>19</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Saturated and colon ideals</td>
      <td>20</td>
    </tr>
    <tr>
      <td>3.2.1</td>
      <td>Bayer's algorithm</td>
      <td>21</td>
    </tr>
    <tr>
      <td>3.2.2</td>
      <td>F<sub>4</sub>SAT, an algorithm for saturated ideals</td>
      <td>21</td>
    </tr>
    <tr>
      <td>3.2.3</td>
      <td>SPARSE-FGLM-COLON, an algorithm for colon ideals</td>
      <td>23</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Parametrizations of the radical</td>
      <td>26</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Guessing Gröbner bases</td>
      <td>29</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Uni-indexed sequences and the Berlekamp-Massey algorithm</td>
      <td>29</td>
    </tr>
    <tr>
      <td>4.1.1</td>
      <td>Matrix viewpoints</td>
      <td>29</td>
    </tr>
    <tr>
      <td>4.1.2</td>
      <td>Polynomial viewpoint</td>
      <td>30</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Extensions to multi-indexed sequences</td>
      <td>32</td>
    </tr>
    <tr>
      <td>4.2.1</td>
      <td>The BMS algorithm</td>
      <td>32</td>
    </tr>
    <tr>
      <td>4.2.2</td>
      <td>The SCALAR-FGLM algorithm</td>
      <td>33</td>
    </tr>
    <tr>
      <td>4.2.3</td>
      <td>Adaptive algorithms</td>
      <td>34</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Benchmarks</td>
      <td>35</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Quasi-commutative Gröbner bases</td>
      <td>37</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>P-relations and quasi-commutative polynomials</td>
      <td>37</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Guessing P-relations</td>
      <td>38</td>
    </tr>
    <tr>
      <td>5.2.1</td>
      <td>A linear-algebra-based algorithm</td>
      <td>38</td>
    </tr>
    <tr>
      <td>5.2.2</td>
      <td>An hybrid approach</td>
      <td>38</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Structured sequences</td>
      <td>39</td>
    </tr>
    <tr>
      <td>6</td>
      <td>Research project</td>
      <td>41</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>MSOLVE and fundamental algorithms</td>
      <td>41</td>
    </tr>
    <tr>
      <td>6.1.1</td>
      <td>Types for the coefficients</td>
      <td>41</td>
    </tr>
  </tbody>
</table>
# Contents

## Contents

<table><tbody><tr><td>6.1.2</td><td>Gröbner bases for a degree order</td><td>41</td></tr><tr><td>6.1.3</td><td>Change of order</td><td>42</td></tr><tr><td>6.1.4</td><td>Saturation and colon ideals</td><td>42</td></tr><tr><td>6.1.5</td><td>Exact solving on a GPU</td><td>42</td></tr><tr><td>6.2</td><td>Faster change of order and guessing C-relations</td><td>43</td></tr><tr><td>6.2.1</td><td>Guessing faster</td><td>43</td></tr><tr><td>6.2.2</td><td>Guessing radical ideals</td><td>44</td></tr><tr><td>6.2.3</td><td>Guessing with multiplicities</td><td>44</td></tr><tr><td>6.2.4</td><td>Polynomial matrices with multiplicities</td><td>44</td></tr><tr><td>6.2.5</td><td>Applications to guessing</td><td>45</td></tr><tr><td>6.3</td><td>Algorithms for quasi-commutative Gröbner bases</td><td>45</td></tr><tr><td>6.3.1</td><td>Moment approach</td><td>45</td></tr><tr><td>6.3.2</td><td>Efficient Gröbner bases computation</td><td>45</td></tr><tr><td>A</td><td>Ideals of C-relations of small degrees</td><td>47</td></tr><tr><td>A.1</td><td>Uni-indexed sequences</td><td>47</td></tr><tr><td>A.2</td><td>Bi-indexed sequences</td><td>47</td></tr><tr><td>A.2.1</td><td>Single roots</td><td>47</td></tr><tr><td>A.2.2</td><td>Root of multiplicity 2</td><td>48</td></tr><tr><td>A.2.3</td><td>Root of multiplicity 3</td><td>48</td></tr></tbody></table>
This habilitation thesis deals with the research that I have done since I have been with the POЛSYS group in Sorbonne Université.

On the one hand, polynomial systems arise in a wide range of areas of scientific domains such as biology [39], chemistry [62], quantum mechanics [85], robotics [91], and computing sciences, including coding theory [94], computer vision [72] and cryptography [51] to cite a few. On the other hand, polynomial system solving is NP-hard, even when the ground field is finite [61, Appendix A7.2]. Moreover, the non-linearity of such systems make reliability issues topical, in particular when complete and exhaustive outputs are required, in the context of numerical algorithms.

Likewise, sequences are a classical mathematical object and computing linear recurrence relations of a multi-indexed sequence or determining the nature of this sequence based on these relations is a fundamental problem in coding theory [65, 94], computer algebra [52, 98, 100] and enumerative combinatorics [28, 30, 31].

Whether it be for solving polynomial systems or for computing or *guessing* linear recurrence relations, one aims to obtain nice generators of an ideal that are able to answer the following questions. Is the number of solutions finite in an algebraic closure of the field of coefficients? How many initial terms and linear recurrence relations do one needs to compute any term of the sequence?

These questions are easily answered when we have a *Gröbner basis* of the ideal at hand. This is why, the main focus of my research is Gröbner bases computations with two main goals: solving polynomial systems exactly and determining all the linear recurrence relations of a multi-indexed sequence.

## 1.1 Gröbner bases: a fundamental object

Buchberger developed the theory of Gröbner bases and designed a first algorithm to compute them in [35]. Since then, many efficient Gröbner basis algorithms were developed.

While lexicographic Gröbner bases are the tool of choice to represent the solution set of a polynomial system, often, they are the hardest Gröbner bases to compute. For *n* generic polynomials of degree *d* in *n* variables, computing the <sub>LEX</sub>-Gröbner basis of the ideal they spanned is bounded by $C_1 d^{C_2 n^3}$, see [37], whereas computing the <sub>DRL</sub>-Gröbner basis of the same ideal is bounded by $C_1 d^{C_2 n^2}$, see [74].

As a caveat to that, Gröbner bases change of orders algorithms have been introduced. They take as an input a Gröbner basis $G_1$ for a monomial order <1 and another monomial order <2 and they return $G_2$, a Gröbner basis of $\langle G_1 \rangle$ for <2. This yields the following framework, in the zero-dimensional case, where $f_1, \dots, f_s$ are the original polynomials that are given as an input to Buchberger's algorithm [35] or to Faugère's $F_4$ [47] or $F_5$ [48] algorithms to compute the <sub>DRL</sub>-Gröbner basis $G_{DRL}$. Then, $G_{DRL}$ is converted into the <sub>LEX</sub>-Gröbner basis $G_{LEX}$ using the so-called FGLM algorithm [50] or faster variants like the SPARSE-FGLM algorithm [52, 53] or more recently [83] and [20]. This framework allows one to compute the <sub>LEX</sub>-Gröbner basis in $C_1 d^{C_2 n^2}$ operations as well.

<figure><img src="image_2.png" /><figcaption>Figure 1.1: Classical Gröbner bases framework for 0-dimensional ideals</figcaption></figure>

Nowadays, many computer algebra software propose implementations of algorithms for computing Gröbner bases, including $F_4$, FGLM and/or SPARSE-FGLM. We can cite MAGMA [24], MAPLE [9],
SINGULAR [42], as general software, but also some efficient C library like FGB [49] and MSOLVE [13].

### 1.2 Solving polynomial systems

Let $\mathbb{K}$ be an algebraic closed field and $n \in \mathbb{N}^*$. We let $\mathbf{x} = (x_1, \dots, x_n)$. We recall that an ideal $I \subseteq \mathbb{K}[\mathbf{x}]$ is zero-dimensional if the solution set $\{\mathbf{a} \in \mathbb{K}^n | \forall f \in I, f(\mathbf{a}) = 0\}$ is finite.

With the convention that $x_n <_{\text{LEX}} \dots <_{\text{LEX}} x_1$, a $<_{\text{LEX}}$-Gröbner basis $\mathcal{G}_{\text{LEX}}$ of $I$ allows one to compute the coordinates of the above solution set. Indeed, for all $1 \le k \le n$, $\mathcal{G}_{\text{LEX}} \cap \mathbb{K}[x_k, \dots, x_n]$ is a $<_{\text{LEX}}$-Gröbner basis of $I \cap \mathbb{K}[x_k, \dots, x_n]$. That is, $<_{\text{LEX}}$-Gröbner bases extend the notion of triangular basis for linear systems.

Assuming $\mathcal{G}_{\text{LEX}}$ is reduced, it suffices then to compute the roots of the (unique) polynomial in $\mathcal{G}_{\text{LEX}} \cap \mathbb{K}[x_n]$ to obtain the last coordinate of each solution. Then, to replace $x_n$ by each of this last coordinate in all the polynomials in $\mathcal{G}_{\text{LEX}} \cap \mathbb{K}[x_{n-1}, x_n]$ to find the corresponding second-to-last coordinate of each solution, and so on and so forth.

Furthermore, if $I$ is radical, then after a generic change of coordinates, its reduced $<_{\text{LEX}}$-Gröbner basis is said to be in *shape position*, i.e.

$$ \mathcal{G}_{\text{LEX}} = \{g_n(x_n), x_{n-1} - g_{n-1}(x_n), \dots, x_1 - g_1(x_n)\}, $$

with $\deg g_n = \deg I$ and for all $1 \le k \le n-1$, $\deg g_k < \deg g_n$. In other words, the last coordinate parameterizes the other coordinates of the solutions.

Thanks to this, computing a $<_{\text{LEX}}$-Gröbner basis of a zero-dimensional ideal $\langle f_1, \dots, f_s \rangle$ is a fundamental step for solving exactly a polynomial system $f_1 = \dots = f_s = 0$. Using the framework of Figure 1.1, we are led to first compute a $<_{\text{DRL}}$-Gröbner basis of this ideal.

Under generic assumptions, the complexity of Gröbner bases change of order from $<_{\text{DRL}}$ to $<_{\text{LEX}}$ is in $\tilde{O}(D^\omega)$ [83] and $O(tD^2)$ [53], where $D$ is the degree of the ideal and $t \le D$ is the number of monomials not in $x_n$ in the $<_{\text{DRL}}$-monomial basis. This parameter $t$ is well-understood when the ideal is spanned by $n$ generic polynomials of degree $d$ [53, Cor. 5.10].

In many applications, the polynomial system at hands is structured. For instance, in polynomial optimization, critical values of a polynomial map restricted to an algebraic sets are solutions of determinantal systems deriving from maximal minors.

**Problem 1 (Complexity of Sparse-FGLM for generic determinantal systems):** Is it possible to give asymptotics for $t$ in generic cases for some of these structured systems?

Furthermore, depending on the asymptotics on $t$ and the theoretical or practical values of $\omega$, it is not clear which complexity is the better one between $\tilde{O}(D^\omega)$ [83] and $O(tD^2)$ [53].

**Problem 2 (Complexity of Gröbner bases change of order):** Under the same generic assumptions, is there an algorithm whose complexity is better than [53, 83] independently on the ratio $t/D$ and on the chosen bound for $\omega$?

In real algebraic geometry, computing sample points in singular real algebraic sets [92] or computing their real dimension [73] require to compute limits of critical points of the restriction of a polynomial map to some algebraic set depending on a parameter. This boils down to the computation of saturated ideals. State-of-the-art algorithms for a saturated ideal $I : \langle\varphi\rangle^\infty$ [4, 90] rely on the Gröbner basis computation of an augmented ideal containing $I$, introducing an extra variable that must be eliminated later. Furthermore, running the $F_4$ algorithm on this kind of ideal, we can observe that the computations are not regular: there a lot of degree falls.

**Problem 3 (Computation of Gröbner bases of saturated ideals):** Is it possible to compute a Gröbner basis of $I : \langle\varphi\rangle^\infty$ on the fly from the generators of $I$ without introducing an extra variable?

### 1.3 Sequences

Let $\mathbf{u} \in \mathbb{K}^{\mathbb{N}^n}$ be a $n$-indexed sequence with values in $\mathbb{K}$, that is $\mathbf{u} = (u_{i_1, \dots, i_n})_{(i_1, \dots, i_n) \in \mathbb{N}^n}$. There is a natural correspondence between finite linear combinations of terms of $\mathbf{u}$ and polynomials in $\mathbb{K}[x_1, \dots, x_n] = \mathbb{K}[\mathbf{x}]$. For $g = \sum_{s \in \mathcal{S}} \gamma_s \mathbf{x}^s$, with $\mathcal{S}$ a finite subset of $\mathbb{N}^n$, we write $\text{Eval}(g, \mathbf{u}) :=$
$\sum_{s \in \delta} \gamma_s u_s$. Since

$$ \mathrm{Eval}(gx^i, u) = \sum_{s \in \delta} \gamma_s u_{s+i}, $$

shifting a linear combination of terms by an index $i$ comes down to multiplying the associated polynomial by $x^i$. Observe that Eval is a bilinear operator from $K[x] \times K^{\mathbb{N}^n}$ to $K$.

**Definition-Proposition 1.1 (e.g. [12, Def. 2 and Prop. 4]):** Let $u \in K^{\mathbb{N}^n}$ be a sequence. A polynomial $g \in K[x]$ defines a linear recurrence relation with constant coefficients, or C-relation for short, on $u$ if, and only if, for all $i \in \mathbb{N}^n$, Eval $(gx^i, u) = 0$.

The set of all such polynomials is an ideal of $K[x]$ called the *ideal of C-relations of u* and denoted $I_C(u)$.

**Definition-Proposition 1.2 (e.g. [11, Def. 2 and Prop. 3]):** A sequence $u \in K^{\mathbb{N}^n}$ is said to be *C-finite* if together with a finite number of terms of $u$ and a finite number of C-relations, one can recover all the terms of $u$. This is equivalent to requiring $I_C(u)$ be 0-dimensional.

In some applications, like in combinatorics, error correcting code, but also in computer algebra through polynomial system solving, one must deal with the ideal of C-relations of a $n$-indexed sequence. In the last two applications, the <LEX>-Gröbner basis of this ideal is the target and its roots the ultimate goal.

It is unclear, at this stage, if given an ideal $I \in K[x]$, we can build a sequence $u \in K^{\mathbb{N}^n}$ such that $I_C(u) = I$. However, we have this following easy result.

**Definition-Proposition 1.3:** If $u$ and $v$ in $K^{\mathbb{N}^n}$ are such that $I_C(u) \subseteq I$ and $I_C(v) \subseteq I$, then for any $\lambda \in K$, $I_C(u + \lambda v) \subseteq I$, where $u + \lambda v = (u_i + \lambda v_i)_{i \in \mathbb{N}^n}$. Hence, for an ideal $I$ of $K[x]$, there exists a whole vector subspace of $K^{\mathbb{N}^n}$, denoted $S_C(I)$ such that for $u \in S_C(I)$, $I_C(u) \subseteq I$.

**Remark 1.4:** By convention, the *zero sequence* is the only sequence whose ideal of C-relations is $\langle 1 \rangle$.

**Theorem 1.5:** Let $u$ and $v$ in $K^{\mathbb{N}^n}$. Let $\lambda, \mu \in K$. Then,

$$ I_C(u) \cap I_C(v) \subseteq I_C(\lambda u + \mu v). $$

*Proof.* Let $g = \sum_{s \in \delta} \gamma_s x^s \in I_C(u) \cap I_C(v)$. By linearity, for all $i \in \mathbb{N}^n$, Eval $(gx^i, \lambda u + \mu v) = \lambda$ Eval $(gx^i, u) + \mu$ Eval $(gx^i, v) = 0$. Hence, $g \in I_C(\lambda u + \mu v)$. $\square$

In this manuscript, most examples will come from sequences such as the following one.

**Example 1.6:** Let $I = \langle (x_2^4 + x_2^3)(x_1^2 + x_2^2 - 1), (x_1 - x_2^2 - x_2)(x_1^2 + x_2^2 - 1) \rangle$. It defines the unit circle with multiplicity 1, through the factor $x_1^2 + x_2^2 - 1$ appearing in both spanning polynomials and two points: the origin $(0,0)$, with multiplicity 3, and $(0,-1)$, with multiplicity 1, as the common vanishing points of $x_2^4 + x_2^3 = x_1 - x_2^2 - x_2 = 0$. Notice that $(0,-1)$ is also on the unit circle and thus it has multiplicity 2 in total.

Let $\varphi = x_1^4 - x_2^4 - 2x_1^2 + 1 = (x_1^2 + x_2^2 - 1)(x_1^2 - x_2^2 - 1)$ defining the union of the unit circle and the unit hyperbola.

The ideal $J = I : \langle \varphi \rangle^{\infty}$ defines the Zariski closure of the set difference of the variety defined by $I$ and that defined by $\varphi$. The resulting ideal defines the origin, with multiplicity 3, and one can prove that $J = \langle x_2^3, x_1 - x_2^2 - x_2 \rangle$.

Now, let $v \in S_C(J)$. By the first polynomial, we have Eval $(x_1^i x_2^{j+3}, v) = v_{i,j+3} = 0$ for all $(i, j) \in \mathbb{N}^2$. Moreover, by the second polynomial, we have $v_{i+1,j} = \mathrm{Eval}(x_1^{i+1} x_2^j, v) = \mathrm{Eval}(x_1^i x_2^{j+2} + x_1^i x_2^{j+1}, v) = v_{i,j+2} + v_{i,j+1}$ for all $(i, j) \in \mathbb{N}^2$. Hence, there exist $\alpha_{0,0}, \alpha_{0,1}$ and $\alpha_{0,2}$ such that

<table><tbody><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>3</td><td>0</td><td>0</td><td>0</td><td>0</td><td>...</td></tr><tr><td>2</td><td>α<sub>0,2</sub></td><td>0</td><td>0</td><td>0</td><td>...</td></tr><tr><td>1</td><td>α<sub>0,1</sub></td><td>α<sub>0,2</sub></td><td>0</td><td>0</td><td>...</td></tr><tr><td>0</td><td>α<sub>0,0</sub></td><td>α<sub>0,1</sub> + α<sub>0,2</sub></td><td>α<sub>0,2</sub></td><td>0</td><td>...</td></tr><tr><td>v</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>
Let us recall that not all ideals can be ideals of C-relations. For instance, any sequence $u$, built from $K = \langle x_1^2, x_1x_2, x_2^2 \rangle$, satisfies the relation Eval $((u_{0,1}x_1 - u_{1,0}x_2)x_1^i x_2^j, u) = u_{0,1}u_{i+1,j} - u_{1,0}u_{i,j+1} = 0$ for all $(i, j) \in \mathbb{N}^2$. We refer to Appendix A for more details. Hence $K \subsetneq I_C(u)$, and $K$ is not an ideal of C-relations. More generally, [32, Prop. 3.3] proves that ideals of C-relations are exactly Gorenstein ideals [63, 76] and problems occur only if the ideal has a zero of multiplicity at least 2. The following theorem can also be found in [45, Th. 8.3].

**Theorem 1.7:** Let $I \subseteq \mathbb{K}[x]$ be a 0-dimensional ideal. The ideal $I$ (resp. ring $\mathcal{R} = \mathbb{K}[x]/I$) is Gorenstein if, equivalently,

(1) $\mathcal{R}$ and its dual are isomorphic as $\mathcal{R}$-modules;

(2) there exists a $\mathbb{K}$-linear form $\tau$ on $\mathcal{R}$ such that the following bilinear form is non degenerate

$$ \mathcal{R} \times \mathcal{R} \to \mathbb{K} \\ (a, b) \mapsto \tau(ab). $$

In order to be self-contained, we detail in Appendix A the link between an ideal of C-relations $I$ and the general term of a sequence in $S_C(I)$.

Guessing the ideal of C-relations of a sequence $u$ is the task of computing candidate polynomials in $I_C(u)$ using only finitely many terms of $u$. In particular, we aim at guessing a $\langle-$Gröbner basis of $I_C(u)$ for a given monomial order $<$. In the uni-indexed sequence case, guessing the essentially unique C-relation can be modeled through the solving a Hankel system. Yet, the best algorithm, which runs in time quasi-linear in the number of input sequence terms, uses fast polynomial arithmetic and fast extended Euclidean algorithm [34].

In the multi-indexed setting, one can still guess the C-relations using linear algebra techniques on a multi-Hankel matrix [11, 12] or a Gram–Schmidt process [81]. The so-called Berlekamp–Massey–Sakata algorithm, due to Sakata [93, 95, 97], guesses the relations using polynomial additions and shifts by a monomial.

**Problem 4 (Polynomial arithmetic for guessing C-relations):** Can we model the guessing of C-relations of multi-indexed sequences with multivariate polynomial arithmetic, including polynomial multiplications and divisions?

In some applications, for instance the polynomial system solving one through the SPARSE-FGLM algorithm, the bottleneck of the algorithm is the computation of the sequence terms, and not the guessing of the C-relations. It is thus essential that only terms strictly needed for the guessing be computed.

**Problem 5 (Minimization of number of queries):** In the case where queries to the sequence are expensive, can we relate the number of queries to the geometry of the staircase of the output Gröbner basis?

In enumerative combinatorics, like 2D/3D-walks, the sequence comes with a structure: the terms are invariant by the action of a group on the indices, the nonzero terms lie on a cone. In particular, when many sequence terms are, this makes the previously mentioned multi-Hankel matrix have a lot of zero rows. This translates into guessed relations that might be false as we would not guess them with more terms at hand. Thus, we need to build a much larger matrix, for instance with many more rows than columns, to prevent guessing these false relations.

**Problem 6 (Structured guessing):** Can we exploit the sparsity of the nonzero terms in order to guess fewer fake relations with fewer queries? What properties do the guess relations have?

## 1.4 Contributions

### 1.4.1 Problem 1: Complexity of SPARSE-FGLM for generic determinantal systems

In polynomial optimization, critical points methods require to compute the vanishing set of the restriction of a linear map to an algebraic set defined by $p$ generic polynomials by means of the simultaneous vanishing of maximal minors of a truncated Jacobian matrix. More generally, we
define below the class of generic determinantal sum ideals that contains this kind of ideals and we provide an asymptotic on the parameter $t$ for them.

**Definition 1.8 (see also [10, Def. 8]):** With $K$ an infinite field, let $I \subset K[x]$ be an ideal which is the sum of $m$ polynomials of degree at most $d$ and the maximal minors of a matrix with polynomial entries also of degree at most $d$. We say that $I$ is a generic determinantal sum ideal if the following three conditions hold:

* the ideal $I$ is zero-dimensional and in shape position;

* the Hilbert series of $K[x]/I$ is

$$H(\tau) = \frac{\det(M(\tau^{d-1})) (1 - \tau^d)^m (1 - \tau^{d-1})^{n-m}}{\tau^{(d-1)(m-1)} (1 - \tau)^n}$$

where $M(\tau)$ is the $(m-1) \times (m-1)$-matrix whose $(i, j)$th entry is $\sum_k \binom{m-i}{k} \binom{n-1-j}{k} \tau^k$;

* for all $e \ge 1$, the Hilbert series of $(K[x]/I) / \langle x_n^e \rangle$ is equal to the series $(1 - \tau)H$ truncated at the first non-positive coefficient.

**Theorem 1.9 (see also [10, Th. 2]):** Let $I$ be a generic determinantal sum ideal of $K[x]$ as in Definition 1.8. Let $\mathcal{G}_{\text{DRL}}$ be the reduced $<_{\text{DRL}}$-Gröbner basis of $I$. Let $t$ be the number of polynomials in $\mathcal{G}_{\text{DRL}}$ whose leading monomial is divisible by $x_n$.

Then, for $d=2$ and $n \ge m$,

$$t = \sum_{k=0}^{m-1} \binom{n-m-1+k}{k} \binom{m}{\lfloor 3m/2 \rfloor - 1 - j}$$

Moreover, for $d \ge 3$ and $n \to \infty$,

$$t \approx \frac{1}{\sqrt{(n-m)\pi}} \sqrt{\frac{6}{(d-1)^2-1}} d^m (d-1)^{n-m} \binom{n-2}{m-1}$$

### 1.4.2 Problem 2: Complexity of Gröbner bases change of order

State-of-the-art FGLM algorithms uses multiplication matrices to compute the polynomial in the reduced Gröbner basis for the target monomial order. These are matrices of size the degree of the ideal and with scalar entries. Using a change of paradigm and relying on polynomial matrix arithmetic, we design a new change of order algorithm whose complexity is in $\tilde{O}(t^{\omega-1}D)$ under the same generic assumption as [53]. As a consequence, it is always asymptotically faster than both [53] and [83]. Furthermore, for $\omega$ close to 2, the complexity becomes quasi-linear in the size of the input.

**Theorem 1.10 (see also [20, Th. 1.1]):** Let $I$ be a zero-dimensional ideal $K[x]$ of degree $D$. Let $\mathcal{G}_{\text{DRL}}$ (resp. $\mathcal{G}_{\text{LEX}}$) be the reduced $<_{\text{DRL}}$ (resp. $<_{\text{LEX}}$) Gröbner basis of $I$ and $S_{\text{DRL}}$ be the $<_{\text{DRL}}$-monomial basis of $K[x]/I$. Assume that $x_1, \dots, x_{n-1}$ are in $S_{\text{DRL}}$, and that for all monomials $\mu \in S_{\text{DRL}}$, either $x_n\mu$ is in $S_{\text{DRL}}$ or it is the $<_{\text{DRL}}$-leading monomial of an element in $\mathcal{G}_{\text{DRL}}$. Assume that $I$ is in shape position. Then, one can compute $\mathcal{G}_{\text{LEX}}$ using $\tilde{O}(t^{\omega-1}D)$ operations in $K$, where $t$ is the number of elements of $\mathcal{G}_{\text{DRL}}$ whose $<_{\text{DRL}}$-leading monomial is divisible by $x_n$.

### 1.4.3 Problem 3: Computation of Gröbner bases of saturated ideals

Given generators $f_1, \dots, f_s$ of an ideal $I$ and a polynomial $\varphi$, we propose an algorithm extending Faugère's $F_4$ algorithm [47] that computes a partial $<_{\text{DRL}}$-Gröbner basis of $I$ and uses this partial information to find polynomials in $I: \langle \varphi \rangle$ not in $I$ using linear algebra. Then, adding these new generators, it resumes the partial $<_{\text{DRL}}$-Gröbner basis computation of $I: \langle \varphi \rangle$. Repeating this process, it finds polynomials in $(I: \langle \varphi \rangle): \langle \varphi \rangle$ and so on and so forth until it stabilizes and has computed the $<_{\text{DRL}}$-Gröbner basis of $I: \langle \varphi \rangle^\infty$. This is the $F_4$SAT algorithm, Algorithm 3.2.

Furthermore, we modify the SPARSE-FGLM algorithm so that it takes as an input the reduced $<_{\text{DRL}}$-Gröbner basis of an ideal $I$ and perform linear algebra operations in the set of multiples of $\varphi$ in $K[x]/I$ in order to compute the reduced $<_{\text{LEX}}$-Gröbner basis of the zero-dimensional ideal $I: \langle \varphi \rangle$. This allows us to perform in one step a colon ideal computation and a change of order. Moreover,
## 1.4. Contributions

this new Algorithm 3.3, SPARSE-FGLM-COLON, does not assume that $I$ is zero-dimensional, which finds application in polynomial optimization.

**Theorem 1.11 (see also [15, Th. 4.12]):** Let $I$ be a positive-dimensional ideal of $\mathbb{K}[\mathbf{x}]$, let $\mathscr{G}_{\text{DRL}}$ be its reduced $<_{\text{DRL}}$-Gröbner basis and $S_{\text{DRL}}$ be the associated staircase. Let $\varphi \in \mathbb{K}[\mathbf{x}] \setminus I$ such that $I : \langle \varphi \rangle$ is zero-dimensional of degree $D'$ and in shape position.

Let $\Sigma$ be a finite staircase of size $N$ containing supp NF $(x_n^i \varphi, \mathscr{G}_{\text{DRL}}, <_{\text{DRL}})$ for all $i \in \mathbb{N}$, where NF $(f, G, <)$ is the normal form of $f$ w.r.t. $G$ and $<$.

Let $t$ be the number of monomials $\sigma$ in $\Sigma$ such that $x_n\sigma \in \text{LM}_{<_{\text{DRL}}}(\mathscr{G}_{\text{DRL}})$ and let $u$ be the number of monomials $\sigma$ in $\Sigma$ such that $x_n\sigma \in (\text{LM}_{<_{\text{DRL}}}(I)) \setminus \text{LM}_{<_{\text{DRL}}}(\mathscr{G}_{\text{DRL}})$.

Then, for a generic choice of vector $r \in \mathbb{K}^N$, Algorithm 3.3 terminates and returns the reduced $<_{\text{LEX}}$-Gröbner basis of $I : \langle \varphi \rangle$. To do so, it requires at most $u+n$ normal form computations w.r.t. $\mathscr{G}_{\text{DRL}}$ and $<_{\text{DRL}}$ plus $O((t+u+n)ND')$ operations in $\mathbb{K}$.

As a byproduct, I designed and implemented in MSOLVE an efficient probabilistic test for the computation of the parametrizations of $x_1, \dots, x_{n-1}$ w.r.t. $x_n$ in $\sqrt{I}$. It relies on the algorithm of [29, 66] using $x_n$ as the linear form. Since their algorithm assumes that the linear form is generic, we cannot assert that $x_n$ parameterizes the other variables. If it does not, the computed parametrizations are meaningless. Hence the necessity for an efficient validity check of the output. Furthermore, the algorithm of [29, 66] builds a sequence before guessing the parametrizations and this sequence building is the actual bottleneck. Thus, the main advantage of this validity check is that it avoids building a brand new sequence. The following lemma gives the main idea.

**Lemma 1.12 (see also [15, Lem. 4.15]):** Let $I$ be a zero-dimensional ideal of $\mathbb{K}[\mathbf{x}]$. Let $\lambda \in \bar{\mathbb{K}}$ be generic. Then, for $1 \le k \le n$, $I = I : \langle x_k + \lambda \rangle$.

It suffices then to compute the parametrizations of such $I : \langle x_k + \lambda \rangle$ using the same sequence $u$ as for $I$. Indeed, $I : \langle x_k + \lambda \rangle$ is the ideal of C-relations of $v = (\text{Eval}((x_k + \lambda)\mathbf{x}^i, u))_{i \in \mathbb{N}^n}$.

### 1.4.4 Problem 4: Polynomial arithmetic for guessing C-relations

While the Berlekamp–Massey–Sakata algorithm [93, 95, 97] performs multivariate polynomials additions and multiplications by a monomial to compute a Gröbner basis of the ideals of C-relations, it fails to generalize the extended Euclidean algorithm. Using the mirror of the truncated generating series of the input multi-index sequence, we design an algorithm performing polynomial divisions and normal forms to compute potential C-relations. Furthermore, we exhibit a criterion based on leading monomials to ensure which of these potential C-relations are correct or not. Finally, when called on a uni-index sequence, the algorithm comes down to performing a truncated extended Euclidean algorithm. All in all, we have the following result.

**Theorem 1.13 (see also [19, Th. 1]):** Let $u$ be a sequence, $<^c$ be a weighted degree monomial order and $M$ be a monomial. Let us assume that the reduced $<^c$-Gröbner basis $\mathscr{G}$ of $I_C(u)$ and its associated staircase $S$ satisfy $\max(S \cup \text{LM}_<^c(\mathscr{G})) \le M$ and for all $m \le M$, $s = \max_{\sigma \le M} \{\sigma \mid \sigma m \le M\}$, we have $\max(S) \le s$. Then, the variant of the Berlekamp–Massey–Sakata algorithm using polynomial arithmetic terminates and computes $\mathscr{G}$ in $O(\#S(\#S + \#\mathscr{G})\#\mathcal{T}_{\le M})$ operations in the base field, where $\mathcal{T}_{\le M}$ is the set of monomials less or equal to $M$.

### 1.4.5 Problem 5: Minimization of number of queries

Assuming the first terms of the sequence are generic enough so that no fake C-relations are wrongly guessed, we design an adaptive algorithm that follows the shape of the staircase of the output Gröbner basis. It follows a strategy similar to the FGLM algorithm, which allows us to minimize the number of queries to the sequence.

**Theorem 1.14 (see also [19, Th. 26]):** Let $u$ be a sequence whose ideal of C-relations is zero-dimensional. Let $<^c$ be a monomial order, $\mathscr{G}$ be the reduced $<^c$-Gröbner basis of $I_C(u)$ and $S$ be the associated staircase.

Assuming the ADAPTIVE SCALAR-FGLM algorithm called on $u$ returns a Gröbner basis $\mathscr{G}'$ with staircase $S'$ and $\#S' = \#S$, then $S' = S$ and $\mathscr{G}' = \mathscr{G}$. Furthermore, the algorithm does not need more that $\#2(S \cup \text{LM}_<^c(\mathscr{G}))$ sequence queries and $O((\#S + \#\mathscr{G})^2\#2S)$ operations to recover $\mathscr{G}$, where $2S$ is the Minkowski sum of $S$ with itself.
### 1.4.6 Problem 6: Structured guessing

P-relations, linear recurrence relations with *polynomial*, in the indices, coefficients can be represented by quasi-commutative polynomials. In this setting, we extend the notion of sparse Gröbner basis of [7, 56]. This corresponds to P-relations between terms of the sequence that are all lying in a predefined cone. This finds application for instance in enumerative combinatorics where the indices of the nonzero terms are far from random but are actually all in a given cone. On the one hand, this allows us to reduce the number of sequence queries to guess the relations and on the other hand, this allows us to guess more relations that are correct and fewer that are fake. For instance, for a subsequence of the Gessel walk, using 3491 sequence terms in the full orthant, we can guess 142 relations amongst which 136 are fake and only 6 are correct. Whereas, taking only sequence terms in a cone allows us to consider sequence terms much further, which in turn allow us to guess more relations. Indeed, with 3010 terms in a cone of the same sequence, we guess 21 relations and all of them are correct. We refer to Table 5.1 for more details. Let us also notice that these fake relations may hide correct ones as their leading monomials could divide the leading monomials of correct relations.

### 1.4.7 MSOLVE

MSOLVE [13, 14] is a new efficient open-source C library, developed by Christian Eder<fnref n="1" />, Mohab Safey El Din<fnref n="2" /> and myself, for solving polynomial systems exactly. It provides implementations of state-of-the-art algorithms for computing Gröbner bases such as $F_4$ [47] for <DRL-Gröbner basis and SPARSE-FGLM [52, 53] for <LEX-Gröbner basis of zero-dimensional ideals.

I implemented in MSOLVE the parametrization of the radical [29, 66], its probabilistic test, see Lemma 1.12 and the SPARSE-FGLM-COLON algorithm for computing the <LEX-Gröbner basis of $I : \langle\varphi\rangle$ whether $I$ is zero- or positive-dimensional.

## 1.5 Organization

First, I start by presenting MSOLVE in Chapter 2. I detail the state-of-the-art algorithms that are available in MSOLVE and compare it with other available implementations of Gröbner bases algorithms. This is a joint development with Christian Eder<fnref n="1" /> and Mohab Safey El Din<fnref n="2" />.

Then, Chapter 3 is dedicated to the design of new algorithms for computing Gröbner bases in order to solve polynomial systems and their complexity analyses. Some of these algorithms are already available in MSOLVE. This is based on joint collaborations with Alin Bostan<fnref n="3" />, Christian Eder<fnref n="1" />, Andrew Ferguson<fnref n="4" />, Vincent Neiger<fnref n="2" /> and Mohab Safey El Din<fnref n="2" />.

In Chapter 4, I give an overview on different algorithms for guessing Gröbner bases of ideal of C-relations of sequences and their efficiency. This is based on joint work with Brice Boyer<fnref n="5" />, Jean-Charles Faugère<fnref n="3" />,<fnref n="6" /> and Mohab Safey El Din<fnref n="2" />.

In Chapter 5, I present extensions of works for guessing ideals of relations, but in the case of relations with *polynomial* coefficients. This is again related to polynomial ideals but in a ring where variables *quasi-commute*. This is based on joint work with Jean-Charles Faugère<fnref n="3" />,<fnref n="6" /> and Mohab Safey El Din<fnref n="2" />.

Finally, in Chapter 6, I propose my research project for the forthcoming years based on increasing the functionalities of MSOLVE by implementing existing state-of-the-art algorithms and also on designing algorithms for computing Gröbner bases in the commutative and quasi-commutative towards polynomial system solving and applications to sequences.

<footnote n="1">¹Rheinland-Pfälzische Technische Universität Kaiserslautern-Landau</footnote>

<footnote n="2">²Sorbonne Université</footnote>

<footnote n="3">³INRIA</footnote>

<footnote n="4">⁴former Ph.D. student at Sorbonne Université</footnote>

<footnote n="5">⁵former Post-doc at Sorbonne Université</footnote>

<footnote n="6">⁶CryptoNext Security</footnote>
# Chapter 1. Introduction

## 1.5. Organization
MSOLVE [13, 14] is an open source C library developed by Christian Eder<fnref n="1" />, Mohab Safey El Din<fnref n="2" /> and myself, for solving multivariate polynomial systems, with a focus on those which have dimension at most 0. It relies on computer algebra methods yielding algebraic parametrizations of their solutions. This allows us to bypass the commonly encountered issues related to accuracy and exhaustively met by numerical methods because of the non-linearity of the input.

This chapter is based on joint work with Christian Eder<fnref n="1" /> and Mohab Safey El Din<fnref n="2" /> and in particular [13].

## 2.1 State-of-the-art algorithms

For positive-dimensional ideals, MSOLVE computes their <DRL-Gröbner basis using the $F_4$ algorithm [47]. For 0-dimensional ideals, it returns a parametrization of their solutions. In generic coordinates, and for a radical ideal, this parametrization is close to the <LEX-Gröbner basis of the ideal. To compute it, it follows the framework of Figure 1.1. The change of order from <DRL to <LEX is performed thanks to the SPARSE-FGLM algorithm [52, 53]. If the ideal is not radical, but its radical is in shape position, then it computes the parametrizations thanks to [29, 66]. In this latter case, we cannot ensure that the ideal is in generic coordinates, hence we designed and implemented into MSOLVE an efficient probabilistic test to ensure that the computed parametrizations are parametrizations of the radical ideal. The presentation of this test is postponed to Section 3.3, as it relies on algorithms presented in Chapter 3.

Since MSOLVE is also now our tool of choice to implement and validate the new algorithms that we design, some algorithms presented in the following chapters, for instance, $F_4$SAT, see Section 3.2.2, and SPARSE-FGLM-COLON, see Section 3.2.3, have been implemented and are also available in MSOLVE.

### 2.1.1 The $F_4$ algorithm

In [35], Buchberger's algorithm introduced the concept of *critical pairs* for computing Gröbner bases. For two polynomials $f_1$ and $f_2$ in a set of generators of an ideal, the critical pair $(f_1, f_2)$ leads to a normal form computation of the *S-polynomial*

$$ \mathrm{sp}_<(f_1, f_2) = \frac{\mathrm{LCM}(\mathrm{LM}_<(f_1), \mathrm{LM}_<(f_2))}{\mathrm{LT}_<(f_1)} f_1 - \frac{\mathrm{LCM}(\mathrm{LM}_<(f_1), \mathrm{LM}_<(f_2))}{\mathrm{LT}_<(f_2)} f_2 $$

w.r.t. the current intermediate basis. The *degree* of such a critical pair is defined as $\deg \mathrm{LCM}(\mathrm{LM}_<(f_1), \mathrm{LM}_<(f_2))$. Notice that this bounds from above $\deg \mathrm{sp}_<(f_1, f_2)$.

In Algorithm 2.1 we state the pseudocode of Faugère's $F_4$ algorithm, highlighting (in red) the main differences to Buchberger's algorithm.

Observe that the termination of the $F_4$ algorithm only relies on Buchberger's first criterion: $\mathcal{G} = \{g_1, \dots, g_i\}$ is a <Gröbner basis of an ideal $I$ if for all $1 \le i, j \le t$, $\mathrm{NF}(\mathrm{sp}_<(g_i, g_j), \mathcal{G}, <) = 0$, see [41, Chap. 2, Sec. 6, Th. 6].

We detail the differences to Buchberger's algorithm.

(1) One can choose several critical pairs at a time, stored in a subset $L \subseteq P$. The so-called *degree strategy* chooses $L$ to be the set of all critical pairs of minimal degree for a total degree monomial order, typically <DRL.

(2) For all terms of all the generators of the S-polynomials, one searches in the current intermediate Gröbner basis $\mathcal{G}$ for possible reducers. One adds those to $L$ and again search for all of their terms for reducers in $\mathcal{G}$. This is the SymbolicPreprocessing function.

(3) Once all reduction data is collected from the last step, one generates a Macaulay-like matrix

<footnote n="1">¹Rheinland-Pfälzische Technische Universität Kaiserslautern-Landau</footnote>

<footnote n="2">²Sorbonne Université</footnote>
**Input:** A list of polynomials $f_1, \dots, f_s$ spanning an ideal $I \subseteq \mathbb{K}[\mathbf{x}]$ and a total degree monomial order $<$.

**Output:** A $<$-Gröbner basis $\mathcal{G}$ of $I$.

$$
1 \mathcal{G} := \{f_1, \dots, f_s\}.
$$

$$
2 P := \{(f_i, f_j) \mid 1 \le i < j \le s\}.
$$

3 While $P \neq \emptyset$ do

$$
\begin{array}{r@{\quad}l}
4 & \text{Choose a subset } L \text{ of } P. \\
5 & \mathcal{P} := \mathcal{P} \setminus L. \\
6 & \mathcal{L} := \text{SymbolicPreprocessing}(L, \mathcal{G}). \\
7 & \mathcal{L} := \text{LinearAlgebra}(L). \\
8 & \textbf{For } h \in L \text{ with } \mathit{LM}_<^<(h) \notin \langle \mathit{LM}_<^<(\mathcal{G}) \rangle \textbf{ do} \\
9 & \quad \mathcal{P} := \mathcal{P} \cup \{(g, h) \mid g \in \mathcal{G}\}. \\
10 & \quad \mathcal{G} := \mathcal{G} \cup \{h\}. \\
11 & \textbf{Return } \mathcal{G}.
\end{array}
$$

<caption>**Algorithm 2.1:** Faugère's $F_4$</caption>

with columns corresponding to the monomials appearing in $L$ (sorted by $<$) and rows corre-
sponding to the coefficients of each polynomial in $L$. Gaussian Elimination is then applied to
reduce now all chosen S-polynomials at once. This is the LinearAlgebra function.

(q) Finally, one adds those polynomials associated to rows in the updated matrix to $\mathcal{G}$ whose leading monomials are not already in $\mathrm{LM}_<^<(\mathcal{G})$.

In order to optimize the algorithm one can now apply Buchberger’s product and chain criteria, see [36, 70]. Thus many useless critical pairs are removed before even being added to P and fewer zero rows are computed during the linear algebra part of F₄. Still, in general, there are many zero reductions left.

Different selection strategies yield different behavior of the algorithm. The degree strategy allows one to compute *truncated Gröbner bases* of ideals in case of early terminations.

**Definition 2.1:** Let $f_1, \dots, f_s$ be polynomials in $\mathbb{K}[\mathbf{x}]$ and $<$ be a monomial order. Let $\mu$ be a monomial and $F_\mu$ be the $\mathbb{K}$-vector subspace of $\langle f_1, \dots, f_s \rangle$ defined as

$$
F_{\mu} = \left\{ \sum_{i=1}^{s} h_i f_i \left| \forall 1 \le i \le s, \text{LT}_{<}(h_i f_i) \le \mu \right. \right\}.
$$

Then, $\mathcal{G} \subset F_\mu$ is a $\mu$-truncated $<$-Gröbner basis of $\langle f_1, \dots, f_s \rangle$ if for all $p \in F_\mu$, there exists $g \in \mathcal{G}$ such that $\mathrm{LM}_<(g) | \mathrm{LM}_<(p)$ and $p - \frac{\mathrm{LT}_<(p)}{\mathrm{LT}_<(g)}g$ is in $F_\mu$.

Observe that taking a triangular basis of $F_\mu$ ordered increasingly w.r.t. $<$ naturally yields a
$\mu$-truncated $<$-Gröbner basis thereof.

**Proposition 2.2:** Let $f_1, \dots, f_s$ be polynomials in $\mathbb{K}[\mathbf{x}]$ and $<$ be a monomial order. Let $\mu$ be a monomial and $F_\mu$ be the $\mathbb{K}$-vector subspace of $\langle f_1, \dots, f_s \rangle$

$$
F_{\mu} = \left\{ \sum_{i=1}^{s} h_i f_i \left| \forall 1 \le i \le s, \mathrm{LM}_{<}(h_i f_i) \le \mu \right. \right\}.
$$

A subset $\mathcal{G} = \{g_1, \dots, g_t\} \subset F_\mu$ is a $\mu$-truncated $<$-Gröbner basis of $\langle f_1, \dots, f_s \rangle$ if, and only if,

$$
F_{\mu} \subseteq G_{\mu} = \left\{ \sum_{j=1}^{t} h_j g_j \left| \forall 1 \le j \le t, \mathrm{LM}_{<}(h_j g_j) \le \mu \right. \right\}.
$$

and for all $(g_i, g_j) \in \mathcal{G}^2$ with $i \neq j$, if $\mathrm{LCM}(\mathrm{LM}_<(g_i), \mathrm{LM}_<(g_j)) \le \mu$, then

$$
\mathrm{NF}(\mathrm{sp}_<(g_i, g_j), \mathcal{G}, <) = 0.
$$

**Remark 2.3:** (1) If $<$ is a total degree monomial order, then for $d \in \mathbb{N}$, we can also define a $d$-truncated $<$-Gröbner basis as a $\mu$-truncated $<$-Gröbner basis for $\mu$ the largest monomial of
degree $d$.

(2) If $\mathcal{G} = \{g_1, \dots, g_l\}$ is a $\mu$-truncated $<\text{Gröbner}$ basis of $\langle f_1, \dots, f_s \rangle$ and

$$ \mu \geq \max_{1 \leq i < j \leq l} \text{LCM}(\text{LM}_< (g_i), \text{LM}_< (g_j)), $$

then $\mathcal{G}$ is a $<\text{Gröbner}$ basis of $\langle f_1, \dots, f_s \rangle$. Indeed, it spans the ideal and by Proposition 2.2, all the S-polynomials reduce to 0 w.r.t. $\mathcal{G}$ and $<$. Hence, by Buchberger's first criterion [41, Chap. 2, Sec. 6, Th. 6], it is a $<\text{Gröbner}$ basis of $\langle f_1, \dots, f_s \rangle$.

(3) Definition 2.1 depends greatly on the set of generators of the ideal. Consider $f_1 = x^n$, $f_2 = (y-1)^n$ and $f_3 = xy-y-1$ for $n \in \mathbb{N} \setminus \{0,1\}$. By Proposition 2.2, $\mathcal{G} = \{f_1, f_2, f_3\}$ is a $n$-truncated $<\text{DRL}$-Gröbner basis of $\langle f_1, f_2, f_3 \rangle$. Yet, this ideal is $\langle 1 \rangle$ hence $\{1\}$ is a $m$-truncated $<\text{DRL}$-Gröbner basis of $\langle 1 \rangle$ for all $m \in \mathbb{N}$.

**Lemma 2.4:** Let $f_1, \dots, f_s \in \mathbb{K}[x]$ be the input polynomials of the $F_4$ algorithm. Let $d \in \mathbb{N}$. Assume that the $F_4$ algorithm uses the degree selection strategy and that, on line 4, $L$ consists in all the critical pairs of degree $d$.

If no new polynomial is added to $\mathcal{G}$ on line 10, then $\mathcal{G}$ is a $d$-truncated $<\text{DRL}$-Gröbner basis of $\langle f_1, \dots, f_s \rangle$.

## 2.1.2 The SPARSE-FGLM algorithm

In this subsection, the input Gröbner basis, $\mathcal{G}_{\text{DRL}}$, is the reduced $<\text{DRL}$-Gröbner basis of a zero-dimensional ideal $I$ of degree $D$. The output is the reduced $<\text{LEX}$-Gröbner basis, $\mathcal{G}_{\text{LEX}}$, of $I$. In [52] and [53, Algo. 3], using [80], the authors observe that the map

$$ \mathbb{K}[x]/I \to \mathbb{K}[x]/I \\ f \mapsto \text{NF}(x_n f, \mathcal{G}_{\text{DRL}}, <_{\text{DRL}}) $$

given in the basis $S_{\text{DRL}}$, the staircase associated to $\mathcal{G}_{\text{DRL}}$, is represented by a matrix, $M_{x_n}$, with a special structure given in the following two lemmas.

**Lemma 2.5:** Let $I$ be a zero-dimensional ideal of $\mathbb{K}[x]$ of degree $D$, $\mathcal{G}_{\text{DRL}}$ be its reduced $<_{\text{DRL}}$-Gröbner basis and $S_{\text{DRL}} = \{\sigma_0, \dots, \sigma_{D-1}\}$ be its associated staircase. Let $M_{x_n}$ be the matrix of the linear map $f \in \mathbb{K}[x]/I \mapsto \text{NF}(x_n f, \mathcal{G}_{\text{DRL}}, <_{\text{DRL}}) \in \mathbb{K}[x]/I$.

Then, one can build the matrix $M_{x_n} = (m_{i,j})_{0 \le i,j < D}$ with the following procedure:

* if $x_n \sigma_j = \sigma_k$, then $m_{k,j} = 1$ and for all $0 \le i < D$, $i \neq k$, $m_{i,j} = 0$;

* otherwise for all $0 \le i < D$, $m_{i,j}$ is the coefficient of $\sigma_i$ in $\text{NF}(x_n \sigma_j, \mathcal{G}_{\text{DRL}}, <_{\text{DRL}})$.

**Lemma 2.6 ([52, 53, 80]):** Let $f_1, \dots, f_n$ be generic polynomials of $\mathbb{K}[x]$ of degrees at most $d$. Let $\mathcal{G}_{\text{DRL}}$ be the reduced $<_{\text{DRL}}$-Gröbner basis of $\langle f_1, \dots, f_n \rangle$. Then, the latter case of Lemma 2.5 only happens if there exists $g \in \mathcal{G}_{\text{DRL}}$ such that $\text{LM}_{<_{\text{DRL}}}(g) = x_n \sigma_j$. As a consequence, one has $\text{NF}(x_n \sigma_j, \mathcal{G}_{\text{DRL}}, <_{\text{DRL}}) = x_n \sigma_j - g$.

Following, we can use Wiedemann algorithm [106] on $M_{x_n}$ to recover its minimal polynomial. Furthermore, whenever the reduced $<_{\text{LEX}}$-Gröbner basis $\mathcal{G}_{\text{LEX}}$ is in shape position, i.e. there exist $g_n, g_{n-1}, \dots, g_1 \in \mathbb{K}[x_n]$ such that

$$ \mathcal{G}_{\text{LEX}} = \{g_n(x_n), x_{n-1} - g_{n-1}(x_n), \dots, x_1 - g_1(x_n)\}, $$

and for all $1 \le k \le n-1$, deg $g_k < \text{deg} g_n$, then $g_1, \dots, g_{n-1}$ can be computed by solving Hankel systems of size $D$. This can be done using the following two algorithms, Algorithms 2.2 and 2.3.

**Proposition 2.7:** Let $M \in \mathbb{K}^{D \times D}$ be a matrix with $s$ nonzero coefficients, $r \in \mathbb{K}^D$ be a row-vector and $c_0, c_1, \dots, c_{n-1} \in \mathbb{K}^D$ be $n$ column-vectors. Then, Algorithm 2.2 is correct and computes the sequences $(rM^i c_0)_{0 \le i < 2D}$ and $(rM^i c_k)_{0 \le i < D}$ for $1 \le k \le n-1$ in $O(sD + nD^2)$ operations in $\mathbb{K}$.

Furthermore, if the vectors $c_0, c_1, \dots, c_{n-1}$ are vectors of the canonical basis, then this complexity drops to $O((s+n)D)$.

**Proposition 2.8:** Let $(v_i^{(0)})_{0 \le i < 2D-1}, (v_i^{(1)})_{0 \le i < D}, \dots, (v_i^{(n-1)})_{0 \le i < D}$ be the first terms of $n$ sequences. Assume that $v^{(0)}$ is linear recurrent of order $D$. Then, Algorithm 2.3 is correct and computes, for
**Input:** A matrix $M \in \mathbb{K}^{D \times D}$, a row-vector $r \in \mathbb{K}^D$ and $n$ column-vectors $c_0, c_1, \dots, c_{n-1} \in \mathbb{K}^D$

**Output:** $(rMc_0)_{0\le i<2D}, (rMc_1)_{0\le i<2D}, \dots, (rMc_{n-1})_{0\le i<2D}$, with $1 = (1, 0, \dots, 0)^T$.

1. $v_0^{(0)} := rc_0, v_0^{(1)} := rc_1, \dots, v_0^{(n-1)} := rc_{n-1}$.

2. For $i$ from 1 to $D-1$ do

3. $r := rM$.

4. $v_i^{(0)} := rc_0, v_i^{(1)} := rc_1, \dots, v_i^{(n-1)} := rc_{n-1}$.

5. For $i$ from $D$ to $2D-1$ do

6. $r := rM$.

7. $v_i^{(0)} := rc_0$.

8. **Return** $(v_i^{(0)})_{0\le i<2D}, (v_i^{(1)})_{0\le i<2D}, \dots, (v_i^{(n-1)})_{0\le i<2D}$.

<caption>**Algorithm 2.2:** Sequences for SPARSE-FGLM</caption>

**Input:** Sequences $(v_i^{(0)})_{0\le i<2D-1}$ and $(v_i^{(k)})_{0\le i<2D}$ for $1 \le k \le n-1$ with coefficients in $\mathbb{K}$.

**Output:** $\gamma_{0,k}, \dots, \gamma_{D-1,k}$ for $1 \le k \le n-1$ such that for all $0 \le i < D$,

$$v_i^{(k)} = \gamma_{D-1,k} v_{D-1+i}^{(0)} + \dots + \gamma_{0,k} v_i^{(0)}.$$

1. For $k$ from 1 to $n-1$ do

2. Solve the Hankel linear system

$$\begin{pmatrix} v_0^{(0)} & v_1^{(1)} & \cdots & v_{D-1}^{(0)} \\ v_1^{(0)} & v_2^{(1)} & \cdots & v_D^{(0)} \\ \vdots & \vdots & \ddots & \vdots \\ v_{D-1}^{(0)} & v_D^{(1)} & \cdots & v_{2D-2}^{(0)} \end{pmatrix} \begin{pmatrix} \gamma_{0,k} \\ \gamma_{1,k} \\ \vdots \\ \gamma_{D-1,k} \end{pmatrix} = \begin{pmatrix} v_0^{(k)} \\ v_1^{(k)} \\ \vdots \\ v_{D-1}^{(k)} \end{pmatrix}.$$

3. **Return** $\gamma_{i,k}$ for $0 \le i < D$ and $1 \le k \le n-1$.

<caption>**Algorithm 2.3:** Hankel system solving for SPARSE-FGLM</caption>
all $1 \le k \le n-1$, $\gamma_{0,k}, \dots, \gamma_{D-1,k}$ such that

$$ \forall 0 \le i < D, v_i^{(k)} = \gamma_{D-1,k} v_{D-1+i}^{(0)} + \dots + \gamma_{0,k} v_i^{(0)} $$

in $O(M(D)(n + \log D))$ operations, where $M(D)$ denote a cost function for multiplying univariate polynomials of degree $D$ with coefficients in $\mathbb{K}$.

We are now in a position to present the SPARSE-FGLM algorithm in the shape position case.

**Input:** The reduced <DRL>-Gröbner basis $\mathcal{G}_{\text{DRL}}$ of a zero-dimensional ideal $I$ and its associated staircase $S_{\text{DRL}}$ of size $D$.

**Output:** The reduced <LEX>-Gröbner basis of $I$, if it is in shape position.

1. Build the matrix $M$ as in Lemma 2.5.

2. Pick $r \in \mathbb{K}^D$ a row-vector at random.

3. $\mathbf{1} := (1, 0, \dots, 0)^T$. // the column-vector of coefficients of NF $(1, \mathcal{G}_{\text{DRL}}, <_{\text{DRL}})$

4. **For $k$ from 1 to $n-1$ do**

5. Build $\mathbf{c}_k$ the column-vector of coefficients of NF $(x_k, \mathcal{G}_{\text{DRL}}, <_{\text{DRL}})$.

6. Compute $(v_i^{(0)})_{0 \le i < 2D}, (v_i^{(1)})_{0 \le i < D}, \dots, (v_i^{(n-1)})_{0 \le i < D}$ with Algorithm 2.2 called on $M, r, \mathbf{1}, \mathbf{c}_1, \dots, \mathbf{c}_{n-1}$.

7. $g_n := \text{Berlekamp-Massey}(v_0^{(0)}, \dots, v_{2D-1}^{(0)})$.

8. **If $\deg g_n < D$ then Return** "Not in shape position or bad vector".

9. Compute $g_1 := \gamma_{D-1,1}x_n^{D-1} + \dots + \gamma_{0,1}, \dots, g_{n-1} := \gamma_{D-1,n-1}x_n^{D-1} + \dots + \gamma_{0,n-1}$ with Algorithm 2.3 called on $(v_i^{(0)})_{0 \le i < 2D-1}, (v_i^{(1)})_{0 \le i < D}, \dots, (v_i^{(n-1)})_{0 \le i < D}$.

10. **Return** $\{g_n(x_n), x_{n-1} - g_{n-1}(x_n), \dots, x_1 - g_1(x_n)\}$.

**Algorithm 2.4: SPARSE-FGLM**

**Theorem 2.9:** Let $I$ be a zero-dimensional ideal of $\mathbb{K}[\mathbf{x}]$ of degree $D$, $\mathcal{G}_{\text{DRL}}$ be its reduced <DRL>-Gröbner basis and $S_{\text{DRL}}$ be its associated staircase. Let $M_{x_n}$ be the matrix of the map $f \in \mathbb{K}[\mathbf{x}]/I \mapsto \text{NF}(x_n f, \mathcal{G}_{\text{DRL}}, <_{\text{DRL}}) \in \mathbb{K}[\mathbf{x}]/I$ in the monomial basis $S_{\text{DRL}}$.

Let us assume that there are $t$ monomials $\sigma$ in $S_{\text{DRL}}$ such that $x_n\sigma \in \text{LT}_{<_{\text{DRL}}}(I)$ and that $x_1, \dots, x_{n-1} \in S_{\text{DRL}}$, that $M_{x_n}$ is known and that the reduced <LEX>-Gröbner basis $\mathcal{G}_{\text{LEX}}$ of $I$ is in shape position. Then, Algorithm 2.4 computes $\mathcal{G}_{\text{LEX}}$ in $O(tD^2 + nM(D))$ operations, where $M(D)$ denote a cost function for multiplying univariate polynomials of degree $D$ with coefficients in $\mathbb{K}$.

Note that the Berlekamp-Massey algorithm and its faster variants return a factor of $g_n$, so if the computed polynomial has degree $D$, i.e. it is the characteristic polynomial of $M_{x_n}$, then it is also its minimal polynomial. Furthermore, based on a deterministic variant of Wiedemann's algorithm, one can also provide a deterministic variant of this algorithm to recover $g_n$ [53, Algo. 4].

**Remark 2.10:** In [29, 66], the authors consider the case where an ideal $J$ is not in shape position but its radical $\sqrt{J}$ is. Let us recall that $\sqrt{J} = \{f \in \mathbb{K}[\mathbf{x}] | \exists k \in \mathbb{N}, f^k \in J\}$, see [41, Chap. 4, Sec. 2, Def. 4]. In that case, the <LEX>-Gröbner basis of $\sqrt{J}$ can be computed in a similar fashion, it suffices to replace the call to Algorithm 2.3 on line 9 by a call to [66, Algo. 2].

**Example 2.11:** Let $\mathcal{H}_{\text{DRL}} = \{x_2^2 - x_1 + x_2, x_1x_2 - x_1 + x_2, x_1^2 - x_1 + x_2\}$ be the <DRL>-Gröbner basis of the ideal $J$. Its associated staircase is $T_{\text{DRL}} = \{1, x_2, x_1\}$ and the matrix of the multiplication by $x_2$ is

$$ \begin{pmatrix} 1 & x_2 & x_1 \\ 0 & 0 & 0 \\ 1 & -1 & -1 \\ 0 & 1 & 1 \end{pmatrix} \begin{pmatrix} 1 \\ x_2 \\ x_1 \end{pmatrix} $$

We build a sequence $u$ with random initial coefficients, e.g. $u_{0,0} = 17$, $u_{1,0} = 5$ and $u_{0,1} = -3$, thanks
to $\mathcal{G}_{\text{DRL}}$:

<table><thead><tr><th>...</th><th>...</th><th>...</th><th>...</th><th>...</th><th>...</th></tr></thead><tbody><tr><td>3</td><td>0</td><td>0</td><td>0</td><td>0</td><td>...</td></tr><tr><td>2</td><td>-8</td><td>0</td><td>0</td><td>0</td><td>...</td></tr><tr><td>1</td><td>5</td><td>-8</td><td>0</td><td>0</td><td>...</td></tr><tr><td>0</td><td>17</td><td>-3</td><td>-8</td><td>0</td><td>...</td></tr><tr><td><i>i</i><sub>2</sub></td><td>0</td><td>1</td><td>2</td><td>3</td><td>...</td></tr></tbody></table>

The subsequence $(u_{0,i_2})_{i_2 \in \mathbb{N}}$ satisfies the relation $u_{0,i_2+3} = 0 = \text{Eval}(x_2^{i_2+3}, u)$ for all $i_2 \in \mathbb{N}$ and no relation of smaller order. Hence $x_2^3 \in J$.

We now solve the Hankel system

$$ \begin{pmatrix} u_{0,0} & u_{0,1} & u_{0,2} \\ u_{0,1} & u_{0,2} & u_{0,3} \\ u_{0,2} & u_{0,3} & u_{0,4} \end{pmatrix} \begin{pmatrix} \gamma_0 \\ \gamma_1 \\ \gamma_2 \end{pmatrix} = \begin{pmatrix} u_{1,0} \\ u_{1,1} \\ u_{1,2} \end{pmatrix} \begin{pmatrix} 17 & 5 & -8 \\ 5 & -8 & 0 \\ -8 & 0 & 0 \end{pmatrix} \begin{pmatrix} \gamma_0 \\ \gamma_1 \\ \gamma_2 \end{pmatrix} = \begin{pmatrix} -3 \\ -8 \\ 0 \end{pmatrix} $$

to find $\gamma_0 = 0$ and $\gamma_1 = \gamma_2 = -1$. Hence, $x_1 - x_2^2 - x_2$ is in $J$. Observe that we already knew this though. Finally, $\mathcal{H}_{\text{LEX}} = \{x_2^3, x_1 - x_2^2 - y\}$.

## 2.2 Experimental results

In order to solve over rational numbers, a multi-modular approach with efficient algorithms is implemented. The $F_4$ part of the computation is performed with a tracer modulo a first prime $p_1$. This first computation allows us to learn all the polynomials and their multiples that are needed at each round and also those that are reduced to 0 or that are only used to reduce other polynomials to 0. Modulo subsequent primes, we apply the tracer so that the matrices we handle are optimal: we do not compute any reduction to 0. Furthermore, we know that this multi-modular approach correctly computes the sought $<_{\text{DRL}}$-Gröbner basis over $\mathbb{Q}$ as long as its projection modulo $p_1$ was correctly computed.

We compare MSOLVE with two other computer algebra systems:

* `magma` -v2.23-6 [24]: using the command `Variety()`.

* `maple` -v2019 [9]: using the command `PolynomialSystem()` from the module `SolveTools` with option `engine=groebner`.

All compared implementations use Faugère's $F_4$ algorithm and variants of the FGLM algorithm and then solve univariate polynomials.

All chosen systems are zero-dimensional with rational coefficients. All computations are done sequentially. First, Table 2.1 compares memory usage of MSOLVE against `maple` and `magma`. It illustrates the low memory usage of MSOLVE. However, we emphasize that MSOLVE is a specialized library while `maple` and `magma` are general purpose computer algebra systems.

Overall, MSOLVE performs very efficiently on a wide range of input systems, using way less memory than its competitors, allowing its users to solve polynomial systems which are not tractable by `maple` and `magma`.

Table 2.2 states various, partly well-known benchmarks, which differ in their specific hardness, like reduction process, pair handling, sparsity of multiplication matrices, etc. It also deals with critical points computations, $CP(d,nv, np)$ describes critical points for a system of $np$ polynomials in $nv$ variables of degree $d$.

For each system we give its degree and if it is radical (all but one are radical). For MSOLVE we give specific timing information, also on the single modular computations: We apply MSOLVE with the tracer option, giving also the timings for the first modular computation learning and generating the tracer ($F_4$ (learn)) and the timings for the further modular computations applying only the
<table><caption>Table 2.1: Maximal memory usage given in MB</caption><thead><tr><th>Examples</th><th>MSOLVE</th><th>maple</th><th>magma</th><th>Examples</th><th>MSOLVE</th><th>maple</th><th>magma</th></tr></thead><tbody><tr><td>Katsura-9</td><td>15</td><td>271</td><td>71</td><td>Henrion-5</td><td>11</td><td>26</td><td>23</td></tr><tr><td>Katsura-10</td><td>25</td><td>276</td><td>223</td><td>Henrion-6</td><td>47</td><td>171</td><td>-</td></tr><tr><td>Katsura-11</td><td>66</td><td>2,279</td><td>-</td><td>Henrion-7</td><td>3,428</td><td>-</td><td>-</td></tr><tr><td>Katsura-12</td><td>229</td><td>1,123</td><td>-</td><td>Noon-7</td><td>209</td><td>419</td><td>-</td></tr><tr><td>Katsura-13</td><td>1,037</td><td>-</td><td>-</td><td>Noon-8</td><td>881</td><td>1,227</td><td>-</td></tr><tr><td>Eco-10</td><td>27</td><td>210</td><td>213</td><td>CP(3,5,2)</td><td>17</td><td>525</td><td>-</td></tr><tr><td>Eco-11</td><td>82</td><td>428</td><td>354</td><td>CP(3,6,2)</td><td>55</td><td>7,885</td><td>-</td></tr><tr><td>Eco-12</td><td>117</td><td>1,027</td><td>-</td><td>CP(3,7,2)</td><td>312</td><td>-</td><td>-</td></tr><tr><td>Eco-13</td><td>318</td><td>8,654</td><td>-</td><td>CP(4,4,3)</td><td>24</td><td>635</td><td>-</td></tr><tr><td>Eco-14</td><td>15,748</td><td>-</td><td>-</td><td>CP(4,5,3)</td><td>2,065</td><td>-</td><td>-</td></tr><tr><td>Phuoc-1</td><td>176</td><td>-</td><td>-</td><td></td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>Table 2.2: Benchmark timings given in seconds.</caption><thead><tr><th rowspan="2">Examples</th><th rowspan="2">System data degree</th><th rowspan="2">radical</th><th colspan="4">MSOLVE single modular computation</th><th colspan="3">MSOLVE overall</th><th colspan="2">maple single modular</th><th colspan="2">Others overall</th></tr><tr><th>F<sub>4</sub> (prob.)</th><th>F<sub>4</sub> (learn)</th><th>F<sub>4</sub> (apply)</th><th>FGLM</th><th># primes</th><th>trace</th><th>independent</th><th>F<sub>4</sub></th><th>FGLM</th><th>maple</th><th>magma</th></tr></thead><tbody><tr><td>Katsura-9</td><td>256</td><td>yes</td><td>0.06</td><td>0.17</td><td>0.03</td><td>0.03</td><td>83</td><td>4.89</td><td>7.49</td><td>0.10</td><td>0.04</td><td>104</td><td>2,522</td></tr><tr><td>Katsura-10</td><td>512</td><td>yes</td><td>0.24</td><td>0.81</td><td>0.09</td><td>0.11</td><td>188</td><td>43.7</td><td>70.5</td><td>0.36</td><td>0.15</td><td>1,278</td><td>82,540</td></tr><tr><td>Katsura-11</td><td>1,024</td><td>yes</td><td>1.34</td><td>6.26</td><td>0.45</td><td>0.49</td><td>388</td><td>424</td><td>814</td><td>1.82</td><td>0.74</td><td>7,812</td><td>-</td></tr><tr><td>Katsura-12</td><td>2,048</td><td>yes</td><td>8.61</td><td>56.1</td><td>3.10</td><td>3.96</td><td>835</td><td>6,262</td><td>11,215</td><td>8.50</td><td>5.40</td><td>120,804</td><td>-</td></tr><tr><td>Katsura-13</td><td>4,096</td><td>yes</td><td>52.8</td><td>425</td><td>18.9</td><td>30.6</td><td>1,772</td><td>89,390</td><td>148,372</td><td>60.9</td><td>35.7</td><td>-</td><td>-</td></tr><tr><td>Katsura-14</td><td>8,192</td><td>yes</td><td>318</td><td>3,336</td><td>128</td><td>210</td><td>3,847</td><td>1,308,602</td><td>2,007,170</td><td>393</td><td>271</td><td>-</td><td>-</td></tr><tr><td>Eco-10</td><td>256</td><td>yes</td><td>0.10</td><td>0.28</td><td>0.05</td><td>0.02</td><td>161</td><td>12.5</td><td>21.2</td><td>0.14</td><td>0.03</td><td>26.3</td><td>6,520</td></tr><tr><td>Eco-11</td><td>512</td><td>yes</td><td>0.39</td><td>1.21</td><td>0.17</td><td>0.07</td><td>327</td><td>90.3</td><td>161</td><td>0.56</td><td>0.12</td><td>312</td><td>214,770</td></tr><tr><td>Eco-12</td><td>1,024</td><td>yes</td><td>2.25</td><td>11,619</td><td>1.07</td><td>0.34</td><td>530</td><td>877</td><td>1,619</td><td>2.97</td><td>0.85</td><td>4,287</td><td>-</td></tr><tr><td>Eco-13</td><td>2,048</td><td>yes</td><td>11.7</td><td>67.3</td><td>6.61</td><td>2.12</td><td>1,225</td><td>12,137</td><td>19,553</td><td>15.1</td><td>6.70</td><td>66,115</td><td>-</td></tr><tr><td>Eco-14</td><td>4,096</td><td>yes</td><td>67.1</td><td>516</td><td>34.8</td><td>25.9</td><td>2,670</td><td>167,798</td><td>254,389</td><td>104.8</td><td>69.1</td><td>-</td><td>-</td></tr><tr><td>Henrion-5</td><td>100</td><td>yes</td><td>0.01</td><td>0.01</td><td>0.004</td><td>0.01</td><td>83</td><td>0.71</td><td>0.83</td><td>0.01</td><td>0.01</td><td>2.7</td><td>93</td></tr><tr><td>Henrion-6</td><td>720</td><td>yes</td><td>0.11</td><td>0.22</td><td>0.07</td><td>0.11</td><td>612</td><td>138</td><td>157</td><td>0.17</td><td>0.16</td><td>1,470</td><td>-</td></tr><tr><td>Henrion-7</td><td>5,040</td><td>yes</td><td>9.55</td><td>27.5</td><td>6.51</td><td>20.46</td><td>4,243</td><td>117,803</td><td>127,456</td><td>12.8</td><td>27.1</td><td>-</td><td>-</td></tr><tr><td>Noon-7</td><td>2,173</td><td>yes</td><td>1.66</td><td>5.3</td><td>0.93</td><td>1.95</td><td>1,305</td><td>4,039</td><td>5,045</td><td>1.97</td><td>3.13</td><td>432</td><td>-</td></tr><tr><td>Noon-8</td><td>6,545</td><td>yes</td><td>26.6</td><td>153</td><td>17.5</td><td>72.3</td><td>6,462</td><td>598,647</td><td>640,177</td><td>32.4</td><td>76.2</td><td>5,997</td><td>-</td></tr><tr><td>Phuoc-1</td><td>1,102</td><td>no</td><td>4.01</td><td>4.65</td><td>3.42</td><td>2.59</td><td>753</td><td>4,467</td><td>5,056</td><td>4.60</td><td>5.91</td><td>-</td><td>-</td></tr><tr><td>CP(3,5,2)</td><td>288</td><td>yes</td><td>0.03</td><td>0.04</td><td>0.01</td><td>0.03</td><td>326</td><td>18.1</td><td>19.2</td><td>0.06</td><td>0.05</td><td>249</td><td>-</td></tr><tr><td>CP(3,6,2)</td><td>720</td><td>yes</td><td>0.22</td><td>0.59</td><td>0.12</td><td>0.16</td><td>1,042</td><td>390</td><td>450</td><td>0.31</td><td>0.22</td><td>23,440</td><td>-</td></tr><tr><td>CP(3,7,2)</td><td>1,728</td><td>yes</td><td>1.97</td><td>8.18</td><td>1.23</td><td>1.20</td><td>3,037</td><td>9,643</td><td>11,511</td><td>2.78</td><td>2.54</td><td>-</td><td>-</td></tr><tr><td>CP(3,8,2)</td><td>4,032</td><td>yes</td><td>18.5</td><td>111.5</td><td>12.2</td><td>19.6</td><td>8,211</td><td>269,766</td><td>323,838</td><td>24.6</td><td>25.3</td><td>-</td><td>-</td></tr><tr><td>CP(4,4,3)</td><td>576</td><td>yes</td><td>0.04</td><td>0.86</td><td>0.03</td><td>0.07</td><td>339</td><td>40.9</td><td>41.8</td><td>0.08</td><td>0.11</td><td>916</td><td>-</td></tr><tr><td>CP(4,5,3)</td><td>3,456</td><td>yes</td><td>3.24</td><td>8.60</td><td>2.23</td><td>4.83</td><td>2,747</td><td>21,528</td><td>23,559</td><td>4.33</td><td>9.21</td><td>-</td><td>-</td></tr><tr><td>CP(3,6,6)</td><td>729</td><td>yes</td><td>0.18</td><td>0.42</td><td>0.11</td><td>0.15</td><td>779</td><td>255</td><td>294</td><td>0.30</td><td>0.23</td><td>-</td><td>-</td></tr><tr><td>CP(4,6,6)</td><td>4,096</td><td>yes</td><td>7.70</td><td>25.6</td><td>5.44</td><td>14.09</td><td>3,476</td><td>71,472</td><td>77,941</td><td>10.2</td><td>16.71</td><td>-</td><td>-</td></tr><tr><td>CP(3,7,7)</td><td>2,187</td><td>yes</td><td>2.49</td><td>8.97</td><td>1.58</td><td>1.86</td><td>2,795</td><td>12,412</td><td>14,375</td><td>3.27</td><td>3.75</td><td>-</td><td>-</td></tr></tbody></table>
tracer (F4 (apply)). We also use MSOLVE with independent modular computations, applying the probabilistic linear algebra in each modular $F_4$ (F4 (prob.)) In any case, we apply the same FGLM implementation. Furthermore, we state the number of primes needed by MSOLVE to solve over $\mathbb{Q}$. For maple and magma we just give the overall timings. Symbol '-' means that the computation was stopped after waiting more than 10 times the runtime of MSOLVE. For all systems, the bottleneck has been the computation of either a `<LEX>`-Gröbner basis in shape position or a rational parametrization of the solution set.

First thing to note is that magma is in all instances slower than MSOLVE or maple. Although, for some examples, magma's modular $F_4$ computation is even a bit faster than the other two, magma's bottleneck is both a not optimized FGLM combined with the fact that magma seems to lift a `<LEX>`-Gröbner basis instead of a rational parametrization (the latter one having in general coefficients of significantly smaller bit size).

For nearly all systems, MSOLVE is faster, sometimes by an order of magnitude, than maple. We report on modular timings of maple for $F_4$ (which is based on a probabilistic linear algebra) and FGLM. It appears that MSOLVE's modular implementations of both $F_4$ and FGLM are faster than the ones in maple (with a speed-up sometimes close to 2, sometimes less). It seems that in the multi-modular process, maple uses its probabilistic variant of $F_4$ while MSOLVE takes advantage of its tracer. Also, maple's documentation indicates that on some examples, an algorithm computing a so-called rational univariate representation (preserving multiplicities), instead of FGLM, may be used. It is likely that on most examples we tried, FGLM is not used. Note also that maple lifts a whole `<DRL>`-Gröbner basis over $\mathbb{Q}$ while MSOLVE avoids this step. Furthermore, our tracer shares some similarities with [86].

There are, of course, few examples, where MSOLVE is not competitive. In particular, for some systems, MSOLVE may need to introduce a generic linear form as previously explained while a rational parametrization can be obtained without it (but up to computing normal forms). Also some systems admit a triangular representation, and/or can be split. It seems that maple can detect and sometimes take advantage of such situations. This is typically the case for the Noon-n examples.
In this chapter, I present some works that I have done on algorithms for computing Gröbner bases.

First, Section 3.1 is dedicated to change of order algorithms, whether my coauthors and I studied its complexity in some special cases or we designed new algorithms. In Section 3.2, we shall look at the problems of saturating or quotienting ideals. From a geometric point of view, saturating comes down to removing a subvariety and then taking the Zariski closure of the remaining set. Then, in Section 3.3, I present an efficient verification of the computation of the parametrization of the radical of an ideal when we cannot ensure that this radical is in shape position.

This chapter is based on joint work with Alin Bostan<fnref n="1" />, Christian Eder<fnref n="2" />, Andrew Ferguson<fnref n="3" />, Vincent Neiger<fnref n="4" /> and Mohab Safey El Din<fnref n="4" />, [10, 13, 15, 20].

## 3.1 Gröbner bases change of order algorithms

Change of order algorithms can be divided into two categories, depending on the dimension of the ideal $I$. In positive dimension, we can mention the Gröbner walk [38] whose goal is to update step by step the Gröbner basis when “walking” along a path from $<1$ to $<2$ in the Gröbner fan of the ideal. In dimension zero, most algorithms rely on linear algebra in the finite-dimensional quotient algebra, following the idea of [50] yielding the so-called FGLM algorithm.

This original algorithm had a complexity $O(nD^3)$, where $D$ is the degree of the ideal, i.e. the dimension of the algebra over the base field. Let us recall that if the ideal is spanned by $n$ generic polynomials of degree $d$, then $D = d^n$. Under the assumption that the change of order is from $<_{\text{DRL}}$ to $<_{\text{LEX}}$ in shape position and some other genericity assumptions, the complexity of the change of order has been improved in [46] to $\tilde{O}(D^\omega)$, where $\omega$ is the matrix multiplication exponent. More recently, Neiger and Schost [83] proposed a general change of order algorithm with complexity $\tilde{O}(D^\omega)$ that only requires a stability assumption on the $<_1$-Gröbner basis and no assumption on the $<_2$-Gröbner basis. This stability assumption is for all monomial $\mu$,

$$ \forall 1 \le i < j \le n, \quad x_j\mu \in \text{LM}_{<1}(I) \implies x_i\mu \in \text{LM}_{<1}(I). $$

### 3.1.1 Complexity of SPARSE-FGLM

Together, Propositions 2.7 and 2.8 and Theorem 2.9 yield the complexity of the SPARSE-FGLM algorithm in the shape position situation. While it is not easy to estimate exactly the number $s$ of nonzero coefficients of the matrix $M_{x_n}$, we can give an upper bound of it. Indeed, from the procedure of Lemma 2.5, $s \in O(tD)$, where $t$ is the number of normal forms to perform. Using now Lemma 2.6, we obtain that for an ideal spanned by $n$ generic polynomials of given degrees, $t$ is exactly the number of polynomials in $\mathcal{G}_{\text{DRL}}$ whose leading monomials are divisible by $x_n$. Therefore, one obtains a bound of the complexity of the SPARSE-FGLM algorithm based on the input Gröbner basis.

This opens a new branch of research to estimate $t$ in order to have a finer complexity estimate of the SPARSE-FGLM algorithm in different situations. In their paper [53, Cor. 5.11], the authors bound $t$ for generic systems of $n$ polynomials of degree $d$ in $\mathbb{K}[x]$ and find that as $d$ tends to $+\infty$, $t$ grows as $\sqrt{\frac{6}{n\pi}}d^{2n-1}$. This is to compare with $D$ that grows as $d^n$. More recently, new estimates for $t$ have been given, especially for polynomial systems that have some of their generators as minors of a multivariate polynomial matrix.

<footnote n="1">¹INRIA</footnote>

<footnote n="2">²Rheinland-Pfälzische Technische Universität Kaiserslautern-Landau</footnote>

<footnote n="3">³former Ph.D. student at Sorbonne Université</footnote>

<footnote n="4">⁴Sorbonne Université</footnote>
In a joint work with Alin Bostan<fnref n="5" />, Andrew Ferguson<fnref n="6" /> and Mohab Safey El Din<fnref n="7" />, we study a case of polynomial ideals that encompass systems defining the critical points of the restriction of a linear map to an algebraic set defined by p generic polynomials by means of the simultaneous vanishing of maximal minors of a truncated Jacobian matrix. For instance, we consider $f_1, \dots, f_m \in K[x]$, the Jacobian matrix

$$ \mathcal{J} = \begin{pmatrix} 1 & 0 & \cdots & 0 \\ \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial f_m}{\partial x_1} & \frac{\partial f_m}{\partial x_2} & \cdots & \frac{\partial f_m}{\partial x_n} \end{pmatrix} $$

of $\varphi_1$, the projection on the $x_1$-axis, $f_1, \dots, f_m$ and the ideal $I = \langle f_1, \dots, f_m \rangle + \langle \text{MaxMinors} \mathcal{J} \rangle$.

**Definition 1.8 (see also [10, Def. 8]):** With $K$ an infinite field, let $I \subset K[x]$ be an ideal which is the sum of $m$ polynomials of degree at most $d$ and the maximal minors of a matrix with polynomial entries also of degree at most $d$. We say that $I$ is a *generic determinantal sum ideal* if the following three conditions hold:

* the ideal $I$ is zero-dimensional and in shape position;

* the Hilbert series of $K[x]/I$ is

$$ H(\tau) = \frac{\det(M(\tau^{d-1})) (1-\tau^d)^m (1-\tau^{d-1})^{n-m}}{\tau^{(d-1)(m-1)} (1-\tau)^n} $$

where $M(\tau)$ is the $(m-1) \times (m-1)$-matrix whose $(i,j)$th entry is $\sum_k \binom{m-i}{k} \binom{n-1-j}{k} \tau^k$;

* for all $e \ge 1$, the Hilbert series of $(K[x]/I) / \langle x_n^e \rangle$ is equal to the series $(1-\tau)H$ truncated at the first non-positive coefficient.

Under some regularity assumptions, the ideal defining the set of critical values of a generic polynomial map restricted to a smooth algebraic set falls in this class, see [10]. Thus, determining an asymptotic for $t$ for such systems allows us to estimate the complexity of the SPARSE-FGLM algorithm in this situation and in particular for computing the critical values of the restriction of a generic polynomial map.

Relying on manipulations of the Hilbert series given in Definition 1.8, we simplified this expression in order to show that this Hilbert series is *unimodal*, which allowed us to prove the following theorem.

**Theorem 3.1 (see also [10, Th. 1]):** Let $I$ be a generic determinantal sum ideal of $K[x]$ as in Definition 1.8. Let $\mathscr{G}_{\text{DRL}}$ be the reduced $<_{\text{DRL}}$-Gröbner basis of $I$. Then, the latter case of Lemma 2.5 only happens if there exists $g \in \mathscr{G}_{\text{DRL}}$ such that $\text{LM}_{<_{\text{DRL}}}(g) = x_n\sigma_j$. As a consequence, one has $\text{NF}(x_n\sigma_j, \mathscr{G}_{\text{DRL}}, <_{\text{DRL}}) = x_n\sigma_j - g$.

As a consequence, the matrix of the multiplication by $x_n$ is obtained for free from the $<_{\text{DRL}}$-Gröbner basis, as in Lemma 2.6 for generic systems. Taking this structure into account, we then prove that this number $t$ is equal to the largest coefficient of the Hilbert series. Finally, using combinatorial techniques, we obtain formulae for $t$ for such generic determinantal sum ideals.

**Theorem 1.9 (see also [10, Th. 2]):** Let $I$ be a generic determinantal sum ideal of $K[x]$ as in Definition 1.8. Let $\mathscr{G}_{\text{DRL}}$ be the reduced $<_{\text{DRL}}$-Gröbner basis of $I$. Let $t$ be the number of polynomials in $\mathscr{G}_{\text{DRL}}$ whose leading monomial is divisible by $x_n$.

Then, for $d=2$ and $n \ge m$,

$$ t = \sum_{k=0}^{m-1} \binom{n-m-1+k}{k} \binom{m}{\lfloor 3m/2 \rfloor - 1 - j} $$

<footnote n="5">⁵INRIA</footnote>

<footnote n="6">⁶former Ph.D. student at Sorbonne Université</footnote>

<footnote n="7">⁷Sorbonne Université</footnote>
Moreover, for $d \ge 3$ and $n \to \infty$,

$$t \approx \frac{1}{\sqrt{(n-m)\pi}} \sqrt{\frac{6}{(d-1)^2-1}} d^m (d-1)^{n-m} \binom{n-2}{m-1}.$$

In Table 3.1, we show the practical accuracy of the formulae in Theorem 1.9, for the parameter $t$. For $d=2$ we use our exact formula, while for $d \ge 3$ we use the asymptotic formula. The column *Actual* is the true density of the matrix of the multiplication by $x_n$ for $<_{\text{DRL}}$. The column *Theoretical* is an approximation given by the ratio $t/D$ and the column *Asymptotic* is another one given by the asymptotic of $t$ divided by $D$. The column *Actual* was already present in [53, Tab. 2], while the other two columns were added in [10]. Let us recall that in this setting, $D = d^m (d-1)^{n-m} \binom{n-1}{m-1}$, see [84, Th. 2.2].

<table><caption>Table 3.1: Density of matrix of the multiplication by $x_n$ for generic critical point systems</caption><thead><tr><th rowspan="2">Parameters<br>(d, m, n)</th><th rowspan="2">Degree<br>D</th><th colspan="3">Matrix Density</th></tr><tr><th>Actual</th><th>Theoretical</th><th>Asymptotic</th></tr></thead><tbody><tr><td>(2, 4, 9)</td><td>896</td><td>30.17%</td><td>30.80%</td><td>30.80%</td></tr><tr><td>(2, 4, 10)</td><td>1344</td><td>31.13%</td><td>31.77%</td><td>31.77%</td></tr><tr><td>(2, 4, 11)</td><td>1920</td><td>31.86%</td><td>32.50%</td><td>32.50%</td></tr><tr><td>(3, 3, 6)</td><td>2160</td><td>17.52%</td><td>18.52%</td><td>27.73%</td></tr><tr><td>(3, 3, 7)</td><td>6480</td><td>17.39%</td><td>18.31%</td><td>26.62%</td></tr><tr><td>(3, 3, 8)</td><td>18144</td><td>17.63%</td><td>18.72%</td><td>25.50%</td></tr><tr><td>(4, 2, 5)</td><td>1728</td><td>14.46%</td><td>15.45%</td><td>21.24%</td></tr><tr><td>(4, 2, 6)</td><td>6480</td><td>14.11%</td><td>15.13%</td><td>19.56%</td></tr><tr><td>(5, 2, 5)</td><td>6400</td><td>11.00%</td><td>11.94%</td><td>15.47%</td></tr><tr><td>(6, 2, 5)</td><td>18000</td><td>8.80%</td><td>9.63%</td><td>12.22%</td></tr></tbody></table>

Let us also mention that during his Ph.D., Andrew Ferguson<fnref n="8" />, with Le, extended this analysis of the asymptotic of $t$ for ideals spanned by $(r+1)$-minors of a polynomial symmetric matrix $(f_{i,j})_{1 \le i,j \le m}$, where the $f_{i,j}$'s are in $n = \binom{m-r+1}{2}$ variables and of degree $d$, see [58]. As above, their asymptotics are derived from the unimodality of the Hilbert series. As the Hilbert series is only known in the cases $r=m-2$, $r=m-1$ and $r=1$, their asymptotics hold in these cases. For the other cases, they hold if the Hilbert series is rightfully conjectured to be unimodal.

### 3.1.2 A polynomial-matrix algorithm

Under the stability assumption, two algorithms, more efficient than the original FGLM algorithm, are at our disposal. The Faugère and Mou's SPARSE-FGLM algorithm [52, 53], which requires that the reduced $<_{\text{LEX}}$-Gröbner basis is in shape position and whose complexity is in $O(tD)$ and Neiger and Schost's algorithm [83] whose complexity is in $\tilde{O}(D^\omega)$.

Despite the practical efficiency of the SPARSE-FGLM algorithm and the asymptotics on $t$ [10, 53], the recent improvements on $\omega$, the best known bounds are $\omega \le 2.3728596$ [1] or even $\omega \le 2.37188$ [43], make Neiger and Schost's algorithm the fastest.

In this subsection, I present a change of order algorithm designed by Vincent Neiger<fnref n="9" />, Mohab Safey El Din<fnref n="9" /> and myself from $<_{\text{DRL}}$ to $<_{\text{LEX}}$ with complexity $\tilde{O}(t^{\omega-1}D)$, assuming that the $<_{\text{DRL}}$-Gröbner basis satisfies the stability assumption and that the $<_{\text{LEX}}$-Gröbner basis is in shape position. Since $t \le D$, this complexity makes it asymptotically faster than the two aforementioned algorithms, independently on the value of $\omega$ and the ratio $t/D$.

To do so, we push forward the study of the properties of the multiplication matrix $M_{x_n}$ by $x_n$. Let $S_{\text{DRL}}$ be the monomial basis of $K[x]/I$ obtained from the reduced $<_{\text{DRL}}$-Gröbner basis $\mathcal{G}_{\text{DRL}}$ of $I$. Under the stability assumption, its columns are either unit vectors (for those $\mu \in S_{\text{DRL}}$ such that $x_n\mu \in S_{\text{DRL}}$) or vectors of coefficients of polynomials in $\mathcal{G}_{\text{DRL}}$ (for those $\mu \in S_{\text{DRL}}$ such that $x_n\mu$ is a leading monomial of one element in $\mathcal{G}_{\text{DRL}}$). The latter ones are usually referred to as a “dense”

<footnote n="8">⁸former Ph.D. student at Sorbonne Université</footnote>

<footnote n="9">⁹Sorbonne Université</footnote>
columns, [53, Sec. 5] and [101, Sec. 4]. This is a well-known matrix structure in $K$-linear algebra, called a *shifted form* in [88, 89], and studied in particular in the context of the computation of the characteristic polynomial or the Frobenius normal form of a matrix over $K$, see [102, Sec. 9.1].

Exploiting the algebraic structure itself, we relate it to operations in a $K[x_n]$-submodule of $I$. Following a classical construction in [102, Sec. 9.1], instead of the multiplication matrix $M_{x_n}$ which is in $K^{D \times D}$, we consider a univariate polynomial matrix $P_{x_n}$ in $K[x_n]^{l \times l}$ whose average column degree is $D/t$. This polynomial matrix can be seen as a “compression” of $M_{x_n}$, or more precisely of the characteristic matrix $x_n Id_D - M_{x_n}$, with smaller matrix dimension but larger degrees.

Under a less restrictive assumption than the stability one, and which is also included in the assumptions of the SPARSE-FGLM algorithm of [52, 53] and of [46], we can build $P_{x_n}$ for free from $\mathscr{G}_{\text{DRL}}$. Then, we retrieve $\mathscr{G}_{\text{LEX}}$ using $K[x_n]$-linear algebra operations on $P_{x_n}$.

**Theorem 1.10 (see also [20, Th. 1.1]):** Let $I$ be a zero-dimensional ideal $K[x]$ of degree $D$. Let $\mathscr{G}_{\text{DRL}}$ (resp. $\mathscr{G}_{\text{LEX}}$) be the reduced $<_{\text{DRL}}$- (resp. $<_{\text{LEX}}$-)Gröbner basis of $I$ and $S_{\text{DRL}}$ be the $<_{\text{DRL}}$-monomial basis of $K[x]/I$. Assume that $x_1, \dots, x_{n-1}$ are in $S_{\text{DRL}}$, and that for all monomials $\mu \in S_{\text{DRL}}$, either $x_n \mu$ is in $S_{\text{DRL}}$ or it is the $<_{\text{DRL}}$-leading monomial of an element in $\mathscr{G}_{\text{DRL}}$. Assume that $I$ is in shape position. Then, one can compute $\mathscr{G}_{\text{LEX}}$ using $\tilde{O}(t^{\omega^{-1}D})$ operations in $K$, where $t$ is the number of elements of $\mathscr{G}_{\text{DRL}}$ whose $<_{\text{DRL}}$-leading monomial is divisible by $x_n$.

**Example 3.2:** Let $\mathscr{H}_{\text{DRL}} = \{x_2^2 - x_1 + x_2, x_1x_2 - x_1 + x_2, x_1^2 - x_1 + x_2\}$ be the $<_{\text{DRL}}$-Gröbner basis of $J$. Then, the associated staircase is $\{1, x_2, x_1\}$ and $D = 3$. One can read the matrix of the multiplication by $x_2$ in $K[x_1, x_2]/J$ in this monomial basis directly from $\mathscr{H}_{\text{DRL}}$. It is

$$M_{x_2} = \begin{pmatrix} 1 & x_2 & x_1 \\ 0 & 0 & 0 \\ \frac{1}{0} & \frac{-1}{1} & \frac{-1}{1} \end{pmatrix}_{x_2 \in K^{D \times D}}$$

where the $i$th column is the image of the $i$th monomial of the staircase by its multiplication by $x_2$.

The univariate polynomial matrix $P_{x_2}$ is obtained from the polynomials in $\mathscr{H}_{\text{DRL}}$ whose leading monomials are multiples of $x_2$, i.e. the first $t = 2$ polynomials of $\mathscr{H}_{\text{DRL}}$, and decomposed in the $K[x_2]$-module $K[x_2] + x_1 K[x_2]$. It is

$$P_{x_2} = \begin{pmatrix} x_2^2 + x_2 & x_2 \\ -1 & x_2 - 1 \end{pmatrix}_{x_1}$$

Its Hermite normal form is

$$\begin{pmatrix} x_2^3 & -x_2^2 - x_2 \\ 0 & 1 \end{pmatrix}_{x_1}$$

which allows us to read that the polynomials $x_2^3$ and $x_1 - x_2^2 - x_2$ are both in $J$ and $K[x_2] + x_1 K[x_2]$. Thus, they form the $<_{\text{LEX}}$-Gröbner basis of $J$.

## 3.2 Saturated and colon ideals

Let $I$ be an ideal of $K[x]$ and $\varphi \in K[x]$. The colon and saturation ideals of $I$ w.r.t. $\varphi$ are defined as

$$I : \langle \varphi \rangle = \{h \in K[x] \mid h\varphi \in I\}, \quad I : \langle \varphi \rangle^{\infty} = \{h \in K[x] \mid \exists k \in \mathbb{N}, h\varphi^k \in I\}.$$

Let us recall that by [41, Chap. 4], the algebraic set $V(I : \langle \varphi \rangle^{\infty}) \subset \overline{K}^n$ is the Zariski closure of the set difference $V(I) \setminus V(\varphi)$.

Computing algebraic representations of saturated ideals arises in many applications ranging from experimental mathematics to engineering sciences (see [33, 60, 87]) since some natural algebraic modelings come with parasite solutions which one excludes through some saturation process. For instance, modeling that some $(p \times q)$-matrix with polynomial entries has rank $r$ through the simultaneous vanishing of its $(r+1)$-minors will include those points at which the matrix has rank less than $r$.
Since

$$
\begin{align*}
I : \langle \varphi \rangle &= \frac{1}{\varphi} (I \cap \langle \varphi \rangle) = (tI + \langle (1-t)\varphi \rangle) \cap \mathbb{K}[\mathbf{x}] \\
I : \langle \varphi \rangle^{\infty} &= (I + \langle 1-t\varphi \rangle) \cap \mathbb{K}[\mathbf{x}],
\end{align*}
$$

using Rabinowitsch's trick [90] and [41, Chap. 4, Sec. 4, Th. 14, (ii)], one can compute Gröbner
bases of these ideals choosing a monomial order eliminating *t* and keeping all polynomials not
involving *t*, see also [41, Chap. 3, Sec. 1, Th. 2 and Ex. 6].

This section is about the computational problem of computing a Gröbner basis associated to
$I : \langle \varphi \rangle^{\infty}$ or $I : \langle \varphi \rangle$ without introducing an extra variable. This is joint work with Christian Eder<fnref n="10" />
and Mohab Safey El Din<fnref n="11" />.

**3.2.1 Bayer’s algorithm**

If *I* is homogeneous, i.e. it is spanned by a set of homogeneous polynomials, Bayer’s algorithm [4]
allows one to compute *I*: <*x*<sub>*n*</sub>>*<sup>∞</sup>*. If it is not, then one can still recover a Gröbner basis of *I*: <*x*<sub>*n*</sub>>*<sup>∞</sup>*
using Algorithm 3.1, still called Bayer’s algorithm:

Input: A list of polynomials $f_1, \dots, f_s$ spanning an ideal $I \subseteq \mathbb{K}[\mathbf{x}]$.

Output: A <DRL>-Gröbner basis of I : ⟨x_n⟩_∞.

1. Homogenize $f_1, \dots, f_s$ with variable $x_0$ to obtain $f_1^h, \dots, f_s^h$.

2 Compute $G$, the <DRL>-Gröbner basis of $\langle f_1^h, \dots, f_s^h \rangle$ with $x_n < \text{DRL} \dots < \text{DRL} x_1 < \text{DRL} x_0$.

3 For each $g$ in $G^h$ do

4. While $x_n \mid g$ do $g := g/x_n$.

5. Set $x_0$ to 1 in $g$.

6 Return $G$.

**Algorithm 3.1:** Bayer’s algorithm

When $\varphi \neq x_n$, one introduces a slack variable $x_{n+1}$, computes the saturation of $I + \langle x_{n+1} - \varphi \rangle$
w.r.t. $x_{n+1}$ and eliminate $x_{n+1}$.

Rabinowitsch's trick and Bayer's algorithm constitute the state-of-the-art algorithms for computing saturations of ideals. Note that they do not take advantage of intermediate data obtained during the Gröbner basis computations since these are used as black boxes.

**3.2.2 F<sub>4</sub>SAT, an algorithm for saturated ideals**

After the first step of the F<sub>4</sub> algorithm in degree *d*, if no new polynomial of degree at most *d* is discovered, then the current Gröbner basis *G* is a *d*-truncated *<*-Gröbner basis of *I*. Therefore, we have a partial information on the staircase of *I*, and thus of *I*: ⟨*φ*⟩<sup>∞</sup>, for *<* since we know monomials that are outside of this staircase. The F<sub>4</sub>SAT algorithm searches for polynomials in *I*: ⟨*φ*⟩<sup>∞</sup> whose supports are entirely included in the given staircase using the fact that (*I*: ⟨*φ*⟩<sup>∞</sup>) : ⟨*φ*⟩ = *I*: ⟨*φ*⟩<sup>∞</sup>. If new polynomials are found, they are added to *G* and the necessary critical pairs are added to the set of pairs to handle. Then, we resume the F<sub>4</sub> algorithm.

The search of new polynomials is done through linear algebra computations. From a *d*-truncated
Gröbner basis of an ideal *J*, *I* ⊆ *J* ⊆ *I*: ⟨*φ*⟩<sup>∞</sup>, we compute a bound *B* on the degree of the
polynomials in the reduced Gröbner basis of *J*: ⟨*φ*⟩<sup>∞</sup> = *I*: ⟨*φ*⟩<sup>∞</sup> using the ComputeMaxDegree
subroutine based on [79, Sec. 4.5, Cor.]. Then, we compute NF(*σφ*, *G*, *<*) for all monomials *σ* in
the associated staircase *S* of degree at most *B*. Finally, we search for vanishing linear combinations
thereof. Indeed, if

$$
\mathrm{NF}(s\varphi, \mathcal{G}, <) - \sum_{\substack{\sigma \in S_d \\ \sigma < s}} c_{\sigma} \mathrm{NF}(\sigma\varphi, \mathcal{G}, <) = 0,
$$

then $(s - \sum_{\sigma \in S_d, \sigma < s} c_{\sigma} \sigma) \varphi \in J$. This yields Algorithm 3.2.

**Theorem 3.3 ([15, Th. 3.1]):** Let $f_1, \dots, f_s$ be a generating family of an ideal $I \subseteq \mathbb{K}[\mathbf{x}]$, $\varphi \in \mathbb{K}[\mathbf{x}]$ be a polynomial and $<^*$ be a total degree monomial order. Then, Algorithm 3.2 terminates and returns

<footnote n="10"><sup>10</sup>Rheinland-Pfälzische Technische Universität Kaiserslautern-Landau</footnote>

<footnote n="11">¹¹Sorbonne Université</footnote>
**Input:** A list of polynomials $f_1, \dots, f_s$ spanning an ideal $I \subseteq \mathbb{K}[x]$, a polynomial $\varphi \in \mathbb{K}[x]$ and a total degree monomial order $<$.

**Output:** A $<$-Gröbner basis $\mathcal{G}$ of $I : \langle \varphi \rangle^{\infty}$.

1 $\mathcal{G} := \{f_1, \dots, f_s\}$.

2 $b := \text{true}$

3 $P := \{(f_i, f_j) \mid 1 \le i < j \le s\}$

4 While $P \neq \emptyset$ do

5 Choose a subset $L$ of $P$.

6 $P := P \setminus L$.

7 $L := \text{SymbolicPreprocessing}(L, \mathcal{G})$.

8 $L := \text{LinearAlgebra}(L)$.

9 For $h \in L$ with $\text{LM}_<^t(h) \notin \langle \text{LM}_<^t(\mathcal{G}) \rangle$ do

10 $P := \{(g, h) | g \in \mathcal{G}\}$.

11 $\mathcal{G} := \mathcal{G} \cup \{h\}$.

12 $b := \text{true}$.

13 If $b$ then

14 $b := \text{false}$.

15 $B := \text{ComputeMaxDegree}(\mathcal{G}, \varphi)$.

16 For $\sigma \notin \text{LM}_<^t(\mathcal{G})$ and $\text{deg} \sigma \le B$ do

17 $\lfloor q_{\sigma} := \text{NF}(\sigma\varphi, \mathcal{G}, <) \rfloor$.

18 Build the matrix $M$ whose rows are given by polynomials $q_{\sigma}$ and columns by each monomials in their support in decreasing order.

19 Compute a lower triangular basis $K$ of the left-kernel of $M$.

20 For each $k \in K$ do

21 $h := \sum_{\sigma \in \langle \text{LM}_<^t(\mathcal{G}) \rangle} k_{\sigma} \sigma$.

22 $P := \{(g, h) | g \in \mathcal{G}\}$.

23 $\mathcal{G} := \mathcal{G} \cup \{h\}$.

24 $b := \text{true}$.

25 Return $\mathcal{G}$.

// tracks if $\mathcal{G}$ has changed

// new information on $\langle \text{LM}_<^t(I : \langle \varphi \rangle^{\infty}) \rangle$

// bounds the degrees in the sought Gröbner basis

// the polynomial whose vector of coefficients is $k$

<caption>Algorithm 3.2: F<sub>4</sub>SAT</caption>
a <Gröbner basis of $I : \langle\varphi\rangle^{\infty}$.

**Example 3.4:** Let us call the $F_4$SAT algorithm on the reduced <DRL-Gröbner basis $\mathscr{G}_{\text{DRL}} = \{(x_2^2 - x_1 + x_2)(x_1^2 + x_2^2 - 1), (x_1^2 - x_1x_2 - x_2^2 + x_1 - x_2)(x_1^2 + x_2^2 - 1)\}$ and $\varphi = x_1^4 - x_2^4 - 2x_1^2 + 1$ in order to compute the reduced <DRL-Gröbner basis of $I : \langle\varphi\rangle^{\infty}$ with

$$I = \langle (x_2^4 + x_2^3)(x_1^2 + x_2^2 - 1), (x_1 - x_2^2 - x_2)(x_1^2 + x_2^2 - 1) \rangle.$$

From $\mathscr{G}_{\text{DRL}}$, we have the staircase $S_{\text{DRL}} = \{1, x_2, x_1, x_2^2, x_1x_2, x_1^2, x_2^3, x_1x_2^2, x_1^3, x_2^4, x_1x_2^5, x_2^6, x_1x_2^7, \dots\}$. Since this staircase is infinite, we first need a degree bound for the considered monomials in $S_{\text{DRL}}$. Let us choose the maximal degree in $\mathscr{G}$ minus one, i.e. 3.

The normal form of $h\varphi$ w.r.t. $\mathscr{G}$ and <DRL, where $h$ of degree 3 has its support in $S_{\text{DRL}}$ is 0, if, and only if,

$$h = c_3(x_2^2 - x_1 + x_2) + c_5(x_1^2 - x_1x_2) + c_6(x_2^3 - x_1x_2 + x_1 - x_2) + c_7x_1x_2^2 + c_8x_1^2x_2 + c_9x_1^3.$$

Thus, we have found new polynomials in $I : \langle\varphi\rangle^{\infty}$. We build the ideal $I'$, spanned by $I$ and these polynomials, and compute its reduced <DRL-Gröbner basis $\mathscr{G}'_{\text{DRL}} = \{x_2^2 - x_1 + x_2, x_1^2 - x_1x_2\}$. The new associated staircase is $S'_{\text{DRL}} = \{1, x_2, x_1, x_1x_2\}$, which is now finite.

The normal form of $h\varphi$ w.r.t. $\mathscr{G}'_{\text{DRL}}$ and <DRL, where $h$ has its support in $S'_{\text{DRL}}$, is 0 if, and only if,

$$h = c_3(x_1x_2 - x_1 + x_2).$$

We have found a new polynomial $x_1x_2 - x_1 + x_2$ and compute the reduced <DRL-Gröbner basis $\mathscr{H}_{\text{DRL}} = \{x_2^2 - x_1 + x_2, x_1x_2 - x_1 + x_2, x_1^2 - x_1 + x_2\}$ of $J = I' + \langle x_1x_2 - x_1 + x_2 \rangle$. The new associated staircase is $T_{\text{DRL}} = \{1, x_2, x_1\}$ and we cannot find a nonzero polynomial $h$ with support in $T_{\text{DRL}}$ whose normal form w.r.t. $\mathscr{H}_{\text{DRL}}$ and <DRL is 0. Therefore $\mathscr{H}_{\text{DRL}}$ is the reduced <DRL-Gröbner basis of $I : \langle\varphi\rangle^{\infty} = J = \langle x_2^3, x_1 - x_2^2 - x_2 \rangle$.

In Table 3.2, we present timings for different systems with coefficients in $\mathbb{Q}$. The $F_4$SAT algorithm is implemented in `MSOLVE` with a multi-modular approach using a tracer. This tracer requires two learning phases modulo two different primes before being able to run its apply phase modulo subsequent primes. We compare it with Rabinowitsch's trick implemented in `MSOLVE` either with probabilistic linear algebra or with the tracer, and with `MAPLE` and `SINGULAR`. Most of the time in the first learning phase is due to the last saturation step, i.e. checking in line 19 that no new polynomial are found in the saturated ideal. As we can see, the `apply` phase of $F_4$SAT is the fastest implementation by an order of magnitude. The `SOS` systems are obtained by making a polynomial $f$ as the sum of the squares of $p$ random polynomials of degree $d$. Then, we consider the ideals $\langle f, \frac{\partial f}{\partial x_1}, \dots, \frac{\partial f}{\partial x_s} \rangle$ or $\langle \frac{\partial f}{\partial x_1}, \dots, \frac{\partial f}{\partial x_s} \rangle$, for $s \le n-1$, saturated by $\frac{\partial f}{\partial x_n}$.

### 3.2.3 SPARSE-FGLM-COLON, an algorithm for colon ideals

For an ideal $I = \langle f_1, \dots, f_s \rangle$ and a polynomial $\varphi$, the colon ideal of $I$ by $\varphi$ is

$$I : \langle\varphi\rangle = \{h | h\varphi \in I\}.$$

This ideal contains $I$ and is included in $I : \langle\varphi\rangle^{\infty}$ defined in Section 3.2.2. In fact, there exists an integer $N$ such that $I : \langle\varphi^N\rangle = I : \langle\varphi^{N+1}\rangle = \dots = I : \langle\varphi\rangle^{\infty}$, see [41, Chap. 4].

Let $\mathscr{G}_{\text{DRL}}$ be the reduced <DRL-Gröbner basis. We assume that the colon ideal $I : \langle\varphi\rangle$ is zero-dimensional, thus $\varphi \notin I$, and that its reduced <LEX-Gröbner basis $\mathscr{H}_{\text{LEX}}$ is in shape position: There exist $h_1, \dots, h_n \in \mathbb{K}[x_n]$, with $\deg h_k < \deg h_n$ for $1 \le k \le n-1$, such that

$$\mathscr{H}_{\text{LEX}} = \{h_n(x_n), x_{n-1} - h_{n-1}(x_n), \dots, x_1 - h_1(x_n)\}.$$

With my co-authors, Christian Eder<fnref n="12" /> and Mohab Safey El Din<fnref n="13" />, we designed a new algorithm in [15] for computing $\mathscr{H}_{\text{LEX}}$ from $\mathscr{G}_{\text{DRL}}$, even when $I$ is positive-dimensional. Our approach is to build a matrix $\tilde{M}_{x_n}$ so that applying Wiedemann's algorithm allows us to recover $\mathscr{H}_{\text{LEX}}$, similarly to the SPARSE-FGLM algorithm [52, 53], see also Algorithms 2.2 and 2.3.

<footnote n="12">¹²Rheinland-Pfälzische Technische Universität Kaiserslautern-Landau</footnote>

<footnote n="13">¹³Sorbonne Université</footnote>
Chapter 3. Computing Gröbner bases

3.2. Saturated and colon ideals

<table><caption>Table 3.2: Timings in seconds, $<_{\text{DRL}}$-Gröbner bases of saturated ideals.</caption><thead><tr><th rowspan="2">Sys-SOS</th><th colspan="3">F<sub>4</sub>SAT</th><th colspan="2">MSOLVE</th><th colspan="2">MSOLVE</th><th rowspan="2">Maple (prob.)</th><th rowspan="2">Singular</th></tr><tr><th>(learn1)</th><th>(learn2)</th><th>(apply)</th><th>(prob.)</th><th>(learn)</th><th>(apply)</th></tr></thead><tbody><tr><td colspan="10" style="text-align:center; font-weight:bold;">positive-dimensional to 0-dimensional</td></tr><tr><td>Steiner</td><td>115</td><td>134</td><td>67.2</td><td>204</td><td>614</td><td>153</td><td>239</td><td>3.642</td><td>8.226</td></tr><tr><td>d<sub>3</sub>-n<sub>6</sub>-p<sub>3</sub></td><td>82.4</td><td>127</td><td>56.7</td><td>51.5</td><td>191</td><td>32.6</td><td>67.4</td><td>-</td><td>-</td></tr><tr><td>d<sub>3</sub>-n<sub>6</sub>-p<sub>4</sub></td><td>1.502</td><td>1.776</td><td>810</td><td>2.123</td><td>5.284</td><td>1.720</td><td>3.585</td><td>-</td><td>-</td></tr><tr><td>d<sub>3</sub>-n<sub>6</sub>-p<sub>5</sub></td><td>9.046</td><td>7.032</td><td>3.321</td><td>7.485</td><td>10.711</td><td>6.466</td><td>7.256</td><td>-</td><td>-</td></tr><tr><td>d<sub>4</sub>-n<sub>6</sub>-p<sub>2</sub></td><td>720</td><td>1.581</td><td>536</td><td>120</td><td>520</td><td>60.6</td><td>135</td><td>24.532</td><td>-</td></tr><tr><td>d<sub>3</sub>-n<sub>6</sub>-p<sub>3</sub></td><td>45.749</td><td>38.657</td><td>18.123</td><td>40.646</td><td>190.009</td><td>35.835</td><td>38.666</td><td>-</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>7</sub>-p<sub>6</sub></td><td>18.9</td><td>41.85</td><td>10.8</td><td>31.8</td><td>101</td><td>19.5</td><td>41.4</td><td>1.773</td><td>-</td></tr><tr><td>d<sub>3</sub>-n<sub>7</sub>-p<sub>2</sub></td><td>28.2</td><td>45.2</td><td>23.9</td><td>5.02</td><td>11.6</td><td>2.63</td><td>8.09</td><td>9.61</td><td>-</td></tr><tr><td>d<sub>3</sub>-n<sub>7</sub>-p<sub>3</sub></td><td>1.462</td><td>2.688</td><td>937</td><td>953</td><td>5.851</td><td>875</td><td>1.108</td><td>-</td><td>-</td></tr><tr><td>d<sub>3</sub>-n<sub>7</sub>-p<sub>4</sub></td><td>48.907</td><td>65.035</td><td>22.844</td><td>40.383</td><td>24.8889</td><td>34.670</td><td>39.729</td><td>-</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>8</sub>-p<sub>4</sub></td><td>2.68</td><td>5.04</td><td>1.89</td><td>3.55</td><td>10.1</td><td>2.02</td><td>4.12</td><td>500</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>8</sub>-p<sub>5</sub></td><td>47.7</td><td>171.9</td><td>37.1</td><td>62.9</td><td>270</td><td>45.3</td><td>48.8</td><td>8.333</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>8</sub>-p<sub>6</sub></td><td>287</td><td>820</td><td>169</td><td>420</td><td>1.599</td><td>391</td><td>350</td><td>54.597</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>8</sub>-p<sub>7</sub></td><td>1.018</td><td>1.841</td><td>442</td><td>997</td><td>3.198</td><td>683</td><td>871</td><td>-</td><td>-</td></tr><tr><td>d<sub>3</sub>-n<sub>8</sub>-p<sub>2</sub></td><td>300</td><td>585</td><td>266</td><td>324</td><td>105</td><td>20.4</td><td>50.7</td><td>9.812</td><td>-</td></tr><tr><td>d<sub>3</sub>-n<sub>8</sub>-p<sub>3</sub></td><td>18.152</td><td>42.436</td><td>11.285</td><td>15.592</td><td>71.595</td><td>8.478</td><td>15.182</td><td>-</td><td>-</td></tr><tr><td colspan="10" style="text-align:center; font-weight:bold;">positive-dimensional to 1-dimensional</td></tr><tr><td>d<sub>3</sub>-n<sub>6</sub>-p<sub>2</sub></td><td>1.31</td><td>0.41</td><td>0.31</td><td>0.77</td><td>2.40</td><td>0.40</td><td>1.12</td><td>52.2</td><td>-</td></tr><tr><td>d<sub>3</sub>-n<sub>6</sub>-p<sub>3</sub></td><td>43.7</td><td>5.55</td><td>1.84</td><td>25.2</td><td>142</td><td>16.6</td><td>35.4</td><td>2.902</td><td>-</td></tr><tr><td>d<sub>3</sub>-n<sub>6</sub>-p<sub>4</sub></td><td>533</td><td>53.1</td><td>19.7</td><td>171</td><td>882</td><td>126</td><td>223</td><td>39.591</td><td>-</td></tr><tr><td>d<sub>3</sub>-n<sub>6</sub>-p<sub>5</sub></td><td>1.863</td><td>184</td><td>104</td><td>276</td><td>1.145</td><td>183</td><td>394</td><td>42.854</td><td>-</td></tr><tr><td>d<sub>4</sub>-n<sub>6</sub>-p<sub>2</sub></td><td>972</td><td>107</td><td>77</td><td>253</td><td>1.176</td><td>191</td><td>394</td><td>28.043</td><td>-</td></tr><tr><td>d<sub>4</sub>-n<sub>6</sub>-p<sub>3</sub></td><td>31.101</td><td>1.316</td><td>596</td><td>7.444</td><td>43.803</td><td>6.336</td><td>8.817</td><td>-</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>7</sub>-p<sub>6</sub></td><td>5.13</td><td>1.82</td><td>0.77</td><td>3.01</td><td>15.3</td><td>1.84</td><td>4.95</td><td>443</td><td>-</td></tr><tr><td>d<sub>3</sub>-n<sub>7</sub>-p<sub>2</sub></td><td>13.4</td><td>3.61</td><td>2.23</td><td>9.59</td><td>54.1</td><td>5.29</td><td>12.5</td><td>872</td><td>-</td></tr><tr><td>d<sub>3</sub>-n<sub>7</sub>-p<sub>3</sub></td><td>1.263</td><td>164</td><td>32.4</td><td>533</td><td>3.647</td><td>406</td><td>984</td><td>-</td><td>-</td></tr><tr><td>d<sub>3</sub>-n<sub>7</sub>-p<sub>4</sub></td><td>22.296</td><td>2.235</td><td>469</td><td>6.605</td><td>47.286</td><td>5.348</td><td>10.001</td><td>-</td><td>-</td></tr><tr><td>d<sub>3</sub>-n<sub>7</sub>-p<sub>5</sub></td><td>126.006</td><td>137.724</td><td>2.881</td><td>29.740</td><td>204.718</td><td>22.925</td><td>33.635</td><td>-</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>8</sub>-p<sub>5</sub></td><td>11.7</td><td>8.37</td><td>1.79</td><td>15.1</td><td>99.9</td><td>7.92</td><td>20.4</td><td>3.972</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>8</sub>-p<sub>6</sub></td><td>95.7</td><td>63.7</td><td>10.5</td><td>54.3</td><td>387</td><td>33.8</td><td>63.1</td><td>15.950</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>8</sub>-p<sub>7</sub></td><td>265</td><td>79.6</td><td>22.2</td><td>81.9</td><td>556</td><td>47.2</td><td>122</td><td>15.125</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>8</sub>-p<sub>2</sub></td><td>288</td><td>276</td><td>18.1</td><td>98.3</td><td>787</td><td>71.7</td><td>135</td><td>15.252</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>8</sub>-p<sub>3</sub></td><td>25.503</td><td>3.710</td><td>471</td><td>11.059</td><td>107.744</td><td>8.984</td><td>13.705</td><td>-</td><td>-</td></tr><tr><td colspan="10" style="text-align:center; font-weight:bold;">positive-dimensional to 2-dimensional</td></tr><tr><td>d<sub>2</sub>-n<sub>8</sub>-p<sub>5</sub></td><td>620</td><td>157</td><td>43</td><td>147</td><td>663</td><td>87</td><td>166</td><td>26.746</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>9</sub>-p<sub>4</sub></td><td>82</td><td>119</td><td>28</td><td>35</td><td>139</td><td>20</td><td>50</td><td>4.563</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>9</sub>-p<sub>5</sub></td><td>2.889</td><td>1.617</td><td>448</td><td>924</td><td>4574</td><td>726</td><td>735</td><td>-</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>9</sub>-p<sub>6</sub></td><td>40.352</td><td>8.290</td><td>2.155</td><td>8.775</td><td>40.921</td><td>7.010</td><td>6.969</td><td>-</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>10</sub>-p<sub>4</sub></td><td>202</td><td>472</td><td>69</td><td>101</td><td>530</td><td>64</td><td>142</td><td>14.237</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>10</sub>-p<sub>5</sub></td><td>11.347</td><td>11.512</td><td>2.801</td><td>3.373</td><td>19.465</td><td>2.710</td><td>3.148</td><td>-</td><td>-</td></tr><tr><td colspan="10" style="text-align:center; font-weight:bold;">positive-dimensional to 3-dimensional</td></tr><tr><td>d<sub>2</sub>-n<sub>9</sub>-p<sub>4</sub></td><td>284</td><td>100</td><td>22</td><td>75</td><td>325</td><td>46</td><td>77</td><td>10.014</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>9</sub>-p<sub>5</sub></td><td>5.575</td><td>758</td><td>201</td><td>747</td><td>3.422</td><td>333</td><td>845</td><td>-</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>10</sub>-p<sub>4</sub></td><td>738</td><td>481</td><td>108</td><td>180</td><td>981</td><td>109</td><td>185</td><td>28.688</td><td>-</td></tr><tr><td>d<sub>2</sub>-n<sub>10</sub>-p<sub>5</sub></td><td>66.845</td><td>12.082</td><td>2.141</td><td>25.707</td><td>60.054</td><td>23.100</td><td>6.885</td><td>-</td><td>-</td></tr></tbody></table>

The main ideal is twofold. First, the SPARSE-FGLM algorithm builds a sequence $\mathbf{u}$ from $\mathcal{G}_{\text{DRL}}$, hence unless $I$ is not Gorenstein and the initial terms of $\mathbf{u}$ are not random enough, $I_C(\mathbf{u}) = I$. Then, we create a new sequence $\mathbf{v} = (\text{Eval}(\varphi\mathbf{x}^i, \mathbf{u}))_{i \in \mathbb{N}^n}$ and prove that it satisfies

$$I_C(\mathbf{v}) = I_C(\mathbf{u}) : \langle \varphi \rangle = I : \langle \varphi \rangle.$$

Finally, it remains to recover the relations satisfied by $\mathbf{v}$. This is where the zero-dimensionality assumption appears, we know that only finitely many monomials $m$ are linearly independent in $\mathbb{K}[\mathbf{x}]/(I : \langle \varphi \rangle)$, which is equivalent to saying that finitely many polynomials NF ($m\varphi, \mathcal{G}_{\text{DRL}}, <_{\text{DRL}}$) are linearly independent in $\mathbb{K}[\mathbf{x}]/I$. Thus, we can work in a finite-dimensional subspace $W$ of $\mathbb{K}[\mathbf{x}]/I$ spanned by a set of monomials that allows us to determine $\mathcal{H}_{\text{LEX}}$.

The shape position assumption on $\mathcal{H}_{\text{LEX}}$ is of upmost importance to mimic the SPARSE-FGLM algorithm by building only one matrix, related to the multiplication by $x_n$. Since $W$ need not be stable by the multiplication by $x_n$, we consider the map multiplication by $x_n$ composed with the projection on $W$. Assuming $W$ is large enough, this new map, and its associated matrix, allow us to determine $\mathcal{H}_{\text{LEX}}$ by tweaking Algorithm 2.2 as follows in Algorithm 3.3.

**Theorem 1.11 (see also [15, Th. 4.12]):** Let $I$ be a positive-dimensional ideal of $\mathbb{K}[\mathbf{x}]$, let $\mathcal{G}_{\text{DRL}}$ be its reduced $<_{\text{DRL}}$-Gröbner basis and $S_{\text{DRL}}$ be the associated staircase. Let $\varphi \in \mathbb{K}[\mathbf{x}] \setminus I$ such that $I : \langle \varphi \rangle$ is zero-dimensional of degree $D'$ and in shape position.

Let $\Sigma$ be a finite staircase of size $N$ containing supp NF ($x_n^i\varphi, \mathcal{G}_{\text{DRL}}, <_{\text{DRL}}$) for all $i \in \mathbb{N}$, where NF ($f, G, <$) is the normal form of $f$ w.r.t. $G$ and $<$.

Let $t$ be the number of monomials $\sigma$ in $\Sigma$ such that $x_n\sigma \in \text{LM}_{<_{\text{DRL}}} (\mathcal{G}_{\text{DRL}})$ and let $\mathbf{u}$ be the number of monomials $\sigma$ in $\Sigma$ such that $x_n\sigma \in \langle \text{LM}_{<_{\text{DRL}}}(I) \rangle \setminus \text{LM}_{<_{\text{DRL}}}(\mathcal{G}_{\text{DRL}})$.

Then, for a generic choice of vector $\mathbf{r} \in \mathbb{K}^N$, Algorithm 3.3 terminates and returns the reduced $<_{\text{LEX}}$-Gröbner basis of $I : \langle \varphi \rangle$. To do so, it requires at most $\mathbf{u} + \mathbf{n}$ normal form computations w.r.t. $\mathcal{G}_{\text{DRL}}$ and $<_{\text{DRL}}$ plus $O((t + u + n)ND')$ operations in $\mathbb{K}$.

As a practical optimization, we show that any zero column of the matrix $\tilde{M}$ can be removed, together with its associated row. Repeating this process can lead to considering a much smaller
**Input:** The reduced <DRL-Gröbner basis $\mathscr{G}_{\text{DRL}}$ of a generic ideal, a polynomial $\varphi \in \mathbb{K}[\mathbf{x}]$, and a finite staircase $\Sigma$ of size $N$ containing supp NF$(x_n^k\varphi, \mathscr{G}_{\text{DRL}}, <_{\text{DRL}})$ for all $k \in \mathbb{N}$.

**Output:** The reduced <LEX-Gröbner basis of $I$: $\langle\varphi\rangle$, if it is in shape position.

1. Build the matrix $\tilde{M}$ of the multiplication by $x_n$ followed by the projection on the space spanned by $\Sigma$.

2. Pick $r \in \mathbb{K}^N$ a row-vector at random.

3. Build $\varphi$ the column-vector of coefficients of $\varphi$ restricted to $\Sigma$.

4. **For $k$ from 1 to $n-1$ do**

5. **$\lfloor$** Build $\psi_k$ the column-vector of coefficients of NF $(x_k\varphi, \mathscr{G}_{\text{DRL}}, <_{\text{DRL}})$ restricted to $\Sigma$.

6. Compute $(w_i^{(0)})_{0 \le i < 2N}, (w_i^{(1)})_{0 \le i < N}, \dots, (w_i^{(n-1)})_{0 \le i < N}$ with Algorithm 2.2 called on $\tilde{M}, r, \varphi, \psi_1, \dots, \psi_{n-1}$.

7. $h_n \leftarrow$ Berlekamp-Massey$(w_0^{(0)}, \dots, w_{2D-1}^{(0)})$, $D' := \deg h_n$.

8. **If** NF $(h_n\varphi, \mathscr{G}_{\text{DRL}}, <_{\text{DRL}}) \neq 0$ **then Return** “Bad vector”.

9. Compute $h_1 := \gamma_{D'-1,1}x_n^{N-1} + \dots + \gamma_{0,1}, \dots, h_{n-1} := \gamma_{D'-1,n-1}x_n^{N-1} + \dots + \gamma_{0,n-1}$ with Algorithm 2.3 called on $(w_i^{(0)})_{0 \le i < 2D'-1}, (w_i^{(1)})_{0 \le i < D'}, \dots, (w_i^{(n-1)})_{0 \le i < D'}$.

10. **For $k$ from 1 to $n-1$ do**

11. **$\lfloor$** If NF $((x_k - h_k(x_n))\varphi, \mathscr{G}_{\text{DRL}}, <_{\text{DRL}}) \neq 0$ **then Return** “Not in shape position”.

12. **Return** $\{h_n(x_n), x_{n-1} - h_{n-1}(x_n), \dots, x_1 - h_1(x_n)\}$.

**Algorithm 3.3: SPARSE-FGLM-COLON**

submatrix. In Table 3.3, we give an example where a 81 068 × 81 068-matrix is replaced by a 32 184 × 32 184-submatrix for instance.

**Example 3.5:** Let $\mathscr{G}'_{\text{DRL}} = \{x_2^2 - x_1 + x_2, x_1^2 - x_1x_2\}$ be the <DRL-Gröbner basis of the ideal $I'$ and $\varphi = x_1^4 - x_2^4 - 2x_1^2 + 1$. We build a sequence $\mathbf{u}$ with random initial coefficients, e.g. $u_{0,0} = -4$, $u_{1,0} = 10$, $u_{0,1} = -3$ and $u_{1,1} = -8$, thanks to $\mathscr{G}'_{\text{DRL}}$:

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

We now consider the auxiliary sequence $\mathbf{v} = (\text{Eval}(x_1^{t_1} x_2^{t_2} \varphi, \mathbf{u}))_{(t_1, t_2) \in \mathbb{N}^2}$, i.e. $\mathbf{v}_{t_1, t_2} = \mathbf{u}_{t_1+4, t_2} - \mathbf{u}_{t_1, t_2+4} - 2\mathbf{u}_{t_1+2, j} + \mathbf{u}_{t_1, t_2}$ for all $(t_1, t_2) \in \mathbb{N}^2$, to find the <LEX-Gröbner basis $\mathscr{H}_{\text{LEX}}$ of $I' : \langle\varphi\rangle = J$,

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

By Example 2.11, we have $\mathscr{H}_{\text{LEX}} = \{x_2^3, x_1 - x_2^2 - x_2\}$.

In Table 3.3, we compare Algorithm 3.3 for computing a <LEX-Gröbner basis of the zero-dimensional colon ideal $I : \langle\varphi\rangle = I : \langle\varphi\rangle^\infty$ for <LEX with MAPLE [9] using the Groebner:-Basis command followed by the Groebner:-FGLM command. The considered input generators are $I = \langle f, \frac{\partial f}{\partial x_1}, \dots, \frac{\partial f}{\partial x_n} \rangle$ and $\varphi = (\frac{\partial f}{\partial x_n})^M$ for $M$ large enough, with $f$ the sum of $p$ squares of polynomials of degree $d$ in $n$ variables. The columns #Σ and #Σ' correspond to the size of the set Σ before and
after reductions by removing the zero columns, while column $D'$ gives the degree of the saturated ideal. Whether it is between $\#\Sigma$ and $\#\Sigma'$ or between $\#\Sigma'$ and $D'$, we can observe ratios going up to around 5. Therefore, it is clear that the algorithm would not be as efficient if one were to work with $\Sigma$ directly. Still, it would be even more beneficial to reduce further the size of $\Sigma'$ to be as close as possible to $D'$.

For `MSOLVE`, we give proportions of time spent to compute the $<_{\text{DRL}}$-Gröbner basis of $I$ with $F_4$, to compute the normal form of correct power of $\varphi$, to build the multiplication matrix, including the computation of extra normal forms, and then to compute the $<_{\text{LEX}}$-Gröbner basis. We also give the time for `MAPLE` in seconds using Rabinowitsch's trick [90] and give the time for computing a $<_{\text{DRL}}$-Gröbner basis and then for the change of order step to obtain a $<_{\text{LEX}}$-Gröbner basis.

We can notice that the `SPARSE-FGLM-COLON` algorithm approach is most efficient when either the change of order step is the most time-consuming or when the ratios between $\#\Sigma$, $\#\Sigma'$ and $D'$ are the smallest. In the former case, the algorithm benefits from the regularity of the computation of the reduced $<_{\text{DRL}}$-Gröbner basis of $I$ compared to the one of $I + (1 - t\varphi)$ in the Rabinowitsch's trick approach. In the latter case, when $\Sigma$ or $\Sigma'$ are large compared to $D'$, the overhead in the linear algebra part becomes overwhelming. Clearly, in a multi-modular approach, one would want to consider an even smaller subset of $\Sigma'$ to perform the computations, once $D'$ is known. All in all, we can see speed-ups that are significant and sometimes higher than 10.

<table><caption>Table 3.3: Timings in seconds, $<_{\text{LEX}}$-Gröbner bases, positive-to-0-dimensional case.</caption><thead><tr><th rowspan="2"></th><th colspan="3">sizes</th><th colspan="4">MSOLVE</th><th colspan="3">MAPLE Groebner</th></tr><tr><th>#Σ</th><th>#Σ'</th><th>D'</th><th>F₄</th><th>Sat. order</th><th>Mat.</th><th>FGLM</th><th>Total</th><th>Basis</th><th>FGLM</th><th>Total</th></tr></thead><tbody><tr><td>d₂-a₈-p₅</td><td>5746</td><td>2636</td><td>1516</td><td>60%</td><td>20%</td><td>7%</td><td>13%</td><td>21</td><td>96%</td><td>4%</td><td>56</td></tr><tr><td>d₂-a₈-p₆</td><td>7901</td><td>5100</td><td>3756</td><td>35%</td><td>9%</td><td>6%</td><td>50%</td><td>140</td><td>88%</td><td>12%</td><td>350</td></tr><tr><td>d₂-a₈-p₇</td><td>8841</td><td>7340</td><td>6444</td><td>33%</td><td>9%</td><td>6%</td><td>52%</td><td>320</td><td>79%</td><td>21%</td><td>890</td></tr><tr><td>d₂-a₉-p₅</td><td>11748</td><td>4548</td><td>2308</td><td>56%</td><td>14%</td><td>8%</td><td>22%</td><td>150</td><td>68%</td><td>32%</td><td>410</td></tr><tr><td>d₂-a₉-p₆</td><td>18829</td><td>10372</td><td>6788</td><td>40%</td><td>7%</td><td>8%</td><td>45%</td><td>1200</td><td>91%</td><td>9%</td><td>3200</td></tr><tr><td>d₂-a₉-p₇</td><td>24332</td><td>17540</td><td>13956</td><td>33%</td><td>5%</td><td>7%</td><td>55%</td><td>4400</td><td>83%</td><td>17%</td><td>12000</td></tr><tr><td>d₂-n₁₀-p₄</td><td>9724</td><td>1996</td><td>652</td><td>67%</td><td>27%</td><td>5%</td><td>1%</td><td>42</td><td>99%</td><td>1%</td><td>68</td></tr><tr><td>d₂-n₁₀-p₅</td><td>22408</td><td>7372</td><td>3340</td><td>52%</td><td>11%</td><td>9%</td><td>28%</td><td>900</td><td>97%</td><td>3%</td><td>1200</td></tr><tr><td>d₂-n₁₀-p₆</td><td>40946</td><td>19468</td><td>11404</td><td>42%</td><td>7%</td><td>9%</td><td>42%</td><td>9900</td><td>92%</td><td>8%</td><td>17000</td></tr><tr><td>d₃-n₅-p₃</td><td>3034</td><td>1320</td><td>672</td><td>39%</td><td>33%</td><td>9%</td><td>19%</td><td>1.1</td><td>97%</td><td>3%</td><td>4</td></tr><tr><td>d₃-n₅-p₄</td><td>3750</td><td>2616</td><td>1968</td><td>27%</td><td>27%</td><td>11%</td><td>35%</td><td>4.3</td><td>95%</td><td>5%</td><td>43</td></tr><tr><td>d₃-n₆-p₃</td><td>10773</td><td>3792</td><td>1632</td><td>60%</td><td>11%</td><td>8%</td><td>21%</td><td>50</td><td>96%</td><td>4%</td><td>77</td></tr><tr><td>d₃-n₆-p₄</td><td>16271</td><td>9192</td><td>5952</td><td>37%</td><td>5%</td><td>6%</td><td>52%</td><td>500</td><td>89%</td><td>11%</td><td>1300</td></tr><tr><td>d₃-n₆-p₅</td><td>18897</td><td>14862</td><td>12432</td><td>12%</td><td>5%</td><td>6%</td><td>77%</td><td>1200</td><td>82%</td><td>18%</td><td>5600</td></tr><tr><td>d₃-n₇-p₃</td><td>35117</td><td>10320</td><td>3840</td><td>52%</td><td>9%</td><td>8%</td><td>31%</td><td>1300</td><td>95%</td><td>5%</td><td>1100</td></tr><tr><td>d₃-n₇-p₄</td><td>62104</td><td>29760</td><td>16800</td><td>19%</td><td>4%</td><td>11%</td><td>66%</td><td>12000</td><td>91%</td><td>9%</td><td>31000</td></tr><tr><td>d₄-n₅-p₃</td><td>15881</td><td>7560</td><td>4104</td><td>41%</td><td>6%</td><td>5%</td><td>48%</td><td>200</td><td>94%</td><td>6%</td><td>620</td></tr><tr><td>d₄-n₅-p₄</td><td>19274</td><td>14088</td><td>11016</td><td>32%</td><td>4%</td><td>4%</td><td>60%</td><td>1000</td><td>86%</td><td>14%</td><td>4700</td></tr><tr><td>d₄-n₆-p₂</td><td>41189</td><td>8424</td><td>1944</td><td>74%</td><td>5%</td><td>9%</td><td>12%</td><td>590</td><td>97%</td><td>3%</td><td>224</td></tr><tr><td>d₄-n₆-p₃</td><td>81068</td><td>32184</td><td>14904</td><td>16%</td><td>3%</td><td>12%</td><td>69%</td><td>11000</td><td>92%</td><td>8%</td><td>27000</td></tr><tr><td>d₃-n₄-p₃</td><td>7235</td><td>4540</td><td>3040</td><td>13%</td><td>24%</td><td>11%</td><td>52%</td><td>9</td><td>93%</td><td>7%</td><td>180</td></tr><tr><td>d₃-n₅-p₃</td><td>54787</td><td>27360</td><td>15360</td><td>8%</td><td>4%</td><td>7%</td><td>81%</td><td>5100</td><td>92%</td><td>8%</td><td>22000</td></tr></tbody></table>

## 3.3 Parametrizations of the radical

**Proposition 3.6 (see also [15, Prop. 4.17]):** Let $I$ be a positive-dimensional ideal of $\mathbb{K}[x]$, let $\mathcal{G}_{\text{DRL}}$ be its reduced $<_{\text{DRL}}$-Gröbner basis and $S_{\text{DRL}}$ be the associated staircase. Let $\varphi \in \mathbb{K}[x] \setminus I$ such that $I : \langle\varphi\rangle$ is zero-dimensional, let $\mathcal{H}_{\text{LEX}}$ be its reduced $<_{\text{LEX}}$-Gröbner basis and let $h_n \in \mathcal{H}_{\text{LEX}} \cap \mathbb{K}[x_n]$.

If $\sqrt{I} : \langle\varphi\rangle$ is in shape position, then one can compute its reduced $<_{\text{LEX}}$-Gröbner basis calling the `SPARSE-FGLM-COLON` algorithm with the following modifications:

(1) On line 7, $h_n$ is the squarefree part of the polynomial returned by the Berlekamp-Massey algorithm.

(2) On line 9, $h_k$ is obtained thanks to [66, Algo. 2], see also [29].

In practice, when an ideal $J$ is not in shape position, it is not easy to check that $\sqrt{J}$ is. Therefore, using the following lemma is the cornerstone of our probabilistic verification algorithm in `MSOLVE` [13, 14] when $J$ is not in shape position but its radical might be, see [13, Sec. 4.4].

**Lemma 1.12 (see also [15, Lem. 4.15]):** Let $I$ be a zero-dimensional ideal of $\mathbb{K}[x]$. Let $\lambda \in \bar{\mathbb{K}}$ be generic. Then, for $1 \le k \le n$, $I = I : \langle x_k + \lambda \rangle$.

We proceed as follows.

(1) Compute the polynomials $x_k - g_k(x_n)$ in $\sqrt{J}$ for $1 \le k \le n-1$, with $\deg g_k$ minimal.

(2) Compute the polynomials $x_k - g'_k(x_n)$ in $\sqrt{J} : \langle x_k + \lambda \rangle$ for $\lambda$ picked at random and $1 \le k \le n-1$,
**Chapter 3. Computing Gröbner bases**

**3.3. Parametrizations of the radical**

with deg $g'_k$ minimal.

(3) Check whether $g_k = g'_k$ for $1 \le k \le n-1$.

By Lemma 1.12, for a generic $\lambda$, $J = J : \langle x_k + \lambda \rangle$, hence both radical ideals are the same. Furthermore, if they are in shape position, then $g_k = g'_k$ for $1 \le k \le n-1$. Therefore, any discrepancy must come from the fact that $\sqrt{J}$ is not in shape position and the polynomials $x_k - g_k(x_n)$ and $x_k - g'_k(x_n)$ are meaningless.
This chapter is devoted to guessing Gröbner bases of ideals of C-relations of sequences. That is, given a finite number of terms of a sequence, the goal is to discover the linear recurrence relations, given by polynomials, satisfied by these few terms.

In some applications, like the Gröbner bases change of order one, see Sections 2.1.2 and 3.1, or in coding theory [23, 65, 96], a bound on the sizes of the relations is known. Therefore, testing the relations with sufficiently enough terms ensures that these are satisfied by the whole sequence. In other cases, like in combinatorics or number theory [3, 25, 30, 31] no a priori bound is known. In which case, one must prove afterwards the validity of the relations.

We first recall how to guess C-relations for uni-indexed sequences using the seminal Berlekamp–Massey algorithm [8, 77] and its faster variants [34] in Section 4.1. Then, we describe extensions of this algorithm to multi-indexed sequences, where my contributions lie, in Section 4.2.

This chapter is based on joint work with Brice Boyer<fnref n="1" />, Jean-Charles Faugère<fnref n="2" />,<fnref n="3" /> and Mohab Safey El Din<fnref n="4" /> [11, 12, 17, 19, 21].

In all this chapter, *multi-Hankel* matrices will be used. These generalize the notion of *Hankel* matrices to multi-indexed as follows: given two sets of monomials $T_1$ and $T_2$ in $x$ and a sequence $u$,

$$H_{T_1, T_2} = x^i \in T_1 \left( \begin{array}{cccc} \dots & & x^j \in T_2 & \dots \\ & \ddots & & \vdots \\ \dots & \cdots & \text{Eval} (x^i x^j, u) & \dots \\ & & \vdots & \ddots \end{array} \right)$$

## 4.1 Uni-indexed sequences and the Berlekamp–Massey algorithm

Let $u = (u_i)_{i \in \mathbb{N}}$ be a one-dimensional sequence that we assume C-finite. The goal is to find the unique polynomial $G_{x^d} = x^d + \gamma_{x^{d-1}}x^{d-1} + \dots + \gamma_1$ with $d$ minimal such that $I_C(u) = \langle G_{x^d} \rangle$. Though, in general $d$ is not known and only the first $N$ terms of $u$: $u_0, \dots, u_{N-1}$ are known. Therefore, we aim to guess $G_{x^d}$ using the fact that it must satisfies

$$\text{Eval}(G_{x^d}, u) = \dots = \text{Eval}(G_{x^d} x^{N-d-1}, u) = 0. \quad (4.1)$$

### 4.1.1 Matrix viewpoints

Requiring that $d$ is minimal and $\gamma_1, \dots, \gamma_{x^{d-1}}$ are such that Equation (4.1) is satisfied is equivalent to looking for the largest collection of vectors

$$\begin{pmatrix} \gamma_1 \\ \vdots \\ \gamma_{x^{d-1}} \\ 1 \\ 0 \\ 0 \\ \vdots \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \gamma_1 \\ \vdots \\ \gamma_{x^{d-1}} \\ 1 \\ 0 \\ \vdots \\ 0 \end{pmatrix}, \dots, \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ \gamma_1 \\ \vdots \\ \gamma_{x^{d-1}} \\ 1 \end{pmatrix} \text{ in the kernel of}$$

$$H_{\{1\}, \{1, \dots, x^{d-1}, x^d, x^{d+1}, \dots, x^{N-1}\}} = 1(u_0 \quad \dots \quad x^{d-1} \quad x^d \quad x^{d+1} \quad \dots \quad x^{N-1} \\ \qquad \dots \quad u_{d-1} \quad u_d \quad u_{d+1} \quad \dots \quad u_{N-1}).$$

<footnote n="1">¹former Post-doc at Sorbonne Université</footnote>

<footnote n="2">²INRIA</footnote>

<footnote n="3">³CryptoNext Security</footnote>

<footnote n="4">⁴Sorbonne Université</footnote>
Equivalently, one also looks for the smallest integer $d$ such that

$$H_{\mathcal{T}_{\le x}^{N-d-1}, \mathcal{T}_{\le x^d}} \begin{pmatrix} \gamma_1 \\ \vdots \\ \gamma_{x^{d-1}} \\ 1 \end{pmatrix} = \begin{pmatrix} 1 & \cdots & x^{d-1} & x^d \\ u_0 & \cdots & u_{d-1} & u_d \\ u_1 & \cdots & u_d & u_{d+1} \\ \vdots & & \vdots & \vdots \\ u_{N-d-1} & \cdots & u_{N-2} & u_{N-1} \\ u_{N-d} & \cdots & u_{N-1} & 0 \\ \vdots & \ddots & \vdots & \vdots \\ u_{N-1} & 0 & \cdots & 0 \end{pmatrix} \begin{pmatrix} \gamma_1 \\ \vdots \\ \gamma_{x^{d-1}} \\ 1 \end{pmatrix} = 0,$$

where $\mathcal{T}_{\le m}$ denotes the set of all monomials less or equal to $m$ for $<$, the only monomial order on $\mathbb{K}[x]$.

With this modeling of the problem, we are led to a polynomial viewpoint.

### 4.1.2 Polynomial viewpoint

The Hankel matrix-vector product can be extended into

$$ \begin{pmatrix} u_0 & \cdots & u_{d-1} & u_d \\ u_1 & \cdots & u_d & u_{d+1} \\ \vdots & & \vdots & \vdots \\ u_{N-d-1} & \cdots & u_{N-2} & u_{N-1} \\ u_{N-d} & \cdots & u_{N-1} & 0 \\ \vdots & \ddots & \vdots & \vdots \\ u_{N-1} & 0 & \cdots & 0 \end{pmatrix} \begin{pmatrix} \gamma_1 \\ \vdots \\ \gamma_{x^{d-1}} \\ 1 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ r_{x^{d-1}} \\ \vdots \\ r_1 \end{pmatrix}, \quad (4.2) $$

representing the product of the two polynomials $P_{\mathcal{T}_{\le x}^{N-1}} = \sum_{i=0}^{N-1} u_i x^{N-i-1}$, given by the sequence terms, and $G_{x^d} = x^d + \sum_{k=0}^{d-1} \gamma_{x^k} x^k$, giving the C-relation, modulo $B = x^N$. This can rewritten as follows, there exists $Q \in \mathbb{K}[x]$ such that

$$ Q \cdot x^N + G_{x^d} \sum_{i=0}^{N-1} u_i x^{N-i-1} = R, \quad \deg R < d, \quad (4.3) $$

which is a Bézout relation between $x^N$ and $P_{\mathcal{T}_{\le x}^{N-1}}$, where we ask that the degree of the right-hand side member, or equivalently its leading term, is less than that of the cofactor of $P_{\mathcal{T}_{\le x}^{N-1}}$.

The main advantage of this viewpoint is that we can now compute $d$ and $G_{x^d}$ using the extended Euclidean algorithm on $x^N$ and $P_{\mathcal{T}_{\le x}^{N-1}}$ and stop the algorithm when a Bézout relation satisfies Equation (4.3).

This viewpoint gives rise to Algorithm 4.1, which is a variant of the BM algorithm and a modified extended Euclidean algorithm. While, as written, its complexity is quadratic in $N$, one can use fast Euclidean algorithm [34] to improve the time complexity to $O(M(n) \log n)$.

**Input:** $u_0, \dots, u_{N-1}$, the first $N$ terms of a sequence $u$.

**Output:** A polynomial $C$ of degree $d$, with $d$ minimal such that $\text{Eval}(Cx^i, u) = 0$ for all $0 \le i \le N-d-1$.

$$ 1 \ m := 1/x, \ L_m := [x^N, 0]. $$

$$ 2 \ m' := 1, \ L_{m'} := \left[ \sum_{i=0}^{N-1} u_i x^{N-i-1}, 1 \right]. $$

3 **While** $\deg L_{m',2} \le \deg L_{m',1}$ **do**

$$ 4 \quad Q := \text{Quo}(L_{m,1}, L_{m',1}). $$

$$ 5 \quad L_{m' \text{ LM}<}(Q) := L_m - Q L_{m'}. $$

$$ 6 \quad m, m' := m', m' \text{ LM<}(Q). $$

7 **Return** $L_{m',2}$.

<caption>**Algorithm 4.1:** Berlekamp–Massey/ Extended Euclidean algorithm</caption>
# Chapter 4. Guessing Gröbner bases

## 4.1. Uni-indexed sequences and the Berlekamp--Massey algorithm

**Example 4.1:** Let us consider the sequence $w = (w_i)_{i \in \mathbb{N}}$,

<table><thead><tr><th>i</th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>...</th></tr></thead><tbody><tr><td></td><td>-4</td><td>10</td><td>-13</td><td>5</td><td>-5</td><td>5</td><td>-5</td><td>5</td><td>-5</td><td>...</td></tr></tbody></table>

and assume $N = 9$.

With the matrix viewpoint, the kernel of

$$H_{\{1\},\{1,x,x^2,x^3,x^4,x^5,x^6,x^7,x^8\}} = 1 \begin{pmatrix} 1 & x & x^2 & x^3 & x^4 & x^5 & x^6 & x^7 & x^8 \\ -4 & 10 & -13 & 5 & -5 & 5 & -5 & 5 & -5 \end{pmatrix},$$

gives us the collection of vectors $\begin{pmatrix} 0 \\ 0 \\ 1 \\ 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 1 \\ 1 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 1 \\ 1 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 1 \end{pmatrix}$, hence $I_C(w) = \langle x^4 + x^3 \rangle$.

Equivalently, we have that

$$H_{\{1,x,x^2,x^3,x^4\},\{1,x,x^2,x^3,x^4\}} \begin{pmatrix} 0 \\ 0 \\ 0 \\ 1 \\ 1 \end{pmatrix} = \begin{pmatrix} 1 & x & x^2 & x^3 & x^4 \\ x & -4 & 10 & -13 & 5 \\ x^2 & 10 & -13 & 5 & -5 \\ x^3 & -13 & 5 & -5 & 5 \\ x^4 & 5 & -5 & 5 & -5 \end{pmatrix} \begin{pmatrix} 0 \\ 0 \\ 0 \\ 1 \\ 1 \end{pmatrix} = 0$$

but

$$H_{\{1,x,x^2,x^3,x^4\},\{1,x,x^2,x^3\}} = \begin{pmatrix} 1 & x & x^2 & x^3 \\ x & -4 & 10 & -13 \\ x^2 & 10 & -13 & 5 \\ x^3 & -13 & 5 & -5 \\ x^4 & 5 & -5 & 5 \\ x^5 & -5 & 5 & -5 \end{pmatrix}$$

has full rank, ensuring that no relation given by a polynomial of degree 3 exists.

With the polynomial viewpoint, we perform the extended Euclidean algorithm on $x^9$ and $P_{\mathbb{F}_{\le x^8}}$ and find

$$L_{1/x} = [x^9, 0]$$

$$L_1 = [-4x^8 + 10x^7 - 13x^6 + 5x^5 - 5x^4 + 5x^3 - 5x^2 + 5x - 5, 1]$$

$$L_x = [24x^7 - 55x^6 + 15x^5 - 15x^4 + 15x^3 - 15x^2 + 15x - 25, 2x + 5]$$

$$L_{x^2} = [-1237x^6 + 285x^5 - 285x^4 + 285x^3 - 285x^2 + 45x - 595, 48x^2 + 110x + 119]$$

$$L_{x^3} = [-2560x^5 + 2560x^4 - 2560x^3 - 3625x^2 + 2200x - 1600, \\ 1237x^3 + 285x^2 - 120x + 320]$$

$$L_{x^4} = [-5x^3, x^4 + x^3] = [R_{x^4}, G_{x^4}],$$

and $G_{x^4} = x^4 + x^3$ is a valid relation.

**Remark 4.2:** The BM algorithm always returns a non-zero relation. If no pair $L_{x^\delta} = [R_{x^\delta}, G_{x^\delta}]$ satisfies the requirements, then it will return a pair $L_{x^d}$ with $\mathrm{LM}_<(G_{x^d}) > x^{N-1}$. From a matrix viewpoint, it returns an element of the kernel of the empty matrix $H_{\emptyset, \mathbb{F}_{\le x^d}}$.
4.2 Extensions to multi-indexed sequences

This section is devoted to several extensions of the uni-indexed case to the multi-indexed one.
It is based on several joint works with Brice Boyer<fnref n="5" />, Jean-Charles Faugère<fnref n="6" />,<fnref n="7" /> and Mohab Safey El
Din<fnref n="8" />.

A first extension to the BM algorithm was proposed by Sakata for bi-indexed sequences [93] and then was extended to *n*-indexed sequences [95, 97]. This is the so-called Berlekamp–Massey–Sakata (BMS) algorithm. Then, Brice Boyer<fnref n="5" />, Jean-Charles Faugère<fnref n="6" />,<fnref n="7" /> and I designed a linear-algebra-based one, called SCALAR-FGLM, in [11, 12] and an adaptive variant thereof. Mourrain also proposed another algorithm using a Gram–Schmidt process in [81], called AGBB. Finally Jean-Charles Faugère<fnref n="6" />,<fnref n="7" /> and I on one hand and Mohab Safey El Din<fnref n="8" /> and I on another hand proposed several extensions and improvements of the SCALAR-FGLM algorithm in [17, 19, 21].

We assume now that < is a monomial order such that for any monomial $x^i \in \mathbb{N}^n$, the set $\mathcal{T}_{\le x^i}$ is finite. Such orders compare first monomials by (weighted) degrees, like $<_{\text{DRL}}$. While $<_{\text{LEX}}$ is not such a monomial order, one can choose a monomial order $<$ such that the $<_{\text{LEX}}$-Gröbner basis and the $<_{\text{LEX}}$-Gröbner basis coincide by putting large weights on the first variables. Though, these weights will depend on the ideal of C-relations one is dealing with.

4.2.1 The BMS algorithm

Given a multi-indexed sequence $\mathbf{u} = (u_i)_{i \in \mathbb{N}^n}$, a monomial order $<_{\mathcal{T}}$ and a monomial $M = x_1^{a_1} \cdots x_n^{a_n}$, the BMS algorithm extends the BM algorithm by computing vectors in the kernel of a multi-Hankel matrix

$$
H_{\{1\}, \mathcal{T} \leq M} = 1(u_{0, \dots, 0}, \dots, u_{a_1, \dots, a_n}),
$$

corresponding to having relations Eval($G_m x^i$, $\mathbf{u}$) = 0, with $\mathit{LM}_<(G_m) = m$ minimal for the division and for all $\mathbf{x}^i$ such that $m\mathbf{x}^i \le M$.

While the BMS algorithm computes the C-relations using polynomial multiplications by monomials and additions, Jean-Charles Faugère<fnref n="6" />,<fnref n="7" /> and I proposed in [17, 19] the following polynomial interpretation. This is a key argument towards the design of our algorithm that guesses the C-relations using polynomial arithmetic, including multivariate polynomial reductions, which extends the univariate Euclidean division.

Guessing a <-Gröbner basis of I<sub>C</sub>(u) comes down to finding the least (for the partial order |) monomials m<sub>1</sub>,...,m<sub>r</sub> &le; M such that dim ker H<sub>T<sub>&le;M</sub> T<sub>&le;m<sub>k</sub></sub></sub> > 0 with s<sub>k</sub> the greatest monomial such that s<sub>k</sub>m<sub>k</sub> &le; M for all k, 1 &le; k &le; r. Indeed, by the minimality of the leading monomials, if sufficiently many terms of u are known, these polynomials form a minimal <-Gröbner basis of I<sub>C</sub>(u).

Then, each multi-Hankel matrix-vector product can be extended further as in equation (4.2), tak-
ing the multi-Hankel matrix $H_{\mathcal{T}_{\leq M}, \mathcal{T}_{\leq m_k}}$ and replacing by zero any sequence term $u_{i+j} = \text{Eval}(\mathbf{x}^i \mathbf{x}^j, \mathbf{u})$
that we do not know, i.e. whenever $\mathbf{x}^i \mathbf{x}^j > M$. This yields the linear system

$$
\left(
\begin{array}{ccc}
\text{Eval}(1, \mathbf{u}) & \cdots & \text{Eval}(m_k^{-}, \mathbf{u}) & \text{Eval}(m_k, \mathbf{u}) \\
\text{Eval}(1^{+}, \mathbf{u}) & \cdots & \text{Eval}(1^{+}m_k^{-}, \mathbf{u}) & \text{Eval}(1^{+}m_k, \mathbf{u}) \\
\vdots & & \vdots & \vdots \\
\text{Eval}(s_k, \mathbf{u}) & \cdots & \text{Eval}(s_k m_k^{-}, \mathbf{u}) & \text{Eval}(s_k m_k, \mathbf{u}) \\
\text{Eval}(s_k^{+}, \mathbf{u}) & \cdots & \text{Eval}(s_k^{+} m_k^{-}, \mathbf{u}) & \mathbf{0} \\
\vdots & & \vdots & \vdots \\
\text{Eval}(M, \mathbf{u}) & \mathbf{0} & \cdots & \mathbf{0}
\end{array}
\right)
\left(
\begin{array}{c}
0 \\
0 \\
\vdots \\
0 \\
f_{M/s_k^{+}} \\
\vdots \\
f_{M/a}
\end{array}
\right)
,
\quad (4.4)
$$

where $\mathcal{M} = \text{LCM}(\mathcal{T}_{\leq M}) = x_1^{N_1} \cdots x_n^{N_n}$ and where, for a monomial $\mu$, $\mu^-$ (resp. $\mu^+$) is its preceding

<footnote n="5"><sup>5</sup>former Post-doc at Sorbonne Université</footnote>

<footnote n="6">⁶INRIA</footnote>

<footnote n="7">⁷CryptoNext Security</footnote>

<footnote n="8">⁸Sorbonne Université</footnote>
(resp. following) monomial for $<$. Let us notice that Eval$(s_k^+ m_k^-, \mathbf{u})$ can also be a $0$ if $s_k^+ m_k^- > M$ and that, more generally, the gray zeroes need not be diagonally aligned like they are in the univariate case.

Such product represents the product of the polynomials $P_{\mathcal{T}_{\leq M}} = \sum x^i \leq M u_i \frac{\mathcal{M}}{x^i}$ and $G_{m_k} = m_k + \sum x^i < m_k \gamma_{x^i} x^i$ modulo $B = (x_1^{N_1}, \dots, x_n^{N_n})$. The requirement for $G_{m_k}$ to encode a valid relation is now that $\text{LM}_<(R_{m_k}) < \frac{\mathcal{M}}{s_k}$ with $R_{m_k} = P_{\mathcal{T}_{\leq M}} G_{m_k} \bmod B$. To compute these relations, one start with the pairs $L_{1/2} = [x_i^{N_i}, 0]$ for all $1 \le i \le n$ and $L_1 = [P_{\mathcal{T}_{\leq M}}, 1]$. Then, one computes S-polynomials and polynomial reductions on the first coordinates of the pairs, and report the operations on the second coordinates of the pairs until finding all pairs $L_m = [R_m, G_m]$, where the leading monomial of $R_m$ is small enough for $G_m$ to encode a valid relation, for all $m$ minimal for the partial order $|$.

**Theorem 1.13 (see also [19, Th. 1]):** Let **u** be a sequence, < be a weighted degree monomial order and **M** be a monomial. Let us assume that the reduced <-Gröbner basis $\mathcal{G}$ of $I_C(\mathbf{u})$ and its associated staircase **S** satisfy $\max(S \cup \text{LM}_<(\mathcal{G})) \le \mathcal{M}$ and for all $m \le \mathcal{M}$, $s = \max_{\sigma \le \mathcal{M}} \{\sigma \mid \sigma m \le \mathcal{M}\}$, we have $\max(S) \le s$. Then, the variant of the Berlekamp–Massey–Sakata algorithm using polynomial arithmetic terminates and computes $\mathcal{G}$ in $O(\#S(\#S + \#\mathcal{G})\#\mathcal{T}_{\le M})$ operations in the base field, where $\mathcal{T}_{\le M}$ is the set of monomials less or equal to $\mathcal{M}$.

**Remark 4.3:** As for the BM algorithm, the BMS algorithm will always return a relation $G_m$ with $\text{LM}_<(G_m) = m$ a pure power in each variable. Therefore, it can return $G_m$ with $m > \mathcal{M}$, corresponding to a vector in the kernel of the empty matrix $H_{\emptyset, \mathcal{T}_{\le m}}$.

4.2.2 The SCALAR-FGLM algorithm

In [11, 12], Brice Boyer<fnref n="9" />, Jean-Charles Faugère<fnref n="10" />,<fnref n="11" /> and I designed an algorithm, SCALAR-FGLM,
that computes the reduced <-Gröbner basis of I_C(u) by means of linear algebra. Together with u,
it takes as an input a set of monomials T, ordered for <, and it computes the right kernel of the
multi-Hankel matrix H_{T,T}. Vectors in this kernel can be seen as polynomials in K[x] and these
polynomials with a leading term minimal for the partial order induced by the division are the ones
returned by the algorithm. Furthermore, if T contains the staircase and the leading monomials
of the reduced <-Gröbner basis of I_C(u), then the SCALAR-FGLM algorithm returns this Gröbner
basis.

**Input:** A sequence $\mathbf{u} = (u_i)_{i \in \mathbb{N}^n}$ with coefficients in $\mathbb{K}$, a monomial ordering $<$, a sufficiently large staircase $T$ and ordered for $<$.

**Output:** A reduced Gröbner basis of the ideal of $C$-relations of $\mathbf{u}$.

1. Build the matrix $H_{T,T}$.

2 Compute the set $S \subseteq T$ of smallest monomials, for <, such that $\text{rank } H_{S,S} = \text{rank } H_{T,T}$.

3 **For all** $m \in T \setminus S$ **do** // stabilize $S$ for the division

4 $\lfloor$ If $\exists s \in S$ such that $m \mid s$ **then** $S := S \cup \{m\}$.

5 $L := T \setminus S$ sorted for <.

6 $G := \emptyset$.

7 **While** $L \neq \emptyset$ **do**

8 $m := \min < L$

9 Solve the linear system $H_{S,S}\gamma + H_{S,\{g\}} = 0$.

10 $G := G \cup \{m + \sum_{s \in S} \gamma_s s\}$.

11 $\lfloor$ Remove $m$ and any of its multiples from $L$.

12 **Return** $G$.

<caption>**Algorithm 4.1: SCALAR-FGLM**</caption>

The algorithm computes the column rank profile of the matrix $H_{T,T}$, that is the set of leftmost linearly independent columns of the matrix. Since these columns are independent from the previous

<footnote n="9">⁹former Post-doc at Sorbonne Université</footnote>

<footnote n="10">¹⁰INRIA</footnote>

<footnote n="11">¹¹CryptoNext Security</footnote>
ones, their labels cannot be the leading monomial, for <, of any polynomial in the ideal of C-
relations, thus they are in the associated staircase of the reduced <Gröbner basis of this ideal.

**Example 4.4:** Let

$$
v =
\begin{array}{r@{\;}c@{\;}l@{\quad}l}
\multicolumn{3}{c}{\vdots} & \vdots & \vdots & \vdots & \vdots & \ddots \\
4 &   & 0 &   & 0 &   & 0 & \cdots \\
3 &   & 0 &   & 0 &   & 0 & \cdots \\
2 &   & -8 &   & 0 &   & 0 & \cdots \\
1 &   & 5 &   & -8 &   & 0 & \cdots \\
0 &   & 17 &   & -3 &   & -8 & \cdots \\
\cline{2-7}
\multicolumn{2}{c}{v_2} &   &   & 0 &   & 1 &   \\
\cline{2-7}
\multicolumn{2}{c}{v_1} &   &   &   &   & 2 &   \\
\cline{2-7}
\multicolumn{3}{c}{} &   &   &   & 3 &   \\
\cline{2-7}
\multicolumn{3}{c}{} &   &   &   & 4 &   \\
\cline{2-7}
\multicolumn{3}{c}{} &   &   &   &   & \cdots
\end{array}
$$

and assume we want to guess its reduced <sub>DRL</sub>-Gröbner basis knowing that it contains no polynomial
of degree higher than 3. We build the multi-Hankel $H_{T,T}$ where $T$ contains all the monomials of
degree at most 3 and then determine its column-rank-profile, i.e. the leftmost independent columns,

$$
H_{T,T} =
\begin{array}{@{}c@{\quad}r@{\quad}r@{\quad}r@{\quad}r@{\quad}r@{\quad}r@{\quad}r@{\quad}r@{\quad}r@{}}
  & \multicolumn{1}{c}{1} & y & x & y^2 & xy & x^2 & y^3 & xy^2 & xy^2 & x^3 \\
\cline{2-10}
  & \multicolumn{1}{c}{1} & 17 & 5 & -3 & -8 & -8 & -8 & 0 & 0 & 0 \\
  & y & 5 & -8 & -8 & 0 & 0 & 0 & 0 & 0 & 0 \\
  & x & -3 & -8 & -8 & 0 & 0 & 0 & 0 & 0 & 0 \\
  & y^2 & -8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
  & xy & -8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
  & x^2 & -8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
  & y^3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
  & xy^2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
  & x^2y & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
  & x^3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\end{array}
$$

Clearly, the first 3 columns are independent, then the other columns, labeled with $y^2$, $xy$, ... linearly depend on them. Thus, $\{1, y, x\} = S_{\text{DRL}}$. The leading monomials of the reduced <DRL-Gröbner basis are $y^2$, $xy$ and $x^2$ and we find the polynomials $y^2 - x + y$, $xy - x + y$ and $x^2 - x + y$.

The polynomial viewpoint presented in Section 4.2.1 for the BMS algorithm is in fact more general. Taking the same notation as in Section 4.2.1, consider the more subtle condition for a relation to be valid: the leading monomial of $R_m$, after discarding some monomials, is sufficiently small. Then, the polynomial viewpoint allows one to compute vectors in the kernel of a more general multi-Hankel, as those encountered in the SCALAR-FGLM algorithm. We refer the reader to [17, 19] for more details.

When the nonzero terms of the sequence all lie in a cone $\mathcal{C}$, we can decide to restrict the construction of the multi-Hankel matrix to monomials in the associated algebra, i.e. the set of monomials $\mathcal{T}(\mathcal{C}) = \{\mathbf{x}^i | i \in \mathcal{C}\}$. In that setting, instead of guessing classical <Gröbner bases, the guessed polynomials form a <-sparse Gröbner basis, as defined in [7, 56].

**Theorem 4.5 (see also [21, Th. 3.2]):** Let $\mathcal{C}$ be cone which is also a submonoid of $\mathbb{N}^n$. Let $<$ be a monomial order on the monomials in $\mathbf{x}$. Let $T \subseteq \mathcal{T}(\mathcal{C})$ be a set of monomials stable by division.

Let **u** be a C-finite sequence such that $I_C(\mathbf{u})$ has a reduced <-sparse Gröbner basis $\mathcal{G}$ with support in $T$. Then, the SCALAR-FGLM algorithm called on **u**, $T$ and < computes $\mathcal{G}$.

**4.2.3 Adaptive algorithms**

In many applications, including the change of order one when the ideal is not in shape posi-
tion, computing sequence terms is expensive so that building such a large multi-Hankel matrix is
prohibitive. In [11, 12], my co-authors and I proposed an adaptive algorithm that starts with a
1 × 1-matrix given by $T = \{1\}$. Then, it makes it grow by adding a monomial in $T$. If the rank of
the matrix grows, then we know that this new monomial is in the staircase of the sought Gröbner
basis. Otherwise, we have guessed a polynomial in the target Gröbner basis.

This approach is similar to the original FGLM algorithm, though we might find *fake* C-relations
when the first terms of the sequence are not generic.
<figure><img src="image_5.png" /><figcaption>Figure 4.1: Numbers of arithmetic operations (2D/3D)</figcaption></figure>

Later, Jean-Charles Faugère<fnref n="12" />,<fnref n="13" /> and I designed an adaptive variant of the BMS algorithm in [18] which can avoid testing some candidate relations, and thus requiring some sequence terms, when a bound of the degree of the zero-dimensional ideal of C-relations is known. Furthermore, in [19], we used the polynomial viewpoint of the SCALAR-FGLM algorithm of [17] and Section 4.2.2 to have an adaptive algorithm using only multivariate polynomial arithmetic. In particular, we prove the following theorem.

**Theorem 1.14 (see also [19, Th. 26]):** Let $\mathbf{u}$ be a sequence whose ideal of C-relations is zero-dimensional. Let $\prec$ be a monomial order, $\mathcal{G}$ be the reduced $\prec$-Gröbner basis of $I_C(\mathbf{u})$ and $S$ be the associated staircase.

Assuming the ADAPTIVE SCALAR-FGLM algorithm called on $\mathbf{u}$ returns a Gröbner basis $\mathcal{G}'$ with staircase $S'$ and $\#S' = \#S$, then $S' = S$ and $\mathcal{G}' = \mathcal{G}$. Furthermore, the algorithm does not need more that $\#2(S \cup \text{LM}_\prec(\mathcal{G}))$ sequence queries and $O((\#S+\#\mathcal{G})^2\#2S)$ operations to recover $\mathcal{G}$, where $2S$ is the Minkowski sum of $S$ with itself.

## 4.3 Benchmarks

In Figure 4.1, we consider three families of $\prec_{\text{DRL}}$-Gröbner basis based on the shape of their staircase:

**Rectangle:** $\text{LM}_\prec(\mathcal{G}) = \{y^{|d/2|}, x^d\}$ in dimension 2 and $\text{LM}_\prec(\mathcal{G}) = \{z^{|d/3|}, y^{|d/2|}, x^d\}$ in dimension 3. This case is the best for the size of the Gröbner basis compared to the size of the staircase.

**L shape:** $\text{LM}_\prec(\mathcal{G}) = \{x y, y^d, x^d\}$ in dimension 2 and $\text{LM}_\prec(\mathcal{G}) = \{y z, x z, x y, z^d, y^d, x^d\}$ in dimension 3. This case is the worst for the number of sequence queries compared to the sizes of the staircase and the Gröbner basis.

**Simplex:** $\text{LM}_\prec(\mathcal{G}) = \{y^d, x y^{d-1}, \dots, x^d\}$ in dimension 2 and $\text{LM}_\prec(\mathcal{G}) = \{z^d, y z^{d-1}, x z^{d-1}, \dots, y^d, x y^{d-1}, \dots, x^d\}$ in dimension 3, i.e. all the monomials of degree $d$. This case is the best for the number of sequence queries and the worst for the size of the Gröbner basis, both compared to the size of the staircase.

The polynomial variant of the BMS and SCALAR-FGLM algorithms performs fewer arithmetic

<footnote n="12">¹²INRIA</footnote>

<footnote n="13">¹³CryptoNext Security</footnote>
# Chapter 4. Guessing Gröbner bases

## 4.3. Benchmarks

operations than the others, for large *d*. More precisely, its number of operations appears to be linear in $(\#S)^2 = O(\#S(\#S + \#\mathscr{G}))$ in fixed dimension suggesting that the complexity estimate in Theorem 1.13 is pessimistic. This can be related to the uni-indexed case: the naive BM algorithm, based on naive extended Euclidean algorithm, is only quadratic in the number of sequence terms that one considers. This is the starting point of a work-in-progress of a new collaboration between Romain Lebreton<fnref n="14" /> and me that I detail in Section 6.2.

<footnote n="14">¹⁴Université de Montpellier</footnote>
This chapter is dedicated to quasi-commutative polynomials, their link with linear recurrence relations with polynomial coefficients, or Prelations of a sequence and how to compute or guess them.

5.1 P-relations and quasi-commutative polynomials

Linear recurrence relations with polynomial, in the indices, coefficients is a large class of lin-
ear recurrence relations containing the class of C-relations. As a classical example, the binomial
sequence $\binom{i_1}{i_2}$, $i_1, i_2 \in \mathbb{N}^2$ satisfies, for all $i_1, i_2 \in \mathbb{N}^2$,

$$
(i_1 - i_2 + 1) \binom{i_1 + 1}{i_2} = (i_1 + 1) \binom{i_1}{i_2}, \quad (i_2 + 1) \binom{i_1}{i_2 + 1} = (i_1 - i_2) \binom{i_1}{i_2}.
$$

The Eval operator, see also Section 1.3, can be extended to deal with these relations by adding a
new set of variables δ = (δ₁,...,δₙ) with the property that δℓ = xℓ ∂ℓ ∂xℓ for all 1 ≤ ℓ ≤ n. The 2n
variables satisfy the following quasi-commutative rules for all 1 ≤ k, ℓ ≤ n, k ≠ ℓ,

$x_k x_\ell = x_\ell x_k,$

$\delta_k \delta_\ell = \delta_\ell \delta_k,$

$x_k \delta_\ell = \delta_\ell x_k,$

$\delta_k x_k = x_k (\delta_k + 1).$

As a consequence, the polynomials in $\delta$ and $x$ with coefficients in $\mathbb{K}$ form a quasi-commutative ring, denoted $\mathbb{K}[\delta] \langle x \rangle$.

These polynomials allow us to represent P-relations through the Eval operator

$$
\text{Eval} \left( \sum_{r \in \mathbb{R}, s \in \mathcal{S}} \gamma_{r,s} \delta^r x^s, \mathbf{u} \right) = \sum_{r \in \mathbb{R}, s \in \mathcal{S}} \gamma_{r,s} s^r u_s.
$$

Observe that $\delta^r x^s \delta^k x^i = (\delta - s)^k \delta^r x^{s+i}$. Therefore, for any $g \in \mathbb{K}[\delta] \langle x \rangle$,

$$
\mathrm{Eval} (g \delta^k x^i, \mathbf{u}) = i^k \mathrm{Eval} (g x^i, \mathbf{u}).
$$

Thus, multiplying $g$ on the right by $\delta^k x^i$ shifts the recurrence relation by $i$ and multiplies it by $i^k$.
We deduce that if $g$ is such that Eval$(gx^i, u) = 0$ for all $i \in \mathbb{N}^n$, then Eval$(g\delta^k x^i, u) = 0$ for all
$k \in \mathbb{N}^n$ and $i \in \mathbb{N}^n$. Hence, the set of polynomials that represent Prelations satisfied by a sequence
$u$ is a right ideal of $\mathbb{K}[\delta] \langle x \rangle$ called the ideal of Prelations of $u$ and denoted $\mathrm{IP}(u)$.

From $\delta^r x^s \delta^k x^i = (\delta - s)^k \delta^r x^{s+i}$, we also notice that for $<$, a monomial order on the variables $\delta$ and $x$, we have $\mathrm{LM}_<(gh) = \mathrm{LM}_<(g) \mathrm{LM}_<(h)$ for any two polynomials $g$ and $h$. Furthermore, the notion of Gröbner basis can be defined for right ideals with the property that a finite subset $\mathcal{G}$ of an ideal $I$ is a Gröbner basis $I$ w.r.t. $<$ if for any $f \in I$, there exist $g \in \mathcal{G}$ and a monomial $m$ such that $\mathrm{LM}_<(f) = \mathrm{LM}_<(g)m$.

In particular, Buchberger's algorithm and criteria can be extended to this setting, see [75].

**Example 5.1:** From $(i_1 + 1 - i_2) \binom{i_1+1}{i_2} - (i_1 + 1) \binom{i_1}{i_2} = 0$ and $(i_2 + 1) \binom{i_1}{i_2+1} - (i_1 - i_2) \binom{i_1}{i_2} = 0$, we deduce that $g_1 = (\delta_1 - \delta_2)x_1 - (\delta_1 + 1)$ and $g_2 = \delta_2 x_2 - (\delta_1 - \delta_2)$ are in the ideal of P-relations of the binomial sequence.

Observe that in $\mathbb{K}(\delta)\langle x \rangle = \mathbb{K}(\delta) \otimes \mathbb{K}[\delta]\langle x \rangle$, that is $\delta$ and $x$ still quasi-commute but we allow rational fractions in $\delta$, these two polynomials are linear in $x$ and linearly independant. They thus span a zero-dimensional ideal that contains $g_1x_2 + g_2 = (\delta_1 - \delta_2)x_1x_2 - (\delta_1 - \delta_2 + 1)x_2 - (\delta_1 - \delta_2) = (x_1x_2 - x_2 - 1)(\delta_1 - \delta_2)$ and thus $x_1x_2 - x_2 - 1$ representing the, so-called Pascal's rule, C-relation:

$$
\forall (i_1, i_2) \in \mathbb{N}^2, \binom{i_1+1}{i_2+1} = \binom{i_1}{i_2+1} + \binom{i_1}{i_2}.
$$
### 5.2 Guessing P-relations

In the uni-indexed case and a sequence $u = (u_i)_{i \in \mathbb{N}}$, one can compute a basis of $I_P(u)$ using the Beckermann-Labahn algorithm [6] and its recent improvements [68], based on a divide-and-conquer approach. Many computer algebra systems propose an implementation of guessing P-relations algorithms, whether be it in the univariate or the multivariate case. We can cite for instance [99] in MAPLE [9], [71] in MATHEMATICA [67]. Some also have dedicated arithmetic for quasi-commutative polynomials like MAPLE or SAGE Math [103] through a package [69].

#### 5.2.1 A linear-algebra-based algorithm

With Jean-Charles Faugère<fnref n="1" />,<fnref n="2" />, I designed an extension to the SCALAR-FGLM algorithm in order to guess P-relations or C-relations between several sequences. For P-relations satisfied by $u$, the idea is to build a generalized multi-Hankel matrix $H_{T,\Theta}$, where $T$ and $\Theta$ are sets of monomials in respectively $\mathbb{K}[x]$ and $\mathbb{K}[\delta]\langle x \rangle$. The entry at the intersection of the column labeled with $\delta^r x^s$ and the row labeled with $x^i$ is Eval $(\delta^r x^s x^i, u) = (s + i)^r u_{s+i}$. Therefore, a vector in the right kernel of this matrix represents a polynomial $G \in \mathbb{K}[\delta]\langle x \rangle$ such that Eval $(Gx^i, u) = 0$ for all $x^i \in T$.

**Example 5.2:** Consider the binomial sequence $\binom{\binom{i_1}{i_2}}{\binom{i_1}{i_2}}_{i_1, i_2} \in \mathbb{N}^2$, $\Theta = \{1, \delta_2, \delta_1, x_2, x_1, \delta_2 x_2, \delta_1 x_2, \delta_2 x_1, \delta_1 x_1\}$ is the set of all monomials of bidegree at most (1,1) in $\delta$ and $x$ for an order $<$ eliminating $x$ and breaking ties on $\delta$ and on $x$ using $<_{\text{DRL}}$. We take $T$ to be the set of all monomials of degree at most 3 so that $T$ is larger than $\Theta$.

$$ H_{T,\Theta} = \begin{pmatrix} 1 & \delta_2 & \delta_1 & x_2 & x_1 & \delta_2 x_2 & \delta_1 x_2 & \delta_2 x_1 & \delta_1 x_1 \\ x_2 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 1 \\ x_1 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 1 \\ x_2^2 & 1 & 0 & 1 & 1 & 1 & 1 & 0 & 2 \\ x_1 x_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 4 \\ x_1^2 & 1 & 0 & 2 & 2 & 1 & 2 & 4 & 3 \\ x_2^3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ x_1 x_2^2 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 2 \\ x_1^2 x_2 & 2 & 2 & 4 & 1 & 3 & 2 & 2 & 9 \\ x_1^3 & 1 & 0 & 3 & 3 & 1 & 3 & 9 & 4 \end{pmatrix} $$

From this matrix, we can see the relations $\delta_2 x_2 - (\delta_1 - \delta_2)$ and $(\delta_1 - \delta_2)x_1 - (\delta_1 + 1)$.

More recently, a divide-and-conquer algorithm was designed and proposed in [82] which generalizes the Beckermann-Labahn algorithm to multi-indexed sequences.

#### 5.2.2 An hybrid approach

In the combinatorics context, the ultimate goal is not always to guess the whole Gröbner basis of the ideal of P-relations for a specific monomial order, but rather to determine if the sequence is P-finite or not. For instance, a necessary condition is that the ideal of P-relations in $\mathbb{K}(\delta)\langle x \rangle$ is zero-dimensional. Therefore, providing an adaptive algorithm that minimizes the number of sequence terms is of upmost importance.

In a similar fashion as in Section 4.2.3, we can discover, step by step, which monomials are in the staircase of the ideal of P-relations of $\mathbb{K}[\delta]\langle x \rangle$ of the input sequence for the input monomial order $<$. In [16], Jean-Charles Faugère<fnref n="1" />,<fnref n="2" /> and I proposed an hybrid approach based on Gröbner bases computations for quasi-commutative polynomials to discover new P-relations without any extra queries to the sequence. The idea is that if two polynomials $g_1, g_2 \in \mathbb{K}[\delta]\langle x \rangle$ encode P-relations satisfied by the sequence, then any polynomial in $\langle g_1, g_2 \rangle$ also encodes a P-relation. Therefore, as soon as two P-relations $g_1$ and $g_2$ are guessed, the goal is to compute a Gröbner basis $\{g_1, g_2, \dots, g_r\}$ of $\langle g_1, g_2 \rangle$. This will yield polynomials, namely $g_3, \dots, g_r$, whose leading monomials are not in $\langle \text{LM}_<(g_1), \text{LM}_<(g_2) \rangle$. The advantage of this method is twofold. First, since $\text{LM}_<(g_3), \dots, \text{LM}_<(g_r) >$

<footnote n="1">¹INRIA</footnote>

<footnote n="2">²CryptoNext Security</footnote>
$\text{LM}_<(g_1), \text{LM}_<(g_2)$, they require more queries to the sequence to be correctly guessed. Yet, such a Gröbner basis computation does not require any more queries. Then, these P-Relations may help us determine that the ideal of P-relations is 0-dimensional in $\mathbb{K}(\delta) \langle x \rangle$, which is a necessary condition for the table to be P-finite.

## 5.3 Structured sequences

In many applications, like in enumerative combinatorics with 2D/3D-walks, the studied sequences are highly structured. For instance, they have a lot of zero-terms. From a linear algebra viewpoint, these zero terms induce void conditions on the C-relations or P-Relations and thus drops of ranks in the generalized multi-Hankel matrix that one deals with to guess the relations. Therefore, we need to build matrices with many more rows than columns to recover correct relations, i.e. C-relations or P-Relations that do not prove to be fake after further testing.

The goal is thus to only deal with the nonzero terms from the beginning. These nonzero terms lie in general in a cone $\mathcal{C}$ which is a submonoid of $\mathbb{N}^n$, i.e. $0 \in \mathcal{C}$ and for all $i, j \in \mathcal{C}$, $(i + j) \in \mathcal{C}$.

In [21], Mohab Safey El Din<fnref n="3" /> and I extended the SCALAR-FGLM algorithm, and its variants, adaptive and for P-Relations, so that it only considers terms of a sequence lying in a cone. In this paper, we make the connexion with sparse Gröbner bases as defined in [7, 56] and extend them to the context of quasi-commutative polynomials in $\mathbb{K}[\delta] \langle x \rangle$ for P-Relations.

Given such a cone $\mathcal{C}$ and polynomials with support in its associated set of monomials $\mathcal{T}(\mathcal{C}) = \{x^i \in \mathcal{T} | i \in \mathcal{C}\}$, one may want to perform all the polynomial operations with monomials in $\mathcal{T}(\mathcal{C})$ in order to take advantage of the structure of the support when computing a Gröbner basis of the ideal they span. While this is not always possible, one can achieve this goal by considering the ideal the polynomials span in the algebra

$$ \mathbb{K}[\mathcal{C}] = \left\{ f = \sum_{s \in \mathcal{C}} f_s x^s \mid \text{supp } f \text{ is finite} \right\}. $$

The ideal spanned by $f_1, \dots, f_m \in \mathbb{K}[\mathcal{C}]$ is defined as

$$ \langle f_1, \dots, f_m \rangle_{\mathcal{C}} = \left\{ \sum_{k=1}^{m} f_k q_k \mid q_1, \dots, q_m \in \mathbb{K}[\mathcal{C}] \right\}. $$

Given a monomial order $<$, a $<$-sparse Gröbner basis of $\langle f_1, \dots, f_m \rangle_{\mathcal{C}}$ is a finite subset thereof whose properties extend those satisfied by a Gröbner basis, in the context of the algebra $\mathbb{K}[\mathcal{C}]$. In particular, the divisibility property on the leading monomials in $\mathcal{T}(\mathcal{C})$ is still satisfied.

**Definition 5.3 ([56, Def. 3.1] and [7, Def. 3.3]):** Let $\mathcal{C} \subseteq \mathbb{N}^n$ be a submonoid, $f_1, \dots, f_m \subseteq \mathbb{K}[\mathcal{C}]$ be polynomials and $<$ be a monomial order. Let $I = \langle f_1, \dots, f_m \rangle_{\mathcal{C}}$. Then, a $<$-sparse Gröbner basis of $I$ is a generating set $\mathcal{G} = \{g_1, \dots, g_l\} \subseteq \mathbb{K}[\mathcal{C}]$ such that for all $f \in I$, $\text{LM}_<(f) = \text{LM}_<(g)x^i$ for some $g \in \mathcal{G}$ and $x^i \in \mathcal{C}$.

The associated staircase of $\mathcal{G}$ is the set of monomials $s$ in $\mathcal{T}(\mathcal{C})$ such that for any $g \in \mathcal{G}$, there is no monomial $x^i \in \mathcal{T}(\mathcal{C})$ such that $s = \text{LM}_<(g)x^i$.

Let us notice that for $\mathcal{C} = \mathbb{N}^n$, $\mathbb{K}[\mathcal{C}] = \mathbb{K}[\mathbb{N}^n] = \mathbb{K}[x]$ and sparse Gröbner bases are classical Gröbner bases. Furthermore, like classical Gröbner bases, sparse Gröbner bases allow one to solve the ideal membership problem in $\mathbb{K}[\mathcal{C}]$ in an effective way.

**Example 5.4:** The cone $\mathcal{C} = \{(i_1, i_2) \in \mathbb{N}^2 | i_1 \le 2i_2, i_2 \le 2i_1\}$ is spanned by $(1,1), (2,1)$ and $(1,2)$ as a submonoid of $\mathbb{N}^2$. Although $x_1x_2$ divides both $x_1^2x_2$ and $x_1x_2^2$ in $\mathbb{K}[x]$, it does not in $\mathbb{K}[\mathcal{C}]$ as the respective quotients are $x_1$ and $x_2$, which are not in $\mathcal{T}(\mathcal{C})$.

**Theorem 5.5 (see also [21, Th. 3.2]):** Let $\mathcal{C}$ be a submonoid cone of $\mathbb{N}^n$ spanned by a finite minimal set of generators. Let $<$ be a monomial ordering on $\mathcal{T}$, the set of monomials in $n$ variables, and let $T \subset \mathcal{T}(\mathcal{C})$ be a set of monomials ordered for $<$ stable by division.

Then, the SCALAR-FGLM algorithm called on sequence $u$, $T$ and $<$ returns a set of polynomials $G$ with support in $\mathcal{T}(\mathcal{C})$, such that for all $s \in T \setminus \langle \text{LM}_<(G) \rangle$, $s$ is in the associated staircase of a

<footnote n="3">³Sorbonne Université</footnote>
### <-sparse Gröbner basis of $I_C(u)$.

Furthermore, if the ideal of C-relations of $u$ is 0-dimensional and has a reduced <-sparse Gröbner basis with support in $T$, then the output of the SCALAR-FGLM algorithm called on $u$ and $T$ is this reduced <-sparse Gröbner basis.

As an illustration, we consider walks in 2D- and 3D-spaces, namely Gessel planar walk $g$ in the nonnegative quadrant $\mathbb{N}^2$ with steps in $\{(1,0), (1,1), (-1,0), (-1,-1)\}$ and the 3D-space Walk-43 $w$ of [25] in the nonnegative octant $\mathbb{N}^3$ with steps in $\{(-1,-1,-1), (-1,-1,1), (-1,1,0), (1,0,0)\}$. In particular, we restrict ourselves to a subsequence of each where one index is 0. These walks come naturally with a cone structure: for instance whenever $n \neq 2n' + 2j$, then $g_{n,0,j} = 0$. Likewise, whenever $n \neq 8n' + 2j + 4k$, then $w_{n,0,j,k} = 0$. Thus, it makes sense to look for the relations given by the sequence terms $g_{2n'+2j,0,j}$ and $w_{8n'+2j+4k,0,j,k}$.

<table><caption>Table 5.1: Guessing fake and correct Relations</caption><thead><tr><th rowspan="2">Type</th><th colspan="4">Cone</th><th colspan="4">Full Orthant</th></tr><tr><th>Matrix size</th><th>Queries</th><th colspan="2">Relations</th><th>Matrix size</th><th>Queries</th><th colspan="2">Relations</th></tr><tr><th></th><th></th><th></th><th>Fake</th><th>Correct</th><th></th><th></th><th>Fake</th><th>Correct</th></tr></thead><tbody><tr><td>$g_{n,0,j}$</td><td>444 x 441</td><td>866</td><td>11</td><td>0</td><td>496 x 495</td><td>946</td><td>48</td><td>0</td></tr><tr><td>$g_{n,0,j}$</td><td>631 x 564</td><td>1 174</td><td>0</td><td>0</td><td>1 326 x 661</td><td>1 942</td><td>84</td><td>0</td></tr><tr><td>$g_{n,0,j}$</td><td>721 x 711</td><td>1 408</td><td>15</td><td>8</td><td>726 x 715</td><td>1 386</td><td>67</td><td>0</td></tr><tr><td>$g_{n,0,j}$</td><td>1 951 x 1 089</td><td>3 010</td><td>0</td><td>21</td><td>2 556 x 1 001</td><td>3 491</td><td>136</td><td>6</td></tr><tr><td>$w_{n,0,j,k}$</td><td>223 x 211</td><td>430</td><td>7</td><td>1</td><td>220 x 210</td><td>395</td><td>24</td><td>0</td></tr><tr><td>$w_{n,0,j,k}$</td><td>444 x 253</td><td>552</td><td>2</td><td>1</td><td>680 x 267</td><td>912</td><td>37</td><td>0</td></tr><tr><td>$w_{n,0,j,k}$</td><td>406 x 400</td><td>799</td><td>11</td><td>6</td><td>406 x 400</td><td>771</td><td>27</td><td>0</td></tr><tr><td>$w_{n,0,j,k}$</td><td>806 x 522</td><td>1 320</td><td>2</td><td>6</td><td>1 540 x 589</td><td>2 073</td><td>68</td><td>0</td></tr></tbody></table>

In Table 5.1, we can see that considering only terms in the cone allows us to almost only discover correct relations compared to the full orthant case, even when the generalized multi-Hankel matrices are almost square.
I want to develop my research project on computing and guessing Gröbner bases in the com-
mutative and quasi-commutative setting relying on MSOLVE as a validation tool. The main goal
is to develop general implementations on the one hand and specific ones for applications from
combinatorics, cryptography or robotics on the other hand.

I present my research project in three axes. The first one, in Section 6.1, is about MSOLVE and its
development. I see it as a transverse axis. The goal is to implement existing algorithms presented
in the previous chapters but also new ones that shall be designed in the other axes. Furthermore,
I shall exploit high-performance computing hardware to enhance the computation capabilities of
MSOLVE.

As a second axis, Section 6.2, I aim to accelerate guessing algorithms, especially for C-relations,
in order to provide theoretical and practical complexity estimates on the SPARSE-FGLM algorithms
when roots have multiplicities. In particular, we will trade the information on the multiplicities for
the efficiency and vice versa, when a second multiplication matrix is required. A polynomial matrix
extension is also targeted.

The last axis, Section 6.3, is dedicated to Gröbner bases for quasi-commutative polynomials in
order to deal with P-relations. The main objective is to bring the complexity of <DRL>-Gröbner bases
for quasi-commutative polynomials to that of classical polynomials.

**6.1 MSOLVE and fundamental algorithms**

A first goal for MSOLVE is to increase its functionalities by implementing state-of-the-art tech-
niques for general ideals which are not covered yet. Furthermore, at the moment the modular
implementation uses native integer types which allow us to handle large characteristic but with
computations that are slower than those using floating-point types.

**6.1.1 Types for the coefficients**

For the moment, MSOLVE uses native integer types to represent modular integers with primes
of size at most 31 bits. While some fast Central Processing Unit (CPU) instructions (like AVX2
or AVX512) are available for integer types, current CPUs are such that even faster instructions are
available to floating-point types. These types allow us to handle modular arithmetic for lesser prime
(up to 26 bits for double-precision floating-point numbers). Assuming not too many bad primes
are picked during a multi-modular Gröbner basis computation, using these types could speed the
computations up before lifting the result over Q.

Likewise, for cryptographic applications, we shall need a dedicated implementation over $\mathbb{F}_2$. In fact, we may have to introduce ways to compute over non-prime finite fields. This is why, I want to integrate fast univariate polynomial-matrix operations in order to deal with computations over finite extension of finite fields, defined by a single polynomial.

**6.1.2 Gröbner bases for a degree order**

While the classical framework of polynomial system solving relies on computing first a <sub>DRL</sub>-Gröbner basis, see Figure 1.1, in some applications, a monomial order < distinct from <sub>DRL</sub> might be more convenient for the computation. This can be a 2-block <sub>DRL</sub>-order, which is already implemented in MSOLVE, or a weighted degree order, see also [54, 55, 105]. When the input system is quasi-homogeneous for a system of weights, one can take advantage of this quasi-homogeneity to improve the runtime and the complexity estimates of the $F_4$ and $F_5$ algorithms [47, 48] by a factor depending on the products of the weight. We want to benefit from this by implementing this kind of weighted degree monomial orders in MSOLVE in order to tackle systems coming from applications such as robotics where the weights are naturally coming from dimensional homogeneity.

Another approach to speed Gröbner bases computations up for (quasi-)homogeneous systems
is to use *Hilbert-driven* Gröbner bases algorithms. These algorithms, such as *Traverso's* [104], see
also [41, Chap. 10, Sec. 2], assume that the Hilbert series of the ideal is already known and provide a criterion based on the degree of a critical pair to discard it without computing its reduction to 0. Combined with the $F_4$ algorithm, this allows us to handle Macaulay-like matrices with both fewer rows and fewer columns. Note, however, that the expected speed-ups may only be observed in positive characteristic or, in the case of characteristic 0, during the learning phase modulo the first prime.

### 6.1.3 Change of order

For the moment, the implemented change of order step relies on the sparse-FGLM algorithm [52, 53] in the shape position case. This assumption rises issues when the solutions of the system have multiplicities: even after introducing a generic linear form, the ideal might not be in shape position. Thus we may not compute the $<_{LEX}$-Gröbner basis of the ideal. While we can compute the parametrizations of the solutions, through the radical of the ideal, this makes us lose the information on the multiplicity. To this end, we want to provide more general change of order algorithms in `MSOLVE`. A first step will be to implement the seminal FGLM algorithm [50] before considering the algorithms I developed based on linear algebra and in particular, the adaptive algorithm of [11, 12], see also Section 4.2.2.

In order to provide an efficient implementation of the guessing algorithms, I want to use the structure of the multi-Hankel matrix and in particular, its *quasi-Hankel* structure. Indeed, if the univariate polynomial in $x_n$ in the reduced $<_{LEX}$-Gröbner basis has a large degree, i.e. of the same magnitude order as the degree of the ideal $D$, the multi-Hankel matrix will be quasi-Hankel with a small displacement rank. We can thus compute its vector in its kernel, and thus C-relations, fast using [26, 27]. Furthermore, depending on the geometry of the $<_{DRL}$- and $<_{LEX}$-staircases, computing the coefficients of this multi-Hankel matrix can just require Algorithm 2.2 with only one matrix, the one of the multiplication by $x_n$, but extra column-vectors. All-in-all, when $D$ is exponential in $n$, the expected complexity will still be dominated by the computation of the vectors of Algorithm 2.2 in $O(tD^2)$ operations.

In addition to this, I want to implement a multi-modular approach for sequences with coefficients in $\mathbb{Q}$. The idea is that the purpose of the adaptive algorithm is to find, at a low cost, the support of the reduced $<_{LEX}$-Gröbner basis. Once this has been computed modulo a first prime, we can switch to the general `SCALAR-FGLM` with one big multi-Hankel matrix whose rows and columns are labeled by the monomials in this support. This shall allow us to compute exactly the number of sequence terms required to recover the sought Gröbner basis modulo subsequent primes.

### 6.1.4 Saturation and colon ideals

With my coauthors, we want to push forward the efficiency of $F_4$SAT by investigating computational tricks to avoid “computing zero”. For instance, we want to provide an early termination criterion to reduce the time spent in the last saturation step.

Likewise, we want to study a signature approach, as in Faugère's $F_5$ algorithm [44, 48], to minimize the reductions to 0 of S-polynomials.

I also want to provide a multi-modular approach and a tracer for `SPARSE-FGLM-COLON` for systems over $\mathbb{Q}$ in `MSOLVE`. The goal is to learn how to minimize $\Sigma$ and $\Sigma'$ to the actual set of required monomials to write the $<_{LEX}$-Gröbner basis of the colon ideal but also to compute the exact number of necessary sequence terms in the Wiedemann part of the algorithm. This will make computations in the apply phase modulo all the primes but the first one optimal and thus the fastest implementation.

### 6.1.5 Exact solving on a GPU

Graphics Processing Units (GPUs) are, by design, well-suited to process large blocks of data in parallel and thus to perform linear algebra routines, more so than CPUs. Furthermore, they natively handle double-precision floating-point number arithmetic but only simulate long integer ones through short integer arithmetic, which comes with an overhead. For instance, NVIDIA CUDA
and TENSOR cores, natively, only have 8-bit integer types, whereas they, natively, have 64-bit floating-point types<fnref n="1" />. Furthermore, it has been shown that GPUs are very efficient for correctly rounding functions evaluations [59], reinforcing even more the advantage of using floating-point arithmetic in order to simulate a modular one. Thus, we aim to take advantage of their computational power to transpose the linear algebra code of our polynomial system solver MSOLVE to efficiently work with GPUs.

The main objective of this task is the design of fast Gröbner bases computation algorithms, based on high performance linear algebra algorithms, in particular exploiting GPUs, and their integration into MSOLVE in order to tackle applications challenges such as multivariate cryptography or robotics. This is Dimitri Lesnoff<fnref n="2" />'s Ph.D. subject, who I have been supervising jointly with Stef Graillat<fnref n="3" /> and Théo Mary<fnref n="4" /> since October 2022. We shall first revisit modular arithmetic at the core of exact algorithms relying on fast low-precision arithmetics such as fp16, bfloat16, fp8, ... for floating-points numbers and int8, int4 for integers. This part will be done in parallel for CPUs, as in Section 6.1.1 and for GPUs. Then, we will adapt the block-Wiedemann approach [40, 66] of the SPARSE-FGLM algorithm. Since the matrix at hand is very particular, it can be seen as the concatenation of a permutation matrix and a dense matrix after reordering the columns, we will study how to balance more efficiently the CPU load and the GPU load to iterate the product of this matrix with some vectors or very thin matrices. We also target the design of a sparse matrix arithmetic which is both efficient for handling Macaulay matrices for the $F_4$ algorithm and dedicated to a GPU or a CPU + GPU architecture.

In order to tackle very large problems whose solving will unlock new advances in critical applications, we have to handle large matrices that cannot even be stored on the GPU. We plan to devise new algorithms that exploit a memory representation of these Macaulay matrices that suits our computations but also the limited RAM, a few tens of gigabytes, of a graphic card. Moreover, computing Gröbner bases at scale will require the use of multiple GPUs and CPUs in parallel. We will therefore work on making the algorithms scalable in a parallel context. Notably, we will minimize memory communication between the CPU and the GPU, by adapting cache-oblivious storing and algorithms [2] to a larger scale such as the RAM of the graphic card.

## 6.2 Faster change of order and guessing C-relations

### 6.2.1 Guessing faster

While the polynomial point of view allows us to bring the guessing of C-relations for multi-indexed sequences closer to the fast guessing for uni-indexed sequences, there is still a complexity gap. Indeed, the stated complexity in Theorem 1.13 becomes cubic, instead of quasi-linear, for uni-indexed sequences. This is due to the complexity analysis that relies, first, on a linear number of polynomial reductions and, second, on naive polynomial reductions. In fact, this cubic complexity stated in Theorem 1.13 comes from an overestimation on the number of needed polynomials to reduce a new one, while in practice, this is not observed, especially on sequences of prescribed ideal of C-relations but with *random* initial terms. Thus, as an on-going collaboration with Romain Lebreton<fnref n="5" />, I want to improve this algorithm by handling univariate polynomials instead of multi-variate ones, in order to take advantage or their fast arithmetic while reestimating the number of reductions to target a complexity which is only quadratic in the number of input sequence terms or the size of the output Gröbner basis.

As a second goal, I want to guess C-relations using polynomial arithmetic instead of linear algebra [16]. From the generating series of a sequence and its derivatives, or more precisely the mirror polynomials of a truncation of these series, the goal is to find algebraic combinations thereof which are small modulo the monomial ideal $\langle x_1^{D_1}, \dots, x_n^{D_n} \rangle$, where $D_1, \dots, D_n$ depend on the sequence terms we allow ourselves to use.

<footnote n="1">¹https://www.nvidia.com/en-us/data-center/tensor-cores/</footnote>

<footnote n="2">²Ph.D. student at Sorbonne Université</footnote>

<footnote n="3">³Sorbonne Université</footnote>

<footnote n="4">⁴CNRS</footnote>

<footnote n="5">⁵Université de Montpellier</footnote>
Efficiency of guessing algorithms is based on two aspects: the number of performed operations and the number of sequence terms that are needed. Indeed, in many applications, computing the sequence is the bottleneck. Thus, to make this approach the most efficient, we shall closely look at the number of different sequences terms that are needed to correctly guess the relations. Furthermore, to optimize the number of operations, we shall rely on efficient algorithms for univariate polynomials and uni-indexed sequences. The main goal is to reach a complexity at most quadratic in the size of the output instead of only cubic.

These guessing algorithms may find fake relations, this happens in general when too few terms are used or when most of the terms are 0. This can be circumvented by *structured* guessing relations using mainly the nonzero sequence terms as in [21]. We will pay attention to these bad sequences so that our new algorithm avoids these fake relations as much as possible.

## 6.2.2 Guessing radical ideals

Radical ideals are important in the polynomial system solving applications as the solutions have all multiplicities 1. Therefore, numerical methods to approximate the real or complex solutions behave much better with this kind of ideals.

In collaboration with Alin Bostan<sup>6</sup>, Manuel Kauers<sup>7</sup> and Christoph Koutschan<sup>8</sup>, I want to investigate an extension of [29] to more general radical ideals. Indeed, in this paper, given a sequence *u*, the authors propose an algorithm for guessing the <sub>LEX</sub>-Gröbner basis of $\sqrt{I_C(\mathbf{u})}$ if it is in shape position. Yet, the output is meaningless if this radical is not in shape position. While BMS and SCALAR-FGLM allow us to guess $I_C(\mathbf{u})$, this may require many sequence terms, especially if the ideal has a large degree. Thus, we want to propose a trade-off between the number of queries to the sequence and the quality of the output by allowing to guess polynomials only in $\sqrt{I_C(\mathbf{u})}$.

## 6.2.3 Guessing with multiplicities

In some situations where the system has roots with multiplicities, the sought <sub>LEX</sub>-Gröbner basis *cannot* be in shape position, even after a generic linear change of variables. Such a system is called 2-thick in [5] and it requires a second structured $D \times D$-matrix, of a similar kind, to be computed: the matrix of the multiplication by $x_{n-1}$, which has $O(\tau D)$ nonzero coefficients. Furthermore, in most situations the stability property is still satisfied which means that the computation of this second matrix is cheap. However, it does not ensure (actually it *almost never can*) that this second matrix is computed for free.

A first goal is to derive a sharp complexity estimate on the computation of this second matrix based on the <sub>DRL</sub>-Gröbner basis, exploiting the stability property and the work of Moreno-Socías [80], for the first matrix. For instance, whenever $\langle \mathrm{LM}_{<_{LEX}}(I) \rangle = \langle x_n^{D_n}, x_{n-1}^{D_{n-1}}, x_{n-2}, \dots, x_1 \rangle$ or $\langle \mathrm{LM}_{<_{LEX}}(I) \rangle = \langle x_n^{D_n}, x_n x_{n-1}, x_{n-1}^2, x_{n-2}, \dots, x_1 \rangle$ which seem to be simplest cases.

As a by-product, we will obtain complexity bounds on the computation of the sequence terms that appear in the multi-Hankel matrix built by SCALAR-FGLM [11, 12] to recover the <sub>LEX</sub>-Gröbner basis. Then, as a second goal, we will rely on the quasi-Hankel structure of this multi-Hankel matrix, and fast algorithms for quasi-Hankel matrices [26, 27], to analyze the complexity on the computation of the sought <sub>LEX</sub>-Gröbner basis. All in all, we will have a complete description and complexity estimate of Faugère and Mou's [52, 53] SPARSE-FGLM algorithm for generic 2-thick systems.

## 6.2.4 Polynomial matrices with multiplicities

A first goal is to adapt the polynomial-matrix algorithm for the change of order, see [20] and Section 3.1.2, to take into account structures. For instance, in [57] and [21], the authors tweak respectively the FGLM and the SPARSE-FGLM algorithms for ideals that are globally invariant under the action of a finite abelian group. The main idea is to split the multiplication matrices into block-matrices, improving the complexity by a factor that depends on the size of the group.

Another objective is to be able to compute the <sub>LEX</sub>-Gröbner basis of a zero-dimensional ideal

<footnote n="6"><sup>6</sup>INRIA</footnote>

<footnote n="7"><sup>7</sup>Johannes Kepler Universität Linz</footnote>

<footnote n="8"><sup>8</sup>Osterreichische Akademie der Wissenschaften</footnote>
$I : \langle \varphi \rangle$ from the <DRL-Gröbner basis of a, potentially positive-dimensional, ideal $I$ using also the polynomial-matrix algorithm for change of order. We also want to investigate the computation of a parametrization of the radical ideal when $I$ is not radical, in order to remove multiplicities. How can we take advantage of the polynomials to add to $I$ to span $\sqrt{I}$ in this polynomial matrix in order to compute the <LEX-Gröbner basis of $\sqrt{I}$.

Then, a last objective, as a follow-up to Section 6.2.3, is to deal with the matrix of the multiplication by $x_{n-1}$ in this polynomial matrix change of order algorithm. That is, we will develop this approach by considering a polynomial matrix $P_{x_{n-1}}$ instead of $M_{x_{n-1}}$. However, this matrix should lie in $\mathbb{K}[x_{n-1}, x_n]$ and thus, efficient bivariate-polynomial matrices operations would have to be developed in order to make this algorithm faster, theoretically and practically, than the SPARSE-FGLM algorithm.

**6.2.5 Applications to guessing**

I also want to target applications such as coding theory and optimization problems. In the former case, I want to extend the existing algorithms for guessing recurrences in order to decode extensions of cyclic codes such as those defined on an algebraic variety. As a first step, I will consider algorithms relying on linear algebra in order to understand the structure of the matrix at hand. In the latter case, I want to adapt our algorithms, in particular SCALAR-FGLM, to handle numerical sequences coming from the moment approach [64]. The main difficulty will be to return non-trivial relations and we shall take inspiration from approximate GCD algorithms to overcome it, see [22].

**6.3 Algorithms for quasi-commutative Gröbner bases**

**6.3.1 Moment approach**

As a starting collaboration with Lorenzo Baldi<fnref n="9" /> and Pierre Lairez<fnref n="10" />, we want to use the theory of holonomic functions, or sequences satisfying P-relations, to provide a general and efficient method for computing the moments of measures that are useful for the applications. To do so, we want to exploit the fact that when computing the volume defined by several polynomial inequalities $f_1, \dots, f_r \ge 0$ in the hypercube $[0,1]^n$, we are led to computing the integrals $\int_{[0,1]^n} f_1^{k_1} \cdots f_r^{k_r} dx$ for many indices $k_1, \dots, k_r$. Yet, the sequence of integrals satisfies P-relations. These relations are fundamental and we want to guess them by computing first a few moments and then exploit them to compute more moments in order to find larger relations.

**6.3.2 Efficient Gröbner bases computation**

Another topic I want to investigate is the computation of Gröbner bases for a total degree monomial order of ideals in a quasi-commutative setting and its theoretical and practical efficiency in relation with Section 5.2.2. A longer term goal is then to implement Gröbner bases of quasi-commutative polynomials into MSOLVE in order to solve applications from combinatorics and physics.

Following [75] and the generalization of Buchberger's criteria, the goal will be to dive into the understanding of how Faugère's F₄ algorithm [47] behaves or can be extended from the commutative setting to this one.

To do so, I want to study the module of trivial syzygies in order to get information on the sizes of the matrices that are built in F₄. What kind of information the commutation rules provide on the syzygies?

Hilbert series and Hilbert polynomials are powerful tools that allow one to understand the complexity of computing Gröbner bases. In the commutative case, one can derive a bound on the degree of the polynomials in a reduced Gröbner basis for a total degree order, thanks to them, see [79, Section 4.5, Corollary]. We will investigate how knowing in advance the Hilbert series can speed the Gröbner bases computations up, or how together with the Hilbert polynomials, they can give us a bound on the degrees of the polynomials in the reduced <DRL-Gröbner basis.

<footnote n="9">⁹Post-doc at Sorbonne Université</footnote>

<footnote n="10">¹⁰INRIA</footnote>
**Chapter 6. Research project**

**6.3. Algorithms for quasi-commutative Gröbner bases**

In the particle physics and algebraic statistics application, computing representants of the quotient $K(\delta)\langle x \rangle/J$, where $J$ is a 0-dimensional ideal allows one to determine the twisted cohomology defined by a likelihood function, see [78]. In this paper, the authors need to compute the contiguity matrices of $J$, which actually correspond to the matrices of multiplication by $x_1, \dots, x_n$ in $K(\delta)\langle x \rangle/J$. Therefore, I want to study how to efficiently build these matrices from a $<_{\text{DRL}}$-Gröbner basis of $J$. Furthermore, depending on the shape of the $<_{\text{DRL}}$-staircase, some of these matrices are obtained from free, for instance the one for $x_n$, like in the commutative case, see [80]. As a longer term goal, I would like to study these conditions in this setting.
## A.1 Uni-indexed sequences

For $d \in \mathbb{N}^*$, it is clear that $S_C(\langle x^d \rangle) = \{(\alpha_0, \dots, \alpha_{d-1}, 0, 0, \dots) | \alpha_0, \dots, \alpha_{d-1} \in \mathbb{K}\}$, whereas for $a \in \mathbb{K}^*$, it is well-known that $S_C(\langle (x-a)^d \rangle) = \{((\alpha_0 + \dots + \alpha_{d-1}i^{d-1})a^i)_{i \in \mathbb{N}} | \alpha_0, \dots, \alpha_{d-1} \in \mathbb{K}\}$. These sequences $(\alpha(i)a^i)_{i \in \mathbb{N}}$ generalize the class of geometric sequences with initial term $\alpha_0$ and common ratio $a$, which are the case $d=1$, and the class of arithmetic sequences with initial term $\alpha_0$ and common difference $\alpha_1$, which are the case $d=2$ and $a=1$. For the sake of the completeness, we recall that if $u \in S_C(\langle (x-a)^d \rangle)$ for $a \in \mathbb{K}$ and $u \notin S_C(\langle (x-a)^{d-1} \rangle)$, then $I_C(u) = \langle (x-a)^d \rangle$. Such a sequence will be denoted $u_{a,d}$.

Using Theorem 1.5, we can deduce the general case.

**Theorem A.1:** Let $\mathbb{K}$ be algebraically closed and $g \in \mathbb{K}[x]$. Assume $g$ factors as

$$g = \prod_{k=1}^{m} (x - a_k)^{d_k},$$

then $u \in S_C(\langle g \rangle)$ if, and only if, there exist unique $u_{a_1,e_1}, \dots, u_{a_m,e_m}$ with $e_1 \le d_1, \dots, e_m \le d_m$ such that

$$u = \sum_{k=1}^{m} u_{a_k, d_k}.$$

As we can see, the roots of the polynomial, and their multiplicities, completely determine the form of the sequence terms.

## A.2 Bi-indexed sequences

From the uni-indexed case and small computations, we can see that $S_C(\langle (x_1 - a_1)^{d_1}, (x_2 - a_2)^{d_2} \rangle)$ is the vector space of sequences $u_{a_1,e_1} \otimes u_{a_2,e_2} = (\alpha_1(i_1)\alpha_2(i_2)a_1^{i_1}a_2^{i_2})_{(i_1,i_2) \in \mathbb{N}^2}$ where $e_1 \le d_1, e_2 \le d_2$ and $\deg \alpha_1 = e_1, \deg \alpha_2 = e_2$. While such a polynomial ideal has $(a_1, a_2)$ as its unique root with multiplicity $d_1 d_2$, this does not encompass all the cases. Before diving into multiplicities, we shall deal with the case of distinct single roots.

To simplify the presentation, we mostly consider the case $a_1 \neq 0$ and $a_2 \neq 0$, but the results on the ideals still hold whenever $a_1, a_2$ or both are 0.

### A.2.1 Single roots

Let us consider two distinct points $(a_1, a_2), (b_1, b_2) \in \mathbb{K}^2$, then we can look at the ideal of polynomials vanishing on these points. This is the ideal $I = \langle (x_2 - a_2)(x_2 - b_2), (a_2 - b_2)(x_1 - a_1) - (a_1 - b_1)(x_2 - a_2), (x_1 - a_1)(x_1 - b_1) \rangle$. Observe how, if $a_2 \neq b_2$, the first two polynomials form a Gröbner basis of $I$ for $x_2 < x_1$, otherwise the last two do. Hence, $S_C(I)$ is a two-dimensional vector space and $u \in S_C(I)$ is uniquely determined by $u_{0,0}$ and, $u_{0,1}$ in the former case or $u_{1,0}$ in the latter one.

More generally, the sequence obtained as a linear combination, with nonzero coefficients, of the sequences $(a_{1,1}^{i_1}a_{1,2}^{i_2})_{(i_1,i_2) \in \mathbb{N}^2}, \dots, (a_{k,1}^{i_1}a_{k,2}^{i_2})_{(i_1,i_2) \in \mathbb{N}^2}$ has its ideal of C-relations which is exactly the radical ideal of all polynomials vanishing on $(a_{1,1}, a_{1,2}), \dots, (a_{k,1}, a_{k,2})$.

It remains to deal with multiplicities. As we shall see, even the multiplicity-2 case is broader than for uni-indexed sequences.
### A.2.2 Root of multiplicity 2

The sequences whose ideal of C-relations only admits $(a_1, a_2)$ as a root, and with multiplicity 2, are all of the type $u = ((\alpha_{0,0} + \alpha_{1,0}i_1 + \alpha_{0,1}i_2)a_1^{i_1}a_2^{i_2})_{(i_1,i_2)\in\mathbb{N}^2}$ with $\alpha_{1,0}$ and $\alpha_{0,1}$ not both 0. The ideal of C-relations of such a sequence $u$ is

$$J_{(a_1,a_2)_2, (\alpha_{1,0},\alpha_{0,1})} = \langle (x_1 - a_1)^2, \alpha_{0,1} \hat{a}_2(x_1 - a_1) - \alpha_{1,0} \hat{a}_1(x_2 - a_2), (x_2 - a_2)^2 \rangle,$$

where $\hat{b}$ equals 1 if $b=0$ and $b$ otherwise. Notice that if $\alpha_{1,0} \neq 0$, the first two polynomials span the ideal, otherwise $\alpha_{0,1} \neq 0$ and the last two polynomials span it.

### A.2.3 Root of multiplicity 3

The sequences whose ideal of C-relations only admits $(a_1, a_2)$ as a root, and with multiplicity 3, are all of the type $u = ((\alpha_{0,0} + \alpha_{1,0}i_1 + \alpha_{0,1}i_2 + \alpha_{2,0}i_1^2 + \alpha_{1,1}i_1i_2 + \alpha_{0,2}i_2^2)a_1^{i_1}a_2^{i_2})_{(i_1,i_2)\in\mathbb{N}^2}$ with $\alpha_{2,0}\alpha_{0,2} = 4\alpha_{1,1}^2$ and $\alpha_{2,0}$ and $\alpha_{0,2}$ not both 0. The ideal of C-relations of such a sequence $u$ is

$$J_{(a_1,a_2)_3, (\alpha_{2,0},\alpha_{0,2})} = \langle (x_1 - a_1)^3, \\ \hat{a}_1^2(x_2 - a_2) - \beta_{1,2}\hat{a}_2(x_1 - a_1)^2 - \beta_{1,1}\hat{a}_1\hat{a}_2(x_1 - a_1), \\ \hat{a}_2^2(x_1 - a_1) - \beta_{2,2}\hat{a}_1(x_2 - a_2)^2 - \beta_{2,1}\hat{a}_1\hat{a}_2(x_2 - a_2), \\ (x_2 - a_2)^3 \rangle,$$

where the $\beta_{i,j}$'s are polynomials in the $\alpha_{i,j}$'s with $1 \le i + j \le 2$. Furthermore, the $\beta_{i,j}$'s are such that either the first two or the last two polynomials span the ideal.

Geometrically, there exists another ideal whose only root is $(a_1, a_2)$, and with multiplicity 3. This is $K_{(a_1,a_2)_3} = \langle (x_1 - a_1)^2, (x_1 - a_1)(x_2 - a_2), (x_2 - a_2)^2 \rangle$.

**Proposition A.2:** Let $(a_1, a_2) \in \mathbb{N}^2$ and $K_{(a_1,a_2)_3} = \langle (x_1 - a_1)^2, (x_1 - a_1)(x_2 - a_2), (x_2 - a_2)^2 \rangle$. Then, there is no sequence $u$ such that $I_C(u) = K_{(a_1,a_2)_3}$.

*Proof.* Let $u$ be such that $K_{(a_1,a_2)_3} \subseteq I_C(u)$. Since $K_{(a_1,a_2)_3}$ contains $\langle (x_1 - a_1)^2, (x_2 - a_2)^2 \rangle$, $u = u_{a_1,\varepsilon_1} \otimes u_{a_2,\varepsilon_2}$ with $\varepsilon_1 \le 2$ and $\varepsilon_2 \le 2$. Assuming $a_1, a_2 \neq 0$, then there exist $\alpha, \beta, \gamma, \delta \in \mathbb{K}$ such that

$$\forall (i_1, i_2) \in \mathbb{N}^2, \quad u_{i_1, i_2} = (\alpha + \beta i_1 + \gamma i_2 + \delta i_1 i_2) a_1^{i_1} a_2^{i_2}.$$

Since $x_1x_2 - a_2x_1 - a_1x_2 + a_1a_2 = (x_1 - a_1)(x_2 - a_2) \in K_{(a_1,a_2)_3}$,

$$\begin{align*} u_{1,1} &= a_2 u_{1,0} + a_1 u_{0,1} - a_1 a_2 u_{0,0} \\ (\alpha + \beta + \gamma + \delta)a_1 a_2 &= (\alpha + \beta)a_1 a_2 + (\alpha + \gamma)a_1 a_2 - \alpha a_1 a_2 \\ \delta a_1 a_2 &= 0 \end{align*}$$

and $\delta = 0$. Now, we saw that if $\gamma \neq 0$, then $I_C(u) \supseteq J_{(a_1,a_2)_2, \beta/\gamma} = \langle x_1 - \frac{\beta}{\gamma}x_2 + \left(\frac{\beta}{\gamma} - 1\right), (x_2 - a_2)^2 \rangle \supseteq K_{(a_1,a_2)_3}$ and otherwise $I_C(u) \supseteq J_{(a_1,a_2)_2, \infty} = \langle (x_1 - a_1)^2, x_2 - a_2 \rangle \supseteq K_{(a_1,a_2)_3}$.

If either $a_1 = 0$ or $a_2 = 0$, then this proof can be adapted using the appropriate sequence $u_{a_1,\varepsilon_1} \otimes u_{a_2,\varepsilon_2}$. This concludes the proof that no sequence has $K_{(a_1,a_2)_3}$ as its ideal of C-relations. $\square$
[1] J. Alman and V. Vassilevska Williams. A Refined Laser Method and Faster Matrix Multiplication. In *Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA)*, pages 522–539, 2021.

[2] M. Bader and C. Zenger. Cache oblivious matrix multiplication using an element ordering based on a peano curve. *Linear Algebra and its Applications*, 417(2):301–313, 2006. Special Issue in honor of Friedrich Ludwig Bauer.

[3] C. Banderier and P. Flajolet. Basic analytic combinatorics of directed lattice paths. *Theoret. Comput. Sci.*, 281(1-2):37–80, 2002. Selected Papers in honour of Maurice Nivat.

[4] D. A. Bayer. *The Division Algorithm and the Hilbert Scheme*. PhD thesis, Harvard University, USA, 1982. AAI8222588.

[5] E. Becker, T. Mora, M. G. Marinari, and C. Traverso. The shape of the shape lemma. In *Proceedings of the International Symposium on Symbolic and Algebraic Computation*, ISSAC '94, pages 129–133, New York, NY, USA, 1994. ACM.

[6] B. Beckermann and G. Labahn. A uniform approach for the fast computation of matrix-type pade approximants. *SIAM J. Matrix Anal. Appl.*, 15(3):804–823, 1994.

[7] M. R. Bender, J.-Ch. Faugère, and E. Tsigaridas. Towards mixed Gröbner basis algorithms: The multihomogeneous and sparse case. In *Proceedings of the 2018 ACM International Symposium on Symbolic and Algebraic Computation*, ISSAC '18, pages 71–78, New York, NY, USA, 2018. ACM.

[8] E. Berlekamp. Nonbinary BCH decoding. *IEEE Trans. Inform. Theory*, 14(2):242–242, 1968.

[9] L. Bernardin, P. Chin, P. Demarco, K. O. Geddes, D. E. G. Hare, K. M. Heal, G. Labahn, J. Mccarron, M. B. Monagan, D. Ohashi, and S. M. Vorkoetter. Maple programming guide, 1996.

[10] J. Berthomieu, A. Bostan, A. Ferguson, and M. Safey El Din. Gröbner bases and critical values: The asymptotic combinatorics of determinantal systems. *J. Algebra*, 602:154–180, 2022.

[11] J. Berthomieu, B. Boyer, and J.-Ch. Faugère. Linear algebra for computing gröbner bases of linear recursive multidimensional sequences. In *Proceedings of the 2015 ACM on International Symposium on Symbolic and Algebraic Computation*, ISSAC '15, pages 61–68, New York, NY, USA, 2015. ACM.

[12] J. Berthomieu, B. Boyer, and J.-Ch. Faugère. Linear Algebra for Computing Gröbner Bases of Linear Recursive Multidimensional Sequences. *Journal of Symbolic Computation*, 83(Supplement C):36–67, 2017. Special issue on the conference ISSAC 2015: Symbolic computation and computer algebra.

[13] J. Berthomieu, Ch. Eder, and M. Safey El Din. Msolve: A library for solving polynomial systems. In *Proceedings of the 2021 on International Symposium on Symbolic and Algebraic Computation*, ISSAC '21, pages 51–58, New York, NY, USA, 2021. Association for Computing Machinery.

[14] J. Berthomieu, Ch. Eder, and M. Safey El Din. msolve: A library for solving polynomial systems, 2021. https://msolve.lip6.fr/.

[15] J. Berthomieu, Ch. Eder, and M. Safey El Din. New efficient algorithms for computing Gröbner bases of saturation ideals (F4SAT) and colon ideals (Sparse-FGLM-colon). preprint, 2022.