Deciding the Consistency of Non-Linear Real Arithmetic Constraints with a Conflict Driven Search Using Cylindrical Algebraic Coverings

We present a new algorithm for determining the satisfiability of conjunctions of non-linear polynomial constraints over the reals, which can be used as a theory solver for satisfiability modulo theory (SMT) solving for non-linear real arithmetic. The algorithm is a variant of Cylindrical Algebraic Decomposition (CAD) adapted for satisfiability, where solution candidates (sample points) are constructed incrementally, either until a satisfying sample is found or sufficient samples have been sampled to conclude unsatisfiability. The choice of samples is guided by the input constraints and previous conflicts. The key idea behind our new approach is to start with a partial sample; demonstrate that it cannot be extended to a full sample; and from the reasons for that rule out a larger space around the partial sample, which build up incrementally into a cylindrical algebraic covering of the space. There are similarities with the incremental variant of CAD, the NLSAT method of Jovanović and de Moura, and the NuCAD algorithm of Brown; but we present worked examples and experimental results on a preliminary implementation to demonstrate the differences to these, and the benefits of the new approach.


Real Algebra
Formulae in Real Algebra are Boolean combinations of polynomial constraints with rational coefficients over real-valued variables, possibly quantified. Real algebra is a powerful logic suitable to express a wide variety of problems throughout science and engineering. The 2006 survey [1] gives an overview of the scope here. Recent new applications include bio-chemical network analysis [2], economic reasoning [3], and artificial intelligence [4]. Having methods to analyse real algebraic formulae allows us to better understand those problems, for example, to find one solution, or symbolically describe all possible solutions for them.
In this paper we restrict ourselves to formulae in which every variable is existentially quantified. This falls into the field of Satisfiability Modulo Theories (SMT ) solving, which grew from the study of the Boolean SAT problem to encompass other domains for logical formulae. In traditional SMT solving, the search for a solution follows two parallel threads: a Boolean search tries to satisfy the Boolean structure of ϕ, accompanied by a theory search that tries to satisfy the polynomial constraints that are assumed to be True in the current Boolean search. To implement such a technique, we need a decision procedure for the theory search that is able to check the satisfiability of conjunctions of polynomial constraints, in other words, the consistency of polynomial constraint sets. The development of such methods has been highly challenging.
Tarski showed that the Quantifier Elimination (QE) problem is decidable for real algebra [5]. That means, for each real-algebraic formula ∀x.ϕ or ∃x.ϕ it is possible to construct another, semantically equivalent formula using the same variables as used to express ϕ but without referring to x. For conjunctions of polynomial constraints it means that it is possible to decide their satisfiability, and for any satisfiable instances provide satisfying variable values. Tarski's results were ground breaking, but his constructive solution was non-elementary (with a time complexity greater than all finite towers of powers of 2) and thus not applicable.

Cylindrical Algebraic Decomposition
An alternative solution was proposed by Collins in 1975 [6]. Since its invention, Collins' Cylindrical Algebraic Decomposition (CAD) method was the target of numerous improvements and has been implemented in many Computer Algebra Systems. Its theoretical complexity is doubly exponential 1 . The fragment of our interest, which excludes quantifier alternation, has lower theoretical complexity of singly exponential time (see for example [11]), however, currently no algorithms are implemented that realise this bound in general (see [12] for an analysis as to why). Current alternatives to the CAD method are efficient but incomplete methods using, e.g., linearisation [13,14], interval constraint propagation [15,16], virtual substitution [17,18], subtropical satisfiability [19] and Gröbner bases [20].
Given a formula in real algebra, a CAD may be produced which decomposes real space into a finite number of disjoint cells so that the truth of the formula is invariant within each cell. Collin's CAD achieved this via a decomposition on which each polynomial in the formula has constant sign. Querying one sample point from each cell can then allow us to determine satisfiability, or perform quantifier elimination.
A full decomposition is often not required and so savings can be made by adapting the algorithm to terminate early. For example, once a single satisfying sample point is found we can conclude satisfiability of the formulae and finish 2 . This was first proposed as part of the partial CAD method for QE [22]. The natural implementation of CAD for SMT performs the decomposition incrementally by polynomial, refining CAD cells by subdivision and querying a sample from each new unsampled cell before performing the next subdivision.
Even if we want to prove unsatisfiability, the decomposition of the state space into sign-invariant cells is not necessary. There has been a long development in how simplifications can be made to achieve truthinvariance without going as far as sign-invariance such as [23], [24], [25].

New Real Algebra Methods Inspired by CAD
The present paper takes this idea of reducing the work performed by CAD further, proposing that the decomposition need not even be disjoint, and neither sign-nor truth-invariant as a whole for the set of constraints. We will produce cells incrementally, each time starting with a sample point, which if unsatisfying we generalise to a larger cylindrical cell using CAD technology. We continue, selecting new samples from outside the existing cells until we find either a satisfying sample, or the entire space is covered by our collection of overlapping cells which we call a Cylindrical Algebraic Covering.
Our method shares and combines ideas from two other recent CAD inspired methods: (1) the NLSAT approach by Jovanović and de Moura, an instance of the model constructing satisfiability calculus (mcSAT) framework [26]; and (2) the Non-uniform CAD (NuCAD) of Brown [27], which is a decomposition of the state space into disjoint cells that are truth-invariant for the conjunction of a set of polynomial constraints, but with a weaker relationship between cells (the decomposition is not cylindrical).
We demonstrate later with worked examples how our new approach outperforms a traditional CAD while still differing from the two methods above: with more effective learning from conflicts than NuCAD and, unlike NLSAT, an SMT-compliant approach which keeps theory reasoning separate from the SAT solver.

Paper Structure
We continue in Section 2 with preliminary definitions and descriptions of the alternative approaches. We then present our new algorithm, first the intuition in Section 3 and then formally in Section 4. Section 1 Doubly exponential in the number of variables (quantified or free). The double exponent does reduce by the number of equational constraints in the input [7], [8], [9]. However, the doubly-exponential behaviour is intrinsic: in the sense that classes of examples have been found where the solution requires a doubly exponential number of symbols to write down [10]. 2 There are other approaches to avoiding a full decomposition, for example, [21] suggests how segments of the decomposition of interest can be computed alone based on dimension or presence of a variety. 5 contains illustrative worked examples while Section 6 describes how our implementation performed on a large dataset from the SMT-LIB [28]. We conclude and discuss further research directions in Section 7.

Formulae in Real Algebra
The general problem of solving, in the sense of eliminating quantifiers from, a quantified logical expression which involves polynomial equalities and inequalities is an old one. Definition 1. Consider a quantified proposition in prenex normal form 3 : where each Q i is either ∃ or ∀ and F is a semi-algebraic proposition, i.e. a Boolean combination of constraints p j (x 1 , . . . , x m , x m+1 , . . . , x n )σ j 0, where p j are polynomials with rational coefficients and each σ j is an element of {<, ≤, >, ≥, =, =}.
The Quantifier Elimination (QE) Problem is then to determine an equivalent quantifier-free semi-algebraic proposition G(x m+1 , . . . , x n ).
Example 1. The formula ∃y. x · y > 0 is equivalent to x = 0; the formula ∃y. x · y 2 > 0 is equivalent to x > 0; whereas the formula ∀y. x · y 2 > 0 is equivalent to False.
The existence of a quantifier-free equivalent is known as the Tarski-Seidenberg Principle [29,5]. The first constructive solution was given by Tarski [5], but the complexity of his solution was indescribable (in the sense that no elementary function could describe it).

Cylindrical Algebraic Decomposition
A better (but still doubly exponential) solution had to await the concept of Cylindrical Algebraic Decomposition (CAD) in [6]. We start by defining what is meant by a decomposition here. Definition 2.

1.
A cell from R n is a non-empty connected subset of R n .

2.
A decomposition of R n is a collection D = {C 1 , . . . , C n } of finitely many pairwise-disjoint cells from R n with R n = ∪ n i=1 C i . 3. A cell C from R n is algebraic if it can be described as the solution set of a formula of the form p 1 (x 1 , . . . , x n )σ 1 0 ∧ · · · ∧ p k (x 1 , . . . , x n )σ k 0, where the p i are polynomials with rational coefficients and variables from x 1 , . . . , x n , and where the σ i are one of {=, >, <} 4 . Equation (2) is the defining formula of C, denoted as Def(C). A decomposition of R n is algebraic if each of its cells is algebraic. 4. A decomposition D of R n is sampled if it is equipped with a function assigning an explicit point Sample(C) ∈ C to each cell C ∈ D.
We are interested in decompositions with certain important properties relative to a set of polynomials as formalised in the next definition.
Definition 3. A cell C from R n is said to be sign-invariant for a polynomial p(x 1 , . . . , x n ) if and only if precisely one of the following is true: A cell C is sign-invariant for a set of polynomials if and only if it is sign-invariant for each polynomial in the set individually. A decomposition of R n is sign-invariant for a polynomial (a set of polynomials) if each of its cells is sign-invariant for the polynomial (the set of polynomials).
For example, the decomposition in Example 2 is sign-invariant for x 2 1 − 1. All of the techniques discussed in this paper to produce decompositions are defined with respect to an ordering on the variables in the formulae.
Definition 4. For positive natural numbers m < n, we see R n as an extension of R m by further dimensions, and denote the coordinates of R m as x 1 , . . . , x m and the coordinates of R n as x 1 , . . . , x n . Unless specified otherwise we assume polynomials in this paper are defined with these variables under the ordering corresponding to their labels, i.e., x 1 ≺ x 2 ≺ · · · ≺ x n . We refer to the main variable of a polynomial to be the highest one in the ordering that is present in the polynomial. If we say that a set of polynomials are in R i then we mean that they are defined with variables x 1 , . . . , x i and at least one such polynomial has main variable x i .
We can now define the structure of the cells in our decomposition. The cells of D which project to C ∈ D are said to form the cylinder over C. For a sampled decomposition, we also insist that the sample point in C be the projection of the sample points of all the cells in the cylinder over C. A (sampled) decomposition D of R n is cylindrical if and only if for each positive natural number m < n there exists a (sampled) decomposition D of R m over which D is cylindrical. So combining the above definitions we have a (Sampled) Cylindrical Algebraic Decomposition (CAD). Note that cylindricity implies the following structure on a decomposition cell.
Here i , u i , s i are functions in i − 1 variables (constants when i = 1). These functions could be indexed root expressions, involving algebraic functions (whose values for given x i , . . . , x n might be algebraic numbers).
Example 3. Fig. 1 shows two decompositions of R 2 in which, each dot, each line segment, and each hatched region are cells. The decomposition of R 2 on the left of the figure is cylindrical over R 1 (horizontal axis), but the decomposition on the right is not (two cells have overlapping non-identical projections onto x 1 ). All cells are locally cylindrical. Note that we are assuming x 1 ≺ x 2 ; for x 2 ≺ x 1 neither of the decompositions would be cylindrical, but still each cell would be locally cylindrical.
Collins' solution to the QE problem in Definition 1 was an algorithm to produce a CAD [6] of R n , sign-invariant for all the p j , and thus truth-invariant for F . Hence it is then only necessary to consider the truth or falsity of F at the finite number of sample points and query their algebraic definitions to form G. Furthermore, if Q i is ∀, we require F to be True at all sample points, whereas ∃ requires the truth of F for at least one sample point. As has been pointed out [23,25], such a CAD is finer than required for QE since as well as answering the QE question asked, it could answer another that involved the same p j ; but possibly different σ j , and even different Q i as long as the variables are quantified in the same order.

Our Setting: SMT for Real Algebra
In this paper we are interested in the subset of QE problems coming from Satisfiability Modulo Theories (SMT) solving [30], namely sentences (formulae whose variables are all bound by quantifiers) which use only the existential quantifier. I.e. (1) with m = n and all Q i = ∃.
Traditional SMT-solving aims to decide the truth of such formulae by searching for solutions satisfying the Boolean structure of the problem using a SAT-solver, and concurrently checking the consistency of the corresponding theory constraints. It additionally restricts to the case that F is a pure conjunction of the p j σ j 0; and makes the following additional requirements on any solution: 1. If the answer is True (generally referred to as SATisfiable), we need only an explicit point at which it is True (rather than a description of all such points); 2. If the answer is False (generally referred to as UNSATisfiable), we want some kind of justification that can be used to guide the SAT-solver's search.
So, in SMT solving we are less interested in complete algebraic structures but rather in deciding satisfiability and computing solutions if they exist. Reflecting this interest, Real Algebra is often called Non-Linear Real Arithmetic (NRA) in the SMT community. The use of CAD (and computer algebra tools more generally) in SMT has seen increased interest lately [31,32] and several SMT-solvers make use of these [33,34].

Relevant Prior Work
Our contribution is a new approach to adapt CAD technology and theory for this particular problem class. There are three previous works that have also sought to adapt CAD ideas to the SMT context: • Incremental CAD is an adaptation of traditional CAD so that it works incrementally, checking the consistency of increasing sets of polynomial constraints, as implemented in the SMT-RAT solver [33,35]. Traditional CAD as formulated by Collins and implemented in Computer Algebra Systems consists of two sequentially combined (projection and construction) phases: it will first perform all projection computations, and then use these to construct all the cells (with sample points). This incremental adaptation instead processes one projection operation at a time, and then uses that to derive any additional cells (with their additional samples), before making the next projection computation. For the task of satisfiability checking, this allows for early termination not only for satisfiable problems (if we find a sample that satisfies all sign conditions then we do not need to compute any remaining projections) but also for unsatisfiable ones (if the CAD for a subset of the input constraints is computed and no sample from its cells satisfies all those sign conditions). Although the implementation is technically involved, the underlying theory is traditional CAD and in the case of UNSAT that is exactly what is produced.
• The NLSAT Algorithm by Jovanović and de Moura [36] lifts the theory search to the level of the Boolean search, where the search at the Boolean and theory levels are mutually guided by each other away from unsatisfiable regions when it can be determined by some kind of propagation or lookahead. Partial solution candidates for the Boolean structure and for the corresponding theory constraints are constructed incrementally in parallel until either a full solution is found or a conflict is detected, meaning that the candidates cannot be extended to full solutions. Boolean conflicts are generalised using propositional resolution. For theory conflicts, CAD technology is used to generalise the conflict to an unsatisfying region (a CAD cell). These generalisations are learnt by adding new clauses that exclude similar situations from further search by the above-mentioned propagation mechanisms. In the case of a theory conflict, the incremental construction of theory solution candidates allows to identify a minimal set of constraints that are inconsistent under the current partial assignment. This minimal constraint set induces a coarser state space decomposition and thus in general larger cells that can be used to generalise conflicts and exclude from further search by learning. The exclusion of such a cell is learnt by adding a new clause expressing that the constraints are not compatible with the (algebraic description of the) unsatisfying cell. This clause will lead away from the given cell not only locally but for the whole future search when the constraints are all True.
• The Non-Uniform CAD (NuCAD) Algorithm by Brown [27] takes as input a set of polynomial constraints and computes a decomposition whose cells are truth-invariant for the conjunction of all input constraints. It starts with the coarsest decomposition, having the whole space as one cell, which does not guarantee any truth invariance yet. We split it to smaller cells that are sign-invariant for the polynomial of one of the input constraints, and mark all refined cells that violate the chosen constraint as final: they violate the conjunction and are thus truth-invariant for it. Each refinement is made by choosing a sample in a non-final cell and generalising to a locally cylindrical cell. At any time all cells are locally cylindrical, but there is no global cylindricity condition. The algorithm proceeds with iterative refinement, until all cells are marked as final. Each refinement will be made with respect to one constraint for which the given cell is not yet sign-invariant. There are two kinds of cells in a final NuCAD: cells that violate one of the input constraints (but these cells are not necessarily sign-nor truth-invariant for all other constraints), and cells that satisfy all constraints (and are sign-invariant for all of them). The resulting decomposition is in general coarser (i.e. it has less cells) than a regular CAD. Its cells are neither sign-nor truth-invariant for individual constraints, but they are truth-invariant for their conjunction and that is sufficient for consistency checking.

From Sample to Cell in Increasing Dimensions
We want to check the consistency of a set of input polynomial sign constraints, i.e., the satisfiability of their conjunction. Traditional CAD first generates algebraic information on the formula we study (the projection polynomials) and then uses these to construct a set of sample points, each of which represents a cell in the decomposition 5 . Our new approach works the other way around: we will select (or guess) a sample point by fixing it dimension-wise, starting with the lowest dimension and iteratively moving up in order to extend lower-dimensional samples to higher-dimensional ones. Thus we start with a zero-dimensional sample and iteratively explore new dimensions by recursively executing the following: • Given an (i−1)-dimensional sample s = (s 1 , . . . , s i−1 ) that does not evaluate any input constraint to False, we extend it to a sample (s 1 , . . . , s i−1 , s i ), which we denote as s × s i .
• If this i-dimensional sample can be extended to a satisfying solution, either by recursing or because it already has full dimension, then we terminate, reporting consistency (and this witness).
• Otherwise we take note of the reason (in the form of a conflict that explains why the sample can not be extended) and exclude from further search not just this particular sample s × s i , but all extensions of s into the ith dimension with any value from a (hopefully large) interval I around s i which will be unsatisfiable for the same reason.
This means that for future exploration, the algorithm is guided to look somewhere away from the reasons of previous conflicts. This should allow us to find a satisfying sample quicker, or in the case of UNSAT build a covering of all R n such that we can conclude UNSAT everywhere with fewer cells than a traditional CAD. In the last item above, the generalisation of an unsatisfying sample s × s i to an unsatisfying interval s × I will use a constraint p(x 1 , . . . , x i )σ0 that is violated by the sample, i.e., p(s × s i )σ0 does not hold. It might be the case that s i is a real zero of p(s) and we then have I = [s i , s i ]. Otherwise we get an interval I = ( , u) with < s i < u where is either the largest real root of p below s i or −∞ if no such root exists (and analogously u is either the smallest real root above s i or ∞). Thus we have p(s × r) = 0 for all r ∈ ( , u).
Since p(s × s i )σ0 is False and the sign of a polynomial does not change between neighbouring zeros, we can safely conclude that p(x 1 , . . . , x i )σ0 is violated by all samples (s × r) with r ∈ I. We continue and check further extensions of s = (s 1 , . . . , s i−1 ), until either we find a solution or the ith dimension is fully covered by excluding intervals. In the latter case, we take the collection of intervals produced and use CAD projection theory to rule out not just the original sample in R i−1 but an interval around it within dimension (i−1), i.e. the same procedure in the lower dimension. Intuitively, each sample s × s i violating a constraint with polynomial p can be generalised to a cell in a p-sign-invariant CAD. So when all extensions of s have been excluded (the ith dimension is fully covered by excluding intervals) then we project all the covering cells to dimension i−1 and exclude their intersection from further search.

Restoring Cylindricity
If we were to generalise violating samples to cells in a CAD that is sign-invariant for all input constraints then in the case of unsatisfiability we would explore a CAD structure. But instead, by guessing samples in not yet explored areas and identifying violated constraints individually per sample, we are able to generalise samples to larger cells that can be excluded from further search: like in the NLSAT approach.
Unlike NLSAT, we then build the intersection creating a cylindrical arrangement at the cost of making the generalisations smaller (but as large as possible while still ordered cylindrically). What is the advantage gained from a cylindrical ordering? In NLSAT the excluded cells are not cylindrically ordered; and so SMTmechanisms are used in the Boolean solver (like clause learning and propagation) to lead the search away from previously excluded cells. In contrast, our aim is to make this book-keeping remain inside the algebraic procedure, which can be done in a depth-first-search approach when the cells are cylindrically ordered.

Cylindrical Algebraic Coverings
So, we maintain cylindricity from CAD, but we relax the disjointness condition on cells in a decomposition, allowing our cells to overlap as long as their cylindrical ordering is still maintained. Instead of decomposition we will use the name covering for these structures. 1. A covering of R n is a collection D = {C 1 , . . . , C n } of finitely many (not necessarily pairwise-disjoint) cells from R n with R n = ∪ n i=1 C i .

2.
A covering of R n is algebraic if each of its cells is algebraic.
3. A covering D of R n is sampled if it is equipped with a function assigning an explicit point Sample(C) ∈ C to each cell C ∈ D.

4.
A cell C from R n is UNSAT for a polynomial constraint (set of constraints) if and only if all points in C evaluate the constraint (at least one of the constraints) to False. A covering of R n is UNSAT for a constraint (set) if each of its cells is UNSAT for the constraint (at least one from the set).

5.
A covering D of R n is said to be cylindrical over a covering D of R m if the projection onto R m of every cell of D is a cell of D . The cells of D which project to C ∈ D form the cylinder over C. For a sampled covering, the sample point in C needs to be the projection of the sample points of all the cells in the cylinder over C. A (sampled) covering D of R n is cylindrical if and only if for each 0 < m < n there exists a (sampled) covering D of R m over which D is cylindrical.
The coverings produced in our algorithm are all UNSAT coverings for constraint sets (i.e. at least one constraint is unsatisfied on every cell).

Differences to Existing Methods
Our new method shares and combines ideas from the related work in Section 2.4 but differs from each.
• A version of the CAD method which proceeds incrementally by constraint or projection factor is implemented in the SMT-RAT solver. Our new method differs as even in the case of unsatisfiability it will not need to compute a full CAD but rather a smaller number of potentially overlapping cells.
• While its learning mechanism made NLSAT the currently most successful SMT solution for Real Algebra, it brings new scalability problems by the large number of learnt clauses. Our method is conflictdriven like NLSAT, but instead of learning clauses, we embed the learning in the algebraic procedure.
Learning clauses allows to exclude unsatisfiable CAD cells for certain combinations of polynomial constraints for the whole future search, but it also brings additional costs for maintaining the truth of these clauses. Our approach remembers the unsatisfying nature of cells only for the current search path, and learns at the Boolean level only unsatisfiable combinations of constraints. To unite the advantages from both worlds, our approach could be extended to a hybrid setting where we learn at both levels (by returning information on selected cells to be learned as clauses).
• Like NuCAD, our algorithm can compute refinements according to different polynomials in different areas of the state space and to stop the refinement if any constraint is violated, however we retain the global cylindricity of the decomposition. Furthermore, driven by model construction, we can identify minimal sets of relevant constraints that we use for cell refinement, instead of the arbitrary choice of these polynomials used by NuCAD. Our expectation is thus that on average this will lead to coarser decompositions, and certainly rule out some unnecessary worst case decompositions.

Interface, I/O, and Data Structure
Our main algorithm, get unsat cover, is presented as Algorithm 2, while Algorithm 1 provides the user interface to it. A user is expected to define the set of constraints whose satisfiability we want to study globally, and then make a call to Algorithm 1 with no input. This will then call the main algorithm with an Algorithm 1: user call() Data: Global set C of polynomial constraints defined over R n . Input : None Output: Either (SAT, S) where S ∈ R n is a full-dimensional satisfying witness; or (UNSAT, C) when no such S exists, where C ⊂ C is also unsatisfiable. 1 (flag, data) := get unsat cover( () ) // Algorithm 2 2 if flag = SAT then 3 return (flag, data) 4 else 5 return (flag, infeasible subset(data)) empty tuple as the initial sample s; the main algorithm is recursive and will later call itself with non-empty input. In these recursive calls the input is a partial sample point s ∈ R i−1 which does not evaluate any global constraint defined over R i−1 to False, and for which we want to explore dimension i and above.
The main algorithm provides two outputs, starting with a flag. When the flag is SAT then the partial sample s was extended to a full sample from R n (the second output) which satisfies the global constraints. When the flag is UNSAT then the method has explored the higher dimensions and determined that the sample cannot be extended to a satisfying solution. It does this by computing an UNSAT cylindrical algebraic covering for the constraints with the partial sample s substituted for the first i−1 variables. Information on the covering, and the projections of these cells, are all stored in the second output here.
More formally, the output I is a set of objects I each of which represent an interval of the real line (specifically s × R). We use I later to mean both the interval, and our data structure encoding it which carries additional algebraic information. In total such a data structure I has six attributes, starting with: • the lower bound ; • the upper bound u; • a set of polynomials L that defined ; • a set of polynomials U that defined u. The bounds are constant, but potentially algebraic, numbers. The polynomials define them in that they are multivariate polynomials which when evaluated at s became univariate with the bound as a real root. The final two attributes are also sets of polynomials: • a set of polynomials P i (multivariate with main variable x i ) • a set of polynomials P ⊥ (multivariate with main variable smaller than x i ). These polynomials have the property that allows for generalisation of s to an interval: the property is that perturbations of s which do not change the signs of these polynomials will result in the interval I (whose numerical bounds may have also perturbed but will still be defined by the same ordered real roots of the same polynomials) remaining a region of unsatisfiability.
Within a covering there must also be special intervals which run to ∞ and −∞. For intervals with these bounds we store the polynomials from the constraints which allowed us to conclude the infinite interval.
In the case of UNSAT, the user algorithm will have to process the data I into an infeasible subset (Algorithm 1 Line 5), i.e. a subset of the constraints that are still unsatisfiable. Ideally this would be minimal (a minimal infeasible subset) although any reduction of the full set would carry benefits. We discuss how we implement this later in Section 4.6. We note that the correctness of Algorithm 1 follows directly from the correctness of its sub-algorithms.

Initial Constraint Processing
The first task in Algorithm 2 is to see what effect the partial sample s has on the global constraints. This is described as Algorithm 3, which will produce those intervals I ⊆ R such that s × I is conflicting with some input constraints (a partial covering). This method resembles a CAD lifting phase where we substitute a sample point into a polynomial to render it univariate, compute the real roots, and decompose the real line into sign-invariant regions. Here we do the same for the truth of our input constraints.
Algorithm 3 Lines 5−9 deal with the case where after substitution the truth value of the constraint may be immediately determined (e.g. the defining polynomial evaluated to a constant). The constraint either provides the entire line as an UNSAT interval, or no portion of it. The rest of the code deals with the case where the substitution rendered a constraint univariate in the ith variable. We use real roots(p, s) to return all real roots of a polynomial p over a partial sample point s that sets all but one of its variables. We do not specify the details of the real root isolation algorithm here 6 but note that it will need to handle potentially algebraic coefficients. We assume roots are returned in ascending numerical order with any multiple roots represented as a single list entry.
The inner for loop queries a sample point in each region of the corresponding decomposition of the line to determine any infeasible regions for the constraint, storing them in the output data structure I. I Algorithm 2: get unsat cover(s) Data: Global set C of polynomial constraints defined over R n . Input : Sample point s = (s 1 , . . . , s i−1 ) ∈ R i−1 such that no global constraint evaluated at s is False. Note that if s = () then i = 1 (i.e. we start from the first dimension). Output: Either (SAT, S) where S ∈ R n is a full-dimensional satisfying witness; or (UNSAT, I) when s cannot be extended to a satisfying sample in R n . I represents a set of intervals which cover s × R and come with algebraic information (see Section 4.1). represents a set of intervals I such that s × I conflicts with some input constraint. The intervals from I may be overlapping, and some may be redundant (i.e. fully contained in others). We discuss this issue of redundancy further in Section 4.5.
It is clear that Algorithm 3 will meet its specification, in that it will define intervals on which constraints are unsatisfiable: the falsity of a constraint caused the inclusion of an interval in the output while the bounds of the interval were defined to ensure that the polynomial defining the failing constraint would not change sign inside. The role and property of the stored algebraic information will be discussed in Section 4.4.
The call to Algorithm 3 initialises I in Algorithm 2 in which we will build our UNSAT covering. It is unlikely but possible that I already covers R after the call to Algorithm 3: it could even contain (−∞, ∞). If I is already a covering then we would skip the main loop of Algorithm 2 and directly return it. Examples of where this may be triggered would be constraints y 2 x < 0 or y 2 + x < 0 at sample s = (x → 0). The former would have been returned early by Line 7 of Algorithm 3 while the latter would have required real root isolation and have been returned in Line 21 of Algorithm 3.

The Main Loop of Algorithm 2
We will iterate through Lines 2 − 12 of Algorithm 2 until the set of intervals represented by I cover all R. In each iteration we collect additional intervals for our UNSAT covering.
To do this we first generate a sample point s i from R \ (∪ I∈I I) using a subroutine sample outside(I) which is left unspecified. This could simply pick the mid-point in any current gap, or perhaps something more sophisticated (a common strategy would be to prefer integers, or at least rationals).
Note that s × s i necessarily satisfies all those constraints with main variable x i , otherwise we would have generated an interval excluding s i at Line 1. This means that: (a) if s × s i has full dimension, we have found a satisfying sample point for the whole set of constraints and can return this along with the SAT flag in Line 5; and (b) if not full dimension then we will meet the input specification for the recursive call on Line 6. The recursion means we will explore s × s i in the next dimension up. The previous check on dimension acts as the base case and thus the recursion clearly terminates.
If the result of the recursive call is SAT we simply pass on the result in Line 8. Otherwise we have an UNSAT covering for s × s i and our next task is to see whether we can generalise this knowledge to rule out Algorithm 3: get unsat intervals(s) Data: Global set C of polynomial constraints defined over R n . Input : Sample point s of dimension i − 1 with i a positive integer (s is empty for i = 1). Output: I which represents a set of unsatisfiable intervals over s × I and comes with some algebraic information (see Section 4.1). Regions Add I to I 21 return I not just s i , but an interval around it. We do this in two steps.
• First in Line 10 we call Algorithm 4 to construct a characterisation from the UNSAT covering: a set of polynomials which were used to determine unsatisfiability and with the property that on the region around s × s i where none of them change sign the reasons for unsatisfiability generalise. I.e., while the exact bounds of the intervals in the coverings may vary: (a) they are still defined by the same ordered roots of the same polynomials (over the sample); and (b) they do not move to the extent that the line is no longer covered.
• Second in Line 11 we call Algorithm 5 to find the interval in dimension i over s in which those characterisation polynomials are sign-invariant.
We describe these two sub-algorithms in detail in the next subsection. The interval produced by Line 11 Algorithm 5 may be (−∞, ∞). In this case all other intervals in I are now redundant and the main loop of Algorithm 2 stops. Otherwise we continue to iterate.
The loop will terminate: although the search space is infinite the combinations of constraints that determine satisfiability is not. Eventually we must either cover all the line or sample in a satisfying region. The correctness of Algorithm 2 thus depends on the correctness of the two crucial sub-algorithms 4 and 5.

Generalising the UNSAT Covering from the Sample
It remains to examine the details of how the UNSAT covering is generalised from the sample s to an interval around it (the calls to Algorithms 4 and 5 on Lines 10 and 11 of Algorithm 2).

Ordering within a covering
The input to Algorithm 4 is an UNSAT covering I whose elements define intervals I which together cover R. For an example of such a covering see Fig. 2. There may be some redundancy here, like the second interval (from 2 to u 2 ) in Fig. 2 which is completely covered by the first and third intervals already. To ensure soundness of our approach we need to remove at least those intervals which are included within a single other interval. We discuss more details on dealing with redundancies in coverings in Section 4.5.
For now we simply make the reasonable assumption of the existence of an algorithm compute cover which computes such a good covering of the real line as a subset of an existing covering I. Since it is not crucial we will not specify the algorithm here, but we note that the ideas in [37] may be useful.
We assume that compute cover orders the intervals in its output according to the following total ordering: i.e. ordered on i with ties broken by u i . We will always have 1 = −∞ and u k = ∞ with the remaining bounds defined as algebraic numbers (possibly not rational). We further require that only possible if we exclude the cases where one interval is a subset of another. Note that this is not just an optimisation but is actually crucial for the correctness of this approach, as explained in Section 4.5.
The intervals we consider are either open (if = u) or point intervals (if = u). Note that it might make sense to extend the presented algorithm to allow for closed (or half-open) intervals as well, for example when built from weak inequalities. This could help to avoid some work on individual sample points like 4 in Fig. 2 and the point intervals we deal with in the examples in Section 5. Such changes are straight-forward to implement and so we do not discuss them here to simplify the presentation.

Constructing the characterisation
The first line of Algorithm 4 ensures the ordering specified above, while the remainder uses CAD projection ideas to collect everything we need to ensure that the UNSAT covering I stays valid when we generalise the underlying sample point s i later on. We include polynomials for a variety of reasons.
• First in Line 5 we pass along any polynomials with a lower main variable that had already been stored in the data structure. These were essentially produced by the following steps at any earlier recursion, by an act of projection that skipped a dimension. For example, the projection of any polynomial f (x 1 , x 3 ) into (x 1 , x 2 )-space will actually give a univariate polynomial in x 1 .
• Next in Lines 6 and 7 we identify polynomials that will ensure the existence of the current lower and upper bounds. We add discriminants, whose zeros indicate where the original polynomial has multiple roots, and leading coefficients (with respect to the main variable), whose zeros indicate asymptotes of the original polynomial. We may also require additional coefficients, as discussed in Section 4.4.4. If we ensure these polynomials do not change sign then we know that the algebraic varieties that defined and u continue to exist (and no other varieties are spawned).
• In Lines 8 and 9 we generate polynomials whose sign-invariance ensures that and u stay the closest bounds. I.e. we avoid the situation where they are undercut by those coming from some other variety.
We need only concern ourselves with those coming from the direction of the bound. For example, when protecting an upper bound we need only worry about roots that are above it and take resultants accordingly on Line 9. This is because any root coming from below would need to first pass through the lower bound and the resultant from Line 8 would block generalisation past that point.
• In Lines 10−11 we finally derive resultants to ensure that the overlapping lower and upper bounds of adjacent intervals do not cross, which would disconnect the UNSAT covering (leaving it not covering some portion of the line). The correctness of this step requires an ordering with lack of redundancy, as specified above and discussed in detail in Section 4.5. This step also has the effect of ensuring the intervals do not overlap further to the extent that one then becomes redundant.
Recall that we may have satisfiability over s refuted by a single constraint in Algorithm 3 (i.e. the defining polynomial cannot change sign over s). Then after Line 1 of Algorithm 4 the data structure I contains a single interval (−∞, ∞) and we have no resultants to compute. Algorithm 4 finishes in Line 12 with some standard CAD simplifications to the polynomials. These all stem from the fact that polynomials matter only so much as where they vanish. E.g.
• Remove any constants, or other polynomials than can easily be concluded to never equal zero.
• Normalise the remaining polynomials to avoid storing multiple polynomials which define the same varieties. E.g. Multiply each polynomial by the constant needed to make it monic (i.e. divide by leading coefficient so for example polynomial 2x − 1 becomes x − 1 2 ). Other normal forms include the primitive positive integral with respect to the main variable.
• Store a square free basis for the factors rather than the polynomials themselves (this much is necessary, else later resultants/discriminants will vanish); or even fully factorising (optional, but generally favoured for the efficiency gains it can bring).

Constructing the generalisation
Now let us discuss how this characterisation (set of polynomials) is used to expand the sample to an interval by Algorithm 5. We first separate the polynomials into P i and P ⊥ where P i are those polynomials that contain x i and P ⊥ the rest. We then identify the crucial points over s beyond which our covering may cease to be. This step essentially evaluates the polynomials with main variable x i over the sample in R i−1 and calculates real roots of the resulting univariate polynomial. There is some additional work within this sub-algorithm call on Line 3 which we discuss in Section 4.4.4.
We then construct the interval around s i from the closest roots and u and collect the polynomials that vanish in and u in the sets L and U , respectively. We supplemented the real roots with ±∞ to ensure that and u always exist, i.e. we can have (−∞, u) or ( , ∞). In this case the corresponding set L or U is empty.
Algorithm 5: interval from characterization(s, s i , P ) Input : Sample point s ∈ R i−1 ; an extension s i to the sample; and set of polynomials P in Q[x 1 , . . . , x i ] that characterize why s × s i cannot be extended to a satisfying witness. Output: An interval I around s i such that on s × I the constraints are unsatisfiable for the same reasons as on s × s i .
In the case where P i has no real roots at all over s then the UNSAT covering is valid unconditionally over s, and the interval (−∞, ∞) is formed and passed back to the main algorithm (where it could be taken to form the next covering in its entirety).
Let us briefly consider some simple examples for Algorithm 5. Suppose we have variable ordering x ≺ y, the partial sample (x → 0) and that P contains only the polynomial whose graph defines the unit circle:

Correctness of the generalisation
The correctness of Algorithms 4 and 5 relies on traditional CAD theory, with Algorithm 4 an analogue of projection and Algorithm 5 an analogue of lifting. However, there are a number of subtleties to discuss.
First, regarding exactly which coefficients are computed in Algorithm 4 Line 7. As explained above, the leading coefficient is essential to include as its vanishing would indicate an asymptote. However, we must also consider what happens within regions of such vanishing: there the polynomial has reduced in degree, and so a lesser coefficient is now leading and must be recorded for similar reasons. We should take successive coefficients until one is taken that can be easily concluded not to vanish (e.g. it is constant) or if it may be concluded that the included coefficients cannot vanish simultaneously. In our context this is simpler: so long as a coefficient evaluates at the sample point to something non-zero then we need not include any subsequent coefficients (since by including that coefficient in the characterisation we ensure we do not generalise to where it is zero). We do need to include preceding coefficients that were zero as we must ensure that the generalisation does not cause them to reappear. This is described by Algorithm 6.
With such calculation of coefficients our characterisations are then close to the McCallum projection [38]: we ensure that individual polynomials are delineable but we cannot claim that for the characterisation as a set since we differ from [38] in which resultants are computed. The inclusion of resultants in CAD projection is to ensure a constant structure of intersections between varieties in each cell. McCallum projection takes all possible resultants between polynomials involved. Instead, we take only those needed to maintain the structure of our covering, as detailed by the arguments in Section 4.4.2.

Completeness of the generalisation
McCallum projection is not complete, i.e. there are some (statistically rare) cases its use is known to be invalid, and these could be inherited in our algorithm. The problem can occur when a polynomial vanishes at a sample of lower dimension (known as nullification) potentially losing critical information. For example, consider a polynomial zy−x which vanishes at the lower dimensional sample (x, y) = (0, 0). The polynomials' behaviour clearly changes with z but that information would be lost.
Such nullifications can be easily identified when they occur. We assume that the sub-algorithm used in Algorithm 5 Line 3 will inform the user of such nullification. What should be done in this case? An extreme option would be to recompute the characterisation to include entire subresultant sequences as in Collins' projection [6], or to use the operator of Hong [39]. A recent breakthrough in CAD theory could offer a better option: [40] proved that Lazard projection [41] is complete 7 . The Lazard operator includes leading and trailing coefficients, and requires more nuanced lifting computations. Instead of evaluating polynomials at a sample and calculating real roots of the resulting univariate polynomials, we must instead perform a Lazard evaluation [40] of the polynomial at the point. This will substitute the sample coordinate by coordinate, and in the event of nullification divide out the vanishing factor to allow the substitution to continue. Thus no roots are lost through nullification.
We have not yet adapted our algorithm to Lazard theory: although SMT-RAT already has an implementation of Lazard evaluation, we are not yet clear on how the polynomial identification in Algorithm 5 should be adapted in cases of nullification. Also, it is not trivial to argue the safe exclusion of trailing coefficients in cases where the leading coefficient is constant, as it is with McCallum, since the underlying Lazard delineability relies heavily on properties of the trailing coefficient. Our current implementation is hence technically incomplete, however, we can produce warnings for such cases. The experiments on the SMT-LIB detailed in Section 6 show that these nullifications are a rare occurrence.

Redundancy of Intervals
We discussed in Section 4.4.1 that the set of intervals in a covering may contain redundancies. We distinguish between two possibilities for how an interval I can be redundant: 1. I is covered by another interval entirely, as interval ( 2 , u 2 ) is by interval ( 1 , u 1 ) in Fig. 3. 2. I is covered but only through multiple other intervals, as the interval defined by the bounds of p 2 would be by those from p 1 and p 3 in Fig. 4.
We now explain why we need to remove redundancies of the first kind, but the second could be kept.

Redundancy of the first kind
The resultants we produce in Lines 10−11 of Algorithm 4 are meant to ensure that consecutive intervals (within the ordering) continue to overlap on the whole region we exclude. Thus we must guarantee that they overlap on our sample point in the first place.
In the examples of Fig. 3 we assume that inside the curves is infeasible, and so the point marked in Fig. 3b is a satisfying witness. The two examples differ only by a small change to the polynomial p 2 : in Fig. 3a p 2 intersects p 3 while in Fig. 3b it does not. For both examples the numbering of the intervals means we will calculate the resultants res(p 1 , p 2 ), res(p 2 , p 3 ) but not res(p 1 , p 3 ). In each case the bounds obtained are indicated by the dashed lines: in Fig. 3a they come from the roots of res(p 2 , p 3 ) while in Fig. 3b this resultant has no real roots and the closest bounds instead come from the discriminant of p 2 .
So, in Fig. 3a the excluded region is all unsatisfiable, i.e. correct, but only by luck! The resultant that bounds the excluded regions has no relation to the true bound (the intersection of p 1 and p 3 ). In Fig. 3b we exclude too much and erroneously exclude the dot which actually marks a satisfying sample.
We do not see an easy way to distinguish the two situations presented in Fig. 3 and we therefore excluded all redundancies of this kind in Section 4.4.1, by requiring that we have the stronger ordering (4) rather than just the weaker version (3).

Redundancy of the second kind
The error described above of excluding a satisfiable point because of a redundancy of the first kind, could not occur in the presence of a redundant interval of the second kind. The algorithm assumes that every adjacently numbered interval overlaps, which is false for a redundancy of the first kind but still true for one of the second.
However, should we still remove redundant intervals of the second kind for efficiency purposes? Removal would mean less intervals in the covering, and thus less projection in Algorithm 4 and less root isolation in Algorithm 5. However, it does not mean the generalisation has to be bigger.
Intuitively, reducing such a redundancy makes the overlap between adjacent intervals smaller. But this may also mean that the interval we can exclude in the lower dimension is smaller than it would have been if we kept the redundant interval, retaining a larger overlap. Consider the two examples from Fig. 4: the polynomials differ but the geometric situation and position of intervals is similar. In both cases an interval defined by the bounds of p 2 will be redundant when combined with those from the bounds of p 1 and p 3 . If we do not reduce the covering by the redundant interval then we exclude I f ull but if we do we exclude I reduced . The excluded region would grow from reduction in Fig. 4a but shrink in Fig. 4b.
We do not see an easy way to check which situation we have (other than completely calculating both and taking the better one). The decision of whether to reduce could be taken heuristically by an implementation. Further investigation into this would be an interesting topic for future work.

Embedding as Theory Solver
Our algorithm can be used as a standalone solver for a set of real arithmetic constraints, but our main motivation is to use it as a theory solver within an SMT solver. Such theory solvers should be SMT compliant: • They should work incrementally, i.e allow for the user to add constraints and solve the resulting problem with minimal recomputation.
• They should similarly allow for backtracking, i.e. the removal of constraints by the user.
• They should construct reasons for unsatisfiability, which usually refers to a subset of the constraints that are already unsatisfiable, often referred to as infeasible subsets.
For infeasible subsets we store, for every interval, a set of constraints that contributed to it, in a similar way to what was called origins in [35]. For intervals created in Algorithm 3 we simply use the respective constraint c as their origin. For intervals that are constructed in Algorithm 5 the set of constraints is computed as the union of all origins used in the corresponding covering in Algorithm 4. The constraints are then gathered together as the final step before returning an UNSAT result in Algorithm 1 (Line 5).
We have yet to implement incrementality and backtracking, but neither pose a theoretical problem. Though possibly involved in the implementation, the core idea for incrementality is straight-forward: after we found a satisfying sample we retain the already computed (partial) coverings for every variable and extend them with more intervals from the newly added constraints. For backtracking, we need to remove intervals from these data structures based on the input constraints they stem from. As we already store these to construct infeasible subsets, this is merely a small addition.

Simple SAT Example in 2D
We start with a simple two-dimensional example to show how the algorithm proceeds in general. Our aim is to determine the satisfiability for the conjunction of the following three constraints: The three defining polynomials are graphed in Fig. 5, with the unsatisfiable regions for each adding a layer of shading (so the white regions are where all three constraints are satisfied). We now simulate how our algorithm would process these constraints, under variable ordering x ≺ y, to find a point in one of those satisfying regions. The user procedure starts by calling the main algorithm get unsat cover() with an empty tuple as the sample.
get unsat cover(s = ()). Since all the constraints are bivariate, the call to get unsat intervals() cannot ascertain any intervals to rule out, initialising I = ∅. We thus enter the main loop of Algorithm 2 and sample s 1 = 0. This sample does not have full dimension and we perform the recursive call in Line 6. get unsat cover(s = (x → 0)). This time, the call of get unsat intervals(s = (x → 0)) sees all three constraints rendered univariate by the sample. Each univariate polynomial has one real root and thus decomposes the real line into three intervals.

Constraint Real Roots Intervals
After analysing each of the 9 regions we identify 6 for which a constraint is infeasible. Thus I is initialised as a set of six intervals as below (where p i refers to the defining polynomial for constraint c i ).
Although (−∞, ∞) is not an interval of I we observe that R is already covered by the two intervals (−∞, 1 2 ) and (−1, ∞). Note how the second of the three constraints is not part of the conflict which simplifies the following calculations.
Since the real line is already covered we do not enter the main loop of Algorithm 2 in this call. Instead we immediately proceed to the final line where we return I to the original call of get unsat cover(s = ()). Since the flag is UNSAT the next step there is constructing a characterisation in Line 10.
construct characterisation(s = (x → 0), I). As already observed, the intervals (−∞, 1 2 ) and (−1, ∞) cover R and thus the call to compute cover() at the start should simplify I to the following: As P i only contains a single polynomial in each interval there are no resultants to calculate in the first loop, and only a single one from the second: res y (4 · y − x − 2, 4 · y − x 2 + 4) = −4 · x 2 + 4 · x + 24 Further, P ⊥ = ∅ and the discriminants and leading coefficients all evaluate to constants: disc y (4 · y − x 2 + 4) = 1, lcoeff y (4 · y − x 2 + 4) = 4.
So the output set from construct characterisation(s = (x → 0), I) consists of a single polynomial: Then in the original call to get unsat cover(s = ()) we continue working in Line 11 with the construction of an interval from the characterisation.
interval from characterisation(s = ∅, s i = (x → 0), P = {x 2 −x−6}). We see that P ⊥ = ∅ and P i = P and obtain the real roots {−2, 3}. Consequently we obtain l = −2, u = 3 with both L and U as P i . We thus return the unsatisfiable interval (−2, 3, P i , P i , P i , ∅). Back in our initial call to Algorithm 2 we add this interval to I and continue with a second iteration of the main loop. This time we must sample some value of x outside of (−2, 3), for example x = −3 or x = 4, both of which can be extended to a full satisfying sample: ((x → −3, y → 0) or (x → −4, y → 2)). In fact, any sample for x outside of this interval can be extended to a full assignment and that would be always discovered during the next recursive call: first Algorithm 5 would rule out any extension that is infeasible and then sample outside would pick a satisfying value in the first iteration of the main loop. This would form a full dimensional sample which is returned on Line 5, and then the initial call would return it to the user via Line 8.

Comparing to Incremental CAD
Let us consider how this would compare with the incremental version of traditional CAD. Recall that this performs one projection step at a time, and then refines the decomposition with respect to the output (producing extra samples). Thus the computation path depends on the order in which projections are performed. The discriminants and coefficients of the input do not contribute anything meaningful to the projections, so it all depends in which order the resultants are computed. If the implementation were unlucky and picked the least helpful resultant it will end up performing more decomposition than is required. For this particular example the effect is mitigated a little by the implementation's preference for picking integer sample points allowing it to find a satisfying sample earlier than guaranteed by the theory: i.e. the cell being formed by the decomposition is not truth invariant for all constraints, but by luck the sample picked is SAT and so the algorithm can terminate early. So for this example our incremental CAD performs no more work than the new algorithm, but that is through luck rather than guidance. In contrast, the superiority of the new algorithm over incremental CAD is certain for the next example.

Comparison to NLSAT
We may also compare to the NLSAT method [36] which also seeks a single satisfying sample point for all the constraints. Like our new method, NLSAT is model driven starting with a partial sample; it explains inconsistent models locally using CAD techniques; and thus exploits the locality to exclude larger regions; combining multiple of these explanations to determine unsatisfiability. The difference is that the conflicts are then converted into a new lemma for the Boolean skeleton and passed to a separate SAT solver to generate a new model. The implementation of NLSAT in SMT-RAT (see Section 6.3 performs the following steps for the first example: theory model explanation clause excluded interval The boundaries of the first excluded interval are the two real roots of the polynomial 2x 2 −2x−7 which is the resultant of the defining polynomials for c 1 and c 2 . Rather than give these as surds or algebraic numbers we use the decimal approximation to shorten the presentation. Throughout, a decimal underlined refers to a full algebraic number that we have computed but choose not to display for brevity. For this example NLSAT took three conflicts to find a satisfying model (compared to two for our new algorithm). However, the difference is due to luck. In the first iteration NLSAT was unlucky to select a subset of constraints that rules out a smaller UNSAT region. It could have instead chosen c 1 and c 3 in the first iteration, essentially yielding the very same computation as our method. Conversely our algorithm could also have chosen c 1 and c 2 as a cover and then possibly have needed another iteration.

Comparison to NuCAD
Our algorithm as stated is designed to solve satisfiability problems, and thus terminates as soon as a satisfying sample is found. Thus it does not compare directly with NuCAD which builds an entire decomposition of R n relative to the truth of a quantifier free logical formula which can then be used to solve more general QE problems stated about that formula 8 . NuCAD constructs the decomposition of one cell (with known truth-value) in R n at a time, so there is a natural variant of the algorithm where we construct cells only until we find one in which the input formula is True. Consider how NuCAD might proceed for the simple example. A natural starting point would be to consider the origin first. NuCAD would recognise that constraints are not satisfied here and choose one of them to process. Let us assume NuCAD chooses in the order the constraints are labelled (i.e. pick c 1 ); then like us it would generalise this knowledge beyond the sample and decompose the plane as in Fig. 6. The shaded cell is known to be UNSAT while the white region is a cell whose truth value is as yet unknown.  Fig. 7 where picking c 2 leads to the decomposition on the left, and picking c 3 to the decomposition on the right (similar to the two possible steps our algorithm had).
It is possible under a natural strategy for NuCAD to find a satisfying point on its second sample, but this is not guaranteed by the theory. In comparison, our algorithm is more guided, if starting at the origin then it could not take more than three iterations for this example.

More Involved UNSAT Example in 2D
In the first example above we saw that the new algorithm was able to learn from conflicts to guide itself towards a satisfying sample. However, due to luck and sensible implementation choices the alternatives could process the example with a similar amount of work to the new algorithm. So we next present a more involved example that will make clearer the savings offered by the new approach. This example will ultimately turn out to be unsatisfiable and thus the incremental CAD of SMT-RAT will in the end perform all projection and produce a regular CAD, as for example would be produced by a traditional CAD implementation like QEPCAD. The advantage of our algorithm in this case is that we can avoid some parts of the projection that are the most expensive.
Our aim for this example is to determine the satisfiability for the conjunction of the following five constraints: They are depicted graphically in Fig. 8, with each constraint again adding a layer of shading where it is unsatisfiable. This time we see there is no white space and so together the constraints are unsatisfiable. Before we proceed let us remark on how the example was constructed. There are two constraints of high degree (c 1 and c 5 ) while the rest are fairly simple. The problem structure was chosen so that c 1 and c 2 conflict on the left part of the figure, c 2 and c 3 in the middle and c 4 , and c 5 on the right. Thus while both c 1 and c 5 are important to determining unsatisfiability, they never need to be considered together to do this. I.e. we have no need of their resultant (with respect to y): an irreducible degree 33 polynomial in x which is dense (34 non-zero terms) with coefficients that are mostly 20 digit integers. In this example c 1 and c 5 are of degree 11 but we could generalise the example by increasing the degree of these polynomials arbitrarily while retaining the underlying problem structure which allows us to avoid working with them together. Hence the advantage over algorithms which do consider them together can be made arbitrarily large.
We now describe how our new algorithm would proceed for this example but skip some details compared to the explanation of the previous example, to avoid unnecessary repetition.
get unsat cover(s = ()). As all constraints contain y we get no unsatisfiable intervals in the first dimension. We enter the loop and sample s 1 = 0 and enter the recursive call.
get unsat cover(s = (x → 0)). All constraints are univariate and get unsat intervals(s = (x → 0)) returns I defining the following 10 unsatisfiable intervals: We can select intervals (−∞, 1 2 ) and (−1, ∞) to cover R and return this to the main call 9 . The main call analyses the covering and constructs the characterisation {x 2 + x − 3} leading to the exclusion of the interval (−2.30, 1.30). A graphical representation is shown in the left image of Fig. 9. We then iterate through the loop again selecting a sample outside of this interval, say s 1 = 2, with which we enter another recursive call. get unsat cover(s = (x → 2)). This call is almost identical to the one above, only we are now forced to use c 4 to build a conflict with c 3 instead of c 2 . The unsatisfiable intervals that initialise I are now as follows: We obtain the (unique) minimal covering from these by selecting (−∞, −1) and (−1.33, ∞) provided by c 5 and c 4 respectively. So we skip the loop and return this to the original parent call.
In the parent call the characterisation (after some simplification) contains two polynomials of degrees nine and eleven from the discriminant for c 5 and the resultant of the pair. They each yield one real root: 1.19 from the resultant (corresponding to the crossing point of c 4 and c 5 ) and 3.18 from the discriminant corresponding to the point of inflection of c 5 . These are visualised in the right image of Fig. 9. Note that the latter point is spurious in the sense that the point of inflection has no bearing on the satisfiability of our problem. But while we could safely ignore it, we have no algorithmic way of doing so. Thus we exclude the interval (1.19, 3.18) and proceed.
We now continue iterating the main loop with first the sample point s 1 = 3.18 and then s 1 = 4, which both yield the same conflict based on c 5 and c 4 and thus the exact same characterisation, excluding [3.18, 3.18] and (3.18, ∞) in turn. Finally we select s 1 = −3 and recurse one last time.
get unsat cover(s = (x → −3)). Once again we compute the unsatisfiable intervals and obtain the following intervals in the initialisation of I: We see that we have covered R and thus we return a final answer of UNSAT. Note that the highest degree of any polynomial we used in the above was eleven: the degree 33 resultant of c 1 and c 5 was never used nor even computed. It is important to recognise the general pattern here that allows our algorithm to gain an advantage over a regular CAD. We have two more complicated constraints involved in creating the conflict, but they are separated in space, or at least the regions where they are needed to construct the conflict are separated. Increasing the degrees within c 1 and c 5 would affect the algorithm only insofar that we have to perform real root isolation on larger polynomials while regular CAD has to deal with the resultant of these polynomials whose degree grows quadratically.

Comparing to Incremental CAD
The full projection contains one discriminant and 10 resultants (plus also some coefficients depending on which projection operator is used). Since no satisfying sample will be found the incremental CAD will actually produce a full sign-invariant CAD for the polynomials in the problem, containing 273 cells.
There is scope for some optimisations here, for example, because the constraints are all strict inequalities we know that they are satisfied if and only if they are satisfied on some full dimensional cell and so we could avoid lifting over any restricted dimension cells. The CAD decomposes the real line into 27 cells and so we could optimise to lift only over the 13 intervals to consider 77 cells in the plane. This avoids some work but we still need to compute and isolate real roots for the full projection set, and perform further lifting and real root isolation over half of the cells in R 1 . This is significantly more real root isolation and even projection than computed by the new algorithm, in particular, the large resultant discussed above.

Comparison to NLSAT
SMT-RAT's implementation of NLSAT performs the following samples and conflicts for this example.
Like the new algorithm it starts with the origin and gradually rules out the whole x-axis by generalising from a sample each time. However, NLSAT required 6 iterations to do this, compared to 5 for the new algorithm.
The main reason for this is that it uses c 2 and c 5 (instead of c 4 and c 5 ) to explain the conflict at x → 5 and thus obtains a slightly smaller interval, leaving a gap between 3.184 and 3.199. However, this poorer choice could have also been made by our new algorithm.

Comparison to NuCAD
This time, because the example is UNSAT, NuCAD will have to complete and produce an entire decomposition of the plane. The order in which cells are constructed will be driven by implementation as the only specification in the algorithm is to choose samples outside of cells with known truth value. Since any reasonable implementation would start with the origin and prefer integer samples it is likely that NuCAD's computation path would follow similarly to our algorithm for this example. But it should be noted that this is not guaranteed. For example, if instead of the origin NuCAD were to start with the model (0, −50) then here c 5 is violated and it may start by constructing the shaded cell in Fig. 10b. It must now pick a model outside of the shaded region. Again, it would be a strange choice but the model (−2, −50) would be acceptable and here c 1 is violated. The necessary splitting would then require the calculation and real root isolation of the large resultant of the defining polynomials of c 1 and c 5 . There is one such real root, indicating a real intersection of the two polynomials as shown in Fig. 10a, and this would necessarily be part of the boundary in the constructed cell.
We accept that the above model choices would be strange but they are a possibility for NuCAD while guaranteed to be avoided by the new algorithm. Of course, it should be possible to shift the coordinate system of the example to make these model choices seem reasonable.

Simple 3D Example
The 2D examples demonstrated how the new algorithm performs less work than a CAD and is more guided by conflict than NuCAD. However, to demonstrate the benefits over NLSAT we need more dimensions.  (7) and on the right the cures of the (x, y)-plane computed as the first characterisation.
Consider the simple problem in 3D of simultaneously satisfying the constraints: These require us to be inside two overlapping spheres as on the left of Fig. 11. Let us consider how the new algorithm would find a witness.
get unsat cover(s = ()). Our algorithm cannot draw any conclusions from the constraints as none are defined only in x, so it samples x = 0 and enters the recursive call.
get unsat cover(s = (x → 0)). Once again no conclusions can be drawn from the constraints yet so we sample y = 0 and recurse again.
get unsat cover(s = (x → 0, y → 0)). This time the call of get unsat intervals(s = (x → 0, y → 0)) produces three unsat intervals: (−∞, −1), and (1, ∞) from the first constraint and (− 5 4 , ∞) from the second. The first and third are together already a covering so we skip the main loop and return to the previous call.
get unsat cover(s = (x → 0)) continued. For the characterisation we calculate three polynomials: the discriminants of each of the two defining polynomials and their resultant which after simplification are y 2 + x 2 − 1, 4y 2 + 4x 2 − 12y + 5 and 12y − 9 respectively. They are graphed on the right of Fig. 11.
When forming the interval around y = 0 we obtain the set Z = {−∞, −1, 1 2 , 3 4 , 1, 5 2 , ∞} and so we generalise to (−1, 1 2 ). Sampling for y anywhere outside ( 1 2 , 1) would lead to a full covering of the z-dimension after the initial querying of constraints in the recursive call. While any sample from within that interval would then be simply extended with z = 0 to a full satisfying witness.

Comparison to NLSAT
NLSAT would proceed similarly in sampling (x, y) = (0, 0) and discovering the conflict. However, it would then immediately build a cell around (0, 0) requiring the computation of the projection of those polynomials in the conflict (e.g. two discriminants and a resultant). Only then would it move to a new partial sample. In contrast, the new algorithm only computes projections with respect to the second variable once it has determined there is no possible y to extend the x-value. Since in this example there is such a y for the first x-value those projection need never be computed.

More Involved 3D Example
We finish with a larger 3D example to demonstrate some of the facets of the algorithm not yet observed in the smaller examples. We define the three polynomials: and seek to determine the satisfiability of We use variable ordering x ≺ y ≺ z, and note that unlike the previous example we have here not only a higher dimension, but also a non-trivial leading coefficient for g and a constraint that is not in the main variable formed by h. Finally please note that z appears only as z 2 in the constraints and thus facts drawn for positive z may be applied similarly for negative 10 . The surfaces defined by the polynomials are visualised in Fig. 12 where the red surface is for f , the blue for g and the green for h. We see that f is a hyperboloid, while g would have been a paraboloid if it were not for the leading coefficient in z. The final constraint simply bounds y to (−10, 10). We see that there are many non-trivial intersections (and self intersections) between the surfaces. A full sign-invariant CAD for the three polynomials may be produced with the Regular Chains Library for Maple [43] in about 20 seconds, and contains 3509 cells.
For use in the discussion later, let us examine and put a label on all the CAD projection polynomials in (x, y) (i.e. those obtained after z is projected out) 11 : Here p 1 is the discriminant of f with respect to z (at least up to a constant factor); p 2 and p 3 are factors of the discriminant of g and also the trailing and leading coefficients of g respectively. The resultant of f and g with respect to z evaluated to p 2 4 . Of course, the projection also includes h itself and is shown in Fig. 13. Also for use in the following discussion, we label the six real roots in y that these four polynomials have when evaluated at x = 0:  get unsat cover(s = (x → 0, y → 6)). This time a complete cover is provided from the initial constraint defined by g so no need for the main loop.
get unsat cover(s = (x → 0)) continued. Thus the characterisation is similar to that obtained from the sample y = 0. We actually have to do a little extra work here: because at (x, y) = (0, 6) the leading coefficient of g is zero we include also the trailing coefficient, but this was actually already a factor of the discriminant so we do have the same characterisation as for y = 0. Because this characterisation produced a Z containing y = 6 we cannot generalise beyond this point. Nevertheless, we have a full cover over x = 0.
get unsat cover(s = ()) continued. The cover is formed from five intervals: two from the initial constraints and three from the recursions as summarised below.
The first of these intervals is redundant: the interval (−∞, −10) from h is entirely inside (−∞, 6) from the sample y = 0. Hence it will be removed by the call to construct cover, where we also order the remaining four intervals as in the table above.
We now form the characterisation at x = 0: it contains the discriminants of p 1 , p 2 and p 4 ; those of p 3 and h were constant and so removed. Similarly, all the leading coefficients were constant and so not required. The only within interval resultants required come from I 3 . We take the resultant of p 3 with both p 1 and p 4 as both the latter produce smaller real roots. The only between interval resultant required is that between p 4 and h which protects the overlap of the final two intervals. After the standard simplifications (including factorisation) we have the following characterisation. Hence we can generalise around x = 0 to the interval (−1, 1). Fig. 14a visualises this generalisation of the cover. We see that I 1 , which was originally the y axis below 6, is generalised to the cylinder bounded at the top by p 3 ; while I 2 , which was originally just the point (x, y) = (0, 6), is generalised to the line segment of y − x − 6 that is above x ∈ (−1, 1). Recall that I 2 was a cell where the leading coefficient of g vanished and observe that this remains true throughout the generalisation. I 3 and I 4 were overlapping segments of the y-axis and now form overlapping cylinders. These different cells are indicated by the different direction of shading in the figure.
Let us pause to compare this cylindrical algebraic covering with an actual CAD. A CAD of the projection {p 1 , p 2 , p 3 , p 4 , h} would also form a cylinder over x ∈ (−1, 1). However, this cylinder would then be split into 17 cells according to the 8 line segments in the image. We would have to work over 17 samples to determine UNSAT over the cylinder, while we obtained the same conclusion using only 4 in with the new algorithm.
The algorithm continues in the call to get unsat cover(s = ()) with a sample outside of (−1, 1). A natural choice here is s 1 = 2, for which an UNSAT cover would be found consisting of the following six intervals (after removal of redundant ones): x y (a) The cover over x ∈ (−1, 1) generalised from x = 0.  The characterisation then allows a generalisation from x = 2 to x ∈ (1, 4.75), promoting a third sample in the main call of x = 5. In the recursion a sample of y = 0 is unsat for the first constraint but cannot be generalised. A second sample of y = 1 leads to the initial constraints providing a cover for all but z ∈ (−1, 1). Hence the main loop is entered and we sample z = 0 to find a satisfying witness: (x, y, z) = (5, 1, 0) which evaluates f, g, h to 1, 15 and −99 respectively.
We note that the covering produces far less cells than a traditional CAD. In particular, not every projection factor is used in every covering. This was demonstrated visually by Fig. 14a where the purple resultant was a cell boundary for positive y but not for negative y.

Dataset
The worked examples show that the new algorithm does have advantages over existing techniques. We will now check whether this translates into a method that is also competitive in practice. This section describes experiments on the benchmarks for quantifier-free nonlinear real arithmetic (QF NRA) from the SMT-LIB [28] initiative. This dataset contains 11489 problem instances from 11 different families. Note that these examples are not presented as sets of constraints but can come with more involved Boolean structure, and so the new algorithm is applied as the theory solver for the larger SMT system. Throughout we apply a time limit of 60 seconds for tackling a problem instance.

New Solver CDCAC
We produced an implementation in the course of a Master thesis within the SMT solver SMT-RAT [33]. The solver CDCAC consists of the regular SAT solver from SMT-RAT, combined with this implementation of the presented method as the only theory solver. We consider this implementation to be a proof of concept, as it does not implement any form of incrementality and there is scope to optimise several subroutines. In the interests of reproducibility it is available from the following URL: https://github.com/smtrat/smtrat/tree/JLAMP-CDCAC With regards to the completeness issue discussed in Section 4.4.5: we observed nullification upon substitution of the sample (Algorithm 5 Line 3) in a small, but significant number of examples: 82 from all 11489 instances (0.7%) involved a nullification within the time limit. This means that for these instances the theory does not guarantee UNSAT results of solver CDCAC as correct. However, we note that in all cases the results provided were still correct (i.e. they matched the declared answer stored in the SMT-LIB), which is consistent with our experience of similar CAD completeness issues in [35].

Competing Solvers
We compare the new algorithm to a theory solver based on traditional CAD as described in [35]. We show this in two configurations: CAD-naive uses no incrementality (by projection factor − see Section 2.4) while CAD-full employs what was called full incrementality in [35]. Solver CDCAC uses no incrementality, as already discussed in Section 4.6, but the difference between the incremental and non-incremental versions of the full CAD gives a glimpse at what further improvements to CDCAC might be possible. Both of these use Lazard's projection operator.
We also compare against the NLSAT method [36]. For better comparability of the underlying algorithms, we used our own implementation of the NLSAT approach which uses the same underlying data structures as the aforementioned solvers. NLSAT uses only the CAD-based explanation backend and the NLSAT variable ordering. This uses the same projection operator as CDCAC and so has similar completeness limitations. We call this solver NLSAT and refer to [44] for a more detailed discussion of the implementation and the relation to the original NLSAT method from [36].
Throughout the paper we also highlighted NuCAD as a relevant competing algorithm. This has been implemented in the Tarski system [42]. However, Tarski only implements open NuCAD (full dimensional cells only) which reduces the comparability, and there are also some technical issues at the time of writing with Tarski's parsing of certain SMT-LIB benchmarks. So we do not report on the use of Tarski here.   The overall performance of the solvers is summarised in Fig. 15. For each solver we show the number of solved problem instances (within the 60 second time limit), both separated into satisfiable and unsatisfiable instances and in total. Additionally, average runtimes for the solved instances are given and the final column shows the overall percentage of the dataset that was solved. The table shows that our new algorithm has led to a competitive solver. It significantly outperforms a solver based on regular CAD computation, even when it has been made incremental. Fig. 16 shows how the solvers tackled instances by their runtime, demonstrating that on trivial examples the regular CAD actually outperforms the others, but for more complex examples the savings of the adapted solvers outweigh their overheads.

Solver
We note that although we use the same SAT solver for CAD-naive, CAD-full and CDCAC, the overall SMT solver is not guaranteed to behave the same, i.e. make the same sequence of theory calls. This is because the theory solvers may construct different infeasible subsets which guide the SAT solver in different ways. Hence the difference seen in Fig. 15 may not be purely due to improved performance within the theory solver. Our prior experience with this issue from [45,35] indicates that the different subsets usually have only a negligible influence on the overall solver performance.

Comparison with NLSAT in Experimental Results
We acknowledge that our new algorithm is outperformed by the NLSAT-based solver. We note that this out-performance is greater on SAT instances than UNSAT, suggesting that the new approach may be particularly useful for proving UNSAT. We are optimistic that NLSAT could be beaten following improvements to our implementation. First, implementing incrementality increased the number of problems the CAD solver could tackle by 7% of the dataset, so gains to CDCAC from incrementality may be particularly fruitful. Second, we have yet to perform optimisation on things like the variable ordering choice for CDCAC, which is known to have a great effect on CAD-based technology (see e.g. [46]) and has been already performed for this NLSAT-based solver to give significant improvements as detailed in [44].   A more detailed comparison of CDCAC and NLSAT is presented in Fig. 17. In Fig. 17a we show a direct comparison of the two solvers, where every dot represents one problem instance. The cluster of points in the lower left corner represent the many easy problem instances that are solved by both solvers. However, more interesting is that the solvers perform very differently on a substantial amount of harder problem instances. Note that almost no problem instance is solved by both solvers after more than 20 seconds, but both solvers can solve such problem instances that the other cannot. There are 555 problems on which CDCAC times out but NLSAT completes, and 358 problems for which NLSAT times out and CDCAC completes.
Figs. 17b and 17c separates the number of problem instances according to how much slower or faster CDCAC is when compared to NLSAT. The data in Fig. 17b are for problems where at least one solver took more than 10ms, and the data in Fig. 17c where at least one took above 100ms (so these exclude instances that the solvers agree are trivial). Note that the rows are cumulative (i.e. instances which are 2× faster are also 1.1× faster). Fig. 17b states that NLSAT was actually slower than CDCAC more often than it was faster. This does not contradict the data in Fig. 15 however, as NLSAT still completed in many of the cases it was slower while CDCAC would time out more often when slower.
Altogether, the data from Fig. 17 indicates that CDCAC and NLSAT have substantially different strengths and weaknesses. There hence may be benefit in combining them, for example in a portfolio solver [47].

Conclusions
We have presented a new algorithm for solving SMT problems in non-linear real arithmetic based on the construction of cylindrical algebraic coverings to learn from conflicts. We demonstrated the approach with detailed worked examples and the performance of our implementation on a large dataset.

Comparison with Traditional CAD
There are clear advantages over a more traditional CAD approach (even one that is made incremental for the SMT context). The new algorithm requires fewer projection polynomials, fewer cells, and fewer sample points to determine unsatisfiability in a given cylinder. It can even sometimes avoid the calculation of projection polynomials entirely, when they represent combinations of polynomials not required to determine unsatisfiability. These advantages manifested in substantially more examples solved in the experiments.

Comparison with NuCAD
The worked examples made clear how the new algorithm more effectively learns from conflicts than the NuCAD approach [27] applied to our setting. Although NuCAD also saves in projection and cell construction when compared to traditional CAD, its savings are not maximised as the search is less directed. The conflicts of separate cells in NuCAD are never combined and so the search within NuCAD is not guided to the same extent. However, NuCAD is able to work in a larger QE context, while this new algorithm is applicable only in a restricted setting for checking the consistency of polynomial constraint sets, as outlined in Section 2.3.

Comparison with NLSAT
Our ideas are closest to those of the NLSAT method [36] which likewise builds partial samples and when these cannot be extended explains the conflict locally and generalises from the point. For several of the worked examples the new algorithm and NLSAT performed similarly, with any differences due to luck or heuristic choices not fundamental to the theory. However, the methods are very different in their presentation and computational flow. There are a number of potential advantages of the new approach compared to NLSAT.
First, on a practical note, our method is SMT compliant and thus allows for a (relatively) easy integration with any existing SMT solver. NLSAT on the other hand is usually separate from the regular SMT solving engine which makes it hard to combine with other approaches, e.g. for theory combination. Second, our approach can avoid some projections performed by NLSAT in cases where a change in the sample in one dimension will find a satisfying witness (see the example in Section 5.3).
Most fundamentally, our approach keeps the theory reasoning away from the SAT solving engine while NLSAT burdens the core solving engine with many lemmas for the (intermediate) theory conflicts. For a problem instance that is fully conjunctive our method may determine the satisfiability alone while NLSAT will still require a SAT solver. It is possible that NLSAT may cause the Boolean reasoning to slow down or use excessive amounts of memory by forcing it to deal with a large number of clauses. Further, once we move on to a different area of the solution space most of the previously generated lemmas are usually irrelevant, but their literals may still leak into the theory reasoning in NLSAT. An interesting claim about NLSAT is that explanations for a certain theory assignment can sometimes be reused in another context [48]. While this may be true, our own experiments suggest that this impact is rather limited. Our new algorithm can not reuse unsatisfiable intervals in such a way, although we can for example retain projection factors to avoid their recomputation.
The experimental results had the NLSAT based solver performing best overall, but our new algorithm was already competitive with only a limited implementation 12 , and outperformed NLSAT on a substantial quantity of problems instances as detailed in Fig. 17.

Future Work
Our next step will be to optimise the implementation of the new algorithm, in particular to have it work incrementally. We will also consider how to optimise some of the heuristic choices required for the algorithm. This most notably includes the variable ordering, critical for the performance of all CAD technology 13 . The new algorithm also introduces its own choices which could be made heuristically, such as whether to remove redundant intervals (Section 4.5) and more generally which intervals to select for a covering.
There are interesting questions regarding the extension of the theory. First of all, there is the desire to use Lazard projection theory in place of McCallum, to avoid the completeness issues of Section 4.4.5; but ideally we would avoid the additional trailing coefficients of Lazard where possible. Second, there is the question of the complexity of our new algorithm. Finally, there are questions of more general extensions: Can the approach outlined here be adapted to more general quantifier elimination? Can we go from "merge intervals for existential quantifier" to "merge intervals for universal quantifier"? Can we use cylindrical coverings to study all first order formulae (instead of just sets of constraints representing conjunction) and thus remove the need for the SAT solver entirely?