The cylindrical algebraic decomposition algorithm is a quanti er elimination method for real-algebraic formulae. We use it as a theory solver in the context of satis ability-modulo-theories (SMT) solving to solve sequences of related real-algebraic satis ability questions. In this paper, we consider some optimizations for handling equational constraints. We review some previously published ideas, in particular Brown's projection operator and some improvements suggested by McCallum. Then we discuss di erent variants of the restricted projection operator to implement them in our SMT solver SMT-RAT and provide experimental results on the SMT-LIB benchmark set QF NRA. We show that the performance improves especially for unsatis able inputs.
Solving the satis ability (SAT) problem is an important sub-problem in many applications, ranging from program analysis [GAB+17] to industrial con guration management [Vol15] or the dependency management in Linux package managers [ADCTZ11]. An extension to the regular SAT problem is the satis abilitymodulo-theories (SMT) problem that allows for richer logics. While SAT only allows for propositional formulae, SMT deals with quanti er-free rst-order logic formulae over one or more theories. One popular theory is the theory of (nonlinear) real arithmetic.
An SMT solver [BHvMW09,KS08] is traditionally separated into a SAT solving part { which makes use of a regular SAT solving engine { and a theory solving part. The SAT solving part deals with the logical structure of the formula and issues theory queries to the theory solving part. The theory solving module checks sets of theory constraints for consistency.
A number of theory solving modules for non-linear real arithmetic have been proposed, including incomplete algorithms like interval constraint propagation [GGI+10,HR97] or virtual substitution [Wei97]. In this paper, we deal with the cylindrical algebraic decomposition (CAD) method [Col75] which is, to the best of our knowledge, the only complete decision procedure for non-linear real arithmetic implemented in any SMT solver.
Given a set (conjunction) of real-algebraic input constraints and an ordering of the variables, the CAD method proceeds in two phases. It rst uses a projection operator to produce a sequence of sets of polynomials of decreasing dimension. The sets of real roots of these polynomials can be used to construct the borders of nitely many semi-algebraic sets that we call cells. The points in each of the cells are equivalent regarding the satisfaction of the input formula. In the second phase, samples are constructed dimensionwise for each of the cells. Thus to decide satis ability it is su cient to check whether any of the generated samples satis es the formula.
A number of di erent projection operators exist, most notably the ones due to Collins [Col75], Hong [Hon90], Lazard [Laz94,MPP17], McCallum [McC85], and Brown [Bro01]. They use the same ingredients and mainly di er in the size of the sets of polynomials, essentially representing the continuing research in this area. It should be noted that the projection operators due to McCallum and Brown are incomplete in the sense that they may not generate some polynomials that are needed to nd satisfying solutions. Past studies [VKA17] have shown, however, that this incompleteness does not surface on the SMT-LIB [BFT16] benchmarks.
In addition to the di erent projection operators, several modi cations have been proposed to reduce the number of polynomials in certain cases. One particular modi cation that we analyse in this paper makes use of equational constraints. If equations are present in the input, they can be used to reduce the work for the projection.
In this paper, we present several possibilities to exploit the fundamental idea, in particular di erent variants of the restricted projection operator suggested by McCallum [McC01] in combination with Brown's projection operator [Bro01]. Then we discuss their impact on a CAD implementation used in the context of SMT solving in the SMT solver SMT-RAT and review some experimental results on the benchmark set QF NRA provided by SMT-LIB. Some of the presented optimizations are proven to be sound while others might lead to incompleteness. However, similar to the projection operators that are incomplete in general, none of our optimizations yielded incorrect satis ability results in our experiments.
After presenting some preliminaries in Section 2 we discuss the optimizations that we have implemented for the CAD projection phase in Section 3. In Section 4 we present and discuss our experimental results and conclude the paper in Section 5. 2 2.1
SMT solvers [BHvMW09,KS08] are used to decide the satis ability of rst-order logic formulae over one or more underlying theories. Traditionally, an SMT solver consists of a SAT solver and one or more theory solving modules. The SAT solver deals with the logical structure of the input formula by iteratively searching for solutions for its Boolean skeleton, which is the propositional logic formula obtained by replacing the theory constraints in the formula with fresh Boolean propositions. During this search, the theory solving modules are used to check
input
constraints
The CAD method [Col75] can be used as a theory solving module for the theory of (non-linear ) real arithmetic. It works with respect to a xed, static variable ordering that we assume to be given.
Polynomials p 2 Z[x1; :::; xn] with integer coe cients over variables x1; : : : ; xn are expressions of the form p = Pm j=1 xjeij with coe cients ai 2 Z and i=1 ai Qn j=1 xeij exponents ei;j 2 N0 for all i = 1; : : : ; m and j = 1; : : : ; n. The products Qn j are called monomials. If n = 1 then p is called univariate, otherwise multivariate. Multivariate polynomials in Z[x1; :::; xn] are usually interpreted as univariate polynomials in xn with polynomial coe cients from Z[x1; :::; xn 1], i.e., as polynomials from Z[x1; :::; xn 1][xn].
Polynomial constraints have the form p 0, where p is a polynomial in Z[x1; : : : ; xn] and where is one of the comparison predicates <; ; =; 6=; ; >. The sign of a polynomial p 2 Z[x1; :::; xn] in a point r 2 Rn is de ned as 1 if p evaluates in r to a strictly negative value, 1 for a strictly positive value, and 0 otherwise. To decide whether p 0 evaluates to true for a point r, it is su cient to determine the sign of p under r, since p is compared to zero. If the sign of p is the same under all points from a set C Rn then C is called p-sign-invariant ; in this case it su ces to determine the sign of p for a single point in the set to evaluate the constraint for all points in the set. Thus the satis ability of p can be determined if we can construct a nite partition C = fC1; : : : ; Cmg of Rn into nitely many p-sign-invariant cells Ci, which we call a decomposition.
A decomposition C = fC1; : : : ; Cmg is called Pn-sign-invariant for a set Pn of polynomials if it is p-sign-invariant for each polynomial from Pn. It is furthermore called algebraic if each Ci 2 C is a connected semi-algebraic1 set. Such a decomposition can be determined by the roots of the polynomials in Pn. It is additionally cylindrical when the cells are cylindrically ordered, which means that for each Ci; Cj 2 C and each 1 k < n the projections of Ci and Cj to Rn k (by removing the last k coordinates) are either disjoint or identical.
For a set Pn of polynomials, a Pn-sign-invariant cylindrical algebraic decomposition (CAD) of Rn as described above can be computed using the CAD method. Though in general CAD works on real arithmetic formulae, in this paper we consider only sets of polynomial constraints as input, as we use the CAD method as an SMT theory solving module. The considered set Pn contains all polynomials of the input constraints.
The CAD method computes the CAD in two phases: the projection phase and the lifting phase. In the projection phase the polynomials that describe the boundaries of the CAD cells are computed using a projection operator. The projection operator is applied to the polynomials in Pn Z[x1; : : : ; xn] to compute a set Pn 1 Z[x1; : : : ; xn 1] of polynomials, for which a CAD is computed recursively. The projection operator has the important property that each CAD Cn 1 for Pn 1 can be extended to a CAD Cn for Pn in an easy way by de ning the cells of Cn to be the Pn-sign-invariant regions in the cylinders Ci0 R for each cell Ci0 2 Cn 1.
There are several possible projection operators, in this paper, Brown's projection operator [Bro01] is used since compared to other operators it is faster and fewer polynomials have to be computed [VKA17]. In the following the properties degree deg (p), and leading coe cient lcf (p) of a polynomial p are used with the usual meaning. To de ne Brown's projection operator in addition the Sylvester matrix of univariate polynomials p = Pik=0 aixin and q = Pli=0 bixin in xn with polynomial coe cients ai; bi 2 Z[x1; : : : ; xn 1], deg (p) = k 1, and 1 A set C Rn that can be described by a conjunction of polynomial constraints is called semi-algebraic. deg (q) = l 1, which is the following (k + l) (k + l)-matrix, needs to be de ned. 0 ak B BBB ...
B
B 0 Syl xn (p; q) := B
BB bl B B BB@ ...
0 ak bl . . . . . .
ak bl a0 b0 a0 b0 : : : . . . . . .
0 19 .. >> . C=>
CC l a0 CCCC>>>; 0... CCCC9>>>=
C k C A>> > ; b0 Furthermore is the resultant of p and q de ned as res(p; q) = det (Syl xn (p; q)), the discriminant of p as disc(p) = det (Syl xn (p; p0)), and the content of p as cont(p) = gcd(a0; :::; ak), where gcd is the greatest common divisor. Finally Brown's projection operator P roj itself is de ned as follows: Let Pn0 = fp01; : : : ; p0mg be the nest square-free basis of Pn.
P roj(Pn) = Pn 1 = flcf (pi); disc(pi); res(pi; pj ) j pi; pj 2 Pn0 ; i 6= jg [ fcont(pi)jpi 2 Pn and cont(pi) is non-zero, non-unitg
Instead of polynomial sets, individual polynomials can be projected incrementally [CKJ+15]. The CAD itself is then modi ed accordingly in an also incremental lifting phase. This way the polynomials can be added and removed individually one after another. That enables to reuse large parts of previously computed CADs instead of recomputing them, which is useful since many similar CAD computations are needed when using the CAD method as a theory solving module in a less lazy SMT solver. This is due to the SAT solver that adds and removes only a few constraints while usually most constraints remain the same.
Given a Pn-sign-invariant CAD it su ces to consider one sample point from each cell to check whether any point in the cell satis es the input constraints. These sample points are computed in the lifting phase. Since the CAD method is used to check whether a set of polynomial constraints is consistent, which is the case if we can assign a real value to each of the variables occurring in the set such that each constraint evaluates to true, it is also su cient to compute only a partial CAD and stop the computation if a satisfying point is found [CH91]. 3
When we construct a CAD for the purpose of theory solving in an SMT solver, we do not need the CAD to be sign-invariant on the input polynomials, instead we are content if every cell is invariant with respect to the evaluations of the constraints. We call a CAD truth-invariant if the truth value of every input constraint is constant on each cell.
Additionally, we can exploit the fact that we only use the CAD method as a theory solving module, which means that it is only used to check the consistency of sets of constraints instead or arbitrary Boolean combinations of constraints. This implies that a part of the solution space where a single constraint evaluates to false can be discarded as a whole, even if other constraints are not signinvariant on this area. This part of the solution space may be more than a single cell of the CAD and we can try to avoid constructing individual cells within this part of the solution space altogether.
In order to exploit this, we consider modi cations of the projection operators based on equational constraints [McC99]. These are all constraints of the form p = 0, where p is a polynomial as de ned above, which we call an equational constraint polynomial, since we only consider conjunctions of constraints as input for the CAD method. Technically, McCallum suggested di erent ways to construct CADs that are sign-invariant with respect to equational constraints and sign-invariant with respect to the other constraints only in the cells where the equational constraints evaluate to true. Of course, this is only possible if equational constraints are present in the set of input constraints.
The advantage of these modi cations is that we can use a coarser CAD using fewer polynomials in the projection phase that still answers our question. This CAD may consider fewer polynomials in the projection phase which leads to a smaller number of sample points in the lifting phase and eventually a fewer number of cells. Therefore we expect the modi ed method to scale better on larger inputs in the presence of equational constraints.
We note that we compute partial CADs, i.e., we let the CAD method terminate if a satisfying sample is found, and thus a full CAD is usually only computed on unsatis able inputs. We also observe that many satis able SMT problems do not produce a lot of unsatis able theory calls since most examples in the used benchmark set have little Boolean structure and large Boolean satisfying regions, therefore we expect the impact of this modi cation to be bene cial mainly on unsatis able inputs. Furthermore, the lifting phase may even pro t if more polynomials are present that are comparably easy { for example with small degrees { in the projection as they may lead to a satisfying sample without considering hard polynomials at all. In such a case these modi cations may actually hinder the lifting phase and slow down our solver. 3.1
The rst modi cation suggested by McCallum is to restrict the projection operator that is applied in the rst step of the projection phase [McC99]. Afterwards, we continue with the original projection operator.
McCallum suggested this approach for his projection operator, however, it can be applied in combination with projection operators other than his own as well. We use this restricted projection with the projection operator due to Brown [Bro01].
For a nite set of constraints with polynomials Pn Z[x1; :::; xn] let E Pn a set of one polynomial which appears in an equational constraint and contains the variable xn that is to be eliminated rst. The restricted projection of Pn relative to E is de ned as follows: P rojE (Pn) =
P roj(E0) [ fres(e; p) j e 2 E0; p 2 Pn0 ; p 2= E0g [ fcont(pi) j pi 2 Pn and cont(pi) is non-zero, non-unitg = flcf (e); disc(e) j e 2 E0g [ fres(e; p) j e 2 E0; p 2 Pn0 ; p 2= E0g [ fcont(pi) j pi 2 Pn and cont(pi) is non-zero, non-unitg with the primitive part of a set of polynomials P de ned as prim(P ) = fp=cont(p) jp 2 P and p=cont(p) is not constantg and the nest square-free basis Pn0 for prim(Pn) and the nest square-free basis E0 for prim(E). Note that we may only be able to apply this restricted projection operator under an appropriate variable ordering as xn must be present in the equational constraint.
When using this restricted projection operator instead of the original projection operator less leading coe cients, discriminants, and resultants are added to the projection. This may reduce the size of the projection signi cantly which we hope to be bene cial for the overall performance. As mentioned before, the removal of comparably easy polynomials (in particular leading coe cients) may also be a disadvantage for satis able sets of constraints, though. We nevertheless hope a partial CAD using this modi ed projection operator to be faster on average.
The result when applying this approach is a CAD that is sign-invariant with respect to the used equational constraint and sign-invariant with respect to the other constraints in those cells where the equational constraint is satis ed. McCallum gave a proof that validates the use of P rojE in the rst projection step in [McC99], as well as in both steps for a 3-dimensional CAD.
In our implementation, we apply P rojE on all levels { provided that appropriate equational constraints are part of the input. Though doing so is not formally proven to be sound, we base our application on the following argument: the places where the proof fails are statistically rare so in the context of solving many problems quickly we accept the risk. Another modi cation suggested by McCallum is the semi-restricted projection operator [McC01], for which the repeated application is formally validated, wherever it is applicable throughout the projection phase. For sets of polynomials Pn Z[x1; :::; xn] and E = feg Pn, where e is an equational constraint polynomial that contains the variable xn, the semi-restricted projection of Pn relative to E is de ned as follows:
P rojE (Pn) = P rojE (Pn) [ fdisc(p)jp 2 Pn0 ; p 2= E0g with the nest square-free basis Pn0 for prim(Pn) and the nest square-free basis E0 for prim(E).
McCallum has shown that this operator can be used whenever applicable and that alternatively the restricted projection operator P rojE can be used for the rst step and the last step in the projection phase, and the semi-restricted projection operator P rojE in every other step [McC01]. For the latter combination a detailed complexity analysis can be found in [EBD15,ED16]. This allows using equational constraint polynomials at every projection step where one is present to reduce the size of the projection. Note that if the underlying projection operator is incomplete (like McCallum's or Brown's) then the restricted projection operator is also incomplete. 3.3
McCallum proposed a method called the resultant rule [McC01] to exploit the semi-restricted projection even if no explicit equational constraint is present for a speci c level. If e1 and e2 are both equational constraint polynomials their resultant res(e1; e2) is a propagated equational constraint polynomial, since e1 = 0 ^ e2 = 0 ) res(e1; e2) = 0. Due to this rule more polynomials in the projection are classi ed as equational constraint polynomials.
This method was only proposed for the semi-restricted projection, as the restricted projection was de ned for the top level only. Given that we assume the restricted projection to be usable on all levels as well, we could also use the resultant rule when we apply the restricted projection multiple times. Currently, we do not use this rule in our implementation. The reasons are speci c requirements for the underlying data structures, causing challenges for the implementation. Intuitively, a polynomial can be part of the projection due to several reasons. In the incremental setting, we would need to keep track of all possible reasons and apply postponed projection steps if e.g. an equational constraint is removed from the input set. 3.4
A di erent approach to modify the projection uses bounds. Bounds are polynomial constraints of the form b x + a 0, with a; b 2 Z, 2 f<; ; ; >g and a variable x. They can be used to neglect some polynomials in the projection, namely those that are for all permitted values of their variables either always positive or always negative since these have no roots and are therefore not needed to determine the CAD. For polynomials that can be neglected due to bounds no successors (leading coe cient, discriminant, and resultants) need to be computed. This approach is described in more detail in [LSC+13]. 4
We implemented several of the modi cations described in the previous section as part of the CAD module in our SMT solver SMT-RAT [CKJ+15]. These implementations are included in the publicly available version. To evaluate them we examine di erent combinations of the restricted and semi-restricted projection, as well as the simpli cation using bounds. We consider the following strategies: B uses only bounds to simplify the projection while R uses the restricted projection operator for as many as possible consecutive steps, starting at the rst step. BR combines these two modi cations. BRI extends it by allowing for an interruption of the restricted projection in the sense that it may be applied even when it is not applied in the step before. BSI employs the semi-restricted projection operator instead of the restricted one, otherwise, there is no di erence to BRI.
For comparison we use Default that is the standard CAD based strategy in SMT-RAT. It uses a { supposedly less powerful { variant of the simpli cation based on bounds, but no optimizations using equational constraints. The modications added to the regular CAD computations in the strategies B, BSI, and Default are provably sound, while for the other strategies no formal soundness proofs are provided yet. Additionally, all strategies are based on the projection operator due to Brown [Bro01] which is incomplete on certain inputs.
As benchmark problems we use the QF NRA benchmark set from the SMTLIB [BFT16] which consists of 12084 problems from 10 di erent applications. We used a time limit of 30 seconds per problem instance and allowed at most 4 GB of memory for every solver strategy and every input problem. The possible outcomes of a solving run are sat, unsat, timeout or memout. For sat and unsat the problem was solved correctly and is satis able respectively unsatis able. For timeout and memout the solver was unable to nd a solution or determine unsatis ability within the given time and memory limit.
We did not nd incorrect results for any of the solvers on any input problem, despite the incompleteness of the projection operator used. For the projection operator due to Brown inputs that actually lead to incorrect results are known, but past experiments [VKA17] already showed that this benchmark set does not contain such examples.
Default 4743 3939 3259 143 B 4764 3951 3225 144 R 4744 3962 3236 142 BR 4750 3964 3228 142 BRI 4752 4039 3151 142 BSI 4752 4038 3152 142
Table 1: Overall solver performances
The overall performance of each strategy is shown in Table 1. As expected all modi ed strategies solve more problems than Default and the improvements are mostly arising on unsatis able inputs. We assume that the reason is that the solver can only take advantage of the modi cations when a full CAD is computed since in that case fewer polynomials are needed. A full CAD is computed on an unsatis able problem to detect a theory con ict, which should occur more often for unsatis able problems. Satis able instances, on the other hand, are usually solved with only a fraction of the projection computed. We examine the number of solved problems for which theory con icts occurred in Table 2. Most satis able problems can be solved without a single theory con ict. Also, most of the problems that could be solved additionally are problems where a theory con ict occurs.
The strategy B was the best on satis able problems, but also the one with the least improvement on unsatis able problems. The pruning of polynomials due to bounds reduces the size of the projection but essentially does not change the lifting phase as the removed polynomials provide no new samples anyway. The restricted projection operator further decreases the size of the projection but also removes samples from the lifting phase. As the restricted projection tends to remove polynomials of smaller degrees, this may actually be detrimental for the lifting phase and thus for satis able instances. We note that in practice there seems to be no signi cant impact of using the restricted projection on satis able instances.
The strategies BRI and BSI achieve the best results overall. Allow for interruptions in between the application of the (semi-)restricted projection operator improves the solver's performance on unsatis able problems signi cantly compared to the other variants. At least on this set of benchmarks, the di erences between the restricted or semi-restricted projection are negligible.
Next, we take a closer look at the running times of the di erent strategies, based on the 8657 input problems that could be solved by all strategies. The results for the average running times are collected in Table 3. Compared to Default the simpli cations by using bounds and the restricted projection operator do speed up the computations signi cantly. The best average running time has the strategy BSI, closely followed by BRI. This directly re ects the superior overall result of these two strategies.
Last we take a look at the memouts, as we can see that none of the modications signi cantly changes the number of memouts. One would expect that signi cantly reducing the size of the projections would mitigate memory issues. However, the modi ed strategies, in contrast to Default, never actually remove polynomials but merely deactivate them. This causes the solver to consume more memory which leads to more memouts. It is possible to delete these polynomials instead of deactivating them when using the modi ed strategies as well. However, we do not expect that to signi cantly improve the overall performance since this is only relevant for large problems that are often hard to solve anyway and are therefore likely to just result in a timeout instead. This gets even more likely due to the fact that deleted polynomials might have to be recomputed later.
We examined this by means of the BRI strategy and implemented the corresponding strategy with the deletion of polynomials, in the following referred to as BRID. In Table 4 the results for these two strategies are compared and they are indeed relatively similar. As expected more timeouts occur when using BRID, while in BRI one more memout occurs. Thus we observe that deleting polynomials even decreases the overall solver performance though the average running time on the problems solved by all strategies is nearly the same for both variants. We examined the impact of restricted projection operators as proposed by McCallum in the CAD method on the performance of an SMT solver. After presenting several variants on how to use these in an actual implementation, we provided some experimental results for the implemented modi cations. We can show that the performance improved especially for unsatis able inputs when the restricted or semi-restricted projection operators are used instead of the original one. It, however, made no noticeable di erence whether the restricted or the semi-restricted projection operator was used, though the repeated application is currently only validated for the semi-restricted operator.
We further investigated the di erence of either deleting polynomials from the projection or only disabling them in the incremental setting during SMT solving. Though one could hope for a decreased memory consumption when deleting polynomials, this change did not make a signi cant di erence.
Our ideas for future investigation concerning equational constraints mainly deal with making the best possible use of the restricted projection operator by applying this idea in as many levels as possible. One direction would be the use of the resultant rule [McC01] that allows propagating equational constraints. Another option would be the modi cation of the variable ordering heuristic depending on the equations that are present in the input. It would also be interesting to investigate whether a heuristic for the choice of which equational constraint to use for the restricted projection operator can be found. Currently we are using the equational constraint that is added rst, however, the number of cells in the resulting CAD depends on the designated equational constraint as shown in [EBD15]. In that paper is furthermore shown how equational constraints can be used to make additional savings in the lifting phase, which could also be implemented in SMT-RAT. [Col75] [EBD15] [ED16] [GAB+17] [GGI+10] [Hon90] [HR97] [KS08] [Laz94] [LSC+13] [McC85] [McC99] [McC01] [MPP17] George E. Collins, Quanti er elimination for real closed elds by cylindrical algebraic decomposition, Automata Theory and Formal Languages, LNCS, vol. 33, Springer, 1975, pp. 134{183.
the use of equational constraints in cylindrical algebraic decomposition,
cal algebraic decomposition with respect to polynomial degree, Computer
Thiemann, Analyzing program termination and complexity automatically with AProVE, Journal of Automated Reasoning 58 (2017), no. 1, 3{31.
narayanan, and Edmund M. Clarke, Integrating ICP and LRA solvers for deciding nonlinear real arithmetic problems, Proceedings of FMCAD'10, IEEE, 2010, pp. 81{90.
Hoon Hong, An improvement of the projection operator in cylindrical algebraic decomposition, ISSAC '90 Proceedings of the International Symposium on Symbolic and Algebraic Computation (1990), 261{264. Stefan Herbort and Dietmar Ratz, Improving the e ciency of a nonlinearsystem-solver using a componentwise Newton method, Tech. Report
Daniel Kroening and Ofer Strichman, Decision procedures: An algorithmic point of view, Springer, 2008.
Daniel Lazard, An improved projection for cylindrical algebraic decomposition, Algebraic Geometry and its Applications: Collections of Papers from Shreeram S. Abhyankar's 60th Birthday Conference (Chandrajit L.
Bernd Becker, A symbiosis of interval constraint propagation and cylindrical algebraic decomposition, Proceedings of CADE-24, LNCS, vol. 7898, Springer, 2013, pp. 193{207.
Scott McCallum, An improved projection operation for cylindrical algebraic decomposition, Tech. report, University of Wisconsin Madison, 1985.
, On projection in CAD-based quanti er elimination with equational constraint, Proceedings of the 1999 International Symposium on
, On propagation of equational constraints in CAD-based quanti er elimination, Proceedings of the 2001 International Symposium on
of Lazard's method for CAD construction, Journal of Symbolic Computation (2017).