# A Direttissimo Algorithm for Equidimensional Decomposition

Christian Eder, Pierre Lairez, Rafael Mohr, Mohab Safey El Din

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

## To cite this version:

Christian Eder, Pierre Lairez, Rafael Mohr, Mohab Safey El Din. A Direttissimo Algorithm for Equidimensional Decomposition. ISSAC 2023 - 48th International Symposium on Symbolic and Algebraic Computation, Jul 2023, Tromsø, Norway. hal-04000248v2

HAL Id: hal-04000248

https://hal.science/hal-04000248v2

Submitted on 9 Jun 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A Direttissimo Algorithm for Equidimensional Decomposition

Christian Eder * Pierre Lairez † Rafael Mohr *‡ Mohab Safey El Din ‡

June 9, 2023

## Abstract

We describe a recursive algorithm that decomposes an algebraic set into locally closed equidimensional sets, i.e. sets which each have irreducible components of the same dimension. At the core of this algorithm, we combine ideas from the theory of triangular sets, a.k.a. regular chains, with Gröbner bases to encode and work with locally closed algebraic sets. Equipped with this, our algorithm avoids projections of the algebraic sets that are decomposed and certain genericity assumptions frequently made when decomposing polynomial systems, such as assumptions about Noether position. This makes it produce fine decompositions on more structured systems where ensuring genericity assumptions often destroys the structure of the system at hand. Practical experiments demonstrate its efficiency compared to state-of-the-art implementations.

## 1 Introduction

**Problem statement** Let $\mathbb{K}$ be an algebraically closed field, let $R = \mathbb{K}[x_1, \dots, x_n]$ be a polynomial ring and let $f_1, \dots, f_c \in R$ be a polynomial system generating an ideal $I \subseteq R$. The zero set $X$ of the polynomials $f_1, \dots, f_c$ in $\mathbb{K}$, decomposes uniquely as a union of irreducible algebraic sets such that none of them contains another. These are the *irreducible components* of $X$ and correspond to the *minimal associated primes* of $I$. The variety $X$ is *equidimensional* if all its irreducible components have the same dimension. It is clear that $X$ always admits a decomposition $X = Y_1 \cup \dots \cup Y_c$ where the $Y_i$ are equidimensional algebraic sets.

Suppose $X = \cup_{i=1}^s X_i$ is the decomposition of $X$ into irreducible components. For $k = 0, \dots, n$ define $Y_k := \cup_{i \text{ s.t. } \dim X_i = k} X_i$, then each $Y_k$ is equidimensional so that $X = \cup_{k=0}^n Y_k$ is an *equidimensional decomposition* of $X$. Given $f_1, \dots, f_c \in R$, we aim at computing such an equidimensional decomposition of $X$.

It will be clear that our algorithms will only have to do arithmetic over the subfield of $\mathbb{K}$ that the coefficients of $f_1, \dots, f_c$ lie in, the output will also be defined over the same subfield. In the following we will work only over $\mathbb{K}$ for the purpose of simplicity of presentation.

This problem finds natural applications in singularity analysis of sensor-based controllers of mechanism design (e.g. García Fontán et al., 2022; Pascual-Escudero et al., 2021, and references therein), in algorithms of real algebraic geometry (e.g. Aubry et al., 2002; Safey El Din & Schost, 2004) and real algebra (Safey El Din et al., 2018, 2021) as well as automated theorem proving and geometry (e.g. Chen et al., 2013; W. Wu & Gao, 2007; Yang et al., 1998, 2001).

**Prior Works** The importance of this computational problem fostered a vast body of literature often also as an intermediate step towards primary decomposition of ideals or prime decomposition of varieties. Algorithms for equidimensional decomposition of algebraic sets can be classified along the data structures which they employ to represent (equidimensional) algebraic sets.

<footnote>*RPTU Kaiserslautern-Landau, Germany</footnote>

<footnote>†Inria, Uni. Paris-Saclay, Palaiseau, France</footnote>

<footnote>‡Sorbonne Uni., CNRS, Paris, France</footnote>
There are two prominent strategies for equidimensional decomposition using Gröbner bases frequently implemented in computer algebra systems. The first one uses algebraic elimination techniques. It combines the knowledge of the dimension of the ideal generated by the input polynomials with the elimination theorem (Cox et al., 2015, Theorem 3.1.2) to compute a description of the projection of the algebraic set under study on a well-suited affine linear subspace to deduce how to split the corresponding ideal (Caboara et al., 1997; Decker et al., 1999; Gianni et al., 1988; Kalkbrener, 1994; Krick & Logar, 1991). The projection of the equidimensional component of highest dimension (frequently called the *equidimensional hull*) of the algebraic set under question will then be cut out by a hypersurface whose defining polynomial has degree equal to the degree of this equidimensional hull. As a consequence, such algorithms have the disadvantage that they need to manipulate polynomials of degree in the order of the Bezout bound of the input system. To circumvent this drawback, another set of methods, called *direct methods* has been introduced by Eisenbud et al. (1992). They rely on homological algebra to reduce the problem of equidimensional decomposition to the computation of syzygies which are then used to split the polynomial ideal under study, while avoiding projections. These algorithms often provide an intermediate step towards primary decomposition of ideals. For this problem, modular techniques and dedicated algorithms for the case where $K$ is a finite field have been designed (Ishihara, 2022; Noro & Yokoyama, 2004; Yokoyama, 2002).

Another body of work uses *lazy* representations of algebraic sets. Frequently, the core idea is to exploit the fact that any equidimensional algebraic set is locally almost everywhere a complete intersection, i.e. an equidimensional algebraic set of codimension $c$ can be represented by the vanishing of $c$ polynomials on a dense Zariski open subset of itself. Hence, equidimensional algebraic sets can be understood as the Zariski closures of locally closed sets defined by polynomial equations and inequations. Taking this perspective, one additionally enforces these $c$ polynomial equations to form a *triangular set*. These have their origin in the Wu-Ritt characteristic sets (Chou & Gao, 1990; Gallo & Mishra, 1991; Ritt, 1950; Wang, 1993; W.-T. Wu, 1986). Triangularity is therein understood with respect to the variables of the underlying polynomial ring (i.e. in a sense analogous to the notion of triangular matrices in linear algebra). This triangular structure naturally also yields the equations of the algebraic set where the $c$ polynomials fail to define the algebraic set at hand and thus triangular sets have a description of the previously mentioned Zariski open subset attached to them in a natural way. Because of their triangular structure they allow the reduction of certain algorithmic challenges to a univariate problem. Of particular importance, especially in the realm of equidimensional decomposition, are certain special triangular sets called *regular chains*, introduced by Kalkbrener (1993) and Lu and Jingzhong (1994). A regular chain models an *unmixed* dimensional ideal and has good algorithmic properties with respect to the ideal it represents. A related frequently implemented algorithm was also given by Lazard (1991). Algorithms using regular chains are prominently part of the computer algebra system Maple (Chen & Moreno Maza, 2012; Chen et al., 2007). We refer to Hubert (2003) and Wang (2001) for introductions to the subject and to Aubry et al. (1999) for a theoretical account as to how certain different notions of triangular sets relate to each other.

It should be noted that, as for methods based on Gröbner bases combined with algebraic elimination, these triangular encodings make use of polynomials whose degrees, in the worst case, can be as high as is the degree of the equidimensional components they do encode. Nonetheless, algorithms based on triangular representations can be quite well behaved compared to Gröbner basis techniques especially on certain sparse polynomial systems.

Another data structure naturally encoding equidimensional algebraic sets is that of a *geometric resolution* developed by Lecerf (2000, 2003). A geometric resolution is a certain zero-dimensional parametrization of an algebraic set in Noether position. In our setting, these zero-dimensional parametrizations are used to encode generic points in the equidimensional components of the algebraic set under study (the numerical counterpart of this encoding is known as the notion of witness sets (Sommese et al., 2005), a notion that will be utilized in this paper as well). Under certain generically satisfied assumptions on the input, these can be combined with *straight line programs* to obtain the best known complexity bounds for equidimensional decomposition. See also (Jeronimo & Sabia, 2002) for a related approach.

To bypass the “projection-degree” problem, incremental approaches have been investigated in combination with Gröbner bases algorithms. Incremental means here that they feed the decomposition algorithm with one input polynomial after another, in the same way as Lazard (1991) or Lecerf (2000), for example, to identify when some polynomial is a zero divisor in the ring of polynomials quotiented by the ideal generated by the previous polynomials. Moroz (2008) combines Gröbner bases computations with
representations of equidimensional algebraic sets by means of locally closed sets. In a previous work, we
also investigated this approach by exploiting properties of signature-based Gröbner bases algorithms to
enhance the detection and exploitation of zero divisors and compute the so-called nondegenerate locus of
a polynomial system (Eder et al., 2023).

**This Work** In this work, we again take the incremental approach previously mentioned. As in the other incremental algorithms, the foundation of our algorithm is a decomposition algorithm to, given an equidimensional algebraic set X and some $f \in R$, determine the locus of proper intersection of $f$ on X, i.e. the set of points $p \in k^n$ such that $X \cap V(f)$ has dimension one less than X. This is then used to iterate over the input equations $f_1, \dots, f_r$. More precisely, one starts by decomposing $V(f_1, f_2)$ then uses the output to decompose $V(f_1, f_2, f_3)$ and so on.

In contrast to a lot of other algorithms for equidimensional decomposition based on Gröbner bases we
borrow from the theory of triangular sets and work with locally closed sets instead of polynomial ideals
similarly to Moroz (2008). In the iterative approach outlined above this turns out to have two benefits.

First, it naturally removes from the output sets of our iterative algorithm certain embedded components
that appear during the decomposition. To illustrate this consider the following example:

*Example 1.1.* Let $R := \mathbb{Q}[x, y, z]$, $X := V(xy), f := xz$. To decompose $X \cap V(f)$ into equidimensional components one may start by decomposing $X = V(x) \cup V(y)$. Then one intersects these two components with $V(f)$ to obtain the equidimensional decomposition $X \cap V(f) = V(x) \cup V(y, xz)$. The latter set has the irreducible component $V(y, x)$ which is embedded in $V(x)$. If one instead splits into a *disjoint union* $X = V(x) \cup [V(y) \setminus V(x)]$ and again intersects both components with $V(f)$ one obtains $X = V(x) \cup (V(y, z) \setminus V(x))$ and the latter component no longer has the irreducible component $V(y, x)$.

Second, an iterative equidimensional decomposition algorithm may produce redundant components,
which, if they are not deduplicated, may yield an exponential blow-up in the number of components: if
one has decomposed $X = \cup_i X_i$ with the $X_i$ sharing a large number of irreducible components then decom-
posing each $X_i \cap V(f)$ to obtain a decomposition of $X \cap V(f)$ results in an even more redundant decom-
position. Because we use locally closed sets to model our equidimensional sets we are enabled to enforce
that every time we decompose a locally closed set the resulting output sets be pairwise set-theoretically
disjoint. Our experiments indicate that this seems to enforce a sufficiently strong irredundancy between
our components to avoid an exponential blow-up in the number of components.

In this paper we provide two methods to work with the locally closed sets appearing in our algorithm:
One method models them “naively” in the sense that we encode them by storing their defining equations
and inequations and use Gröbner bases of their associated ideals to work with them algorithmically. The
other method tries to avoid having to know a Gröbner basis for the ideal associated to a locally closed set
as much as possible by storing instead a Gröbner basis for a *witness set* of the locally closed set in question.
Using Gröbner bases here with the graded reverse lexicographical ordering has the effect that, compared
to algorithms using triangular sets, we are able to

* avoid computing projections of the algebraic sets to be decomposed and certain frequently made genericity assumptions such as ideals being in Noether position;

• obtain descriptions of these sets with lower degree polynomials.

Borrowing further from the theory of triangular sets we also adopt the heuristic that it is a good
idea to decompose given algebraic sets as often and as finely as possible when working with them. This
philosophy is baked into the recursive structure of our algorithms which exists so as to decompose a given
locally closed set as much as possible given generating sets for certain saturation ideals.

We implemented our algorithm in the computer algebra system OSCAR (The OSCAR team, 2023) using
its interface to the library `msolve` (Berthomieu et al., 2021) for all necessary Gröbner basis computations.
Experimental results indicate that our algorithm is able to tackle polynomial systems which are out of
reach of state-of-the art implementations of algorithms for equidimensional decomposition which are
available in leading computer algebra systems.
# 2 Algorithms

## 2.1 Principles

To illustrate the basic principles behind our equidimensional decomposition algorithm, consider an equidimensional variety X in the affine space $\mathbb{K}^n$. Let $f \in \mathbb{R}$. The variety X is partitioned into:

1. Points $p$ where $f$ is a non zero divisor locally at $p$ (that is in the ring $R_p/I(X)R_p$). The polynomial $f$ takes nonzero values in any open neighborhood of $p$ in X. This defines an open subset $X_{\text{proper}}$ of X.

2. Points $p$ contained in an irreducible component of X on which $f$ vanishes identically. This defines a closed subset $X_{\text{improper}}$ of X.

It is clear that $X = X_{\text{proper}} \cup X_{\text{improper}}$ (where $\cup$ denotes a disjoint union) and that $X_{\text{improper}} \subseteq V(f)$, so that

$$X \cap V(f) = (X_{\text{proper}} \cap V(f)) \cup X_{\text{improper}}.$$

By construction, the $X_{\text{proper}} \cap V(f)$ is a proper intersection: it is equidimensional of dimension $\dim X - 1$, or empty. As a union of irreducible components of X, the closed set $X_{\text{improper}}$ is equidimensional, with the same dimension as the one of X, unless it is empty. So we obtain an equidimensional decomposition of $X \cap V(f)$. Given defining equations for X, this process can be applied iteratively to obtain an equidimensional decomposition of any affine algebraic variety.

In our algorithm we apply the above idea without directly computing $X_{\text{textproper}}$ and $X_{\text{improper}}$. Let $I(X) \subset R$ be an ideal such that $V(I(X)) = X$. Further, we denote by $(I(X) : f^{\infty})$ the saturation ideal of $I(X)$ by $f$. Recall that $V((I(X) : f^{\infty}))$ is the Zariski closure of $X \setminus V(f)$ (Cox et al., 2015, Theorem 4.4.10).

We look for an element $g \in (I(X) : f^{\infty}) \setminus I(X)$. If there is none, this implies that $X_{\text{improper}} = \emptyset$ so $X \cap V(f)$ is equidimensional. If there is such a $g$, then we consider the following partition of X:

1. the closed locus $X_1$ of points $p$ where $g$ has nonzero values in any neighborhood of $p$ in X;

2. the open locus $X_2$ of points $p$ where $g$ is zero in some neighborhood of $p$ in X.

These two sets are equidimensional. By construction, $fg$ vanishes identically on X, so $X_1 \subseteq X_{\text{improper}}$ and this gives the following decomposition of X:

$$X = X_1 \cup (X_2 \cap V(f)). \qquad (1)$$

The ideal of $X_1$ is given by $(I(X) : g^{\infty})$. The term $X_2 \cap V(f)$ may not be equidimensional but we may apply the above idea recursively: We again split $X_2$ along an element in $(I(X_2) : f^{\infty}) \setminus I(X_2)$ if it exists. This leads to Algorithm split.

The set $X_2$ is not closed, this raises the need to deal not only with closed sets of the affine space, but more generally locally closed sets. We do so by partitioning them into special locally closed sets, more precisely into closed sets in the complement of a hypersurface in the affine space, which we call *affine cells*. Concretely, suppose that $I(X_1) = (I(X) : g^{\infty})$ is given by a finite generating set $H \cup \{h\} \subset R$. We then recursively decompose $X_2 = X \setminus V(H \cup \{h\})$ via

$$X_2 = X \setminus V(H \cup \{h\}) = (X \setminus V(h)) \cup ((X \setminus V(H)) \cap V(h)).$$

The intersection with $V(h)$ is computed with split to ensure equidimensionality. Algorithm remove below performs these operations. Findally we obtain an equidimensional decomposition algorithm following an incremental strategy by repeated application of split, see Algorithm equidim.

The primitive operations we use to manipulate affine cells are presented next, while the proof of correctness and termination of the algorithms are in Section 2.3.
## 2.2 Primitives

**Definition 2.1.** An affine cell X is a locally closed set of $\mathbb{K}^n$ of the form $Z \setminus V(g)$ where Z is an algebraic set and $g \in R$. An affine cell X is equidimensional if all the irreducible components of the Zariski closure $\bar{X}$ have the same dimension.

Regardless of the mode of representation of an affine cells, we assume that we can perform the following operations on any affine cell X:

(1) Given $f \in R$, compute the affine cell $X \cap V(f)$;

(2) Given $f \in R$, compute the affine cell $X \setminus V(f)$;

As often in effective algebraic geometry, algebraic sets are defined by ideals that are not always radical so our affine cells come with a distinguished ideal $I(X) \subseteq R$ such that $\bar{X} = V(I(X))$. The radical of $I(X)$ is denoted $\mathrm{rad} I(X)$. We assume that operations (1) and (2) satisfy $I(X) + \langle f \rangle \subseteq I(X \cap V(f))$ and $I(X) \subseteq I(X \setminus V(f))$. We assume further that we can perform the following operations on any affine cell X:

(3) Given $f \in R$, decide if $f \in I(X)$;

(4) Compute a basis of $I(X)$, denoted $\mathrm{basis}(X)$.

For example, we may represent an affine cell X by a pair $(F,g)$, where F is a Gröbner basis of $I(X)$, for some monomial ordering, and g a polynomial such that $X = \bar{X} \setminus V(g)$ (see Becker & Weispfenning, 1993, for an introduction to Gröbner bases). We denote $X = V(F;g)$. For a set $F \subseteq R$ and an element $g \in R$, let $\mathrm{sat}(F,g)$ denote a Gröbner basis of the saturation ideal $(\langle F \rangle : g^\infty)$. Recall that

$$ (I : g^\infty) \stackrel{\mathrm{def}}{=} \{f \in R \mid \exists k \in \mathbb{N}, fg^k \in I\}. $$

Using these primitive $\mathrm{sat}$, we can perform all the four operations above:

(1) $V(F;g) \cap V(f) = V(\mathrm{sat}(F \cup \{f\},g);g)$;

(2) $V(F;g) \setminus V(f) = V(\mathrm{sat}(F,f);fg)$;

(3) $f \in I(X)$ if and only if the normal form of $f$ w.r.t. $F$ is zero;

(4) $\mathrm{basis}(V(F;g)) = F$.

**Remark 2.1.** In Section 3 we explain how to perform the above primitive operations on an affine cell X using a notion called *witness sets*, introduced for the purpose of equidimensional decomposition by Lecerf (2003) under the name lifting fibers. This leads to a lazier representation of X, one where a Gröbner basis for $I(X)$ is not always required.
**Algorithm 1** Equidimensional decompositions

**Input** An equidimensional affine cell $X$, an element $f \in R$

**Output** A partition of $X \cap V(f)$ into equidimensional affine cells

1 **function** split$(X, f)$

2 $G \leftarrow \text{basis}(X \setminus V(f))$

3 **if** $G \subseteq I(X)$ **[can be replaced by $G \subseteq \text{rad } I(X)$]**

4 **return** $\{X \cap V(f)\}$

5 **else**

6 $g \leftarrow$ any element of $G \setminus I(X)$

7 $H \leftarrow \text{basis}(X \setminus V(g))$

8 $\mathcal{D} \leftarrow \{X \cap V(H)\}$

9 **for** $Y \in \text{remove}(X \cap V(g), H)$

10 $\mathcal{D} \leftarrow \mathcal{D} \cup \text{split}(Y, f)$

11 **end**

12 **return** $\mathcal{D}$

13 **end**

14 **end**

**Input** An affine cell $X$, a finite set $H \subset R$

**Precondition** $X \setminus V(H)$ is equidimensional

**Output** A partition of $X \setminus V(H)$ into equidimensional affine cells

1 **function** remove$(X, H)$

2 **if** $H = \emptyset$

3 **return** $\emptyset$

4 **else**

5 $h \leftarrow$ any element of $H$

6 $\mathcal{D} \leftarrow \{X \setminus V(h)\}$

7 **for** $Y \in \text{remove}(X, H \setminus \{h\})$

8 $\mathcal{D} \leftarrow \mathcal{D} \cup \text{split}(Y, h)$

9 **end**

10 **return** $\mathcal{D}$

11 **end**

12 **end**

**Input** a finite set $F \subseteq R$

**Output** A partition of $V(F)$ into equidimensional affine cells

1 **function** equidim$(F)$

2 $\mathcal{D} \leftarrow \{V(\emptyset; 1)\}$ [the full affine space]

3 **for** $f$ **in** $F$

4 $\mathcal{D} \leftarrow \bigcup_{X \in \mathcal{D}} \text{split}(X, f)$

5 **end**

6 **return** $\mathcal{D}$

7 **end**

*Example 2.1.* To illustrate Algorithm *split* we spell out how it behaves on the input $X := V(xy, zw)$ and $f := xz$. Using the notation of Algorithm *split* we find $G = \{y, w\}$. This is not contained in $I(X)$, so we may choose $g := y$ in line 6 of Algorithm *split*. Then we find $H = \{x, zw\}$. Note that $X \setminus V(zw) = \emptyset$ and
so Algorithm split returns

$$ 
\begin{aligned}
& X \cap V(H) \text{ and } \text{split}(\text{remove}(X, H), f) \\
& = V(zw, x) \text{ and } \text{split}(V(y, zw) \setminus V(x), xz).
\end{aligned}
$$

This second call to Algorithm split finds $G = \{y, w\}$, again this set is not contained in $I(V(y, zw) \setminus V(x))$, and so we can choose $g := w$ in line 6. Then we find $H = \{z\}$ which this time yields

$$ 
\begin{aligned}
\text{split}(V(y, zw) \setminus V(x), xz) &= V(y, z) \setminus V(x) \\
&\text{and } \text{split}(V(y, w) \setminus V(xz), xz)
\end{aligned}
$$

The last call to split simply finds the empty set and so all in all we have obtained the decomposition

$$ V(xy, zw, xz) = V(x, zw) \cup V(y, z) \setminus V(x). $$

*Remark 2.2.* Example 2.1 illustrates the fact that Algorithm split may split an algebraic set even if it is equidimensional. Heuristically, the finer the intermediate decomposition in Algorithm equidim is, the computationally easier subsequent steps will be.

## 2.3 Correctness and Termination

When computing an intersection of an equidimensional affine cell $X$ with a hypersurface $V(f)$, we distinguish two cases, depending on whether $V(f)$ intersects $X$ properly or not. Lemma 2.2 deals with the first case, while Lemma 2.3 deals with the second.

**Lemma 2.2.** Let $X$ be an equidimensional affine cell. Let $f \in R$, such that $(I(X) : f^\infty) \subseteq \text{rad } I(X)$. Then $X \cap V(f)$ is empty or equidimensional with dimension $\dim X - 1$.

*Proof.* Let $I = I(X)$. We may assume that $I$ is radical: the assumption $(I : f^\infty) \subseteq \text{rad } I$ also implies $(\text{rad}(I) : f^\infty) \subseteq \text{rad } I$. If $a \in (\text{rad}(I) : f^\infty)$, then $af^r \in \text{rad } I$, for some $r \ge 0$, and so $(af^r)^s \in I$, for some $s \ge 0$. In particular, $a^s \in (I : f^\infty) \subseteq \text{rad } I$. So $a \in \text{rad } I$. Suppose that $X \cap V(f)$ is not empty. By Krull's principal ideal theorem any minimal prime over $I + \langle f \rangle$ has codimension at most $\text{codim } I + 1$. The condition $(I : f^\infty) \subseteq \text{rad } I$ means geometrically that $X \subseteq X \setminus V(f)$, so that $f$ has nonzero values in the neighborhood of any point in $X$. So $f$ is a not a zero divisor in $R/I$. In particular, there is a regular sequence of length $\text{codim } I + 1$ in $I + \langle f \rangle$. Since the polynomial ring $R$ is Cohen-Macaulay it follows that every minimal prime over $I + \langle f \rangle$ has at least codimension $\text{codim } I + 1$. $\square$

**Lemma 2.3.** Let $X$ be an equidimensional affine cell. Let $f \in R$, let $g \in (I(X) : f^\infty)$ and let $I_g = (I(X) : g^\infty)$. Let $X_1 = X \cap V(I_g)$ and $X_2 = (X \cap V(g)) \setminus V(I_g)$. Then:

(i) $X = X_1 \cup X_2$;

(ii) $X \cap V(f) = X_1 \cup (X_2 \cap V(f))$;

(iii) $X_1$ is empty or equidimensional with $\dim X_1 = \dim X$;

(iv) $X_2$ is empty or equidimensional with $\dim X_2 = \dim X$;

*Proof.* Obviously $X = X_1 \cup (X \setminus V(I_g))$. As a set, $X_1$ is the union of the components of $X$ on which $g$ is not identically zero. In particular $X \setminus V(I_g)$ is the set of points of $X$ in a neighborhood of which $g$ is identically zero. Therefore $X \setminus V(I_g) \subseteq V(g)$, so we obtain

$$ X \setminus V(I_g) = (X \cap V(g)) \setminus V(I_g), $$

which gives (i).

Next, we have $I(X_1) = I(X) + I_g = (I(X) : g^\infty)$. Moreover $f \in \text{rad } I(X_1)$. Indeed, $gf^k \in I(X)$ for some $k \ge 0$, by definition of $g$, and therefore $f \in \text{rad}(I(X) : g) \subseteq \text{rad}(I(X) : g^\infty)$. So $X_1 \subseteq V(f)$. It follows that $X \cap V(f) = X_1 \cup (X_2 \cap V(f))$. This proves (ii).

Since $X$ is equidimensional, it follows that $X_1$ (as a union of components of $X$) is also equidimensional of same dimension, unless it is empty. This proves (iii). As for $X_2$, it is open in $X$, so it inherits the equidimensionality and the dimension of $X$, unless it is empty. This proves (iv). $\square$
We now prove correctness and termination of Algorithms split and remove with a mutual induction. On line 3, the test $G \subseteq I(X)$ can be replaced by $G \subseteq \text{rad } I(X)$, or any condition which holds when $G \subseteq I(X)$ and does not hold when $G \not\subseteq \text{rad } I(X)$, this does not affect correctness or termination. We will use this variant in Section 3.

**Theorem 2.4.** For any affine cell X:

(i) If X is equidimensional, then for any $f \in R$, the procedure split terminates on input X and f and outputs a partition of $X \cap V(f)$ into equidimensional affine cells Y with $I(X) \subseteq I(Y)$.

(ii) For any finite set $H \subset R$ such that $X \setminus V(H)$ is equidimensional, the procedure remove terminates on input X and H and outputs a partition of $X \cap V(H)$ into equidimensional affine cells Y with $I(X) \subseteq I(Y)$;

*Proof.* We proceed by Noetherian induction on $I(X)$ and assume the statement holds for any affine cell $X'$ with $I(X) \subseteq I(X')$.

We begin with split. Let $f \in R$ and let $I_f = (I(X) : f^\infty)$. If $I_f \subseteq I(X)$, then Lemma 2.2 applies and $X \cap V(f)$ is equidimensional. So split($X, f$) terminates and is correct in this case.

Assume now that there is some $g \in I_f \setminus I(X)$. Let $I_g = (I : g^\infty)$. Lemma 2.3 applies: an equidimensional decomposition of $X \cap V(f)$ is given by $X \cap V(I_g)$ and an equidimensional decomposition of $((X \cap V(g)) \setminus V(I_g)) \cap V(f)$. Moreover $(X \cap V(g)) \setminus V(I_g)$ is equidimensional. Since $g \notin I(X)$, we have $I(X) \subseteq I(X \cap V(g))$ so remove($X \cap V(g), H$) (using the notations of Algorithm split, where $H$ is a generating set of $I_g$) is correct and terminates, by induction hypothesis. Moreover, it outputs affine cells Y such that $I(X) \subseteq I(X \cap V(g)) \subseteq I(Y)$. So the recursive calls split($Y, f$) are correct and terminate.

As for remove, let $H \subset R$ finite such that $X \setminus V(H)$ is equidimensional. If $H = \emptyset$, then (ii) holds trivially. As for the case $H \neq \emptyset$, let $h \in H$ and $H' = H \setminus h$. Since $V(H) = V(h) \cap V(H')$, we have

$$X \setminus V(H) = (X \setminus V(h)) \cup ((X \setminus V(H')) \cap V(h)). \quad (2)$$

The set $X \setminus V(h)$ and $X \setminus V(H')$ are open in $X \setminus V(H)$ so equidimensional (or empty). By induction on the cardinal of $H$, we assume that remove($X, H'$) is a partition of $X \setminus V(H')$ into equidimensional affine cells, and that every cell Y of this partition satisfies $I(X) \subseteq I(Y)$. By (i), the calls split($Y, h$) terminates and yields a partition of $(X \setminus V(H')) \cap V(h)$ into cells Y with $I(X) \subseteq I(Y)$. Moreover the affine cell $Y = X \setminus V(h)$ also satisfies $I(X) \subseteq I(Y)$. By (2), remove($X, H$) terminates too and is a partition of $X \setminus V(H)$ into cells Y with $I(X) \subseteq I(Y)$. $\square$

**Corollary 2.5.** Algorithm equidim is correct and terminates.

# 3 Implementation and experimental results

## 3.1 Implementation Details

In this section we give some implementation details and alternatives. In particular, we show a lazier data structure for affine cells which is able to delay some Gröbner basis computations at the cost of a Monte Carlo randomization. We have implemented both the method described in Section 2 and the method described in this section.

For either method, we will need an algorithm that, given generators for an ideal $I$ and an element $f \in R$, computes generators for the saturation $(I : p^\infty)$. Even for our lazy representation, this will still sometimes be needed to compute a Gröbner basis for the ideal $I(X)$, where $X$ is an affine cell. In the probabilistic setting, some saturations will be replaced by saturations of zero dimensional ideals.

In our implementation we chose the standard method of performing saturations using Gröbner bases. To compute generators for $(I : p^\infty)$, fix a monomial order $\le$ on $R[t]$ for a new variable $t$ such that $\le$ eliminates $t$. Compute a Gröbner basis $G$ for the ideal $I + \langle tp - 1 \rangle \subset R[t]$ w.r.t. $\le$. Then the elements in $G$ that do not contain the variable $t$ give a Gröbner basis of $(I : p^\infty)$ by the elimination theorem. Other saturation methods also exist such as the methods presented in Eder et al. (2023) or Berthomieu et al. (2022).
Randomization relies on intersecting with random linear subspaces of appropriate dimension to reduce to the zero-dimensional case. This idea is well known in symbolic computation (Lecerf, 2003) and numerical algebraic geometry (Bates et al., 2013, e.g.) wherein these intersections of algebraic sets with random suitable random linear subspaces are known under the name witness sets.

**Proposition 3.1.** Let $X \subseteq \mathbb{K}^n$ be an equidimensional affine cell of dimension $d$ and let $f \in R$. Then, for a generic linear subspace $L \subset \mathbb{K}^n$ of codimension $d$ the following statements hold:

1. $f \in \mathrm{rad} I(X)$ if and only if $f \in \mathrm{rad} I(X \cap L)$.

2. $I(X \setminus V(f)) \subseteq \mathrm{rad} I(X)$ if and only if $X \cap L \cap V(f) = \emptyset$.

*Proof.* We always have $\mathrm{rad} I(X) \subseteq \mathrm{rad} I(X \cap L)$. Conversely, assume that $f \notin \mathrm{rad} I(X)$. Let $U = \{p \in X | f(p) \neq 0\}$. It is an open subset of $X$ and it is non empty, by hypothesis. Since $X$ is equidimensional, $U$ has dimension $d$ and the intersection $U \cap L$ is nonempty (because $L$ is generic). Therefore $f$ is nonzero on a nonempty subset of $X \cap L$. In particular, $f \notin \mathrm{rad} I(X \cap L)$. This proves the first point.

For the second point note that $I(X \setminus V(f)) \subseteq \mathrm{rad} I(X)$ if and only if $X$ and $V(f)$ intersect properly, that is $X \cap V(f)$ is equidimensional of dimension $d-1$. The intersection of $X \cap V(f)$ with the codimension $d$ generic space $L$ is empty if and only if the dimension of $X \cap V(f)$ is less than $d$. The proves the second point. $\square$

In this setting, we represent an equidimensional affine cell $X$ by a triple $(F, G, W, d)$, where $F, G$ and $W$ are subsets of $R$ and $d$ is an integer such that $\dim X = d$, $X = V(F) \setminus V(\prod_{g \in G} G)$ and $W$ (stands for witness set) is a Gröbner basis of $I(X \cap L)$ for some generic linear subspace space $L$ of $\mathbb{K}^n$ of codimension $d$. We denote $X = V(F; G, W, d)$. In practice, $L$ will only be random and sufficient genericity will only hold with high probability (assuming that $\mathbb{K}$ has enough elements). Given only $F, G$ and $d$, we can compute a suitable set $W$ by choosing a set $J \subseteq R$ of $d$ random linear forms and computing a Gröbner basis of $((F : g_1^\infty) : \cdots : g_r^\infty)$, where $G = \{g_1, \cdots, g_r\}$. This procedure is denoted witness$(F, G, d)$.

The four primitive operations are perfomed as follows. For the intersection operation, we need some additional knowledge on the expected dimension of the output. Let $X = V(F; G, W, d)$ be an equidimensional cell.

(1) [Proper intersection] Given $f \in R$ such that $X$ intersects $V(f)$ properly,

$$X \cap V(f) = V(F'; G, \mathrm{witness}(F', G, d-1), d-1),$$

$$\text{with } F' = F \cup \{f\};$$

(1') [Purely improper intersection] Given $H \subset R$ such that $X \cap V(H)$ is a union of components of $X$,

$$X \cap V(H) = V(F \cup H; G, gb(W \cup H), d),$$

where $gb(W \cup H)$ denotes a Gröbner basis of the ideal generated by $X \cup H$;

(2) for $f \in R$, $X \setminus V(f) = V(F; G \cup \{f\}, \mathrm{sat}(H, f), d)$;

(3) $f \in \mathrm{rad} I(X)$ if and only if $1 \in (W : f^\infty)$;

(4) $I(X)$ is computed by saturating $\langle F \rangle$ successively by all the elements of $G$.

In the decomposition algorithm, we always know a priori the kind of each intersection. The intersection on line 4 of split is proper, the intersection on line 8 is purely improper. The one on line 9 is more subtle. Indeed, the decomposition algorithm may produce here a nonequidimensional cell when considering $X \cap V(g)$. With the notations of this algorithm, the cell $X' = X \cap V(g)$ is only equidimensional outside of $V(H)$ (of dimension $\dim X$). This nonequidimensional cell will go through only one operation among the four primitives: $X' \setminus V(h)$ for some $h \in H$. This operation restores equidimensionality. So we can mostly ignore this issue and compute the intersection $X \cap V(g)$ as a purely improper intersection, pretending that $X \cap V(g)$ is equidimensional.
**Algorithm 2 Proper intersection check**

**Input** An equidimensional affine cell $X$, an element $f \in R$

**Output** true if $X \cap V(f)$ is a proper intersection, false otherwise

1 **function** isProper($X, f$)

2 $W \leftarrow$ the withness set of $X$

3 $W' \leftarrow$ a Gröbner basis of $(\langle W \rangle : f^{\infty})$

4 **return** $1 \in W'$

In addition we obtain a fifth operation: a probabilistic algorithm to check if $X \cap V(f)$ is empty or equidimensional of dimension one less than $X$ (or, equivalently $(I(X) : f^{\infty}) \subseteq \text{rad } I(X)$). This is given by Algorithm *isProper*. Equipped with this algorithm, we can replace the if-condition in line 3 of Algorithm *split* with *isProper*($X, f$). Only if this is not satisfied we proceed to compute a Gröbner basis for $I(X \setminus V(f))$.

Lastly we want to note the following: In Algorithm *split*, on input $X$ and $f$, we may have to compute $G := \text{basis}(X \setminus V(f))$ but we use only one element in $G$ in line 6 of Algorithm *split*. This situation can be improved by a simple caching mechanism: Note that in line 10 of Algorithm *split* we call *split*($Y, f$) with affine cells $Y$ satisfying $Y \subset X$. This certainly means $G = \text{basis}(X \setminus V(f)) \subseteq \text{basis}(Y \setminus V(f))$. Hence we may first try to pick an element from the already computed set $G$ in 6 of the call *split*($Y, f$) before computing *basis*($Y \setminus V(f)$).

### 3.1.1 Rationale for the new Data Structure

Always knowing a Gröbner basis for the affine cells appearing in Algorithm *equidim* puts a large penalty on the cost of our algorithms.

This is actually related to well-known observations on the complexity of Gröbner bases under some regularity assumptions. Indeed, for a regular sequence in strong Noether position, the cost of linear algebra steps needed to compute intermediate Gröbner bases in an incremental manner is higher than the final steps (Bardet et al., 2015). Dimension dependent complexity bounds provide another confirmation of this behaviour (Hashemi & Seiler, 2017).

Using witness sets we can potentially avoid a lot of intermediate Gröbner basis computations in our algorithms. In our experience, for a large number of cases, using witness sets greatly improves the efficiency of our algorithm which is theoretically backed up by the previously mentioned complexity results.

Furthermore, in the data structure for affine cells presented in the last subsection we store the defining inequation of our affine cells as a factorization. If one wants to saturate a polynomial ideal $I$ by an element $f \in R$ which is known to have a factorization $f = \prod_{g \in G} g$ given by a finite set $G$ then it is expected to be cheaper to saturate by the elements $g \in G$ one-by-one instead of saturating by $f$ directly using the above elimination method. This lowers the degrees of the polynomials involved.

### 3.1.2 A Better Version of remove

Furthermore, we encountered the following problem when implementing Algorithm *remove* as presented in Section 2. In this algorithm $H$ is a Gröbner basis so it tends to be very redundant (that is very far from being a minimal set of generators). So it often happens that there are two or more elements $h_1, h_2 \in H$ such that for $X_1 := X \setminus V(h_1)$ and $X_2 := X \setminus V(h_2)$ we have $I(X_1) = I(X_2)$. Eventually the sets $X_1$ and $X_2$ become disjoint, since eventually $X_1$ is intersected with $V(h_1)$ or $X_2$ is intersected with $V(h_1)$, but Algorithm *split* may have to split $X_1$ and $X_2$ before that happens. Since splitting an affine cell with Algorithm *split* depends only on the underlying ideals we may then repeat the exact same operations on the level of ideals twice or more. This issue then compounds exponentially due to the recursive nature of our algorithms. We therefore modified Algorithm *remove* to obtain disjoint equidimensional affine cells from $X_1$ and $X_2$ as fast as possible, resulting in Algorithm 3. Note that when we use witness sets this algorithm avoids knowing Gröbner bases for the ideals $I(X_i)$ until potentially line 12.
**Algorithm 3 remove'**

**Input** An affine cell X, a finite set $H \subset R$

**Output** A partition of $X \setminus V(H)$ into equidimensional affine cells

1 **function remove'(X, H)**

2 $D \leftarrow \emptyset$

3 **for** $i$ **from** 1 **to** $r$

4 $X_i \leftarrow X \setminus V(h_i)$

5 $H_i \leftarrow \emptyset$

6 **for** $j$ **from** 1 **to** $i-1$

7 **if** isProper($X_i$, $h_j$)

8 $X_i \leftarrow X_i \cap V(h_j)$

9 **else**

10 $H_i \leftarrow H_i \cup \{h_j\}$

11 **end**

12 **end**

13 $D_i \leftarrow$ decomposition of $X_i \cap V(H_i)$

14 by repeated application of split

15 $D \leftarrow D \cup D_i$

16 **end**

17 **return** $D$

18 **end**

## 3.2 Experimental Results

In this section we give some experimental results. We compare the timings of our implementation of Algorithm *equidim* and methods for equidimensional decomposition of algebraic sets available in various computer algebra systems in a table given below. Some of the timings are discussed in more detail in the next section. We compared to the following implementations:

1. The function *Triangularize* from the *RegularChains* library in Maple (Lemaire et al., 2005), which decomposes a polynomial system into regular chains,

2. The function *equidimensional_decomposition_weak* in *OSCAR* (The OSCAR team, 2023) which is a wrapper around a corresponding *Singular* function (Decker et al., 2021).

3. The *Magma* (Bosma et al., 1997) functions
   *EquidimensionalDecomposition* (corresponding to the column "Magma" in the table) and
   *ProbablePrimeDecomposition* (corresponding to the column "Magma (prime dec.)" in the table)
   and

4. The numerical polynomial systems solver Bertini (Bates et al., 2013) which we ran on each system at hand by requesting a witness set decomposition into irreducible components with fixed precision set to Bertini's default value.

The implementation of our algorithms is itself done in *OSCAR* which is written in the programming language Julia (Bezanson et al., 2017). Its source code is available at

https://github.com/RafaelDavidMohr/Decomp.jl

For all necessary Gröbner basis computations we employ the library *msolve* (Berthomieu et al., 2021) for which *OSCAR* offers an interface. Our suite of example systems is comprised as follows:

1. Cyclic(8), coming from the classical Cyclic(n) benchmark.

2. The systems P4L 1 to 3 come from the perspective-four line problem in robotics, see García Fontán et al. (2022).
3. The systems $C_1$ to $C_3$ are certain jacobian ideals of single multivariate polynomials which define singular hypersurfaces.

4. $Ps(n)$, encoding pseudo-singularities via polynomials

$$f_1, \dots, f_{n-1}, g_1, \dots, g_{n-1}$$

with $f_i \in \mathbb{K}[x_1, \dots, x_{n-2}, z_1, z_2]$, $g_i \in \mathbb{K}[y_1, \dots, y_{n-2}, z_1, z_2]$, the $f_i$ being chosen as a random dense quadrics, and $g_i$ chosen such that $g_i(x_1, \dots, x_{n-2}, z_1, z_2) = f$, i.e. as a copy of $f_i$ in the variables $y_1, \dots, y_{n-2}, z_1, z_2$.

5. $sos(s, n)$, encoding the critical points of the restriction of the projection on the first coordinate to a hypersurface which is a sum of $s$ random dense quadrics in $\mathbb{K}[x_1, \dots, x_n]$.

$$f, \frac{\partial f}{\partial x_2}, \dots, \frac{\partial f}{\partial x_n}, \quad f = \sum_{i=1}^{s} g_i^2.$$

6. $sing(n)$, encoding the critical points of the restriction of the projection on the first coordinate to a (generically singular) hypersurface which is defined by the resultant of two random dense quadrics $A, B$ in $\mathbb{K}[x_1, \dots, x_{n+1}]$:

$$f, \frac{\partial f}{\partial x_2}, \dots, \frac{\partial f}{\partial x_n}, \quad f = \text{resultant}(A, B, x_{n+1}).$$

7. The Steiner polynomial system, coming from Breiding et al. (2020).

8. All remaining examples are part of the BPAS library (Asadi et al., 2021). The BPAS library offers an alternative to the RegularChains library in Maple with special emphasis on paralellism and it will be interesting to compare it to our algorithm in the future.

To obtain the timings in the table below we almost exclusively used the witness set based data structure for affine cells. Every polynomial system was computed with in characteristic 65521 with the exception of Bertini which, as a numerical piece of software, computes over the complex numbers. Due to this difference a comparison between Bertini's and our timings needs to be considered carefully. We tried to indicate this in the table below by coloring the Bertini column in grey. All computations except for Magma were done on an single core of an Intel Xeon Gold 6244 CPU @ 3.60GHz. All Magma computations were done on a single core of an Intel Xeon E5-2690 @ 2.90GHz. We let every algorithm run for at least an hour or 50 times the time it took for the fastest algorithm to complete the system in question, whichever was bigger.

Using the witness sets of our output we also did the following to compare to Bertini: We ran our algorithm in a large random prime charateristic. We then removed the embedded irreducible components from each of our output components and computed the degrees of the output components. This gives us the degree in each dimension of the algebraic set defined by the input. Whenever Bertini reports different degrees, we marked it in the respective column. Due to the randomly chosen large characteristic these degrees should be the same one obtains when considering the algebraic set in question over the complex numbers.

In the second column of this table, we additionally provide the number of affine cells that Algorithm *equidim* decomposed the respective system into. All timings in this table are given in seconds. Due to the way we measured the timings of Bertini we can only report them without any decimal places, rounded up.

## 3.3 Discussion of Experimental Results

We provide here some further information about some of the examples and the behaviour of the different implementations on these examples compared below.

Our algorithm, i.e. Algorithm *equidim*, seems to behave best in comparison to the other implementations when the input system is dense in the sense that each of the input equations of the system in
question involves most, or all, of the variables. This is the case for cyclic 8, the class of the $Ps(\bullet)$ systems, the class of the $Sing(\bullet)$ systems, the class of the $sos(\bullet, \bullet)$ systems and the Steiner polynomial system.

On certain polynomial systems, where each input equation involves only a small subset of the variables, we were able to improve our timings by foregoing the witness set based data structure and instead running a deterministic version of our algorithm akin to the version in Section 2. The improvement we thusly obtained can be explained by the fact that intersecting very sparse systems with random hyperplanes can "destroy their sparsity" and make certain Gröbner basis computations much harder. This was the case for the example Leykin-1: Here running the deterministic version improved our timing to 2.6 seconds.

The Gonnet and dgp6 polynomial systems demonstrate that our algorithm is highly sensitive to the ordering of the input equations: By default we ran our implementation by iterating over the input equations degree by degree in Algorithm *equidim*. With this ordering, our algorithm did not terminate within several hours of computation. When we changed this ordering on these two examples and sorted the input equations instead by length of support, our algorithm terminated in less than one second on these two examples. Our algorithms practical efficiency depends highly on the difficulty of the intermediate polynomial systems which it encounters, these in turn depend on the order of the input equations. On the Gonnet polynomial system, the new ordering resulted in only a subset of the variables being involved in the first few intermediate systems, thus making Gröbner basis for them more tractable. On the dgp6 example the new ordering resulted in an intermediate polynomial system consisting of monomials and binomials which our algorithm decomposes very finely, making the treatment of the remaining equations substantially easier

The system sys2874 can be attacked by both changing the order of the input equations to be ordered by length of support and by using the deterministic version of our algorithm: Doing this, the timing improved by several orders of magnitude to 0.26 seconds.

We also remark that OSCAR's timings improved significantly on the examples sys2449, sys2297 and Leykin-1 (each to less than one second) if one decomposes the radicals of these systems instead of the systems themselves.

For the examples KdV and sys2882 we seem to be bottlenecked by very difficult Gröbner basis computations and less by the inherent structure of our algorithm. Informal experiments where we tried to compute just a Gröbner basis for these systems using msolve suggest that even this is a highly non-trivial computation. For these two systems, techniques involving regular chains seem to be vastly superior over anything that involves Gröbner basis computations.

All in all, these experiments illustrate that on a wide range of examples, our algorithm performs on average better than state-of-the-art implementations and can tackle some problems which were previously unrea

## Acknowledgements

The authors wish to thank Marc Moreno Maza for providing feedback on the BPAS library and its benchmark examples. Additionally the authors thank the two anonymous reviewers for their helpful comments and careful reading of this manuscript.

This work has been supported by the Agence nationale de la recherche (ANR), grant agreements ANR-18-CE33-0011 (SESAME), ANR-19-CE40-0018 (Re Rerum Natura); by the joint ANR-Austrian Science Fund FWF grant agreements ANR-19-CE48-0015 (ECARP) and ANR-22-CE91-0007 (EAGLES); by the EOARD-AFOSR grant agreement FA8665-20-1-7029; by the DFG Sonderforschungsbereich TRR 195; the Forschungsinitiative Rheinland-Pfalz; and by the European Research Council (ERC) under the European Union's Horizon Europe research and innovation programme, grant agreement 101040794 (10000 DIGITS).
<table><thead><tr><th>name</th><th>nb. comp.</th><th>equidim</th><th>Maple</th><th>Oscar</th><th>Magma</th><th>Magma (prime dec.)</th><th>Bertini</th></tr></thead><tbody><tr><td>8-3-config-Li</td><td>23</td><td>1.6</td><td>16 x10</td><td>&gt;1h</td><td>&gt;1h</td><td>65 x40</td><td>4 x2.5</td></tr><tr><td>cyclic8</td><td>6</td><td>381</td><td>&gt;5h</td><td>&gt;5h</td><td>&gt;5h</td><td>&gt;5h</td><td>126 x0.3</td></tr><tr><td>dgp6</td><td>3</td><td>0.2</td><td>53</td><td>2.2</td><td>&gt;1h</td><td>1.2</td><td>75</td></tr><tr><td>Gonnet</td><td>3</td><td>0.2</td><td>2.1</td><td>2.8</td><td>&gt;1h</td><td>1.4</td><td>74</td></tr><tr><td>P4L1</td><td>6</td><td>0.3</td><td>2.4</td><td>1.8</td><td>0.7</td><td>1.5</td><td>21</td></tr><tr><td>P4L3</td><td>8</td><td>0.3</td><td>3.3</td><td>10</td><td>&lt;0.1</td><td>1.5</td><td>11</td></tr><tr><td>KdV</td><td></td><td>&gt;4h</td><td>353</td><td>&gt;4h</td><td>&gt;4h</td><td>7109 x20</td><td>&gt;4h</td></tr><tr><td>Leykin-1</td><td>13</td><td>2.6 x1.9</td><td>4 x3.2</td><td>641 x468</td><td>&gt;1h</td><td>1.4</td><td>6</td></tr><tr><td>C1</td><td>4</td><td>129</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>6</td></tr><tr><td>C2</td><td>4</td><td>0.3</td><td>100</td><td>152</td><td>&gt;1h</td><td>2.0</td><td>6</td></tr><tr><td>C3</td><td>13</td><td>10</td><td>55</td><td>7</td><td>0.3</td><td>1.5</td><td>7</td></tr><tr><td>MontesS16</td><td>6</td><td>1.9 x1.4</td><td>2.7 x1.9</td><td>2.0 x1.4</td><td>1.4</td><td>1.5 x1.1</td><td>9 x5</td></tr><tr><td>Ps(10)</td><td>2</td><td>1.7</td><td>&gt;1h</td><td>30 x17</td><td>&gt;1h</td><td>6 x3.3</td><td>9 x5</td></tr><tr><td>Ps(12)</td><td>2</td><td>51</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>2060 x40</td><td>38 x0.7</td></tr><tr><td>Ps(6)</td><td>2</td><td>&lt;0.1</td><td>0.2</td><td>1.7</td><td>0.5</td><td>0.3</td><td>2.0</td></tr><tr><td>Ps(8)</td><td>2</td><td>&lt;0.1</td><td>4</td><td>1.7</td><td>1.2</td><td>0.8</td><td>6</td></tr><tr><td>Sing(10)</td><td>2</td><td>0.4</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>495</td></tr><tr><td>Sing(4)</td><td>2</td><td>&lt;0.1</td><td>76</td><td>2.2</td><td>8</td><td>5</td><td>5</td></tr><tr><td>Sing(5)</td><td>2</td><td>&lt;0.1</td><td>&gt;1h</td><td>4</td><td>7</td><td>1636</td><td>1.0</td></tr><tr><td>Sing(6)</td><td>2</td><td>&lt;0.1</td><td>&gt;1h</td><td>51</td><td>&gt;1h</td><td>&gt;1h</td><td>8</td></tr><tr><td>Sing(7)</td><td>2</td><td>&lt;0.1</td><td>1704</td><td>399</td><td>&gt;1h</td><td>&gt;1h</td><td>54</td></tr><tr><td>Sing(8)</td><td>2</td><td>0.1</td><td>&gt;1h</td><td>995</td><td>&gt;1h</td><td>&gt;1h</td><td>139</td></tr><tr><td>Sing(9)</td><td>2</td><td>0.2</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>271</td></tr><tr><td>sos(4,2)</td><td>2</td><td>&lt;0.1</td><td>16</td><td>2.2</td><td>1.3</td><td>1.2</td><td>1.0</td></tr><tr><td>sos(4,3)</td><td>2</td><td>&lt;0.1</td><td>694</td><td>2.6</td><td>3.4</td><td>6</td><td>3.0</td></tr><tr><td>sos(5,2)</td><td>2</td><td>&lt;0.1</td><td>&gt;1h</td><td>1.8</td><td>3.7</td><td>1.2</td><td>3.0</td></tr><tr><td>sos(5,3)</td><td>2</td><td>&lt;0.1</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>149</td><td>15</td></tr><tr><td>sos(5,4)</td><td>2</td><td>0.5</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>21</td></tr><tr><td>sos(6,2)</td><td>2</td><td>&lt;0.1</td><td>&gt;1h</td><td>2.0</td><td>5</td><td>1.6</td><td>5</td></tr><tr><td>sos(6,3)</td><td>2</td><td>0.1</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>34</td></tr><tr><td>sos(6,4)</td><td>2</td><td>5</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>69 x14</td></tr><tr><td>sos(6,5)</td><td>2</td><td>14</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>40 x2.9</td></tr><tr><td>steiner</td><td>2</td><td>870</td><td>&gt;12h</td><td>&gt;12h</td><td>&gt;12h</td><td>&gt;12h</td><td>69 x14</td></tr><tr><td>sys2128</td><td>20</td><td>1.1</td><td>9 x9</td><td>5 x4</td><td>&gt;1h</td><td>1.9 x1.8</td><td>829 x790</td></tr><tr><td>sys2161</td><td>33</td><td>8</td><td>29 x3.8</td><td>&gt;1h</td><td>&gt;1h</td><td>8 x1.0</td><td>1196 x159</td></tr><tr><td>sys2297</td><td>11</td><td>0.5</td><td>14</td><td>49</td><td>&gt;1h</td><td>1.7</td><td>497</td></tr><tr><td>sys2353</td><td>13</td><td>1.6 x1.2</td><td>5 x3.8</td><td>2.0 x1.5</td><td>5 x3.7</td><td>1.3</td><td>&gt;1h</td></tr><tr><td>sys2449</td><td>24</td><td>1.3</td><td>28 x21</td><td>60 x46</td><td>&gt;1h</td><td>2.2 x1.7</td><td>1338 x1014</td></tr><tr><td>sys2647</td><td>2</td><td>&lt;0.1</td><td>7</td><td>4</td><td>&gt;1h</td><td>2.0</td><td>9</td></tr><tr><td>sys2874</td><td>5</td><td>0.3</td><td>202</td><td>1.9</td><td>8</td><td>10</td><td>&gt;1h</td></tr><tr><td>sys2880</td><td>50</td><td>4 x2.4</td><td>144 x80</td><td>1.8</td><td>3.4 x1.9</td><td>4 x2.2</td><td>324 x180</td></tr><tr><td>sys2882</td><td></td><td>&gt;1h</td><td>39</td><td>&gt;1h</td><td>&gt;1h</td><td>&gt;1h</td><td>324 x180</td></tr><tr><td>sys2885</td><td>2</td><td>0.3</td><td>6</td><td>2.2</td><td>&gt;1h</td><td>1.2</td><td>76</td></tr><tr><td>sys2945</td><td>5</td><td>0.5</td><td>3.2</td><td>1.8</td><td>0.6</td><td>1.0</td><td>120</td></tr><tr><td>sys2946</td><td>7</td><td>0.2</td><td>0.7</td><td>2.1</td><td>&gt;1h</td><td>1.6</td><td>3.0</td></tr><tr><td>W2</td><td>4</td><td>0.9</td><td>6</td><td>1.9</td><td>7</td><td>1.0</td><td>61</td></tr><tr><td>W44</td><td>3</td><td>0.5</td><td>13</td><td>4</td><td>&gt;1h</td><td>1.5</td><td>66</td></tr><tr><td>Wu-Wang</td><td>3</td><td>3.1</td><td>3.2</td><td>1.8</td><td>0.7</td><td>1.3</td><td>106</td></tr></tbody></table>

Timings are in seconds, except otherwise indicated. The ratio with respect to the best time is given when the latter is over 1 second.

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

We made some minor preparation of the input (like reordering the input equations, or disabling the probabilistic representation of affine cells) to improve the timing.

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

Bertini terminated the computation with an error.

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

The result given by Bertini is not consistent with our result in terms of degree/dimension.
## References

Asadi, M., Brandt, A., Chen, C., Covanov, S., Kazemi, M., Mansouri, F., Mohajerani, D., Moir, R. H. C., Moreno Maza, M., Talaashrafi, D., Wang, L., Xie, N., & Xie, Y. (2021). Basic Polynomial Algebra Subprograms (BPAS) (version 1.791).

Aubry, P., Lazard, D., & Moreno Maza, M. (1999). On the Theories of Triangular Sets. *J. Symb. Comput.*, **28**(1), 105–124.

Aubry, P., Rouillier, F., & Safey El Din, M. (2002). Real solving for positive dimensional systems. *J. Symb. Comput.*, **34**(6), 543–560.

Bardet, M., Faugère, J.-C., & Salvy, B. (2015). On the complexity of the F5 Gröbner basis algorithm. *J. Symb. Comput.*, **70**, 49–70.

Bates, D. J., Hauenstein, J. D., Sommese, A. J., & Wampler, C. W. (2013). *Numerically Solving Polynomial Systems with Bertini*. SIAM.

Becker, T., & Weispfenning, V. (1993). *Gröbner bases* (Vol. 141). Springer-Verlag.

Berthomieu, J., Eder, C., & Safey El Din, M. (2021). Msolve: A library for solving polynomial systems. *Proc. ISSAC 2021*.

Berthomieu, J., Eder, C., & Safey El Din, M. (2022). *New efficient algorithms for computing Gröbner bases of saturation ideals (F4SAT) and colon ideals (Sparse-FGLM-colon)*. arXiv: 2202.13387

Bezanson, J., Edelman, A., Karpinski, S., & Shah, V. B. (2017). Julia: A fresh approach to numerical computing. *SIAM Rev.*, **59**(1), 65–98.

Bosma, W., Cannon, J., & Playoust, C. (1997). The Magma algebra system. I. The user language. *J. Symb. Comput.*, **24**(3-4), 235–265.

Breiding, P., Sturmfels, B., & Timme, S. (2020). 3264 conics in a second. *Not. Am. Math. Soc.*, **67**(1), 30–37.

Caboara, M., Conti, P., & Traverse, C. (1997). Yet another ideal decomposition algorithm. *AAECC 1997*, 39–54.

Chen, C., Corless, R. M., Moreno Maza, M., Yu, P., & Zhang, Y. (2013). An application of regular chain theory to the study of limit cycles. *Int. J. Bifurc. Chaos*, **23**(09), 1350154.

Chen, C., Lemaire, F., Moreno Maza, M., Pan, W., & Xie, Y. (2007). Efficient computations of irredundant triangular decompositions with the regularchains library. *Int. Conf. Comput. Sci.*, 268–271.

Chen, C., & Moreno Maza, M. (2012). Algorithms for computing triangular decomposition of polynomial systems. *J. Symb. Comput.*, **47**(6), 610–642.

Chou, S.-C., & Gao, X.-S. (1990). Ritt-Wu's decomposition algorithm and geometry theorem proving. *10th Int. Conf. Autom. Deduc.*, 207–220.

Cox, D. A., Little, J., & O'Shea, D. (2015). *Ideals, varieties, and algorithms* (4th ed.). Springer.

Decker, W., Greuel, G.-M., & Pfister, G. (1999). Primary Decomposition: Algorithms and Comparisons. *Algorithmic Algebra Number Theory*, 187–220.

Decker, W., Greuel, G.-M., Pfister, G., & Schönemann, H. (2021). *SINGULAR 4-2-1* — A computer algebra system for polynomial computations.

Eder, C., Lairez, P., Mohr, R., & El Din, M. S. (2023). A signature-based algorithm for computing the nondegenerate locus of a polynomial system. *Journal of Symbolic Computation*, 119, 1–21.

Eisenbud, D., Huneke, C., & Vasconcelos, W. (1992). Direct methods for primary decomposition. *Invent. Math.*, **110**(1), 207–235.

Gallo, G., & Mishra, B. (1991). Efficient algorithms and bounds for Wu-Ritt characteristic sets. In T. Mora & C. Traverso (Eds.), *Effective methods in algebraic geometry* (pp. 119–142). Birkhäuser.

García Fontán, J., Nayak, A., Briot, S., & Safey El Din, M. (2022). Singularity analysis for the perspective-four and five-line problems. *Int. J. Comput. Vis.*, **130**, 909–932.

Gianni, P., Trager, B., & Zacharias, G. (1988). Gröbner bases and primary decomposition of polynomial ideals. *J. Symb. Comput.*, **6**(2), 149–167.

Hashemi, A., & Seiler, W. M. (2017). Dimension-dependent upper bounds for gröbner bases. *Proc. ISSAC 2017*, 189–196.

Hubert, E. (2003). Notes on Triangular Sets and Triangulation-Decomposition Algorithms I. *Symb. Numer. Sci. Comput.*, 1–39.

Ishihara, Y. (2022). Modular techniques for intermediate primary decomposition. *Proc. ISSAC 2022*, 479–487.
Jeronimo, G., & Sabia, J. (2002). Effective equidimensional decomposition of affine varieties. *J. Pure Appl. Algebra*, **169**(2), 229-248.

Kalkbrener, M. (1993). A generalized Euclidean algorithm for computing triangular representations of algebraic varieties. *J. Symb. Comput.*, **15**(2), 143–167.

Kalkbrener, M. (1994). Prime decompositions of radicals in polynomial rings. *J. Symb. Comput.*, **18**(4), 365–372.

Krick, T., & Logar, A. (1991). An Algorithm for the Computation of the Radical of an Ideal in the Ring of Polynomials. In *AAECC 1991* (pp. 195–205). Springer.

Lazard, D. (1991). A new method for solving algebraic systems of positive dimension. *Discrete Appl. Math.*, **33**(1), 147–160.

Lecerf, G. (2000). Computing an equidimensional decomposition of an algebraic variety by means of geometric resolutions. *Proc. ISSAC 2000*, 209–216.

Lecerf, G. (2003). Computing the equidimensional decomposition of an algebraic closed set by means of lifting fibers. *J. Complex.*, **19**(4), 564–596.

Lemaire, F., Moreno Maza, M., & Xie, Y. (2005). The RegularChains library in MAPLE. *ACM SIGSAM Bull.*, **39**(3), 96–97.

Lu, Y., & Jingzhong, Z. (1994). Searching dependency between algebraic equations: An algorithm applied to automated reasoning.

Moroz, G. (2008). Regular Decompositions. In D. Kapur (Ed.), *Comput. Math*. (pp. 263–277). Springer.

Noro, M., & Yokoyama, K. (2004). Implementation of prime decomposition of polynomial ideals over small finite fields. *J. Symb. Comput.*, **38**(4), 1227–1246.

Pascual-Escudero, B., Nayak, A., Briot, S., Kermorgant, O., Martinet, P., Safey El Din, M., & Chaumette, F. (2021). Complete singularity analysis for the perspective-four-point problem. *Int. J. Comput. Vis.*, **129**(4), 1217–1237.

Ritt, J. F. (1950). *Differential Algebra*. AMS.

Safey El Din, M., & Schost, É. (2004). Properness defects of projections and computation of at least one point in each connected component of a real algebraic set. *Discrete Comput. Geom.*, **32**(3), 417–430.

Safey El Din, M., Yang, Z.-H., & Zhi, L. (2018). On the complexity of computing real radicals of polynomial systems. *Proc. ISSAC 2018*, 351–358.

Safey El Din, M., Yang, Z.-H., & Zhi, L. (2021). Computing real radicals and S-radicals of polynomial systems. *J. Symb. Comput.*, **102**, 259–278.

Sommese, A. J., Verschelde, J., & Wampler, C. W. (2005). Introduction to numerical algebraic geometry. In M. Bronstein, A. M. Cohen, H. Cohen, D. Eisenbud, B. Sturmfels, A. Dickenstein, & I. Z. Emiris (Eds.), *Solving Polynomial Equations* (pp. 301–337). Springer.

The OSCAR team. (2023). OSCAR – *open source computer algebra research system*.

Wang, D. (1993). An Elimination Method for Polynomial Systems. *J. Symb. Comput.*, **16**(2), 83–114.

Wang, D. (2001). *Elimination Methods*. Springer.

Wu, W., & Gao, X. (2007). Mathematics mechanization and applications after thirty years. *Front. Comput. Sci. China*, **1**(1), 1–8.

Wu, W.-T. (1986). Basic principles of mechanical theorem proving in elementary geometries. 2(3), 221–252.

Yang, L., Hou, X., & Xia, B. (1998). Automated discovering and proving for geometric inequalities. *Int. Workshop Autom. Deduc. Geom.*, 30–46.

Yang, L., Hou, X., & Xia, B. (2001). A complete algorithm for automated discovering of a class of inequality-type theorems. *Sci. China Ser. F Inf. Sci.*, **44**(1), 33–49.

Yokoyama, K. (2002). Prime decomposition of polynomial ideals over finite fields. In *Mathematical software* (pp. 217–227). World Scientific.