Cylindrical Algebraic Decomposition with Equational Constraints

TL;DR

This paper reviews how equational constraints can optimize Cylindrical Algebraic Decomposition (CAD), reducing complexity and computational costs, and introduces new theoretical insights and practical methods for leveraging these constraints.

Contribution

It provides a comprehensive tutorial on CAD with ECs, introduces a new complexity analysis, and demonstrates how ECs can reduce both polynomial count and degree in CAD.

Findings

Double exponential complexity bound reduces with more ECs.

ECs decrease the number of polynomials and their degrees.

Leverages McCallum's reduced projection theory for efficiency.

Abstract

Cylindrical Algebraic Decomposition (CAD) has long been one of the most important algorithms within Symbolic Computation, as a tool to perform quantifier elimination in first order logic over the reals. More recently it is finding prominence in the Satisfiability Checking community as a tool to identify satisfying solutions of problems in nonlinear real arithmetic. The original algorithm produces decompositions according to the signs of polynomials, when what is usually required is a decomposition according to the truth of a formula containing those polynomials. One approach to achieve that coarser (but hopefully cheaper) decomposition is to reduce the polynomials identified in the CAD to reflect a logical structure which reduces the solution space dimension: the presence of Equational Constraints (ECs). This paper may act as a tutorial for the use of CAD with ECs: we describe all necessary background and the current state of the art. In particular, we present recent work on how McCallum's theory of reduced projection may be leveraged to make further savings in the lifting phase: both to the polynomials we lift with and the cells lifted over. We give a new complexity analysis to demonstrate that the double exponent in the worst case complexity bound for CAD reduces in line with the number of ECs. We show that the reduction can apply to both the number of polynomials produced and their degree.

Figures

Tables

Table 1. Table 1: Projection under operator ( 5 ).

Variables	Number	Degree
$n$	$m$	$d$
$n-1$	$2m^2$	$2d^2$
$n-2$	$4m^4$	$8d^4$
⋮	⋮	⋮
$n-r$	$2^{2^{r-1}}{m}^{2^r}$	$2^{2^r-1}{d}^{2^r}$
⋮	⋮	⋮
1	$2^{2^{n-2}}{m}^{2^{n-1}}$	$2^{2^{n-1}-1}{d}^{2^{n-1}}$

Table 2. Table 2: Projection with ( 7 ) ℓ ℓ \ell times and then ( 5 ).

Variables	Number	Degree
$n$	$m$	$d$
$n-1$	$3m$	$2d^2$
⋮	⋮	⋮
$n-\ell$	$3^{\ell }m$	$2^{2^{\ell }-1}{d}^{2^{\ell }}$
$n-\left(\ell +1\right)$	$3^{2\ell }{m}^2$	$2^{2^{\ell +1}-1}{d}^{2^{\ell +1}}$
⋮	⋮	⋮
$n-\left(\ell +r\right)$	$3^{2^r\ell }{m}^{2^r}$	$2^{2^{\ell +r}-1}{d}^{2^{\ell +r}}$
⋮	⋮	⋮
1	$3^{2^{\left(n-1-\ell \right)}\ell }{m}^{2^{n-1-\ell }}$	$2^{2^{n-1}-1}{d}^{2^{n-1}}$

Table 3. Table 3: Maximum degree of projection polynomials produced for CAD when using projection operator ( 7 ) for the first ℓ ℓ \ell projections and then ( 5 ) for the remaining.

Variables	Maximum Degree
	EC	Others
$n$	$d$	$d$
$n-1$	$4d^2$	$d^2$
$n-2$	$27d^3$	$d^4$
$n-3$	$256d^4$	$d^7$
⋮	⋮	⋮
$n-\ell$	${\ell }^{\ell }{d}^{\ell +1}$	$d^{\ell \left(\ell +1\right)\left(1/2\right)+1}$
$n-\left(\ell +1\right)$	$d^{\ell \left(\ell +1\right)+2}$
$n-\left(\ell +2\right)$	$d^{2\ell \left(\ell +1\right)+2^2}$
$n-\left(\ell +3\right)$	$d^{2^2\ell \left(\ell +1\right)+2^3}$
⋮	⋮
$n-\left(\ell +r\right)$	$d^{2^{r-1}\ell \left(\ell +1\right)+2^r}$
⋮	⋮
1	$d^{2^{n-\ell -2}\ell \left(\ell +1\right)+2^{n-\ell -1}}$

Equations

$${\ell }_i(x_1,\dots ,x_{i-1})<$$

$$x_i<u_i(x_1,\dots ,x_{i-1})$$

$$x_i=s_i(x_1,\dots ,x_{i-1})$$

$$\mathrm{cont}(A)$$

$$\mathrm{coeff}(A)$$

$$\mathrm{res}(A)=\left\{\mathrm{res}(f_i,f_j)|f_i\in A,f_j\in A,f_i\ne f_j\right\}.$$

$$P(B)$$

$$P_F(B)$$

$$P_F^{\ast }(B)$$

$$f_1=x+y^2+z,f_2=x-y^2+z,g=x^2+y^2+z^2-1.$$

$$r_1$$

$$R_1$$

$$P_E(A)=\{y,y^4+2x{y}^2+2x^2+y^2-1\}.$$

$$R={\mathrm{res}}_y(y,y^4+2x{y}^2+2x^2+y^2-1)=2x^2-1.$$

$$f_1:=x-y+z^2,f_2:=z^2-u^2+v^2-1,g:=x^2-1,$$

$$f_3:=x+y+z^2,f_4:=z^2+u^2-v^2-1,h:=z.$$

$$\varphi =f_1=0\wedge f_2=0\wedge f_3=0\wedge f_4=0\wedge g\ge 0\wedge h\ge 0.$$

$$\{u=\pm v,x=-1,y=0,z=1\}.$$

$${\mathrm{res}}_z(f_1,f_2)=(v^2-u^2+y-x-1{)}^2$$

$$\displaystyle\mathop{\rm res}\nolimits_{y}\big{(}r_{1},\mathop{\rm res}\nolimits_{z}(f_{1},f_{3})\big{)}=16(u^{2}-v^{2}+x+1)^{4},$$

$$\displaystyle\mathop{\rm res}\nolimits_{y}\big{(}r_{1},\mathop{\rm res}\nolimits_{z}(f_{1},f_{4})\big{)}=4(u^{2}-v^{2})^{2}$$

$$A_5=\{f_1,f_2,f_3,f_4,g,h\}$$

$$A_4:=P_{F_5}(B_5)$$

$$(u^2-v^2-x+y-1{)}^2,4y^2,x-y\}.$$

$$B_4$$

$$A_3:=C\cup P_{F_4}^{\ast }(B_4)=\{x^2-1,u^2-v^2+x+1,u^2-v^2,u^2-v^2+1\}.$$

$$B_3:=\{x^2-1,u^2-v^2+x+1\},C:=\{u^2-v^2,u^2-v^2+1\},$$

$$A_2:=\{u^2-v^2,u^2-v^2+1,u^4-2u^2{v}^2+v^4+2u^2-2v^2\}.$$

$$\{y^2,u^2-v^2+x-y+1,-u^2+v^2+x-y+1,$$

$$u^2-v^2+x+y+1,-u^2+v^2+x+y+1\},$$

$$\{x+1,-u^2+v^2+x+1,u^2-v^2+x+1\},$$

$$\{u^2-v^2\}.$$

$$\prod _{i=1}^n\left[2m_i{d}_i+1\right]=(2md+1)\prod _{r=1}^{n-1}\left[2\left(2^{2^{r-1}}{m}^{2^r}\right)\left(2^{2^r-1}{d}^{2^r}\right)+1\right].$$

$$(2d{)}^{2^n-1}{m}^{2^n-1}{2}^{2^{n-1}-1}.$$

$$\lfloor \frac{1}{2}(3m+1)\rfloor +(m-1)<3m$$

$$\prod _{i=1}^n\left[2m_i{d}_i+1\right]=\prod _{i=1}^{n-(\ell +1)}\left[2m_i{d}_i+1\right]\times \prod _{i=n-\ell }^n\left[2m_i{d}_i+1\right]$$

Full text

The original algorithm produces decompositions according to the signs of polynomials, when what is usually required is a decomposition according to the truth of a formula containing those polynomials. One approach to achieve that coarser (but hopefully cheaper) decomposition is to reduce the polynomials identified in the CAD to reflect a logical structure which reduces the solution space dimension: the presence of Equational Constraints (ECs).

This paper may act as a tutorial for the use of CAD with ECs: we describe all necessary background and the current state of the art. In particular, we present recent work on how McCallum’s theory of reduced projection may be leveraged to make further savings in the lifting phase: both to the polynomials we lift with and the cells lifted over. We give a new complexity analysis to demonstrate that the double exponent in the worst case complexity bound for CAD reduces in line with the number of ECs. We show that the reduction can apply to both the number of polynomials produced and their degree.

1 Introduction

1.1 Cylindrical algebraic decomposition

A Cylindrical Algebraic Decomposition (CAD) splits $\mathbb{R}^{n}$ into cells to maintain an invariance structure relative to an input. Traditionally, the cells are produced to be sign-invariant for a set of input polynomials: meaning that throughout each cell, each of those polynomials has a constant sign. The first CAD algorithm was introduced by Collins [28] to perform Quantifier Elimination (QE) over real closed fields. We describe the necessary details and terminology for CAD in Section 2.1. The invariance property of a CAD means that problems for non-linear polynomial systems such as QE are reduced to testing a finite number of sample points; with the nature of the cells produced allowing for the easy generation of solution descriptions.

However, the use of CAD is often limited by the complexity of computing one. CAD is known to have worst case complexity doubly exponential in the number of variables [33, 19]. Broadly speaking (see Theorem 5 for a precise result) if the input has $m$ polynomials of degree at most $d$ then CAD complexity could be in the order of $(2dm)^{2^{O(n)}}$. For some problems there exist algorithms with better complexity (see the textbook by Basu et al. [7] for example), and there are also many specialised algorithms for restricted inputs; but CAD implementations remain the only general purpose approach for many problems.

This complexity statement is obtained by considering the necessary size of the output for certain examples, and so can only be reduced with changes to the output requirements. Further, unlike some other theoretical results, this one is clearly felt in practice: when increasing dimensionality one will hit the “doubly exponential wall” where progress becomes infeasible. However, extensive investigation into CAD and its sub-algorithms has allowed for the wall to be “pushed back” to the point of allowing many useful computations.

Applications of CAD include: motion planning [71], weight minimisation for truss design [26], epidemic modelling [20], steady state analysis of biological networks [10], economic reasoning [68], artificial intelligence to pass exams [76], parametric optimisation [43], theorem proving [70], derivation of optimal numerical schemes [41], reasoning with multi-valued functions [32], and much more.

1.2 CAD and SC2

CAD has long been important within Symbolic Computation with implementations in multiple computer algebra systems. However, in recent years CAD has been of interest to the separate community of Satisfiability Checking. There, search based algorithms developed for the Boolean SAT problem, make use of heuristics and learning (see the textbook by Biere et al. [8] for details). Success here led to research on domains other than the Booleans, and Satisfiability Module Theory (SMT)-Solvers which use SAT algorithms on the Boolean skeleton of a problem with queries to theory solvers to see if a satisfying Boolean assignment is valid in the domain (learning new clauses if not) [6, 54].

For the SMT domain of non-linear real arithmetic (NRA), CAD and more generally computer algebra systems can play the role of such theory solvers111However, as discussed by Ábrahám et al. [4] a more custom approach is beneficial.. SMT-RAT contains a tailored CAD implementation for use in SMT [57, 53], while Z3 contains an algorithm by Jovanovic and de Moura [52] which uses CAD theory without producing actual CADs. The latter inspired new developments in symbolic computation such as non-uniform CAD by Brown [18]. Further collaboration is informed by the SC2 project: forging interaction between Symbolic Computation and Satisfiability Checking [4].

CADs are produced relative to a problem statement expressed in logic connectives between atoms involving (potentially non-linear) polynomials with integer coefficients. The original CAD algorithm produces decompositions according to the signs of these polynomials, essentially ignoring the logical structure entirely and so producing decompositions fine enough to solve all problems for all logical formulae formed by those polynomials. A Satisfiability Checking approach like that of Jovanovic and de Moura [52] takes an opposite focus, analysing and extending the logical skeleton of the formula until a solution is found with the correct algebraic properties. In order to derive a full solution from a CAD what is truly required is a decomposition on whose cells the truth of the overall logical formula is constant222A sign-invariant decomposition for the polynomials in the formula achieves this, but with far more cells and computation than required..

1.3 Equational constraints

A CAD complexity analysis, such as that in Section 6, does not just conclude the upper bound $(2dm)^{2^{O(n)}}$: it actually shows that in one dimension we must isolate the roots of at most $M$ polynomials of degree $D$, where $D=d^{2^{O(n)}}$ and $M=m^{2^{O(n)}}$. The same orders have been found for lower bounds: for $D$ by Davenport and Heintz [33] and for $M$ by Brown and Davenport [19]. The formulae demonstrating this are not straightforward but the underlying polynomials are surprisingly simple (all bar two linear with each only involving a bounded number of variables, independent of $n$). This demonstrates that the difficulty of CAD resides in the complicated number of ways simple polynomials can interact. Improvement must come from reducing the number of interactions we track.

Definition 1**.**

A Quantifier Free Tarski Formula (QFF) is made up of a finite number of atoms connected by the standard Boolean operators $\land$, $\lor$ and $\neg$. The atoms are statements about the signs of polynomials with integer coefficients: $f\,\sigma\,0$ where $\sigma\in\{=,<,>\}$ (and by combination also $\{\leq,\geq,\neq\}$).

Definition 2**.**

An Equational Constraint (EC) is a polynomial equation logically implied by a QFF. If it is an atom of the formula it is said to be explicit and if not then it is implicit.

Example 1**.**

Let $f$ and $g$ be polynomials. (a) The formula $f=0\land g>0$ has explicit EC $f=0$. (b) The formula $f=0\lor g=0$ has no explicit EC but it does have the implicit EC $fg=0$. (c) The formulae $f^{2}+g^{2}\leq 0$ also has no explicit EC but this one has two implicit ECs $f=0$ and $g=0$.

Collins [29] was the first to suggest that CAD could be simplified in the presence of an EC. He noted that a CAD need only be sign-invariant for the defining polynomial of an EC, and sign-invariant for any others only within those cells where the EC polynomial is zero. He sketched an intuitive approach to produce this by refining the polynomials identified by his CAD algorithm. This approach was formalised and verified by McCallum [63]. A complexity analysis [12, Section 2] showed that making use of a single EC in this way reduces the double exponent of $m$ in the complexity bound for CAD by $1$. Some natural questions arising are: {itemizeshort}

Can savings be made iteratively in the presence of multiple ECs?

Do those savings further reduce the double exponent?

Can corresponding savings be made for the double exponent of $d$?

The first question was answered affirmatively by McCallum [64], although the extension was not trivial333and contained a small mistake as described in Remark 3.. The other questions are answered affirmatively by the present paper. Such questions are of growing importance as CAD finds new application domains with increasing numbers of equations: e.g. in biology by Bradford et al. [10] and England et al. [39]; and in economics by Mulligan et al. [68, 69]. Indeed, many problems that arise in the Satisfiability Checking context contain far more equalities than inequalities (see the NRA benchmarks in the SMT-LIB [5]).

1.4 Contribution and plan

In Section 2 we define the necessary CAD terminology and revise the theory for projection, and reduced projection in the presence of an EC, of McCallum [61, 63, 64]. Then in Section 3 we present recent work on how to leverage this for savings elsewhere in CAD. In Section 4 we propose and verify the corresponding algorithm.

We demonstrate the benefit first with a worked example in Section 5 and then a complexity analysis in Section 6. The latter observes the double exponent in the bound on the number of projection polynomials reducing by the number of ECs. In Section 7 and 8 we explain how a similar reduction can be observed for their degree if we assume CAD input is pre-processed with a Gröbner Basis: a common step in CAD implementations but this is the first theoretical justification for it. Together these show that CAD is doubly exponential in number of variables minus number of ECs.

In Section 9 we examine the main caveat: an assumption of primitivity to the ECs. We demonstrate its presence in the key results of Davenport and Heintz [33], Brown and Davenport [19] proving worst case CAD complexity, and discuss what might be done. We finish by discussing lessons for the Satisfiability Checking community who may call CAD within SMT-solvers in Section 10, and giving a summary in Section 11.

The main contributions in this paper were presented at ISSAC 2015 [37] and CASC 2016 [38]. These conference publications addressed the savings in the number of polynomials and their degrees separately. The present paper unifies the results into a coherent whole providing a single statement of the state of the art. It also expands on some details, such as the need for primitivity in Section 9.

We include all necessary background theory, allowing the paper to act as a tutorial for CAD with ECs, timely given the increased use of CAD outside of computer algebra systems, as part of satisfiability checkers. We further expose the results to the wider SC2 by considering implications for CAD in SMT-solvers in Section 10.

2 Background Material

2.1 CAD computation and terminology

We work under variable ordering $\bm{x}=x_{1}\prec\ldots\prec x_{n}$. The main variable of a polynomial or formula (${\rm mvar}$) is the greatest present under the ordering.

Definition 3**.**

A Cylindrical Algebraic Decomposition (CAD) is a decomposition of $\mathbb{R}^{n}$ into connected cells such that:

each cell is a semi-algebraic set meaning it is defined by a finite sequence of polynomial equations or inequalities; and

• 2.

the cells are arranged cylindrically meaning the projections of any two cells in the decomposition on any lower dimensional space with respect to the ordering are either equal of disjoint.

The latter condition means that each CAD cell is defined by a sequence of conditions: $c_{1}(x_{1}),c_{2}(x_{1},x_{2}),\dots,c_{n}(x_{1},\dots,x_{n})$ where each $c_{i}$ is one of:

[TABLE]

The $\ell_{i},u_{i},s_{i}$ are constants when $i=1$ and otherwise most likely an indexed root expression (a particular root of a polynomial). The former condition tells us that an equivalent semi-algebraic description can be found.

We describe the computation scheme and terminology that the Collins-descended CAD algorithms share. Assume a set of input polynomials (possibly derived from formulae). The first phase of CAD, projection, applies projection operators repeatedly, each time producing another set of polynomials in one less variable (following the variable ordering). Together these are the projection polynomials used in the second phase, lifting.

First $\mathbb{R}$ is decomposed into cells according to the real roots of polynomials univariate in $x_{1}$. Each cell is either a point and therefore its own sample, or an interval inside which we choose a convenient sample (often the dyadic rational with least denominator). We next decompose $\mathbb{R}^{2}$ by repeating a process over each cell in $\mathbb{R}^{1}$. In each case we take the bivariate projection polynomials in $(x_{1},x_{2})$, evaluate them at the sample point of the cell in $\mathbb{R}^{1}$ to give univariate polynomials in $x_{2}$ whose roots we can count and isolate.

The cells identified in $\mathbb{R}^{2}$ fall into two categories. Sections are defined according to the vanishing of a polynomial as in (4), and correspond to the real root of a univariate polynomial. Sectors are usually defined as the regions between two sections as in (2), corresponding to the intervals between real roots of a univariate polynomial. The exceptions are the two infinite sections at either end of the decomposition as in (3) and (1); or if there was no need to decompose at all there may be a single infinite sector with no restrictions on $x_{2}$. In each case the sample point is extended from that of the underlying cell: for sections to include the algebraic number isolated for the real root and for sectors to include any convenient number from the interval (certainly in $\mathbb{Q}$). Together the sections and sectors form a stack over the cell in $\mathbb{R}^{1}$.

Taking the union of these stacks gives the CAD of $\mathbb{R}^{2}$. The process may then be repeated, each time producing a CAD of larger $\mathbb{R}^{i}$, until a CAD of $\mathbb{R}^{n}$ is produced. The subspaces $\mathbb{R}^{i}$ decomposed are those implied by the ordering: $(x_{1})$- space, $(x_{1},x_{2})$-space, $(x_{1},x_{2},x_{3})$-space etc.

Cells are represented at a minimum with a sample point and index. The latter is a list of integers, with the $k$th describing variable $x_{k}$ according to the ordered real roots of the projection polynomials. If the integer is $2i$ the cell is a section corresponds to the $i$th root (counting low to high) and if $2i+1$ it is the adjacent sector444The dimension of a cell is hence easily identified from the index.. The projection operator is chosen so polynomials are delineable in a cell: the portion of their zero set in the cell consists of disjoint sections. A set of polynomials are delineable if each is delineable individually, and the sections of different polynomials are identical or disjoint. If all projection polynomials are delineable then the input polynomials must be sign-invariant: have constant sign in each cell of the CAD.

There have been a great many extensions and optimisations of CAD since its inception, with the survey article by Collins [29] highlighting those from the first 20 years. Since then further highlights have included: symbolic-numeric lifting schemes [73, 50]; local projection approaches [17, 74]; comprehensive Gröbner basis approaches [44] and decompositions via complex space [27, 9].

2.2 Projection operators

A key improvement to CAD has been in the projection operator to reduce the number of projection polynomials computed [47, 61, 63, 64, 14, 11, 46, 12].

The minimal complete projection operator (for sign-invariant CAD) proposed is that of Lazard [55], however the proof of its correctness was shown to be flawed by McCallum and Hong [66]. Shortly before this article went to press a corrected proof was published by McCallum et al. [67] (requiring changes elsewhere in the CAD lifting phase). The theory in the present paper uses the family of projection operators by McCallum [61, 63, 64] but, in due course they may be improved by extending Lazard’s family into the EC theory. We note that the relative savings from the ECs laid out in the present paper would still be maintained if recast into Lazard projection.

Throughout, let $\mathop{\rm cont}\nolimits$, $\mathop{\rm prim}\nolimits$, $\mathop{\rm disc}\nolimits$, $\mathop{\rm ldcf}\nolimits$ and $\mathop{\rm coeff}\nolimits$ denote the content, primitive part, discriminant, leading coefficient, and set of all coefficients of polynomials respectively (in each case taken with respect to a given mvar). When applied to a set of polynomials we interpret these as applying the operation to each polynomial in the set. e.g.

[TABLE]

We let $\mathop{\rm res}\nolimits$ denote the resultant of a pair of polynomials, and for a set

[TABLE]

Recall that a set $B\subset\mathbb{Z}[{\bf x}]$ is an irreducible basis if the elements of $B$ are of positive degree in the mvar, irreducible and pairwise relatively prime. Throughout this section suppose $B$ is an irreducible basis for a set of polynomials, and further that every element of $B$ has mvar $x_{n}$ and that $F\subseteq B$.

We may now define the projection operators introduced respectively by McCallum [61, 63, 64]:

[TABLE]

In the general case with $A$ a set of polynomials and $E\subseteq A$ we proceed with projection by: letting $B$ and $F$ be irreducible bases of the primitive parts of $A$ and $E$ respectively; applying the operators as defined above; and then taking the union of the output with $\mathop{\rm cont}\nolimits(A)$.

Remark 1**.**

We see that (6) is contained in (5) and will usually be smaller. It excludes discriminants and coefficients of $B\setminus F$ which are then reinstated by (7). It also excludes those resultants which involve two polynomials from $B\setminus F$, an exclusion that is maintained by (7). Thus we have $(\ref{eq:ECProj})\subseteq(\ref{eq:ECProjStar})\subseteq(\ref{eq:P})$.

Remark 2**.**

The full set of coefficients of a polynomial is usually unnecessary for CAD [14]. We require knowledge of when the polynomial’s degree drops, and thus the vanishing of the leading coefficient is of most importance. But in the case that it did vanish then the next coefficient becomes leading and thus must also be studied. To guarantee correctness in all cases the full set is taken but most implementations will safely optimise this. E.g. if the leading coefficient is constant then it can never vanish and no further coefficients need to considered.

Remark 3**.**

Operator (7) is different to the $P_{F}^{*}(B)$ of McCallum [64], which excluded the coefficients of $B\setminus F$. The 2001 definition was a mistake, pointed out to us by the anonymous referee of this paper and confirmed in a private communication with McCallum. However, it is not necessary for us to prove a corrected version of the theorems from that paper because, as was noted by McCallum [64] (just after Theorem 2.1) if we allow the additional coefficients then we can assume degree invariance and thus use the existing theorems of McCallum [61] to validate $P_{F}^{*}(B)$.

2.3 Sign invariant CAD Projection

The theorems below use the idea of order-invariance, meaning each polynomial has constant order of vanishing within each cell, which of course implies sign-invariance. We say that a polynomial with mvar $x_{k}$ is nullified over a cell in $\mathbb{R}^{k-1}$ if it vanishes identically throughout.

Theorem 1** (McCallum [61], Thm. 1).**

Let $S$ be a connected submanifold of $\mathbb{R}^{n-1}$ in which each element of $P(B)$ is order-invariant. Then on $S$, each element of $B$ is either nullified or analytic delineable555A variant on delineability defined by McCallum [61].. Further, the sections of elements of $B$ that are not nullified are either identical or pairwise disjoint, and each element of $B$ is order-invariant on such sections.

Suppose we apply $P$ repeatedly to generate projection polynomials. Repeated use of Theorem 1 concludes that a CAD produced by lifting with respect to these projection polynomials is order-invariant so long as no projection polynomial with mvar $x_{k}$ is nullified over a cell in the CAD of $\mathbb{R}^{k-1}$. This condition is known as well-orientedness. It is common for problems to be well oriented and the condition can be easily checked for during lifting666Recall that to lift over a cell we first substitute the cell sample point into the polynomials with main variable one higher: so at this stage we check if any vanish entirely.. In the case that a problem is not well-oriented we cannot conclude sign-invariance (at least over the cell in which nullification occurred). There are some cases where we can rescue the computation (see the report by Brown [16]) but in some cases the only option will be to use an alternative complete projection operator, such as that of Hong [47]. We note that the Lazard operator, recently validated by McCallum et al. [67], removes the well-orientedness condition for sign-invariant CAD777Instead a modified lifting stage checks for nullification and adapts such polynomials to recover the lost information..

2.4 CAD projection for a formula with a single EC

A second theorem allows us to understand how $P_{E}(A)$ is validated.

Theorem 2** (McCallum [63], Thm. 2.2).**

Let $f$ and $g$ be integral polynomials with mvar $x_{n}$ and $r(x_{1},\ldots,x_{n-1})$ be their resultant. Suppose $r\neq 0$. Let $S$ be a connected subset of $\mathbb{R}^{n-1}$ on which $f$ is delineable and $r$ order-invariant. Then $g$ is sign-invariant in every section of $f$ over $S$.

Suppose $A$ was derived from a formula with EC defined by $E=\{f\}$, and that we apply $P_{E}(A)$ once and then $P$ repeatedly to generate projection polynomials. Assuming the input is well-oriented, we can use Theorem 1 to conclude the CAD of $\mathbb{R}^{n-1}$ order invariant for $P_{E}(A)$. The CAD of $\mathbb{R}^{n}$ is then sign-invariant for $E$ using Theorem 1 and sign-invariant for $A$ in the sections of $E$ using Theorem 2. Hence the CAD is truth-invariant for the formula.

2.5 CAD projection with multiple ECs

We now consider the case of multiple ECs. We could of course apply the previous technology by designating one EC for special use and treating the rest as any other constraint (heuristics can help with the choice [13]). But this does not gain any further advantage from the additional ECs, which should be reducing the dimension of our solution space. It is important to note that we cannot simply add multiple EC defining polynomials into $E$ and use $P_{E}(A)$. That would result in a CAD truth-invariant for the disjunction of the ECs, not the conjunction implied by multiple ECs.

Let us first assume that we have two ECs: one whose mvar is that of the system, $x_{n}$, and another whose mvar is $x_{n-1}$. Consider applying first the operator $P_{E}(A)$ where $E$ defines the first EC and then $P_{E^{\prime}}(A^{\prime})$ where $A^{\prime}=P_{E}(A)$ and $E^{\prime}\subseteq A^{\prime}$ contains the second EC. Unfortunately, Theorem 2 does not validate this approach. While it could be applied once for the CAD of $\mathbb{R}^{n-1}$ it cannot then validate the CAD of $\mathbb{R}^{n}$ because the first application of the theorem provided sign-invariance while the second requires the stronger condition of order invariance. The approach is acceptable if $n=3$ (since in two variables the conditions are equivalent for squarefree bases).

Example 2**.**

We consider the formula $\phi=f_{1}=0\land f_{2}=0\land g\geq 0$ with the following polynomials:

[TABLE]

The polynomials are graphed in Figure 1 where $g$ is the sphere, $f_{1}$ the upper surface and $f_{2}$ the lower. We see that $f_{1}$ and $f_{2}$ only meet when $y=0$ and this projection is on the right of Figure 1. It shows that the solution requires $|x|\geq\sqrt{2}/2$ and $z=-x$. How could this be ascertained using CAD?

With variable ordering $z\succ y\succ x$ a sign-invariant CAD for $(f_{1},f_{2},g)$ has 1487 cells using the Qepcad-B program by Brown [15]. We could then test a sample point of each cell to identify the ones where $\phi$ is true. It is preferable to use the presence of ECs. Declaring an EC to Qepcad will ensure it uses the algorithm by McCallum [63] based on a single use of $P_{E}(A)$ followed by $P$. Either choice results in 289 cells. In particular, the solution set is described using 8 cells: all have $y=0,z=-x$ but the $x$-coordinate unnecessarily splits cells at $\frac{1}{2}(1\pm\sqrt{6})$. This is identified due to the projection polynomial $d=\mathop{\rm disc}\nolimits_{y}(\mathop{\rm res}\nolimits_{z}(f_{i},g))$.

For problems with $n>3$ it is still possible to make use of multiple ECs. However, we must include the extra information necessary to provide order-invariance of the non-EC polynomials in the sections of ECs. The following theorem explains that discriminants are required to maintain order invariance, along with degree invariance (and hence coefficients).

Theorem 3** (McCallum [61], Thm. 2).**

Let f be a polynomial with main variable $x_{n}$ with positive degree and $d=\mathop{\rm disc}\nolimits(f)$ which we suppose to be non-zero. Let $S$ be a connected submanifold of $\mathbb{R}^{n-1}$ on which $f$ is degree invariant, and does not vanish identically, and in which $d$ is order invariant. Then $f$ is analytic delineable on $S$ and is order-invariant in each section of $f$ over $S$.

Theorem 3 was used originally as a tool to prove Theorem 1: it gives us the order invariance of polynomials and individual delineability (adding the resultants then extends this to delineability of the set). We can use it again now to show that the strengthening of (6) with the added discriminants and coefficients to form (7) allows for the order invariance conclusion needed for continued application and validation of the operator.

Suppose we have a formula with two ECs, one with mvar $x_{n}$ and the other with mvar $x_{n-1}$. We may now use a reduced operator twice. We first calculate $A^{\prime}=P_{E}(A)$ where $E$ contains the defining polynomial of the first EC, and then $P_{E^{\prime}}^{*}(A^{\prime})$ where $E^{\prime}$ contains the defining polynomial of the other. Subsequent projections simply use $P$. When lifting we use: first Theorem 1 to verify the CAD of $\mathbb{R}^{n-2}$ as order-invariant for $P_{E^{\prime}}^{*}(A^{\prime})$; then Theorem 1 to verify the CAD of $\mathbb{R}^{n-1}$ delineable and order-invariant for $E^{\prime}$, and Theorem 3 to verify it order-invariant and delineable for the individual polynomials of $A^{\prime}$ in the sections of $E^{\prime}$; and finally Theorem 1 and 2 to verify the CAD of $\mathbb{R}^{n}$ order-invariant for $E$ and sign-invariant for $A$ in those cells that are both sections of $E$ and $E^{\prime}$.

2.6 EC propagation

We can now maximise savings in projection when we have ECs in different main variables. Further, if we have two ECs with the same mvar we can usually derive another with a lower mvar by taking the resultant of the two defining polynomials. McCallum [64] defined this as EC propagation. The process requires the two original ECs to be independent, i.e. the satisfaction of one does not imply the satisfaction of the other.

Given additional ECs one can perform multiple rounds of propagation to obtain implicit ECs in a sequence of different main variables. Actually in this case there will be more ECs than we are able to use. For example, given (independent) ECs $f_{i}(x,y,z)=0$ for $i=1,2,3$ in variables $z\prec y\prec x$ then a further three implicit ECs can be found with main variable $y$ and another three with main variable $x$:

[TABLE]

Of course, the latter three will not be independent (the vanishing of one should imply the vanishing of another) but even then there may still be questions of efficiency over which to use. While Bradford et al. [13]h have developed heuristics to help with the choice of which EC to use, there is likely room for improvement888One possibility is the use of machine learning classifiers to make such choices. This is a growing topic within mathematical software, with a recent survey given by England [35]. It has been applied to CAD by Huang et al. [49, 48]..

Example 3**.**

Consider again the example problem from Example 2. We can propagate the two ECs $f_{1}=0$ and $f_{2}=0$ to find implicit EC $r_{1}=0$ as defined above. The resultant of the two defining polynomials is $-2y^{2}$ so we may simplify the EC to $y=0$.

If we declare both ECs in $z$ to Qepcad then it will perform the propagation for us and use reduced projection twice. It will actually apply $P_{E}(A)$ twice (allowed since $n=3$) to produce a CAD with 133 cells. The solution set is now described using only 4 cells (the minimum possible). Note that $d$ (see Example 2) was no longer produced as a projection polynomial.

3 Reductions in the Lifting Stage

The first main contribution of the present paper is to realise that the theorems from the previous section allow for significant savings in the lifting phase (beyond those achieved from reduced projection). To implement these we must discard two embedded principles of CAD:

That the projection polynomials are a fixed set.

• 2.

That the invariance structure of the final CAD can be expressed in terms of sign-invariance of polynomials.

The first was also abandoned by Chen et al. [27], Bradford et al. [9], Jovanovic and de Moura [52], Brown [18], while the second by Brown and McCallum [21], McCallum and Brown [65].

3.1 Minimising the number of polynomials used for lifting

Consider Theorem 2: it allows us to conclude that $g$ is sign-invariant in the sections of $f$ produced over a CAD of $\mathbb{R}^{n-1}$ order-invariant for $P_{\{f\}}(\{f,g\})$. Therefore, it is sufficient to perform the final lift with respect to $f$ only (decompose cylinders according to the roots of $f$ but not $g$). The decomposition imposes sign-invariance for $f$ while Theorem 2 guarantees it for $g$ in the cells where it matters (where those signs could change the truth of the formula).

Example 4**.**

We return to Example 2. Recall that designating either EC and using the algorithm by McCallum [63] produced a CAD with 289 cells. If we follow this approach but lift only with respect to the designated EC at the final step (implemented in the Maple package by England et al. [40]) we obtain a CAD with 141 cells: less than half the output.

This improved lifting follows from the theorems of McCallum [63], but was only noticed 15 years later during the generalisation of that work to the case of multiple formulae by Bradford et al. [11, 12]. Experiments there demonstrated its importance, particularly for problems with many constraints: see Section 8.3 of [12].

When we apply a reduced operator at two levels then we can make such reductions at both the corresponding lifts.

Example 5**.**

We return to the problem from Example 2. Set $A=\{f_{1},f_{2},g\}$ and $E=\{f_{1}\}$. The first projection to eliminate $z$ finds

[TABLE]

These are the resultants of $f_{1}$ with $f_{2}$ and $g$ (after the first is simplified as discussed in Example 3). The discriminant of $f_{1}$ was a constant and so could be discarded, as was its leading coefficient (meaning no further coefficients were required as explained in Remark 2). Set $A^{\prime}=P_{E}(A)$ and $E^{\prime}={y}$ (since $y$ defines an EC for the problem as discussed in Example 3). We have $P_{E^{\prime}}(A^{\prime})=\{R\}$ where

[TABLE]

The other possible entries (the discriminants and coefficients from $E^{\prime}$) are all constants. We hence build a 5 cell CAD of the real line with respect to the two real roots of $R$. We then lift above each cell with respect to $y$ only, in each case splitting the cylinder into three cells about $y=0$, to give a CAD of $\mathbb{R}^{2}$ with 15 cells.

Finally, we lift over each of these 15 cells with respect to $f_{1}$ to give 45 cells of $\mathbb{R}^{3}$. This compares to 133 from Qepcad, which used reduced projection but then lifted with all projection polynomials. No polynomials were nullified, so using Theorems 1 and 2, the output is concluded truth-invariant for $\phi$.

Remark 4**.**

We note that not only is the additional lifting that Qepcad performed unnecessary for the problem at hand, it also provides no further structure on the output. For example, if we had lifted with respect to $f_{2}$ at the final stage in Example 5 then we would be doing so without the knowledge that it is delineable. Hence splitting the cylinder at the sample point offers no guarantee that the cells produced are sign-invariant away from that sample point. So the extra work does not allow us to conclude that $f_{2}$ is sign-invariant (except on sections of $f_{1}$ where we could have concluded it already).

Remark 5**.**

The improvement outlined above not only decreases output size (and computation time) but also the risk of failure999Although we note that if the recent work of McCallum et al. [67] can be extended to ECs then such worries may be unnecessary altogether. from non well-oriented input: we only need worry about nullification of polynomials we lift with.

3.2 Minimising the cells for stack generation

We can achieve more savings by abandoning the aim of producing a CAD sign-invariant with respect to any polynomial, instead insisting only on truth-invariance. We may then lift cells already known to be false trivially to cylinders. The idea of avoiding computations over false cells was presented by Seidl [72], and one could argue that it is the basis of the Partial CAD for QE problems by Collins and Hong [30]. Our contribution here is to explain how such cells can easily be identified in the presence of ECs. We demonstrate with our example.

Example 6**.**

Return to the problem from Examples 2 $-$ 5 and in particular the CAD of $\mathbb{R}^{2}$ produced with 15 cells in Example 5. On 5 of these 15 cells the polynomial $R$ is zero and on the others it is either positive or negative.

Now, $\phi$ can only be satisfied above the 5 cells, as elsewhere the two explicit EC defining polynomials cannot share a root and thus cannot vanish together. We can already conclude the truth value for the 10 cells (false) and thus we do not need to lift over them, except in the trivial sense of extending them to a cylinder in $\mathbb{R}^{3}$. Lifting over the 5 cells where $R=0$ with respect to $f_{1}$ gives 15 cells, which when combined with the 10 cylinders gives a CAD of $\mathbb{R}^{3}$ with 25 cells that is truth-invariant for $\phi$.

Remark 6**.**

The improvements in this subsection are affecting the concluded structure of the output101010In comparison with Remark 4.. The final 25 cell truth-invariant CAD in Example 6 is not sign-invariant for $f_{1}$. The cylinders above the 10 cells in $\mathbb{R}^{2}$ where $R\neq 0$ may have $f_{1}$ varying sign, but since $f_{1}$ can never equal zero at the same time as $f_{2}$ in these cells it does not affect the truth of $\phi$.

Identifying the 5 cells in $\mathbb{R}^{2}$ where $R=0$ was trivial since they are simply the sections of the second lift: those cells with second entry even in the cell index. Similarly, all sections in the third lift are those where $f_{1}$ is zero, however, we cannot conclude that $f_{2}$ is also zero on these as Theorem 2 only guarantees that $f_{2}$ is sign-invariant on them. So we must still finish by evaluating the polynomials at the sample points, but only for the sections.

Reducing the number of cells for stack generation clearly decreases output size, and since the cells can be identified using only an integer parity check computation time decreases (less real root isolation is performed). As described in Remark 5 for the improvements in Section 3.1, this also decreases the risk of non well-oriented input.

4 Algorithm

We present Algorithm 1 (note that it is split into two parts) to build a truth-invariant CAD for a formula through the use of ECs. The input is a quantifier free formula in $x_{1},\dots,x_{n}$ (we assume a fixed variable ordering). The output is a CAD of $\mathbb{R}^{n}$ which we interpret as a set of cells where each cell comes with a cell index and sample point. Each of these is a list of $n$ numbers (recall Section 2.1 for how sample points and indices are represented and extended). This is the minimum information needed, but implementations may choose to store more, such as the formulae for cells.

The first two steps process the input formula into sets of polynomials: $A_{n}$ contains all polynomials in the input formula; while the $E_{k}$ are subsets of $A_{n}$ which each contain the defining polynomial for a primitive EC of main variable $x_{k}$ if one is available, and are empty otherwise. Step 1 is a simple extraction but Step 1 is non-trivial: it must identify suitable ECs for use in projection, and these may not be explicit in the formula (see Definition 2 and Example 1). In practice this will likely require EC propagation (as described in Section 2.6); and EC designation (choosing which of a variety of potential ECs with the same mvar to identify for reduced projection). We discuss the intricacies of this latter step in more detail later.

Steps 1 $-$ 1 run the projection phase of the algorithm. Each projection starts by identifying contents and primitive parts. This is not required for $E_{i}$: since we assume primitive ECs we have $\mathop{\rm cont}\nolimits(E_{i})=\emptyset$, $\mathop{\rm prim}\nolimits(E_{i})=E_{i}$ and the set of factors of $E_{i}$ contained in $B_{i}$ for each $i$.

When there is no declared EC ($E_{i}$ is empty) the projection operator (5) is used (step 1). Otherwise the operator (7) is used (step 1), unless it is the very first or very last projection (step 1) when we use (6). This follows the theory detailed in Section 2. In each case the output of the projection is combined with the contents to form the next layer of projection polynomials.

Steps 1 $-$ 1 construct a CAD for the real line (returning it for univariate input), in what is called the base phase. If there is a declared EC in the smallest variable then the real line is decomposed according to its roots; otherwise according to the roots of all the univariate projection polynomials.

Steps 2 $-$ 2 run the lifting phase, incrementally building CADs of $\mathbb{R}^{k}$ for $k=2,\dots,n$. For each $k$ there are two considerations:

First, whether there is a declared EC with mvar $x_{k}$. If so we lift only with respect to this (step 2) and if not we use all projection polynomials with mvar $x_{k}$ (step 2). See Section 3.1.

• 2.

Second, whether there is a declared EC with mvar $x_{k-1}$. If so we restrict stack generation to those cells where the EC was satisfied. These are simply those with the final entry of the cell index $I_{k-1}$ even (step 2). We lift the other cells trivially to a cylinder in step 2. See Section 3.2.

Algorithm 1 clearly terminates. We will verify that it produces a truth-invariant CAD for the formula if the input is well-oriented, as defined below.

Definition 4**.**

For $k=2,\dots,n$ define sets:

$L_{k}$* $-$ the lifting polynomials: defining polynomial of the declared EC with mvar $x_{k}$ if it exists; else all projection polynomials with mvar $x_{k}$.*

• 2.

$\mathcal{C}_{k}$* $-$ the lifting cells: those cells in the CAD of $\mathbb{R}^{k-1}$ in which the designated EC with mvar $x_{k-1}$ vanishes if it exists, and all cells in that CAD otherwise.*

The input of Algorithm 1 is well-oriented if for $k=2,\dots,n$ no element of $L_{k}$ is nullified over an element of $\mathcal{C}_{k}$.

Theorem 4**.**

Algorithm 1 satisfies its specification.

Proof.

We must show the CAD is truth-invariant for $\phi$, unless the input is not well-oriented when FAIL is returned.

First consider input with $n=1$. The projection phase would not run, with the algorithm jumping to the CAD construction in step 1, returning the output in step 1. If there was no declared EC then the CAD is sign-invariant for all polynomials defining $\phi$ and thus every cell is truth invariant for $\phi$. If there was a declared EC then the output is sign-invariant for its defining polynomial. Cells would either be intervals where the formula must be false; or points, where the EC is satisfied, and the formula either identically true or false depending on the signs of the other polynomials.

Next suppose that the input were not well-oriented (Definition 4). For a fixed $k$, the conditional in steps 2 $-$ 2 sets the lifting polynomials $L_{k}$ to $L$ and the conditional in steps 2 $-$ 2 the lifting cells $\mathcal{C}_{k}$ to $\mathcal{C}_{a}$. Thus it is exactly the conditions of Definition 4 which are checked by step 2, returning FAIL in step 2 when they are not satisfied. Hence if the lifting phase completes then the input is well-oriented.

From now on we suppose $n>1$ and the input is well-oriented. For a fixed $k$ with $2\leq k\leq n$ define admissible cells to be those in the CAD $\mathcal{D}_{k-1}$ of $\mathbb{R}^{k-1}$ produced by Algorithm 1 where all declared ECs with mvar smaller than $x_{k}$ are satisfied, or to be all cells in that CAD if there are no such ECs. Then let $I(k)$ be the following statement in italics.

*Over admissible cells (in $\mathbb{R}^{k-1}$) the CAD $\mathcal{D}_{k}$ of $\mathbb{R}^{k}$ produced by Algorithm 1 is: (a) order-invariant for any EC with mvar $x_{k}$; (b) order- (sign- if $k=n$) invariant for all projection polynomials with mvar $x_{k}$ on sections of the EC over admissible cells, or over all admissible cells if no EC exists. *

We shall prove that, for all k with $1\leq k\leq n$, $I(k)$ is true. We have already proved $I(1)$ (the induction base). Now let $1<k\leq n$ and assume $I(k-1)$ as the induction hypothesis. The truth of $I(k)$, which completes the induction, is then a consequence of the following remarks: {itemizeshort}

When $E_{k}$ is empty we use Theorem 1 to assert all projection polynomials with mvar $x_{k}$ are order-invariant in the stacks over admissible cells giving (a) and (b).

When $E_{k}$ is not empty and $k=2$ we used the projection operator (6). Theorem 2 allows us to conclude (b) and that the EC is sign-invariant in admissible cells. The stronger property of order-invariance follows since the lifting polynomials form a squarefree basis in two variables.

When $E_{k}$ is not empty and $k=n$ we used the projection operator (6). Theorem 2 allows us to conclude (b), but also (a) since in the case $k=n$ the statement requires only sign-invariance.

When $E_{k}$ is not empty and $2<k<n$ we used the projection operator (7). Theorem 3 explains that the additions (7) makes to (6) are sufficient to conclude the statement.

In each case the assumptions of the theorems are met by the inductive hypothesis exactly over admissible cells, according to whether $E_{k-1}$ was empty.

From the definition of admissible cells, we know that $\phi$ is false (and thus trivially truth invariant) upon all cells in the CAD of $\mathbb{R}^{n}$ built over an inadmissible cell of $\mathbb{R}^{k}$, $k<n$. Coupled with the truth of (a) for $k=1,\dots,n$, this implies the CAD of $\mathbb{R}^{n}$ is truth-invariant for the conjunction of ECs (although it may not be truth-invariant for any one individually). The truth of (b) implies that on those cells where all ECs are satisfied, the other polynomials in $\phi$ are sign-invariant and thus $\phi$ is truth-invariant. ∎

5 Worked Example

We consider an example with sufficient variables to show all the features of the algorithm but still small enough to discuss in text. Assume variable ordering $z\succ y\succ x\succ u\succ v$ and define

[TABLE]

We consider the formula

[TABLE]

The solution can be found manually by decomposing the system into blocks. The surfaces $f_{1}$ and $f_{3}$ are graphed in $(x,y,z)$-space on the left of Figure 2. They meet only on the plane $y=0$ and this projection is shown on the right. The surfaces $f_{2}$ and $f_{4}$ are graphed in $(z,u,v)$-space on the left of Figure 3 and meet only when $z=\pm 1$. We consider only $z=+1$ due to $h\geq 0$, with this projection plotted on the right. We thus see that the solution set is

[TABLE]

To ascertain this by Algorithm 1 we must first propagate and designate ECs in Step 1. We choose to use $f_{1}$ first, calculate

[TABLE]

and assign $r_{1}:=v^{2}-u^{2}+y-x-1$. So $r_{1}$ is the defining polynomial for an EC with mvar $y$. Similarly consider

[TABLE]

and assign $r_{2}:=u^{2}-v^{2}+x+1$, $r_{3}:=u^{2}-v^{2}$. These are defining polynomials for ECs with mvar $x$ and $u$ respectively. There is no series of resultants that leads to an EC with mvar $v$ (they all result in constants by that stage). We hence identify $\{E_{j}\}_{k=1}^{n}:=\{f_{1}\},\{r_{1}\},\{r_{2}\},\{r_{3}\},\{\,\}$ in Step 1.

The algorithm continues by extracting the defining polynomials

[TABLE]

and finds $B_{5}=A_{5},F_{5}=E_{5}$ (in fact $F_{i}=E_{i}$ for all $i=1,\dots,5$).

We now start the projection phase. There is a declared EC for the first projection so we use operator (6) to derive

[TABLE]

Hence $C:=\{x^{2}-1\}$ and

[TABLE]

For the next projection we must use operator (7), giving

[TABLE]

For this example the extra discriminants in (7) all evaluated to constants and so could be discarded, while for all polynomials the leading coefficient was constant and so could be discarded with no further coefficients considered (see Remark 2). Then

[TABLE]

and the next projection also uses (7) to produce

[TABLE]

For the final projection there is no EC and so we use operator (5) to find $A_{1}:=\{v^{2}\}.$ The base phase of the algorithm hence produces a 3-cell CAD of the real line isolating [math].

For the first lift we have $L=\{u^{2}-v^{2}\}$ and $C_{a}$ containing all 3 cells. Above the two intervals we split into 5 cells by the curves $u=\pm v$, while above $v=0$ we split into three cells about the origin. From these 13 cells of $\mathbb{R}^{2}$ we select the 5 which were sections of $u^{2}-v^{2}$ for $C_{a}$. These are lifted with respect to $L=\{r_{2}\}$, and the other 8 are simply extended to cylinders in $\mathbb{R}^{3}$. Together this gives a CAD of $\mathbb{R}^{3}$ with 23 cells. The next two lifts are similar, producing first a CAD of $\mathbb{R}^{4}$ with 53 cells and finally a CAD of $\mathbb{R}^{5}$ with 113 cells. The entire calculation takes less than a second in Maple.

5.1 Choice in EC designation

Algorithm 1 could have been initialised with alternative EC designations. There were the 4 explicit ECs with mvar $z$, and by taking repeated resultants we discover the following implicit ECs, in sets with decreasing mvar:

[TABLE]

There are hence $60$ possible permutations of EC designation, but they lead to only 3 different output sizes: 113, 103 and 93 cells. Heuristics for other questions of CAD problem formulation [34, 13, 49, 78] could perhaps be adapted to assist here. None of these are the minimal truth invariant CAD for $\phi$ as all split the CAD of $\mathbb{R}^{1}$ at $v=0$ (from the discriminant $u^{2}-v^{2}$).

5.2 Comparison with previous EC theory

A sign-invariant CAD of $\mathbb{R}^{5}$ for the 6 input polynomials can be produced by Qepcad with 1,118,205 cells. Neither the RegularChains Library in Maple [27] nor our Maple package [40] could do this in under an hour.

Our implementation of the algorithm by McCallum [63], which uses operator (6) once but also performs the final lift with respect to the EC only, can produce a CAD with either 3023, 10,935 or 48,299 (twice) cells depending on which EC is designated. The Qepcad implementation of that algorithm gives 11,961, 30,233, 158,475, or 158,451 cells. Comparing these sets of figures we see the dramatic improvements from just a single reduced lift.

Allowing Qepcad to propagate the four ECs (so a similar projection phase as Algorithm 1 but then a normal CAD lifting phase) produces a CAD with 21,079 cells. By declaring only a subset of the four (which presumably changes the designations of implicit ECs) a CAD with 5,633 cells can be produced, still much more than using Algorithm 1.

The RegularChains Library can also make use of multiple ECs111111but only the latest version from www.regularchains.org, as detailed by Bradford et al. [9], a CAD can be produced instantly. There are choices with analogies to designation [36]), but they all lead to a 137 cell output. In particular, they all have an induced CAD of the real line which splits at $v=\pm 1$ as well as $v=0$.

We note that our Maple implementation is unrefined and unoptimised. We do not claim it as a leading CAD implementation. The purpose of the paper is to illustrate the state of the art in CAD with EC theory, that all CAD implementations should adapt to reproduce. The worked example shows the clear benefits of the improved lifting techniques, which we next generalise with a complexity analysis.

6 Complexity Analysis of CAD with EC

We build on recent work by Bradford et al. [12] to measure the dominant term in bounds on the number of CAD cells produced. Numerous studies have shown this to be closely correlated to the computation time [34, 11, 9]. We assume CAD input with $m$ polynomials of maximum degree $d$ in any one of $n$ variables.

Definition 5**.**

Consider a set of polynomials $p_{j}$. The combined degree of the set is the maximum degree (taken with respect to each variable) of the product of all the polynomials in the set: $\textstyle\max_{i}(\deg_{x_{i}}(\,\prod_{j}p_{j}\,)).$

The set has the (m,d)-property if it may be partitioned into $m$ subsets, each with maximum combined degree $d$.

For example, $\{y^{2}-x,y^{2}+1\}$ has combined degree $4$ and thus the $(1,4)$-property, but also the $(2,2)$-property.

We will measure complexity by keeping track of the number and degree of projection polynomials. Of course, by replacing $\{f,g\}$ with $\{fg\}$ we can reduce the number at the cost of increasing the degree, but since it is much easier to find the roots of $\{f,g\}$ than $\{fg\}$, we do not want to do that. The $(m,d)$-property, introduced in the thesis of McCallum [60], is in some sense the optimal measure of these properties.

Bradford et al. [12] proved that if $A$ has the $(m,d)$-property then $P(A)\cup\mathop{\rm cont}\nolimits(A)$ has the $(M,2d^{2})$-property with $M=\left\lfloor\tfrac{1}{2}(m+1)^{2}\right\rfloor$. When $m>1$, we can bound $M$ by $m^{2}$ (but we need $2m^{2}$ to cover $m=1$).

6.1 Complexity of sign-invariant CAD

If $A$ has the $(m,d)$-property then so does its squarefree basis. Hence applying this result recursively (as in Table 1) measures the growth in $(m,d)$-property during projection under operator (5). After the first projection there are multiple polynomials and so the tighter bound for $M$ is used.

The number of real roots in a set with the $(m,d)$-property is at most $md$ (although in practice many will be in $\mathbb{C}\setminus\mathbb{R}$). The number of cells in the CAD of $\mathbb{R}^{1}$ is thus bounded by twice the product of the final two entries, plus 1. Similarly, if we let $d_{i}$ and $m_{i}$ be the entries in the Number and Degree columns of Table 1 from the row corresponding to $i$ variables, then the total number of cells in the CAD of $\mathbb{R}^{n}$ is bounded by

[TABLE]

Omitting the $+1$s will leave us with the dominant term of the bound, which evaluates to give the following result.

Theorem 5**.**

The dominant term in the bound on the number of CAD cells in $\mathbb{R}^{n}$ produced using (5) is

[TABLE]

6.2 Reduced projection from ECs

From now on assume $\ell$ equational constraints, $0<\ell\leq\min(m,n)$, all with different mvar. For simplicity we assume these variables are $x_{n},\dots,x_{n-\ell+1}$, i.e. the first $\ell$ projections are the reduced ones.

Lemma 6**.**

Suppose $A$ is a set with the $(m,d)$-property and $E\subset A$ has the $(1,d)$-property. Then $\mathop{\rm cont}\nolimits(A)\cup P_{E}^{*}(A)$ has the $(3m,2d^{2})$-property.

Proof.

Bradford et al. [12] proved that applying $P_{E}(A)\cup\mathop{\rm cont}\nolimits(A)$ gives a set of polynomials of size at most $\left\lfloor\tfrac{1}{2}(3m+1)\right\rfloor$ with combined degree $2d^{2}$.

We now have the additional discriminant and coefficients of (7) to take care of. Each polynomial in $A\setminus E$ will generate an additional discriminant of degree at most $d(d-1)$ and $(d+1)$ additional coefficients of degree at most $d$. Multiplying all these polynomials together gives a single polynomial of degree at most $2d^{2}$. There are $m-1$ polynomials in $A\setminus E$ and so in total this projection generates

[TABLE]

polynomials of degree $2d^{2}$. ∎

We apply this recursively in the top part of Table 2, with the bottom derived via the process for $P$, as in Table 1.

Define $d_{i}$ and $m_{i}$ as the entries in the Number and Degree columns of Table 2 from the row corresponding to $i$ variables. We can bound the number of real roots of projection polynomials in $i$ variables by $m_{i}d_{i}$. If we lifted with respect to all projection polynomials the cell count would be bounded by

[TABLE]

Omitting the $+1$ from each product allows us to calculate the dominant term of the bound explicitly as

[TABLE]

6.3 Reduced lifting from ECs

Now we consider the benefit of improved lifting. Start by considering the CAD of $\mathbb{R}^{n-(\ell+1)}$. There can be no reduced lifting until this point and so the cell count bound is given by the second product in (10), which we will denote by $(\dagger)$. The lift to $\mathbb{R}^{n-\ell}$ will involve stack generation over all cells, but only with respect to the EC. This can have at most $d_{n-\ell}$ real roots and so the CAD at most $(2d_{n-\ell}+1)\times(\dagger)$ cells.

The next lift, to $\mathbb{R}^{n-\ell-1}$, will lift the sections with respect to the EC, and the sectors only trivially (to produce the same number of cylinders). Hence the cell count bound is

[TABLE]

with dominant term $2d_{n-(\ell-1)}d_{n-\ell}(\dagger)$. Subsequent lifts follow the same pattern and so $2d_{n}d_{n-1}\dots d_{n-(\ell-1)}d_{n-\ell}(\dagger)$ is the dominant term in the bound for $\mathbb{R}^{n}$. This evaluates to give the following result.

Theorem 7**.**

Consider the CAD of $\mathbb{R}^{n}$ produced using Algorithm 1 in the presence of ECs in the top $\ell$ variables of the ordering. The dominant term in the bound on the number of cells is

[TABLE]

6.4 Summary of complexity analysis

The bound in Theorem 7 is strictly less than the one in Theorem 5. The double exponent of $m$ has decreased by the number of ECs; the result of the improved projection in (11). Then improved lifting has reduced the single exponents in the bound further still in (12).

However, even with this maximal use of ECs, CAD is still doubly exponential in the number of variables due to the first term in (12), the one whose degree is the degree term. This should not be surprising: the theory of ECs is based around reducing the number of polynomials identified in each projection, but not the number of projections which controls the degree growth. Indeed, we can see directly from Tables 1 and 2 that at the end of projection we are dealing with univariate polynomials of degree doubly exponential in $n$ regardless of whether we used ECs or not. Reduced lifting allows us to avoid isolating the real roots of many of these polynomials, but we will always need to consider at least one (the EC defining polynomial). To control degree growth we must show this to be of a lower degree.

7 Controlling Degree Growth

7.1 Degree growth through iterated resultant calculations

The doubly exponential degree comes from the use of iterated resultant calculations during projection: the resultant of two degree $d$ polynomials is the determinant of a $2d\times 2d$ matrix whose entries all have degree at most $d$, and thus a polynomial of degree at most $2d^{2}$. This increase in degree compounded by $(n-1)$ projections gives the first term of the bound (9). Note that the derivation of ECs themselves via EC propagation (Section 2.6) is itself such an iterated resultant calculation. So even though the EC theory of the previous sections allows us to avoid constructing or lifting with many such polynomials, the ECs themselves encode the degree.

The purpose of the resultant in CAD construction is to ensure that the points in lower dimensional space where polynomials vanish together are identified, and thus that the behaviour over a sample point in a lower dimensional cell is indicative of the behaviour over the cell as a whole. The iterated resultant (and discriminant) calculations involved in CAD have been studied previously, for example by McCallum [62] and Lazard and McCallum [56]. We will follow the work of Busé and Mourrain [24] who consider the iterative application of the univariate resultant to multivariate polynomials, demonstrating decompositions into irreducible factors involving the multivariate resultants121212They follow the formalisation of Jouanolou [51] as laid out in [24, §2].. They show that the approach will identify polynomials of higher degree than the true multivariate resultant and thus more than required for the purpose of identifying implicit equational constraints. For example, given 3 polynomials in 3 variables of degree $d$ the true multivariate resultant has degree $\mathcal{O}(d^{3})$ rather than $\mathcal{O}(d^{4})$.

The key result of Busé and Mourrain [24] for our purposes follows. Note that this considers polynomials of a given total degree. However, the CAD complexity analysis discussed above and later is (following previous work on the topic) with regards to polynomials of degree at most $d$ in a given variable. For clarity we use the Fraktur font when discussing total degree and Roman fonts when the maximum degree.

Corollary 8** (Busé and Mourrain [24, Cor. 3.4]).**

Given three polynomials $f_{k}(\emph{{x}},y,z)$ of the form

[TABLE]

where $S$ is any commutative ring, then the iterated univariate resultant

[TABLE]

is of total degree at most $\mathfrak{d}_{1}^{2}\mathfrak{d}_{2}\mathfrak{d}_{3}$ in x, and we may express it in multivariate resultants [51], denoted $\mathop{\rm Res}\nolimits$, as

[TABLE]

Moreover, if the polynomials $f_{1},f_{2},f_{3}$ are sufficiently generic and $n>1$, then this iterated resultant has exactly total degree $\mathfrak{d}^{2}_{1}\mathfrak{d}_{2}\mathfrak{d}_{3}$ in x and both resultants on the right hand side of the above equality are distinct and irreducible.

Remark 7**.**

Although not stated as part of the result by Busé and Mourrain [24], under these generality assumptions, $\mathop{\rm Res}\nolimits_{y,z}(f_{1},f_{2},f_{3})$ has total degree $\mathfrak{d}_{1}\mathfrak{d}_{2}\mathfrak{d}_{3}$ and the second resultant on the right hand side of (13) has total degree $\mathfrak{d}_{1}(\mathfrak{d}_{1}-1)\mathfrak{d}_{2}\mathfrak{d}_{3}$ (see [24, Proposition 3.3] and [62, Theorem 2.6]).

Busé and Mourrain [24] interpret this result in the following quote:131313The quote contains a correction in the description of the second set of roots (removing a dash from $y_{1}$ in the second distinct point). The mistake was identified by the anonymous referees of [38]..

The resultant $r_{12}:=\mathop{\rm res}\nolimits_{z}(f_{1},f_{2})$ defines the projection of the intersection curve between the two surfaces $\{f_{1}=0\}$ and $\{f_{2}=0\}$. Similarly, $r_{13}:=\mathop{\rm res}\nolimits_{z}(f_{1},f_{3})$ defines the projection of the intersection curve between the two surfaces $\{f_{1}=0\}$ and $\{f_{3}=0\}$. Then the roots of $\mathop{\rm res}\nolimits_{y}(r_{12},r_{13})$ can be decomposed into two distinct sets: the set of roots $x_{0}$ such that there exists $y_{0}$ and $z_{0}$ such that

$f_{1}(x_{0},y_{0},z_{0})=f_{2}(x_{0},y_{0},z_{0})=f_{3}(x_{0},y_{0},z_{0}),$

and the set of roots $x_{1}$ such that there exist two distinct points $(x_{1},y_{1},z_{1})$ and $(x_{1},y_{1},z^{\prime}_{1})$ such that

$f_{1}(x_{1},y_{1},z_{1})=f_{2}(x_{1},y_{1},z_{1})\quad\mbox{and}\quad f_{1}(x_{1},y_{1},\allowbreak z^{\prime}_{1})=f_{3}(x_{1},y_{1},z^{\prime}_{1}).$

The first set gives rise to the term $\mathop{\rm Res}\nolimits_{y,z}(f_{1},f_{2},f_{3})$ in the factorization of the iterated resultant $\mathop{\rm res}\nolimits_{y}(\mathop{\rm res}\nolimits_{12},\mathop{\rm res}\nolimits_{13})$, and the second set of roots corresponds to the second factor.

If the $f_{i}$ are all ECs then only the first set are of interest to us as the truth of the formula of interest needs them all to vanish at once. However, for a general CAD construction, the second set of roots may also be necessary as they indicate points where the geometry of the sectors changes.

7.2 How large are these resultants?

Consider three ECs defined by $f_{1},f_{2}$ and $f_{3}$ of degree at most $d$ in each variable separately; and that we wish to eliminate two variables $z=x_{n}$ and $y=x_{n-1}$. We may naïvely set each $\mathfrak{d}_{i}=nd$ to bound the total degree.

The following approach does better. Let $K=S[x_{1},\ldots,x_{n-2},y,z]$ and $L=S[\xi_{1},\ldots,\xi_{N},y,z]$. Only a finite number of monomials in $x_{1},\ldots,x_{n-2}$ occur as coefficients of the powers of $y$, $z$ in $f_{1}$, $f_{2}$ and $f_{3}$. Map each such monomial $x^{\alpha}=\prod_{i=1}^{n-2}x_{i}^{\alpha_{i}}$ to $\widetilde{m_{j}}:=\xi_{j}^{\max\alpha_{i}}$ (using a different $\xi_{j}$ for each monomial141414We could economise: if $x_{1}x_{2}^{2}\mapsto\xi_{1}^{2}$, then we could map $x_{1}^{2}x_{2}^{4}$ to $\xi_{1}^{4}$ rather than a new $\xi_{2}^{4}$. Since this is for the analysis and not in implementation, we ignore such possibilities.) and let $\widetilde{f_{i}}\in L$ be the result of applying this map to the monomials in $f_{i}$. Operation $\widetilde{\ }$ commutes with taking resultants in $y$ and $z$ (but not $x_{i}$).

The total degree in the $\xi_{j}$ of $\widetilde{f_{i}}$ is the same as the maximum degree in all the $x_{1},\ldots,x_{n-2}$ of $f_{i}$, i.e. bounded by $d$, and hence the total degree of the $\widetilde{f_{i}}$ in all variables is bounded by $3d$ ($d$ for the $\xi_{i}$, $d$ for $y$ and $d$ for $z$). If we apply (13) to the $\widetilde{f_{i}}$, we see that \mathop{\rm res}\nolimits_{y}\big{(}\mathop{\rm res}\nolimits_{z}(\widetilde{f_{1}},\widetilde{f_{2}}),\mathop{\rm res}\nolimits_{z}(\widetilde{f_{1}},\widetilde{f_{3}})\big{)} has a factor $\mathop{\rm Res}\nolimits_{y,z}(\widetilde{f_{1}},\widetilde{f_{2}},\widetilde{f_{3}})$ of total degree (in the $\xi_{j}$) $(3d)^{3}$. Hence, by inverting $\widetilde{\ }$, we may conclude $\mathop{\rm Res}\nolimits_{y,z}(f_{1},f_{2},f_{3})$ has maximum degree, in each $x_{i}$, of $(3d)^{3}$.

The results of Jouanolou [51] and Busé and Mourrain [24] apply to any number of eliminations. In particular, if we have eliminated not $2$ but $\ell-1$ variables we will have a polynomial $\mathop{\rm Res}\nolimits_{x_{n-\ell+1}\ldots x_{n}}(f_{n-\ell},\allowbreak\ldots,f_{n})$ of maximum degree $\ell^{\ell}d^{\ell}$ in the remaining variables $x_{1},\ldots,x_{n-\ell}$ as the last implicit EC. Therefore the multivariate resultants we need, $\mathop{\rm Res}\nolimits_{x_{n-\ell+1}\ldots x_{n}}$, only have singly-exponential growth, rather than the doubly-exponential growth of the iterated resultants: can we compute them?

7.3 Gröbner basis instead of iterated resultants

A Gröbner Basis $G$ is a particular generating set of an ideal $I$ (within the ring of polynomials over an algebraically closed field) defined with respect to a monomial ordering. One definition is that the ideal generated by the leading terms of $I$ is generated by the leading terms of $G$. Gröbner Bases (GB) allow properties of the ideal to be deduced such as dimension and number of zeros and so are one of the main practical tools for working with polynomial systems. Their properties and an algorithm to derive a GB for any ideal were introduced in the 1965 PhD thesis of Buchberger [22] (since republished). There has been much research to improve and optimise GB calculation, with the $F_{5}$ algorithm of Faugère [42] perhaps the most used approach currently.

Like CAD the calculation of a GB is necessarily doubly exponential in the worst case [58] (with lexicographic monomial ordering). Recent work by Mayr and Ritscher [59] showed that rather than being doubly exponential with respect to the number of variables present the dependency is in fact on the dimension of the ideal. Despite this bound GB computation can often be done very quickly usually to the point of instantaneous for any problem tractable by CAD, as demonstrated for example by Wilson et al. [77].

A reasonably common CAD technique is to precondition systems with multiple ECs by replacing the ECs by their GB. I.e. let $E=\{e_{1},e_{2},\dots\}$ be a set of polynomials; $G=\{g_{1},g_{2},\dots\}$ a GB for $E$; and $B$ any Boolean combination of constraints, $f_{i}\,\sigma_{i}\,0$, where $\sigma_{i}\in\{<,>,\leq,\geq,\neq,=\}$) and $F=\{f_{1},f_{2},\dots\}$ is another set of polynomials. Then

[TABLE]

are equivalent. A truth-invariant CAD for $\Psi$ is also truth-invariant for $\Phi$.

If we consider GB preconditioning of CAD in the knowledge of the improved projection schemes for ECs then we see an additional benefit. It provides implicit ECs which are not in the main variable of the system removing the need for EC propagation. Since our aim is to produce one EC in each of the last $\ell$ variables, we need to choose an ordering on monomials which is lexicographic with respect to $x_{n}\succ x_{n-1}\succ\cdots\succ x_{n-\ell+1}$: it does not actually matter (in regards to the theory) how we tie-break after that151515Research suggests that ‘total degree reverse lexicographic in the rest’ is most efficient..

Let us suppose that we have $\ell$ ECs $f_{1},\ldots,f_{\ell}$ (at least one of them, say $f_{1}$ must include $x_{n}$, and similarly we can assume $f_{2}$ includes $x_{n-1}$ and so on), such that these imply (even over $\mathbb{C}$) that the last $\ell$ variables are determined (not necessarily uniquely) by the values of $x_{1},\ldots,x_{n-\ell}$. Then the vanishing of polynomials $f_{1}$, $\mathop{\rm Res}\nolimits_{x_{n}}(f_{1},f_{2})$, $\mathop{\rm Res}\nolimits_{x_{n},x_{n-1}}(f_{1},f_{2},f_{3})$ etc. are all implied by the ECs. Hence either they are in the GB, or they are reduced to 0 by the GB, which implies that smaller polynomials are in the GB. Hence our GB will contain polynomials (which are ECs) of degree (in each variable separately) at most

[TABLE]

Note that we are not making, and in the light of the work Mayr and Ritscher [59] cannot make, any similar claim about the polynomials in fewer variables. Also, it is vital that the ECs be in the last variables for our use of the work of Jouanolou [51] and Busé and Mourrain [24] to work. So our results do not directly extend from the case we study, first applying $\ell$ reduced CAD projections in the presence of ECs before reverting to standard projection), to the more general case of having any $\ell$ of the projections be reduced.

7.4 Inclusion in Algorithm 1

There are two routes to include the above suggestion in Algorithm 1.

Directly replace the explicit ECs in $\phi$ by those from the GB as suggested above. This is a pre-processing of the input to Algorithm 1. The identification of ECs in Step 1 involves only a minimal designation choice when there are multiple explicit ECs in $\phi$ with the same mvar. 2. 2.

Encode this process into a sub-algorithm for Step 1. The GB polynomials become additional options for designated ECs along with those from EC propagation and choices are made based on minimal degree or some other criteria [77, 48]. However, in this case, if GB polynomials are designated they must be added to the input set $A_{n}$ (a reinterpreting of Step 1 so it extracts both from $\phi$ and $\{E_{k}\}_{k=1}^{n}$).

The first approach is the one commonly used in implementations (and the one assumed in later discussions). The benefit of the second is that it caps any increase in the number of polynomials from the use of GBs.

It is unlikely that the GB would produce more polynomials in the main variable than explicit ECs (since we are starting with a generating set all in the main variable and deriving another which would mostly not be) but we have yet to rule it out. Of course, the number of polynomials in the input can bear little relation to the number generated by projection. But with the second approach any increase in the initial $m$ is capped to the number of additional designated ECs taken from the GB. The second option may become preferable in the event of development of a good (cheap) heuristic.

8 Evaluating the Use of GBs for ECs

8.1 Worked Example

Let us work with variable ordering $z\succ y\succ x\succ w$; polynomials

[TABLE]

and the following QFF for which we seek a truth-invariant CAD.

[TABLE]

In theory, we could analyse this system with a sign-invariant CAD for the four polynomials however none of the CAD implementations in Maple could do this within 30 minutes. Instead, let us take advantage of the ECs. There are 3 explicit ECs all with mvar $z$ meaning only one can be designated for the first projection. We can propagate to find additional implicit ECs:

[TABLE]

all with main variable $y$. Continuing the propagation with

[TABLE]

gives the three polynomials in the Appendix, each degree 16 in $x$. These are different polynomials161616most easily verified by comparing the final lines of each. but a numerical plot shows them all to have overlapping real part. All possible resultants to eliminate $x$ evaluate to [math].

Step 1 could hence produce $3\times 3\times 3=27$ possible configurations if ECs are identified by propagation. Our implementation could build CADs for only 6 of these configurations171717The common factor of these 6 was the designation of $r_{2}$ for second projection., when using a time limit of 30 minutes. Of the 6 completed there was an average of 423 cells calculated in 113 seconds. The optimal configuration gave 227 cells in 36 seconds using a designation of $f_{2},r_{3}$ and $R_{2}$.

Now consider instead taking a GB of $\{f_{1},f_{2},f_{3}\}$. We use a plex monomial ordering on the same variable ordering as the CAD to achieve a basis:

[TABLE]

This is an alternative generating set for the ideal defined by the explicit ECs and thus all $g_{i}=0$ are ECs for $\phi$. Note that the degrees of the GB polynomials (with respect to any one variable) are on average lower (and never greater) than those of the (corresponding) iterated resultants.

Deriving ECs this way removes the choice for EC with mvar $z$ or $x$ but there are 3 possibilities for the designation with mvar $y$. Designating $g_{2}$ yields 83 cells while either $g_{3}$ or $g_{4}$ result in 55 cells. All 3 configurations took less than 20 seconds to compute (with designating $g_{4}$ the quickest).

8.2 Effect on the complexity bound

We now consider how using a GB to produce the designated ECs will improve the complexity analysis of Section 6. The number of polynomials will be the same as found earlier in Table 2. But we must now track separately the degree of the designated EC and the degree of the main projection polynomials as they are derived differently. For simplicity we will ignore the constant term and focus on the exponents.

As described above, the designated ECs will have degrees $d,4d^{2},27d^{3},\dots$ as tracked in the middle column of Table 3. For the projection polynomials in the top half of the table the reduced projection operator $P_{F}(B)$ will take discriminants and coefficients of the EC polynomial; and resultants of them with the other projection polynomials. Thus the highest degree polynomial produced will have degree that is the sum of the degree of the EC polynomial and the highest degree other polynomial. This generates the right column of Table 3. We see that the degree exponents here form the so called Lazy Caterer’s sequence181818The On-Line Encyclopedia of Integer Sequences (2010), Sequence Number A000124, https://oeis.org/A000124 otherwise known as the Central Polygonal Numbers. The remaining projections recorded in the bottom half of the table use the sign-invariant projection operator and so the degree is squared each time.

Now let us use the the top line of equation (10) derived earlier as the bound when using improved EC projection and lifting before applying the degrees of the projection polynomials. We can substitute here with the degrees from Table 2 as the $d_{i}$. The term with base $d$ may be computed by

[TABLE]

The exponent of $d$ evaluates to

[TABLE]

8.2.1 The ignored constants

Above, we tracked only the degree of the monomial in the bound, and not the constants that multiply it. As well as for simplicity, this was because we could not find a closed form expression for the product of constants generated. However, it is simple to check that the constant factors derived by the GB grow exponentially in $\ell$ while those from iterated resultants grow doubly exponentially. Further, the constant term can be shown to be strictly lower for all but the first few projections. Finally, note that in Section 7.1 we saw that the multivariate resultant was itself a factor of the iterated resultant.

8.2.2 Comparison with base $m$ term

Let us compare the derived exponent (14) with that for the term with base $m$ from (12): $2^{n-\ell}-2$. We see that both show the double exponent of the complexity bound reducing by $\ell$, the number of ECs used. However, the reduction in degree is not quite as clean as the exponential term in the single exponent is multiplied by a quadratic in $\ell$. This is to be expected as the singly exponential dependency on $\ell$ in the Number column of Table 1 was only in the term with constant base while for Table 2 the term with base $d$ is itself single exponential in $\ell$.

8.3 Should one always use GBs?

In Section 8.1 we showed the significant savings available if one derived ECs with GBs and in Section 8.2 we showed this follows through into a theoretical lowering of the worse case complexity bound. The latter offers the first theoretical justification for what is a widely used CAD optimisation.

However, experimental studies by Buchberger and Hong [23], Wilson et al. [77], Huang et al. [48] have shown that it is not always beneficial to pre-process CAD with GB. The most recent experiment by Huang et al. [48] found that 75% of a data set of 1200 randomly generated CAD problems benefited from GB preconditioning. So it is certainly worth giving consideration to how ECs are derived. As noted earlier, the cost of computing the GB itself is usually negligible in comparison to the CAD so it is reasonable to first compute the GB and then decide whether or not to use it. A simple man-made heuristic was presented by Wilson et al. [77] to make the decision while Huang et al. [48] described the training of a machine learning classifier to decide.

There is no contradiction here with the complexity analysis above: the analysis is for the worst case and large input and makes no claim to the average complexity or what happens for smaller input. However, we hypothesise that repeating those studies using the new multiple EC technology would see a reduction in the cases where GB hindered CAD.

9 Caveats and the Need for Primitivity

There are a few caveats to the results presented above. First, Algorithm 1 can fail for non-well oriented input, but as noted earlier, this restriction may be lifted if the new theory for Lazard’s projection operator validated by McCallum et al. [67] can be extended to the EC case. Second, the complexity analysis (both in Section 6 and 8.2) assumes the designated ECs are in strict succession at the start of projection. For the first analysis it was only made to simplify the working, but for the second analysis it was crucial. However, Algorithm 1 itself does not carry this restriction and savings will still clearly be made in this case.

The only substantial restriction in the paper is that the designated ECs191919whether they be explicit in the formula or calculated via iterated results or GBs. be defined by primitive polynomials in the main variable of the projection. The restriction is common in the literature and present in all the underlying theory of McCallum [63, 64].

There are analogies to be made with the well-oriented issue (when a projection factor is nullified). Non-primitive projection factors are not a problem for general CAD because we can factorize prior to projection (order-invariance of factors implies order-invariance of the product). We cannot do the same for a non-primitive EC though, as the next example shows. Also, unlike the well-orientedness issue, the primitivity restriction it is not likely to be removed by developing a Lazard family of EC projection operators.

9.1 Possibilities to use non-primitive ECs?

Example 7**.**

Consider $\phi:=zy=0\land\varphi.$ under ordering $\dots\succ z\succ y\succ\dots$. Polynomial $zy$ is not primitive, so Algorithm 1 cannot use the explicit EC.

We may be tempted to take $E=\{z\}$ as the primitive part, project with operator (6) and include the content $y$ in the first projection. The CAD of $(y,\dots)$-space would be sign-invariant for $y$ and thus the CAD of $(z,y,\dots)$-space truth invariant for the EC (over admissible cells). But we can no longer say only sections are admissible for the next lift as there may be cells with $z\neq 0$ and $y=0$. We must instead lift over all cells of $(y,\dots)$-space, saying: {itemizeshort}

Over sections of $y$: $z$ is no longer an EC (as $zy=0$ is forced by $y=0$), so we lift onto all polynomials.

Over sectors of $y$: $z$ is an EC, so we only lift with respect to this.

In $(z,y,\dots)$-space, some cells are admissible (either $y=0$ or $z=0$) and the rest are not ($zy\neq 0$). There are difficulties in forming a general algorithm:

What would happen if the main variable of the content with respect to $x_{i}$ were not $x_{i-1}$? 2. 2.

What if the content with respect to $x_{i}$ were itself not primitive as a polynomial in $x_{i-1}$? 3. 3.

What if there were another equational constraint in $x_{i-1}$?

The first two can probably (but we have not implemented this yet) be solved by replacing the logic at lines 2– in the algorithm by a dynamic determination of which cells were admissible.

Alternatively in that example we might rewrite $\phi$ as

[TABLE]

so each clause has its own EC. The theory of truth-table invariant CADs (TTICADs) developed by Bradford et al. [11, 12] is designed to deal with such input. More generally, given a formula of the form

[TABLE]

TTICAD allows for an improvement on the standard EC theory. Since $f_{1}f_{2}=0$ is an implicit EC of (16) standard EC theory allows us to avoid studying the $g_{i}$ away from where any $f_{i}$ is zero. By utilising a TTICAD from [11] we can also avoid studying the $g_{i}$ away from where the corresponding $f_{i}$ is zero. This approach was extended in [12] to also consider formulae such as

[TABLE]

where there is no single implicit EC. Although (15) looks closer to (16) it is actually more like (17) since $y$ is not an EC in the main variable.

Although there is the possibility of applying TTICAD for this problem it would first require its own extension to use beyond the first projection (analogous to the present work for standard ECs).

9.2 Classical non primitivity

We can see the importance of the primitivity restriction in the classic complexity results of Brown and Davenport [19], Davenport and Heintz [33]. Both rest on the following construction. Let $P_{k}(x_{k},y_{k})$ be the statement $x_{k}=f(y_{k})$ and then define recursively

[TABLE]

This is $\exists z_{k}\left(z_{k}=f(y_{k-1})\land x_{k-1}=f(z_{k})\right)$, i.e. $x_{k-1}=f(f(y_{k-1}))$. Repeated nesting of this procedure builds the doubly-exponential growth. So

[TABLE]

gives $x_{k-2}=f(f(f(f(y_{k-2}))))$ etc. Rewriting (22) in prenex form gives

[TABLE]

The negation of (23) is therefore

[TABLE]

where the $\overline{\quad}$ operator interchanges $\forall$ and $\exists$. Now, $L_{k}$ can be rewritten as

[TABLE]

and further

[TABLE]

which shows $L_{k}$ to be a conjunction of (non primitive) ECs. This is true for any $L_{i}$, hence the propositional part of (24) is a conjunction of eight ECs, mostly non primitive, and $\neg P_{k}(x_{k},y_{k})$. Hence by induction we have that the whole family of examples $\neg P_{i}$ may be written as complete conjunction of (mostly non primitive) ECs. Furthermore there are equalities whose main variables are the first variables to be projected if we try to produce a quantifier-free form of (24). But that quantifier-free form describes the complement of the semi-algebraic varieties in [19] or [33] (depending which $P_{k}$ we take) and these have doubly-exponential complexity in $n$.

So we observe that the classical results proving the doubly exponential complexity of CAD are not tackled by our EC technology.

10 Lessons for SC2

The SC2 community already appreciates that the logical structure of CAD input is important and should be exploited where ever possible. The main additional lesson from the present paper is that this exploitation can take place not only at the Boolean skeleton level but also in the computer algebra.

10.1 Reasons for optimism

There are high barriers to implementing CAD without the support of a computer algebra system, however, SMT solvers such as SMT-RAT by Loup et al. [57], Kremer and Ábrahám [53] and Z3 by Jovanovic and de Moura [52] show it is possible. Indeed, the developers of SMT-RAT are now beginning to expand their CAD module to include a variety of projection operators [75] and even EC functionality as described by Haehn et al. [45].

One particular barrier is the need multivariate factorization algorithms, which those developing in a computer algebra system can take for granted but represent a significant implementation cost. On this point we highlight to the SMT community the availability of CoCoALib which is a free C++ library that can perform computer algebra computations without the requirement for the accompanying Computer Algebra System (in this case CoCoA) [1]. Further, CoCoA is now actively developing features for use in SMT as described by Abbott and Bigatti [2], Abbott et al. [3].

10.2 Incrementality

A key requirement for the effective use of CAD by SMT-solvers is that the CAD technology be incremental: that polynomials can be added and removed to the input with the data structures of the CAD edited rather than recalculated. Such incremental CAD algorithms are now under development as part of the SC2 by Kremer and Ábrahám [53], Cowen-Rivers and England [31].

An additional advantage from incremental CAD would be with regards to the issue of well-orientedness. I.e. if a particular operator is found to not be well-oriented at the end of a CAD calculation the next step would be to revert to a less efficient operator which is a superset of the original. Refining an existing decomposition should be cheaper than recomputing from scratch. Although on this point, the development of the Lazard projection theory may remove the well-orientedness condition all together.

However, the use of CAD with ECs incrementally requires additional development work. First, it introduces additional decisions to be taken such as EC designation and whether to pre-processing with GB (not to mention whether that can also be done incrementally). Second, this growing number of decisions needs to be taken in tandem, prompting exponential growth in the number of possibilities that overwhelms existing heuristics. Machine learning techniques may be one way forward, as outlined by England [35].

Finally, existing heuristics that guide the Boolean search may not be suitable since the use of ECs could prompt what appears as strange behaviour in the SMT context. For example, removing a constraint that was equational could actually grow the output CAD since it necessitates the use of a larger projection operator. Correspondingly, adding an equational constraint could allow a smaller operator and shrink the output. SMT solver search heuristics will need to be adapted to handle these possibilities.

11 Summary

We have presented much of the state of the art in the theory of CAD with Equational Constraints. This included how ECs may be leveraged for savings in the lifting phase as well as projection. We demonstrated the benefits of the theory with worked examples and complexity analysis. The latter shows that the worst CAD bound has double exponent that reduced from $n$ by the number of ECs. Crucially, this is the global double exponent covering both the number and degree of polynomials, if we allow for Groebner Basis pre-processing.

The main avenues for future work are an exploration of dealing with non-primitive ECs; the extension of the Lazard projection operator to a family of operators for ECs; the development of heuristics for choosing which ECs to designate; and the development of incremental EC technology. We note that the current results and any future progress have benefits not only for Symbolic Computation but the wider SC2 community.

Acknowledgements

This work was originally supported by EPSRC grant EP/J003247/1 and later by EU H2020-FETOPEN-2016-2017-CSA project $\mathcal{SC}^{2}$ (712689).

We are grateful to all the anonymous referees of this paper, and also those of our conference papers at ISSAC 2015 [37] and CASC 2016 [38], for many helpful comments. We are particularly grateful to the referee who pointed out the mistake in the literature on operator (7) and Dr McCallum for discussing this with us.

We also thank Prof. Buchberger for reminding JHD that Gröbner bases were applicable to the problem of degree growth.

Appendix A The Iterated Resultants From Section 8.1

[TABLE]

[TABLE]