Improving the use of equational constraints in cylindrical algebraic decomposition
IImproving the use of equational constraints incylindrical algebraic decomposition
Matthew England
Coventry University
Russell Bradford
University of Bath
James H. Davenport
University of Bath
ABSTRACT
When building a cylindrical algebraic decomposition (CAD)savings can be made in the presence of an equational con-straint (EC): an equation logically implied by a formula.The present paper is concerned with how to use multipleECs, propagating those in the input throughout the pro-jection set. We improve on the approach of McCallum inISSAC 2001 by using the reduced projection theory to makesavings in the lifting phase (both to the polynomials we liftwith and the cells lifted over). We demonstrate the benefitswith worked examples and a complexity analysis.
Categories and Subject Descriptors
I.1.2 [
Symbolic and Algebraic Manipulation ]: Algo-rithms—
Algebraic algorithms, Analysis of algorithms
General Terms
Algorithms, Experimentation, Theory
Keywords cylindrical algebraic decomposition, equational constraint
1. INTRODUCTION A cylindrical algebraic decomposition (CAD) splits R n intocells arranged cylindrically , meaning the projections of anypair are either equal or disjoint, and such that each can bedescribed with a finite sequence of polynomial constraints.Introduced by Collins for quantifier elimination in realclosed fields, applications of CAD include: derivation of op-timal numerical schemes [18], parametric optimisation [19],epidemic modelling [9], theorem proving [27], reasoning withmulti-valued functions [13], and much more.CAD has complexity doubly exponential in the number ofvariables [14]. For some applications there exist algorithmswith better complexity (see [2]), but CAD implementationsremain the best general purpose approach for many. This Permission to make digital or hard copies of part or all of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor profit or commercial advantage, and that copies bear this notice and the full ci-tation on the first page. Copyrights for third-party components of this work must behonored. For all other uses, contact the owner/author(s). Copyright is held by theauthor/owner(s).
ISSAC’15,
July 6–9, 2015, Bath, United Kingdom.ACM 978-1-4503-3435-8/15/07.DOI: http://dx.doi.org/10.1145/2755996.2756678. may be due to the many extensions and optimisations ofCAD since Collins including: partial CAD (to lift only whennecessary for quantifier elimination); symbolic-numeric lift-ing schemes [29, 22]; local projection approaches [8, 30]; anddecompositions via complex space [11, 3]. Collins originalalgorithm is described in [1] while a more detailed summaryof recent developments can be found, for example, in [5].
We describe the computation scheme and terminology thatmost CAD algorithms share. We assume a set of input poly-nomials (possibly derived from formulae) in ordered vari-ables x = x ≺ . . . ≺ x n . The main variable of a polynomial(mvar) is the greatest variable present under the ordering.The first phase of CAD, projection , applies projection op-erators repeatedly, each time producing another set of poly-nomials in one fewer variables. Together these contain the projection polynomials used in the second phase, lifting , tobuild the CAD incrementally. First R is decomposed intocells which are points and intervals according to the realroots of polynomials univariate in x . Then R is decom-posed by repeating the process over each cell with the bi-variate polynomials in ( x , x ) evaluated at a sample point.This produces sections (where a polynomial vanishes) and sectors (the regions between) which together form the stack over the cell. Taking the union of these stacks gives the CADof R and this is repeated until a CAD of R n is produced.At each stage cells are represented by (at least) a samplepoint and an index . The latter is a list of integers, withthe k th describing variable x k according to the ordered realroots of the projection polynomials in ( x , . . . , x k ). If theinteger is 2 i the cell is over the i th root (counting from lowto high) and if 2 i + 1 over the interval between the i th and( i + 1)th (or the unbounded intervals at either end).The projection operator is chosen so polynomials are de-lineable in a cell: the portion of their zero set in the cellconsists of disjoint sections. A set of polynomials are de-lineable if each is individually, and the sections of differentpolynomials are identical or disjoint. If all projection poly-nomials are delineable then the input polynomials must be sign-invariant : have constant sign in each cell of the CAD. Most applications of CAD require truth-invariance for log-ical formulae, meaning each formula has constant booleantruth value on each cell. Sign-invariance for the polynomi-als in a formula gives truth-invariance, but we can obtainthe latter more efficiently by using equational constraints. a r X i v : . [ c s . S C ] M a y efinition We use
QFF to denote a quantifier freeTarski formula: Boolean combinations ( ∧ , ∨ , ¬ ) of state-ments about the signs ( = 0 , > , < ) of integral polynomials.An equational constraint (EC) is a polynomial equationlogically implied by a QFF. If an atom of the formula it issaid to be explicit and is otherwise implicit . Collins first suggested that the projection phase of CADcould be simplified in the presence of an EC [12]. He notedthat a CAD sign-invariant for the defining polynomial of anEC, and sign-invariant for any others only on sections of thatpolynomial, would be sufficient. An intuitive approach toproduce this is to consider resultants of the EC polynomialwith the other polynomials, in place of them. This approachwas first formalised and verified in [24].A recent complexity analysis [5] showed that using an ECin this way reduces the double exponent in the complexitybound for CAD by 1. A natural question is whether this canbe repeated in the presence of multiple ECs. An algorithmfor CAD in the presence of two ECs was detailed in [25].The main idea was to observe that the resultant of the poly-nomials defining two ECs is itself an EC, and so the sameideas could be applied for the second projection as for thefirst. However, this approach was complicated as the keyresult verifying [24] could not be applied recursively.
This paper discusses how we can extend the theory of ECsto produce CADs more efficiently. In Section 2.1 we revisekey components of the theory for reduced projection in thepresence of an EC from [24, 25]. Then in Section 2.2 weexplain how it can also give reductions in the lifting phase,allowing us to propose and verify a new algorithm in Section3 for making use of multiple ECs. This breaks with the tra-dition of producing CADs sign-invariant for EC polynomials,instead guaranteeing only invariance for the truth of theirconjunction. We demonstrate our contributions in Sections4 and 5 with a worked example and complexity analysis.All experiments in
Maple were conducted using
Maple
18. All code and data created for this paper is openly avail-able from http://dx.doi.org/10.15125/BATH-00071 .
2. CAD WITH MULTIPLE EQUATIONALCONSTRAINTS2.1 Key theory from [23, 24, 25]
We recall some of the key theory behind McCallum’s oper-ators. Let cont, prim, disc, coeff and ldcf denote the content,primitive part, discriminant, coefficients and leading coeffi-cient of polynomials respectively (in each case taken withrespect to a given mvar). Let res denote the resultant of apair of polynomials. When applied to a set of polynomialswe interpret these as producing sets of polynomials, e.g.res( A ) = { res( f i , f j ) | f i ∈ A, f j ∈ A, f j (cid:54) = f i } . Recall that a set A ⊂ Z [ x ] is an irreducible basis if the el-ements of A are of positive degree in the mvar, irreducibleand pairwise relatively prime. Throughout this section sup-pose B is an irreducible basis for a set of polynomials, thatevery element of B has mvar x n and that F ⊆ B . Define P ( B ) := coeff( B ) ∪ disc( B ) ∪ res( B ) , (1) P F ( B ) := P ( F ) ∪ { res( f, g ) | f ∈ F, g ∈ B \ F } , (2) P ∗ F ( B ) := P F ( B ) ∪ res( B \ F ) , (3)as the projection operators introduced respectively in [23,24, 25]. In the general case with A a set of polynomialsand E ⊆ A we proceed with projection by: letting B and F be irreducible basis of the primitive parts of A and E respectively; applying the operators as defined above; andthen taking the union of the output with cont( A ).The theorems in this section validate the use of these oper-ators for CAD. They use the condition of order-invariance ,meaning each polynomial has constant order of vanishingwithin each cell, which of course implies sign-invariance. Wesay that a polynomial with mvar x k is nullified over a cellin R k − if it vanishes identically throughout. Theorem 1 ([23]).
Let S be a connected submanifoldof R n − in which each element of P ( B ) is order-invariant.Then on S , each element of B is either nullified or analyticdelineable (a variant on delineability, see [23]). Further, thesections of B not nullified are pairwise disjoint, and eachelement of such B is order-invariant on such sections. Suppose we apply P repeatedly to generate projection poly-nomials. Repeated use of Theorem 1 concludes that a CADproduced by lifting with respect to these projection polyno-mial is order-invariant so long as no projection polynomialwith mvar x k is nullified over a cell in the CAD of R k − (acondition known as well-orientedness which can be checkedduring lifting). If this condition is not satisfied then P can-not be used (and we should restart the CAD constructionusing a different projection operator, such as Hong’s [20]). Theorem 2 ([24]).
Let f and g be integral polynomialswith mvar x n , r ( x , . . . , x n − ) be their resultant, and sup-pose r (cid:54) = 0 . Let S be a connected subset of R n − on which f is delineable and r order-invariant.Then g is sign-invariant in every section of f over S . Suppose A was derived from a formula with EC defined by E = { f } , and that we apply P E ( A ) once and then P repeat-edly to generate a set of projection polynomials. Assumingthe input is well-oriented, we can use Theorem 1 to concludethe CAD of R n − order invariant for P E ( A ). The CAD of R n is then sign-invariant for E using Theorem 1 and sign-invariant for A in the sections of E using Theorem 2. Hencethe CAD is truth-invariant for the formula.What if there are multiple ECs? We could designate onefor special use and treat the rest as any other constraint(heuristics can help with the choice [6]). But this does notgain any more advantage than one EC gives. However, wecannot simply add multiple polynomials into E at the toplevel as this would result in a CAD truth-invariant for thedisjunction of the ECs, not the conjunction.Suppose we have a formula with a second EC. If this hasa lower mvar then we may consider applying the reducedprojection operator again at this lower level. In fact, evenif the second EC is also in the mvar of the system we can propagate it to the lower level by noting that the resultantof the two ECs is itself an EC in one fewer variable.So we consider applying first the operator P E ( A ) where E defines the first EC and then P E (cid:48) ( A (cid:48) ) where A (cid:48) = P E ( A )and E (cid:48) ⊆ A (cid:48) contains the EC in one variable fewer. Unfor-tunately, Theorem 2 does not validate this approach. Whileit could be applied once for the CAD of R n − it cannot then igure 1: The polynomials from Example 1. validate the CAD of R n because the first application of thetheorem provided sign-invariance while the second requiresthe stronger condition of order invariance. Note however,that this approach is acceptable if n = 3 (since in two vari-ables the conditions are equivalent for squarefree bases). Example The following are graphed in Figure 1, with g the sphere, f the upper surface and f the lower: f = x + y + z, f = x − y + z, g = x + y + z − . We consider the formula φ = f = 0 ∧ f = 0 ∧ g ≥ . Thesurfaces f and f only meet on the plane y = 0 and thisprojection is on the right of Figure 1). From this it is clearthe solution requires | x | ≥ √ / and z = − x .How could this be ascertained using CAD? With variableordering z (cid:31) y (cid:31) x a sign-invariant CAD for ( f , f , g ) has1487 cells using Qepcad [7]. We could then test a samplepoint of each cell to identify the ones where φ is true.It is preferable to use the presence of ECs. Declaring anEC to Qepcad will ensure it uses the algorithm in [24] basedon a single use of P E ( A ) followed by P . Either choice re-sults in 289 cells. In particular, the solution set is describedusing 8 cells: all have y = 0 , z = − x but the x -coordinateunnecessarily splits cells at (1 ± √ . This is identified dueto the projection polynomial d = disc y (res z ( f i , g )) .If we declare both ECs to Qepcad then it will use thealgorithm in [25] applying P E ( A ) twice (allowed since n = 3 )to produce a CAD with 133 cells. The solution set is nowdescribed using only 4 cells (the minimum possible). Notethat d was no longer produced as a projection polynomial. For problems with n > P ∗ E ( A ) is appropriate. Theorem 3 ([25]).
Let f and g be integral polynomialswith mvar x n , r = res( f, g ) , d = disc ( g ) , and suppose r, d (cid:54) =0 . Let S be a connected subset of R n − on which f is analyticdelineable, g is not nullified and r and d are order-invariant.Then g is order-invariant in each section of f over S . Suppose we have a formula with two ECs, one with mvar x n and the other with mvar x n − . The second could be ex-plicit in the formula or implicit (a resultant as described ear-lier). Theorem 3 allows us to use a reduced operator twice.We first calculate A (cid:48) = P E ( A ) where E contains the defin-ing polynomial of the first EC, and then P ∗ E (cid:48) ( A (cid:48) ) where E (cid:48) contains the defining polynomial of the other. Subsequentprojections simply use P . When lifting we use Theorem 1 to verify the CAD of R n − as order-invariant for P ∗ E (cid:48) ( A (cid:48) );Theorem 1 to verify the CAD of R n − order-invariant for E (cid:48) everywhere and Theorem 3 to verify it order-invariantfor A (cid:48) in the sections of E (cid:48) ; and Theorem 1 and 2 to verifythe CAD of R n order-invariant for E and sign-invariant for A in those cells that are both sections of E and E (cid:48) . The main contribution of the present paper is to realisethat the theorems above also allow for significant savings inthe lifting phase of CAD. However, to implement these wemust discard two embedded principles of CAD:1. That the projection polynomials are a fixed set.2. That the invariance structure of the final CAD can beexpressed in terms of sign-invariance of polynomials.Abandoning the first is key to recent work in [11, 3], whilethe second was also investigated in [10, 26].
Consider Theorem 2: it allows us to conclude that g issign-invariant in the sections of f produced over a CADof R n − order-invariant for P { f } ( { f, g } ). Therefore, it issufficient to perform the final lift with respect to f only(decompose cylinders according to the roots of f but not g ). The decomposition imposes sign-invariance for f whileTheorem 2 guarantees it for g in the cells where it matters. Example We return to Example 1. Recall that desig-nating either EC and using [24] produced a CAD with 289cells. If we follow this approach but lift only with respectto the designated EC at the final step (implemented in our
Maple package [17]) we obtain a CAD with 141 cells.
This improved lifting follows from the theorems in [24],but was only noticed 15 years later during the generalisa-tion of [24] to the case of multiple formulae in [4, 5]. Exper-iments there demonstrated its importance, particularly forproblems with many constraints (see Section 8.3 of [5]).When we apply a reduced operator at two levels then wecan make such reductions at both the corresponding lifts.
Example We return to the problem from Example 1.Set A = { f , f , g } and E = { f } . Then project out z using P E ( A ) = { y , y + 2 xy + 2 x + y − } . These are the resultants of f with f and g . The discrimi-nant of f was a constant and so could be discarded, as wasits leading coefficient (meaning no further coefficients wererequired). We set A (cid:48) = P E ( A ) , E (cid:48) = res z ( f , f ) = y and R = res y ( y , y + 2 xy + 2 x + y −
1) = (2 x − . We have P E (cid:48) ( A (cid:48) ) = { R } since the other possible entries (thediscriminants and coefficients from E (cid:48) ) are all constants.We hence build a 5 cell CAD of the real line with respect tothe two real roots of R . We then lift above each cell withrespect to y only, in each case splitting the cylinder intothree cells about y = 0 , to give a CAD of R with 15 cells.Finally, we lift over each of these 15 cells with respectto f to give 45 cells of R . This compares to 133 from Qepcad , which used reduced projection but then lifted withall projection polynomials. No polynomials were nullified, sousing Theorems 1 and 2, the output is truth-invariant for φ . he additional lifting that Qepcad performed does notprovide any further structure. For example, if we had liftedwith respect to f at the final stage in Example 3 then wewould be doing so without the knowledge that it is deline-able. Hence splitting the cylinder at the sample point offersno guarantee that the cells produced are sign-invariant awayfrom that point. So the extra work does not allow us to con-clude that f is sign-invariant (except on sections of f ).Note that using fewer projection polynomials for liftingnot only decreases output size (and computation time) butalso the risk of failure from non well-oriented input: we onlyneed worry about nullification of polynomials we lift with. We can achieve still more savings from the theory in Sec-tion 2.1 by abandoning the aim of producing a CAD sign-invariant with respect to any polynomial, instead insistingonly on truth-invariance for the formula. We may then lifttrivially to cylinders over cells already known to be false,only identifying sections of projection polynomials if thereis a possibility the formula may be true. The idea of avoidingcomputations over false cells was presented in [28]. Our con-tribution is to explain how such cells can easily be identifiedin the presence of ECs. We demonstrate with our example.
Example Return to the problem from Examples 1 − R produced with 15 cells inExample 3. On 5 of these 15 cells the polynomial R is zeroand on the others it is either positive or negative throughout.Now, φ can only be satisfied above the 5 cells, as elsewherethe two EC defining polynomials cannot share a root and thusvanish together. We can already conclude the truth value forthe 10 cells (false) and thus we do not need to lift over them,except in the trivial sense of extending them to a cylinder in R . Lifting over the 5 cells where R = 0 with respect to f gives 15 cells, which combined with the 10 cylinders gives aCAD of R with 25 cells that is truth-invariant for φ . This 25 cell CAD is not sign-invariant for f . The cylin-ders above the 10 cells in R where R (cid:54) = 0 may have f varying sign, but since f can never equal zero at the sametime as f in these cells it does not affect the truth of φ .Identifying the 5 cells where R = 0 in the CAD of R was trivial since they are simply the sections of the secondlift, and hence those cells with second entry even in the cellindex. Those sections produced in the third lift are similarlyall cells where f is zero, however, we cannot conclude that f is also zero on these. Theorem 2 only guarantees that f is sign-invariant on such cells, so to determine those signswe must still evaluate the polynomials at the sample point.Reducing the number of cells for stack generation clearlydecreases output size, and since the cells can be identified us-ing only a parity check on an integer, computation time de-creases also. As with the improvements in Section 2.2.1, thisalso decreases the risk of non well-oriented input: we onlyneed worry about nullification over these identified cells.
3. ALGORITHM
We present Algorithm 1 to build a truth-invariant CADfor a formula in the presence of multiple ECs. We assumethat the ECs are already identified as input to the algorithm(they may have been first computed through propagation asdescribed in Section 2). We assume further that each EC
Algorithm 1:
CAD using multiple ECs
Input : A formula φ in variables x , . . . , x n , and asequence of sets { E k } nk =1 ; each either empty orcontaining a single primitive polynomial withmvar x k which defines an EC for φ . Output : Either: D , a truth-invariant CAD of R n for φ (described by lists I and S of cell indices andsample points); or FAIL , if not well-oriented. Extract from φ the set of defining polynomials A n ; for k = n, . . . , do Set B k to the finest squarefree basis for prim( A k ); Set C to cont( A k ); Set F k to the finest squarefree basis for E k ; if F k is empty then Set A k − := C ∪ P ( B k ); else if k = n or k = 2 then Set A k − := C ∪ P F i ( B i ); else Set A k − := C ∪ P ∗ F i ( B i ); If E is not empty then set p to be its element;otherwise set p to the product of polynomials in A ; Build D := ( I , S ) according to the real roots of p ; if n = 1 then return D ; for k = 2 , . . . , n do Initialise D k = ( I k , S k ) with I k and S k empty sets; if F k is empty then Set L := B k ; else Set L := F k ; if E k − is empty then Set C a := D k − and C b empty; else Set C a to be cells in D k − with I k − [ −
1] even; Set C b := D k − \ C a ; for each cell c ∈ C a do if An element of L is nullified over c then return FAIL; Generate a stack over c with respect to thepolynomials in L , adding cell indices and samplepoints to I k and S k ; for each cell c ∈ C b do Extend to a single cell in R k (cylinder over c ),adding index and sample point to I k and S k ; return D n = ( I n , S n ).is primitive, and that all the ECs have different mvar (so inpractice a choice of designation may have been made).Steps 1 −
12 run the projection phase of the algorithm.Each projection starts by identifying contents and primitiveparts. When there is no declared EC ( E i is empty) the pro-jection operator (1) is used (step 7). Otherwise the operator(3) is used (step 12), unless it is the very first or very lastprojection (step 10) when we use (2). In each case the out-put of the projection operator is combined with the contentsto form the next layer of projection polynomials.teps 13 −
16 construct a CAD for the real line (andreturn it if the input was univariate). This is sometimesreferred to in the literature as the base phase . If there isa declared EC in the smallest variable then the real line isdecomposed according to its roots, otherwise according tothe roots of all the univariate projection polynomials.Steps 17 −
33 run the lifting phase, incrementally buildingCADs of R k for k = 2 , . . . , n . For each k there are twoconsiderations. First, whether there is a declared EC withmvar x k . If so we lift only with respect to this (step 22)and if not we use all projection polynomials with mvar x k (step 20). Second, whether there is a declared EC with mvar x k − . If so we restrict stack generation to those cells wherethe EC was satisfied. These are simply those with I k − [ − Definition For k = 2 , . . . , n define sets: • L k − the lifting polynomials : the defining polynomialof the declared EC with mvar x k if one exists, or allprojection polynomials with mvar x k otherwise. • C k − the lifting cells : those cells in the CAD of R k − in which the designated EC with mvar x k − vanishesif it exists, and all cells in that CAD otherwise.The input of Algorithm 1 is well-oriented if for k = 2 , . . . , n no element of L k is nullified over an element of C k . Theorem Algorithm 1 satisfies its specification.
Proof.
We must show the CAD is truth-invariant for φ ,unless the input is not well-oriented when FAIL is returned.First consider the case where n = 1. The projection phasewould not run, with the algorithm jumping to the CADconstruction in step 13, returning the output in step 16. Ifthere was no declared EC then the CAD is sign-invariantfor all polynomials defining φ and thus every cell is truthinvariant for φ . If there was a declared EC then the outputis sign-invariant for its defining polynomial. Cells wouldeither be intervals where the formula must be false; or points,where the EC is satisfied, and the formula either identicallytrue or false depending on the signs of the other polynomials.Next suppose that the input were not well-oriented (Def-inition 2). For a fixed k , the conditional in steps 19 − L k to L and the conditional insteps 23 −
27 the lifting cells C k to C a . Thus it is exactlythe conditions of Definition 2 which are checked by step 29,returning FAIL in step 30 when they are not satisfied. If thelifting phase completes then the input is well-oriented.From now on we suppose n > k define admissible cells to be thosein the induced CAD of R k − where all declared ECs withmvar smaller than x k are satisfied, or to be all cells in thatinduced CAD if there are no such ECs. Then let I ( k ) be thefollowing statement for the CADs produced by Algorithm 1. Over admissible cells (in R k − ) the CAD of R k is:(a) order-invariant for any EC with mvar x k ;(b) order- (sign- if k = n ) invariant for all projection poly-nomials with mvar x k on sections of the EC over admis-sible cells, or over all admissible cells if no EC exists. We have already proved I (1), and I ( n ) may be proved byinduction. To assert the truth of I ( k ) we note the following: • When E k is empty we use Theorem 1 to assert all pro-jection polynomials with mvar x k are order-invariantin the stacks over admissible cells giving (a) and (b). • When E k is not empty and k = 2 we used the projec-tion operator (2). Theorem 2 allows us to conclude(b) and that the EC is sign-invariant in admissiblecells. The stronger property of order-invariance fol-lows automatically since the lifting polynomials forma squarefree basis in two variables. • When E k is not empty and k = n we used the pro-jection operator (2). Theorem 2 allows us to conclude(b), but also (a) since in the case k = n the statementrequires only sign-invariance. • When E k is not empty and 2 < k < n we used theprojection operator (3). Theorem 3 then allows us toconclude the statement.In each case the assumptions of the theorems are met by theinductive hypothesis, exactly over admissible cells as definedaccording to whether E k − was empty or not.From the definition of admissible cells, we know that φ is false (and thus trivially truth invariant) upon all cellsin the CAD of R n built over an inadmissible cell of R k , k < n . Coupled with the truth of (a) for k = 1 , . . . , n , thisimplies the CAD of R n is truth-invariant for the conjunctionof ECs (although it may not be truth-invariant for any oneindividually). The truth of (b) implies that on those cellswhere all ECs are satisfied, the other polynomials in φ aresign-invariant and thus φ is truth-invariant.
4. WORKED EXAMPLE
Assume variable ordering z (cid:31) y (cid:31) x (cid:31) u (cid:31) v and define f := x − y + z , f := z − u + v − , g := x − ,f := x + y + z , f := z + u − v − , h := z. We consider the formula φ = f = 0 ∧ f = 0 ∧ f = 0 ∧ f = 0 ∧ g ≥ ∧ h ≥ . The solution can be found manually by decomposing thesystem into blocks. The surfaces f and f are graphed in( x, y, z )-space on the left of Figure 2. They meet only onthe plane y = 0 and this projection is shown on the right.The surfaces f and f are graphed in ( z, u, v )-space on theleft of Figure 3 and meet only when z = ±
1. We consideronly z = +1 due to h ≥
0, with this projection plotted onthe right. We thus see that the solution set is given by { u = ± v, x = − , y = 0 , z = 1 } . To ascertain this by Algorithm 1 we must first propagateand designate ECs. We choose to use f first, calculateres z ( f , f ) = ( − u + v − x + y − and assign r to be the square root: the defining polynomialfor an EC with mvar y . Similarly considerres y (cid:0) r , res z ( f , f ) (cid:1) = 16( u − v + x + 1) , res y (cid:0) r , res z ( f , f ) (cid:1) = 4( u − v ) and assign r := u − v + x + 1, r := u − v : definingpolynomials for ECs with mvar x and u respectively. Thereis no series of resultants that leads to an EC with mvar u . We hence have { E j } nk =1 := { f } , { r } , { r } , { r } , { } asinput for Algorithm 1, along with φ . igure 2: The polynomials f and f from Section 4.Figure 3: The polynomials f and f from Section 4. The algorithm starts by extracting the defining polynomi-als A = { f , f , f , f , g, h } and finds B = A , F = E (infact F i = E i for all i = 1 , . . . , A := P F ( B ) = { ( x − , ( − u + v − x + y − , ( u − v − x + y − , y , x − y } . Hence C := { x − } and B := { y, y − x, − u + v − x + y − , u − v − x + y − } . For the next projection we must use operator (3), giving A := C ∪ P ∗ F ( B ) = { x − , u − v + x +1 , u − v , u − v +1 } noting that for this example the extra discriminants in (3)all evaluated to constants and so could be discarded. Then B := { x − , u − v + x + 1 } , C := { u − v , u − v + 1 } , and the next projection also uses (3) to produce A := { u − v , u − v + 1 , u − u v + v + 2 u − v } . For the final projection there is no EC and so we use operator(1) to find A := { v } . The base phase of the algorithmhence produces a 3-cell CAD of the real line isolating 0.For the first lift we have L = { u − v } and C a containingall 3 cells. Above the two intervals we split into 5 cells by thecurves u = ± v , while above v = 0 we split into three cellsabout the origin. From these 13 cells of R we select the 5which were sections of u − v for C a . These are lifted withrespect to L = { r } , and the other 8 are simply extended tocylinders in R . Together this gives a CAD of R with 23cells. The next two lifts are similar, producing first a CADof R with 53 cells and finally a CAD of R with 113 cells.The entire calculation takes less than a second in Maple . Choice in EC designation
Algorithm 1 could have been initialised with alternative ECdesignations. There were the 4 explicit ECs with mvar z , and by taking repeated resultants we discover the followingimplicit ECs, organised in sets with decreasing mvar: { y , u − v + x − y + 1 , − u + v + x − y + 1 ,u − v + x + y + 1 , − u + v + x + y + 1 } , { x + 1 , − u + v + x + 1 , u − v + x + 1 } , { u − v } . There are hence 60 possible permutations of EC designation,but they lead to only 3 different outputs, with 113, 103 and93 cells. Heuristics for other questions of CAD problemformulation [15, 6, 21, 31] could likely be adapted to assisthere. We note that 93 cells is not a minimal truth invariantCAD for φ as it splits the CAD of R at v = 0 (identifiedfrom the discriminant of the only EC with mvar u ). Comparison with other CAD implementations
A sign-invariant CAD of R for the 6 polynomials in theexample could be produced by Qepcad with 1,118,205 cells.Neither the
RegularChains
Library in
Maple [11] nor our
Maple package [17] could produce one in under an hour.Our implementation of [24], which uses operator (2) oncebut also performs the final lift with respect to the EC only,can produce a CAD with either 3023, 10935 or 48299 (twice)cells depending on which EC is designated. The
Qepcad implementation of [24] gives 11961, 30233, 158475, or 158451cells. Comparing these sets of figures we see the dramaticimprovements from just a single reduced lift.Allowing
Qepcad to propagate the 4 ECs (so a similarprojection phase as Algorithm 1 but then a normal CADlifting phase) produces a CAD with 21079 cells. By declar-ing only a subset of the 4 (which presumably changes thedesignations of implicit ECs) a CAD with 5633 cells can beproduced, still much more than using Algorithm 1.The
RegularChains
Library can also make use of multipleECs, as detailed in [3]. The version in
Maple
18 times outafter an hour, however, with the development version a CADcan be produced instantly. There are choices (with analogiesto designation [16]) but they all lead to a 137 cell output.In particular, they all have an induced CAD of the real linewhich splits at v = ± v = 0.
5. COMPLEXITY ANALYSIS
We build on recent work in [5] to measure the dominantterm in bounds on the number of CAD cells produced. Nu-merous studies have shown this to be closely correlated tothe computation time. We assume input with m polynomi-als of maximum degree d in any one of n variables. Definition Consider a set of polynomials p j . The combined degree of the set is the maximum degree (takenwith respect to each variable) of the product of all the poly-nomials in the set: max i (deg x i ( (cid:81) j p j )) . The set has the (m,d)-property if it may be partitionedinto m subsets, each with maximum combined degree d . For example, { y − x, y + 1 } has combined degree 4 andthus the (1 , , A has the ( m, d )-property then P ( A ) ∪ cont( A ) has the( M, d )-property with M = (cid:4) ( m + 1) (cid:5) . When m > M by m (but we need 2 m to cover m = 1).f A has the ( m, d )-property then so does its squarefreebasis. Hence applying this result recursively (as in Table 1)measures the growth in ( m, d )-property during projectionunder operator (1). After the first projection there are mul-tiple polynomials and so the tighter bound for M is used.The number of real roots in a set with the ( m, d )-propertyis at most md . The number of cells in the CAD of R is thusbounded by twice the product of the final two entries, plus1. Similarly, the total number of cells in the CAD of R n by(2 md + 1) (cid:81) n − r =1 (cid:104) (cid:16) r − m r (cid:17) (cid:16) r − d r (cid:17) + 1 (cid:105) . (4)Omitting the +1s will leave us with the dominant term ofthe bound, which evaluates to give the following result. Theorem The dominant term in the bound on thenumber of CAD cells in R n produced using (1) is (2 d ) n − m n − n − − . (5)From now on assume (cid:96) ECs, 0 < (cid:96) ≤ min( m, n ), all withdifferent mvar. For simplicity we assume these variables are x n , . . . , x n − (cid:96) +1 (the first (cid:96) projections are reduced). Lemma Suppose A is a set with the ( m, d ) -propertyand E ⊂ A has the (1 , d ) -property. Then cont( A ) ∪ P ∗ E ( A ) has the (2 m, d ) -property. Proof.
In [5] we proved that applying P E ( A ) ∪ cont( A )gives a set of (cid:4) (3 m + 1) (cid:5) polynomials of combined degree2 d . The extra m − d ( d − (cid:100) ( m − (cid:101) sets of combined degree at most 2 d . Then (cid:4) (3 m + 1) (cid:5) + (cid:100) ( m − (cid:101) = m + (cid:4) ( m + 1) (cid:5) + (cid:98) m (cid:99) and since m ∈ Z this always equals 2 m .We apply this recursively in the top half of Table 2, withthe bottom derived via the process for P , as in Table 1.Define d i and m i as the entries in the Number and Degreecolumns of Table 2 from the row with i Variables. We canbound the number of real roots of projection polynomialsin i variables by m i d i . If we lifted with respect to all theseprojection polynomials, the cell count would be bounded by (cid:81) ni =1 [2 m i d i + 1] = (cid:81) (cid:96)s =0 (cid:104) (cid:16) s m s − d s (cid:17) + 1 (cid:105) · (cid:81) n − (cid:96) − r =1 (cid:104) (cid:16) r (cid:96) m r (cid:96) + r − d (cid:96) + r (cid:17) + 1 (cid:105) . (6)Omitting the +1 from each product allows us to calculatethe dominant term of the bound explicitly as(2 d ) n − m n − (cid:96) + (cid:96) − (cid:96) n − (cid:96) + (cid:96) ( (cid:96) − / . (7)Now we consider the benefit of improved lifting. Start byconsidering the CAD of R n − ( (cid:96) +1) . There can be no reducedlifting until this point and so the cell count bound is givenby the second product in (6), which we will denote by † .The lift to R n − (cid:96) will involve stack generation over all cells,but only with respect to the EC. This can have at most d n − (cid:96) real roots and so the CAD at most [2 d n − (cid:96) + 1]( † ) cells.The next lift, to R n − (cid:96) − , will lift the sections with respectto the EC, and the sectors only trivially (to produce thesame number of cylinders). Hence the cell count bound is[2 d n − ( (cid:96) − + 1] d n − (cid:96) ( † ) + ( d n − (cid:96) + 1)( † ) with dominant term2 d n − ( (cid:96) − d n − (cid:96) ( † ). Subsequent lifts follow the same patternand so 2 d n d n − . . . d n − ( (cid:96) − d n − (cid:96) ( † ) is the dominant term inthe bound for R n . This evaluates to give the following result. Table 1: Projection under operator (1).
Variables Number Degree n m dn − m d n − m d ... ... ... n − r r − m r r − d r ... ... ...1 2 n − m n − n − − d n − Table 2: Projection with (3) (cid:96) times and then (1).
Variables Number Degree n m dn − m d ... ... ... n − (cid:96) (cid:96) m (cid:96) − d (cid:96) n − ( (cid:96) + 1) 2 (cid:96) m (cid:96) +1 − d (cid:96) +1 ... ... ... n − ( (cid:96) + r ) 2 r (cid:96) m r (cid:96) + r − d (cid:96) + r ... ... ...1 2 ( n − − (cid:96) ) (cid:96) m n − − (cid:96) n − − d n − Theorem Consider the CAD of R n produced using Al-gorithm 1 in the presence of ECs in the top (cid:96) variables. Thedominant term in the bound on the number of cells is (cid:81) (cid:96)s =0 (cid:104) s − d s (cid:105) (cid:81) n − (cid:96) − r =1 (cid:104) (cid:16) r (cid:96) m r (cid:96) + r − d (cid:96) + r (cid:17)(cid:105) = (2 d ) n − m n − (cid:96) − (cid:96) n − (cid:96) − (cid:96) . (8)The bound in Theorem 7 is strictly less than the one inTheorem 5. The double exponent of m has decreased by thenumber of ECs; the result of the improved projection in (7).Improved lifting reduced the single exponents further still.
6. CONCLUSIONS AND FUTURE WORK
We have explained how the existing theory for CAD pro-jection using ECs can also be leveraged for significant sav-ings in the lifting phase. We can reduce both the projectionpolynomials used for lifting and the cells over which stacksare generated. We have formalised these ideas in Algorithm1, verified their use in Theorem 4, and demonstrated thebenefit with a worked example and complexity analysis.A key question is how to best deal with non-primitiveECs? Consider φ := zy = 0 ∧ ϕ. Under ordering · · · (cid:31) z (cid:31) y (cid:31) . . . the EC zy = 0 is not primitive, so Algorithm1 cannot use it. We may be tempted to take E = { z } asthe primitive part, project with operator (2) and includethe content y in the first projection. The CAD of ( y, . . . )-space would be sign-invariant for y and thus the CAD of( z, y, . . . )-space truth invariant for the EC (over admissiblecells). But we can no longer say only sections are admissiblefor the next lift as there may be cells with z (cid:54) = 0 and y = 0.We could instead lift over all cells. Alternatively we mightrewrite φ as φ := ( z = 0 ∧ ϕ ) ∨ ( y = 0 ∧ ϕ ) , so each clause hasits own EC. The theory of truth-table invariant CADs [4, 5]is designed to deal with such input, but would require itsown extension to use beyond the first projection. Of course,this extension would also be valuable in its own right. cknowledgements Thanks to the the referees for their helpful comments. Thiswork was supported by EPSRC grant: EP/J003247/1.
7. REFERENCES [1] D. Arnon, G.E. Collins, and S. McCallum. Cylindricalalgebraic decomposition I: The basic algorithm.
SIAMJ. of Computing , 13:865–877, 1984.[2] S. Basu, R. Pollack, and M.F. Roy. Algorithms in RealAlgebraic Geometry. (Volume 10 of Algorithms andComputations in Mathematics). Springer-Verlag, 2006.[3] R. Bradford, C. Chen, J.H. Davenport, M. England,M. Moreno Maza, and D. Wilson. Truth tableinvariant cylindrical algebraic decomposition byregular chains. In
Computer Algebra in ScientificComputing (LNCS 8660), pages 44–58. Springer, 2014.[4] R. Bradford, J.H. Davenport, M. England,S. McCallum, and D. Wilson. Cylindrical algebraicdecompositions for boolean combinations. In
Proc.ISSAC ’13 , pages 125–132. ACM, 2013.[5] R. Bradford, J.H. Davenport, M. England,S. McCallum, and D. Wilson. Truth table invariantcylindrical algebraic decomposition.
Submitted forPublication . Preprint: arXiv:1401.0645 , 2015.[6] R. Bradford, J.H. Davenport, M. England, andD. Wilson. Optimising problem formulations forcylindrical algebraic decomposition. In
IntelligentComputer Mathematics (LNCS 7961), pages 19–34.Springer Berlin Heidelberg, 2013.[7] C.W. Brown. QEPCAD B: A program for computingwith semi-algebraic sets using CADs.
ACM SIGSAMBulletin , 37(4):97–108, 2003.[8] C.W. Brown. Constructing a single open cell in acylindrical algebraic decomposition. In
Proc. ISSAC’13 , pages 133–140. ACM, 2013.[9] C.W. Brown, M. El Kahoui, D. Novotni, andA. Weber. Algorithmic methods for investigatingequilibria in epidemic modelling.
J. SymbolicComputation , 41:1157–1173, 2006.[10] C.W. Brown and S. McCallum. On using bi-equationalconstraints in CAD construction. In
Proc. ISSAC ’05 ,pages 76–83. ACM, 2005.[11] C. Chen, M. Moreno Maza, B. Xia, and L. Yang.Computing cylindrical algebraic decomposition viatriangular decomposition. In
Proc. ISSAC ’09 , pages95–102. ACM, 2009.[12] G.E. Collins. Quantifier elimination by cylindricalalgebraic decomposition – 20 years of progress. In
Quantifier Elimination and Cylindrical AlgebraicDecomposition , pages 8–23. Springer-Verlag, 1998.[13] J.H. Davenport, R. Bradford, M. England, andD. Wilson. Program verification in the presence ofcomplex numbers, functions with branch cuts etc. In
Proc. SYNASC ’12 , pages 83–88. IEEE, 2012.[14] J.H. Davenport and J. Heintz. Real quantifierelimination is doubly exponential.
J. SymbolicComputation , 5(1-2):29–35, 1988.[15] A. Dolzmann, A. Seidl, and T. Sturm. Efficientprojection orders for CAD. In
Proc. ISSAC ’04 , pages111–118. ACM, 2004. [16] M. England, R. Bradford, C. Chen, J.H. Davenport,M. Moreno Maza, and D. Wilson. Problemformulation for truth-table invariant cylindricalalgebraic decomposition by incremental triangulardecomposition. In
Intelligent Computer Mathematics (LNAI 8543), pages 45–60. Springer, 2014.[17] M. England, D. Wilson, R. Bradford, and J.H.Davenport. Using the Regular Chains Library to buildcylindrical algebraic decompositions by projecting andlifting.
Mathematical Software – ICMS 2014 (LNCS8592), pages 458–465. Springer Heidelberg, 2014.[18] M. Erascu and H. Hong. Synthesis of optimalnumerical algorithms using real quantifier elimination(Case Study: Square root computation). In
Proc.ISSAC ’14 , pages 162–169. ACM, 2014.[19] I.A. Fotiou, P.A. Parrilo, and M. Morari. Nonlinearparametric optimization using cylindrical algebraicdecomposition. In
Proc. CDC-ECC ’05 , pages3735–3740, 2005.[20] H. Hong. An improvement of the projection operatorin cylindrical algebraic decomposition. In
Proc. ISSAC’90 , pages 261–264. ACM, 1990.[21] Z. Huang, M. England, D. Wilson, J.H. Davenport,L. Paulson, and J. Bridge. Applying machine learningto the problem of choosing a heuristic to select thevariable ordering for cylindrical algebraicdecomposition. In
Intelligent Computer Mathematics (LNAI 8543), pages 92–107. Springer, 2014.[22] H. Iwane, H. Yanami, H. Anai, and K. Yokoyama. Aneffective implementation of a symbolic-numericcylindrical algebraic decomposition for quantifierelimination. In
Proc. SNC ’09 , pages 55–64, 2009.[23] S. McCallum. An improved projection operation forcylindrical algebraic decomposition. In
QuantifierElimination and Cylindrical Algebraic Decomposition ,pages 242–268. Springer-Verlag, 1998.[24] S. McCallum. On projection in CAD-based quantifierelimination with equational constraint. In
Proc.ISSAC ’99 , pages 145–149. ACM, 1999.[25] S. McCallum. On propagation of equationalconstraints in CAD-based quantifier elimination. In
Proc. ISSAC ’01 , pages 223–231. ACM, 2001.[26] S. McCallum and C.W. Brown. On delineability ofvarieties in CAD-based quantifier elimination with twoequational constraints. In
Proc. ISSAC ’09 , pages71–78. ACM, 2009.[27] L.C. Paulson. Metitarski: Past and future. In
Interactive Theorem Proving (LNCS 7406), 1–10.Springer, 2012.[28] A. Seidl. Cylindrical decomposition underapplication-oriented paradigms. PhD Thesis(University of Passau, Germany), 2006.[29] A. Strzebo´nski. Cylindrical algebraic decompositionusing validated numerics.
J. Symbolic Computation ,41(9):1021–1038, 2006.[30] A. Strzebo´nski. Cylindrical algebraic decompositionusing local projections. In
Proc. ISSAC ’14 , pages389–396. ACM, 2014.[31] D. Wilson, M. England, J.H. Davenport, andR. Bradford. Using the distribution of cells bydimension in a cylindrical algebraic decomposition.