Cylindrical Algebraic Decomposition Using Validated Numerics

Adam Strzebonski

We present a version of the Cylindrical Algebraic Decomposition algorithm which uses approximate sample points in the lifting phase, whenever the results can be validated.This allows for substantial time savings by avoiding computations with exact algebraic numbers.

Introduction

The Cylindrical Algebraic Decomposition (CAD) algorithm is the best currently known algorithm for solving many classes of problems related to systems of real polynomial equations and inequalities. It not only allows elimination of quantifiers, but also provides an explicit description of the solution set in terms of real algebraic functions. The CAD algorithm has been used in practice in several areas including control system design, stability analysis, multidimensional integration and graphical representation of semialgebraic sets, and global optimization and assumption propagation in a computer algebra system.

Semialgebraic sets

A system of polynomial equations and inequalities in variables (Subscript[x, 1]), …, (Subscript[x, n]) is a formula

\(S \((x\_1, …, 
          x\_n)\) = \(\(∨\)\(\ \)\)\+\(i = 1\)\%l\ \( ∧ \+\(j = 
              1\)\%m\)\ \(f\_\(i, j\)\) \((x\_1, …, 
            x\_n)\) \(ρ\_\(i, j\)\) 0\)

where (Subscript[f, i, j]) ∈ ℝ[
          (Subscript[x, 1]), …, (Subscript[x, n])], and each (Subscript[ρ, i, j]) is one of <, ≤, >, ≥, =, or ≠.

A subset of  (ℝ^n) is semialgebraic if it is a solution set of a system of polynomial equations and inequalities.

A quantified system of real polynomial equations and inequalities in free variables (Subscript[x, 1]), …, (Subscript[x, n]) and quantified variables (Subscript[t, 1]), …, (Subscript[t, m]) is a formula

(Subscript[Q, 1]) (Subscript[t, 1]) ...  (Subscript[Q, m]) 
      (Subscript[t, m]) S ((Subscript[t, 1]), …, (Subscript[t, m]),  (Subscript[x, 1]), …, 
        (Subscript[x, n]))

where (Subscript[Q, i]) is ∃ or ∀, and S is a system of real polynomial equations and inequalities in (Subscript[t, 1]), …, (Subscript[t, m]), (Subscript[x, 1]), …, (Subscript[x, n]).

By Tarski's theorem, solution sets of quantified systems of real polynomial equations and inequalities are semialgebraic.

Cell decomposition

Every semialgebraic set can be represented as a finite union of disjoint cells defined recursively as follows.

- A cell in ℝ is a point or an open interval.

- A cell in (ℝ^k) has one of the two forms

{((Subscript[a, 1]), …, (Subscript[a, k]), 
            (Subscript[a, k + 1])) : ((Subscript[a, 1]), …, (Subscript[a, k])) ∈ 
            (Subscript[C, k])⋀(Subscript[a, k + 1]) = 
        r ((Subscript[a, 1]), …, (Subscript[a, k]))}
{((Subscript[a, 1]), …, (Subscript[a, k]), 
          (Subscript[a, k + 1])) : ((Subscript[a, 1]), …, (Subscript[a, k])) ∈ 
          (Subscript[C, k])⋀(Subscript[r, 1]) ((Subscript[a, 1]), …, (Subscript[a, k])) < 
            (Subscript[a, k + 1]) < (Subscript[r, 2]) ((Subscript[a, 1]), …, 
                (Subscript[a, k]))}

where (Subscript[C, k]) is a cell in (ℝ^k), r  is a continuous algebraic function, and (Subscript[r, 1]) and (Subscript[r, 2]) are continuous algebraic functions, -∞, or ∞, and (Subscript[r, 1]) < (Subscript[r, 2]) on (Subscript[C, k]).

By an algebraic function we mean a function  r : 
        (Subscript[C, k])⟶ℝ for which there is a polynomial

\(f = \(c\_0\) x\_\(k + 1\)\^m + \(c\_1\) 
            x\_\(k + 1\)\^\(m - 1\) + …c\_m ∈ ℝ[x\_1, …, x\_k, 
          x\_\(k + 1\)]\)

such that

(Subscript[c, 0]) ((Subscript[a, 1]), …, (Subscript[a, k])) ≠ 
        0⋀f ((Subscript[a, 1]), …, (Subscript[a, k]), 
            r ((Subscript[a, 1]), …, (Subscript[a, k]))) = 0

The CAD algorithm, introduced by Collins (1975), allows to compute a cell decomposition of any semialgebraic set presented by a quantified system of polynomial equations and inequalities. The objective of the original Collins algorithm was to eliminate quantifiers from a quantified system of polynomial equations and inequalities and to produce an equivalent quantifier-free system of polynomial equations and inequalities. After finding a cell decomposition the algorithm performed an additional step of finding an implicit representation of the semialgebraic set in terms of polynomial equations and inequalities in the free variables. Our objective here is somewhat different. Given a semialgebraic set presented by a system of polynomial equations and inequalities, quantified or not, we find a cell decomposition of the set, explicitly written in terms of algebraic functions.

IneqPlot code

<< Experimental`;
<< Developer`;
<< Utilities`FilterOptions`;
<< Graphics`ImplicitPlot`;
<< Graphics`FilledPlot`;
Off[General::spell1];
Off[General::spell];
Off[$MaxExtraPrecision::meprec];

\(\(\(IneqPlot[ineq_, {x_, x0_, x1_}, {y_, y0_, y1_}, 
        opts___?OptionQ] := 
	
      Module[{cyls, \ plots, \ 
          bds}, 
		cyls = \(GenericCylindricalAlgebraicDecomposition[
              ineq && x0 < x < x1 && y0 < y < y1, {x, 
                y}]\)[\([1]\)]; 
		bds = \(Boundaries /. {opts}\) /. {Boundaries -> True}; 
		plots = 
			If[cyls === False, 
            Plot[y1 + 1, {x, x0, x1}, PlotRange -> {{x0, x1}, {y0, y1}}, 
              Evaluate[FilterOptions[\ Graphics, \ opts\ ]], 
              DisplayFunction -> Identity], 
					If[
              Head[cyls] === Or, cyls = List @@ cyls, 
              cyls = {cyls}]; 
					\(cylinderplot[#, x, 
                  FilterOptions[FilledPlot, \ opts]] &\) /@ cyls]; 
		If[
          bds\  && \ cyls =!= False, \ 
          plots = {plots, 	\(boundaryplot[#, x, 
                    FilterOptions[Graphics, \ opts]] &\) /@ 
                cyls}]; 
	\ \ Show[\ 
          plots, \ 
\ \ \ \ \ \ FilterOptions[\ Graphics, \ 
            opts], \ 
\ \ \ \ \ \ Axes -> True, \ 
          AxesLabel -> {x, 
              y}, \ 
\ \ \ \ \ \ DisplayFunction -> $DisplayFunction, 
          AspectRatio -> Automatic, 
          PlotRange -> {{x0, x1}, {y0, y1}}\ ]]\)\(
\)
    \)\)
\(\(\(cylinderplot[cyl_, x_, opts___?OptionQ] := 
	 (*\ 
        assumes\ form\ a < x < 
            b\  && \ \((\(\(c < y < d\)\(\ \)\(||\)\)\  ... )\)\ *) 
	
      Module[{x0 = cyl[\([1, 1]\)], x1 = cyl[\([1, 5]\)], ys}, 
		ys = 
          If[Head[cyl[\([2]\)]] === Or, 
            List @@ \((cyl[\([2]\)])\), {cyl[\([2]\)]}]; 
		\(FilledPlot[\ {#[\([1]\)], #[\([5]\)]}, \ {x, \ x0, x1}, \ DisplayFunction -> Identity, 
              opts] &\) /@ ys]\)\(		\)\(
\)
    \)\)
\(\(\(boundaryplot[cyl_, x_, opts___?OptionQ] := 
	 (*\ 
        assumes\ form\ a < x < 
            b\  && \ \((\(\(c < y < d\)\(\ \)\(||\)\)\  ... )\)\ *) 
	
      Module[{x0 = cyl[\([1, 1]\)], x1 = cyl[\([1, 5]\)], ys}, 
		ys = 
          If[Head[cyl[\([2]\)]] === Or, 
            List @@ \((cyl[\([2]\)])\), {cyl[\([2]\)]}]; 
		\({Graphics[
                Line[{{x0, #[\([1]\)] /. x -> x0}, {x0, #[\([5]\)] /. 
                        x -> x0}}]], 
			Graphics[
                Line[{{x1, #[\([1]\)] /. x -> x1}, {x1, #[\([5]\)] /. 
                        x -> x1}}]], 
              Plot[\ {#[\([1]\)], #[\([5]\)]}, \ {x, \ x0, x1}, \ 
                DisplayFunction -> Identity, opts]} &\) /@ 
          ys]\)\(		\)\(
\)
    \)\)
xypoints[{xl_, xu_, yl_, yu_}, pp_] := 
      Block[{x, y}, 
        Module[{ε = 10. ^-6, dx, dy, nx}, δ = 
            1. /pp; nx = 
            Max[Floor[(xu - xl - 2  ε)/δ], 
              1]; δx = (xu - xl - 2  ε)/nx; 
          dx = xu - xl; 
          dy = 0; 
Do[
            x = xl + ε + ix δx; 
            dy = Max[dy, yu - yl], {ix, 0, nx}]; 
{Ceiling[
              dx/δ], Ceiling[dy/δ]}]]

    
makePolygons3D[
        Inequality[xl_, Less, x, Less, xu_] && 
          Inequality[yl_, Less, y, Less, yu_] && 
          Inequality[zl_, Less, z, Less, zu_], pp_] := 
      Block[{x, y}, 
        Module[{ε = 10. ^-12, nx, ny}, {nx, ny} = 
            xypoints[{xl, xu, yl, yu}, pp]; 
          If[nx < 1 || ny < 1, 
            Print[{nx, ny}]]; δx = (xu - xl - 2  ε)/
              nx; 
pointsl = 
            Table[x = 
                xl + ε + 
                  ix δx;  δy = (yu - yl - 
                      2  ε)/ny; 
Table[
                y = yl + ε + iy δy;  {x, y, 
                  Re[zl]}, {iy, 0, ny}], {ix, 0, 
                nx}]; 
pointsu = 
            Table[x = 
                xl + ε + 
                  ix δx;  δy = (yu - yl - 
                      2  ε)/ny; 
Table[
                y = yl + ε + iy δy;  {x, y, 
                  Re[zu]}, {iy, 0, ny}], {ix, 0, 
                nx}]; 
{Table[
              Polygon[{pointsl[[i, j]], pointsl[[i + 1, j]], 
                  pointsl[[i + 1, j + 1]], pointsl[[i, j + 1]]}], {i, 
                nx}, {j, ny}], 
            Table[Polygon[{pointsu[[i, j]], pointsu[[i + 1, j]], 
                  pointsu[[i + 1, j + 1]], pointsu[[i, j + 1]]}], {i, 
                nx}, {j, ny}]}]]

    
CADPlot3D[gcad_, opts___] := 
      With[{gc = 
            gcad[[1]] //. 
              a_ && b_Or :> (a && # & /@ 
                    b)}, 
With[{polyg = makePolygons3D[#, 
                    Evaluate[(PlotPoints /. {opts}) /. 
                        PlotPoints -> 15]] & /@ 
                Switch[Head[gc], Or, List @@ gc, And, {gc}, _, {}]}, 
          Print[Length[Flatten[polyg]]]; 
          Graphics3D[{EdgeForm[], polyg}, FilterOptions[Graphics3D, opts], 
            PlotRange -> All, 
            AspectRatio -> Automatic]]]

\(CADPlot2D[cad_, opts___] :=  (*assumes\ nonempty\ and\ bounded*) 
	
      Module[{cyls, \ plots, \ x, y}, 
cyls = 
          cad //. a_ && b_Or :> \((\(a && # &\) /@ b)\); 
        If[Head[cyls] === Or, cyls = List @@ cyls, cyls = {cyls}]; 
        cyls = Select[cyls, FreeQ[#, Equal] &]; x = cyls[\([1, 1, 3]\)]; 
        y = cyls[\([1, 2, 3]\)]; 
        plots = \(cylinderplot[#, x, FilterOptions[FilledPlot, \ opts]] &\) /@ 
            cyls; 
	\ \ Show[\ 
          plots, \ 
\ \ \ \ \ \ FilterOptions[\ Graphics, \ 
            opts], \ 
\ \ \ \ \ \ Axes -> True, \ 
          AxesLabel -> {x, 
              y}, 
\ \ \ \ \ \ DisplayFunction -> $DisplayFunction, 
          AspectRatio -> Automatic, PlotRange -> All\ ]]\)

Example of a cell decomposition of a semialgebraic set

ineqs = (((x + (7/4)))^2) + (((y - (9/4)))^2) < 
          (1/9) || (((x - (9/4)))^2) + (((y - (9/4)))^2) < (1/9) || 
        (x^2) + (y^2) < 25 && (((x - 2))^2) + (((y - 2))^2) > 
            1 && (((x + 2))^2) + (((y - 2))^2) > 
            1 && ((x^2) + (((y + 1))^2) > 5 || 
              (x^2) + (((y - 1))^2) < 9); 
    IneqPlot[ineqs, {x, -5, 5}, {y, -5, 5},  
      Fills :> Hue[Random[]]]

Representation of real algebraic functions

In Mathematica algebraic functions can be represented as Root objects

\(TraditionalForm\`Root[\(c\_0\) #\^m + \(c\_1\) #\^\(m - 1\) + … + c\_m &, \ i]\)

where i is the sequential number of the root. The root numbers of real roots start with 1, and the roots are numbered by the increasing value, counting multiplicities.

Algebraic functions are piecewise continuous.

Plot3D[
      Root[a #(1^5) + ((a^2) + 1) #1 + b &, 1], {a, -3, 
        3}, {b, -3, 3}, PlotPoints -> 30]

We can use the CAD algorithm to find the possible discontinuity sets of algebraic functions, i.e. the sets of points where at least two roots intersect.

cad = 
      CylindricalAlgebraicDecomposition[
        Resultant[a (x^5) + ((a^2) + 1) x + b, 
            D[a (x^5) + ((a^2) + 1) x + b, x], x] == 0, {a, 
          b}]

Plot[{cad[[1, 2, 1, 2]], cad[[1, 2, 2, 2]]}, {a, -3, 0}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}, 
      PlotRange -> {{-3, 3}, {-3, 3}}, 
      AspectRatio -> Automatic]

Here is a plot of all real roots of a quintic.

f[i_] := 
      Root[1 + 3 (a^3) #1 - (a^2) #(1^2) - 5 a #(1^3) + #(1^5) &, 
        i]
real[x_] := If[Im[x] == 0, x, 100]
Plot[Table[real[f[i]], {i, 5}], {a, -4, 4}, 
      PlotStyle -> Table[Hue[(i/5)], {i, 0, 4}], 
      PlotRange -> {-4, 4}]

The projection phase of the CAD algorithm

Finding a cell decomposition of a semialgebraic set using the CAD algorithm consists of two phases, projection and lifting. In the projection phase we start with the set (Subscript[A, n]) of  factors of the polynomials present in the system, and eliminate variables one by one using a projection operator P such that

(Subscript[P, k + 1]) : ℝ[(Subscript[x, 1]), …, (Subscript[x, k]), 
          (Subscript[x, k + 1])] ⊃ 
        (Subscript[A, k + 1])⟶(Subscript[A, k]) ⊂ ℝ[
          (Subscript[x, 1]), …, (Subscript[x, k])]

and, generally speaking, if all polynomials of (Subscript[A, k]) have constant signs on a cell C ⊂ (ℝ^k), then all polynomials of (Subscript[A, k + 1]) are delineable over C,  i.e. each has a fixed number of real roots on C as a polynomial in (Subscript[x, k + 1]), the roots are continuous functions on C, they have constant multiplicities, and two roots of two of the polynomials are equal either everywhere or nowhere in C. This way the roots of polynomials of (Subscript[A, 1]), …, (Subscript[A, n]) are the algebraic functions needed in the construction of the cell decomposition of the semialgebraic set.

There have been several improvements made reducing the size of the original Collins projection.The currently best projection operator applicable in all cases is due to Hong, however in most situations we can use a smaller projection operator given by McCallum, with an improvement by Brown. There are even smaller projection operators that can be applied in some special cases. When equational constraints are present we can use the projection operator suggested by Collins, and developed and proven by McCallum. When there are no equations and only strict inequalities, and there are no free variables, or we are interested only in the full-dimensional part of the semialgebraic set we can use an even smaller projection operator (and we need not use algebraic number sample points, so the  improvement presented here does not apply in this case).

The lifting phase of the CAD algorithm

In the lifting phase we find a cell decomposition of the semialgebraic set. Generally speaking (the actual details depend on the projection operator used), we start with cells in (ℝ^1) consisting of all distinct roots of (Subscript[A, 1]), and the open intervals between the roots. We find a sample point in each of the cells, and remove the cells whose sample points do not satisfy the system describing the semialgebraic set (the system may contain conditions involving only (Subscript[x, 1])). Now we lift the cells to cells in (ℝ^n), one dimension at a time. Suppose we have lifted the cells to (ℝ^k).  To lift a cell C ⊂ (ℝ^k) to (ℝ^(k + 1)) we find the real roots of  (Subscript[A, k + 1]) with (Subscript[x, 1]), …(Subscript[x, k]) replaced with the coordinates of the sample point c in C. Since the polynomials of  (Subscript[A, k + 1])are delineable on \(TraditionalForm\`C\), each root  r  is a value of a continuous algebraic function at c, and the function can be represented as a p-th Root of a polynomial f ∈ (Subscript[A, k + 1])  such that r is the p-th root of f(c, (Subscript[x, k + 1])). Now the lifting of the cell C to (ℝ^(k + 1)) will consist of graphs of these algebraic functions, and of the slices of C×ℝ between the subsequent graphs. The sample points in each of the new cells will be obtained by adding the k + 1-st coordinate to c, equal to one of the roots, or to a number between two subsequent roots. Similarly as in the first step, we remove those lifted cells whose sample points do not satisfy the system describing the semialgebraic set.

The coefficients of sample points computed this way are in general algebraic numbers, and algebraic number computations are often dominating the computation time of CAD. Our improvement is to use Mathematica  arbitrary-precision floating point number approximations of the coefficients, whenever the results can be validated. Arbitrary-precision numbers are floating point numbers with unrestricted number of digits, which contain information about the number of digits considered accurate. The arithmetic automatically keeps track of precision. Note that using approximate arithmetic may be enough to prove that two roots of a polynomial or a pair of polynomials are distinct, and to find a nonzero sign of a polynomial at a sample point. What we cannot prove with approximate arithmetic is that two roots of a polynomial or a pair of polynomials are equal, or that a polynomial is zero at a sample point. However, we can often use information about the origins of the cell to resolve these problems. For instance, if we know that the resultant of two polynomials vanishes on the cell, and these two polynomials have exactly one pair of complex roots that can be equal within the precision bounds, we can conclude that these roots are equal. Similarly, if the last coordinate of a sample point was a root of a factor of the given polynomial, we know that this polynomial is zero at the sample point. If we cannot resolve all the uncertainties using the collected information about the cell we compute the exact algebraic number values of the coordinates.

In the ``generic'' situation we will never encounter  ``unexpected'' singularities, i.e.ones that do not follow from our construction of the cell. And indeed, in experiments with randomly generated systems we very rarely see cells whose sample points need to be recomputed exactly, and we see a dramatic speed improvement compared to the exact version of CAD. We have also tried several examples coming from practical applications of CAD,and those examples tend to have some singularities. However in many examples we still see a substantial speed-up, even though some of the cells need to be computed using exact arithmetic.

The Algorithm

Projection phase

We use the smallest projection appropriate for the given example. Whenever applicable we use the equational constraints, otherwise we attempt to use the McCallum's projection with Brown's improvement. When the system turns out to be not well-oriented, we compute the Hong's projection.

During the projection phase we annotate each projection polynomial p with ``genealogy'' information. The possible annotations may be
- InputPolynomial(f)  when p is a factor of a polynomial f from the input system.
- Coefficient (f, x, k) when p is a factor of a coefficient at (x^k) of  f.
- PSC (f, x, k) when p is a factor of the k-th principal subresultant coefficient of f and ∂(f/∂)x in x.
- PSC(f, g, x, k) when p is a factor of the k-th principal subresultant coefficient of f and g in x.

First step of the lifting phase

In the first step of the lifting phase we use exact real root isolation to find the roots of the univariate projection polynomials, and if there is no equational constraint, we find a rational number between each two subsequent roots. We select those of the roots and rational numbers which do not make the input system false, and we set the roots that are not rational numbers to a fixed finite precision. This way we obtain the sample points in cells in (ℝ^1). We annotate each sample point with all annotations of all projection polynomials that are zero at this root.

k+1-st step of the lifting phase

Suppose we have a cell C ⊂ (ℝ^k) with a sample point c = ((Subscript[c, 1]), …, (Subscript[c, k])). Let  F = {(Subscript[F, 1]), …, (Subscript[F, s])} be the set of k + 1-variate projection polynomials.

For each (Subscript[f, i]) ∈ F we check whether we know from the construction of C that some leading terms of (Subscript[f, i]) as a polynomial in (Subscript[x, k + 1]) are zero on C, i.e. we look for  annotations of c of the form Coefficient((Subscript[f, i]), (Subscript[x, k + 1]), m). If we find out that some  (Subscript[f, i]) ∈ F are identically zero on C, we collect their annotations and later add them to annotations of all sample points extending c (and, if we are using equational constraint or McCallum's projection, we use this information to alter the projection polynomials or to decide to redo the computation with Hong's projection).

Let G = {(Subscript[g, 1]), …, (Subscript[g, t])} consist of elements of F with the identified terms that vanish on C removed, and with (Subscript[f, i]) ∈ F that are identically zero on C removed or replaced with delineating polynomials, and let

H = {(Subscript[g, 1]) (c, (Subscript[x, k + 1])), …, (Subscript[g, t]) (c, 
            (Subscript[x, k + 1]))}

We need to find approximations of all the real roots of elements of H, and identify which roots are equal and which are distinct. We find approximations of all complex roots of elements of H using the Jenkins-Traub method, and compute their guaranteed precision bounds. If the computation fails due to loss of precision (which includes the case when one of the polynomials has an inexact zero leading coefficient), we revert to exact computations. How to revert to exact computations we will describe later on.

We decide which of the roots are real and which pairs of roots are equal using the ``worst case scenario'', i.e. if a root may be real within the precision bounds it is considered real, and if two roots may be equal within the precision bounds they are considered equal. Now, for each pair of elements of  H which has common real roots, we find the number r of common complex roots, and check whether the construction of \(TraditionalForm\`C\) implies that these polynomials have at least  r common complex roots, i.e.whether we can find the necessary principal subresultant coefficients among the annotations of  c. Similarly, for each element of H which has a multiple real root we find the total multiplicity r of its roots, equal to its degree less the number of its distinct complex roots, and check whether the total multiplicity r  is implied by the principal subresultant coefficient annotations of c. If yes, we are guaranteed that all approximations that are equal within the precision bounds indeed represent equal roots, and we can continue the construction of cells extending C. If not, we need to revert to exact computation .(Of course, unless we use Hong's projection, the annotations may contain only resultants, and so each time a polynomial has more than one double root, or a pair of polynomials has more than one pair of common roots we need to revert to exact computations.)

If we found no ``unexpected'' multiple roots or common roots we can continue with the construction of cells in (ℝ^(k + 1))extending C. The cells will consists of graphs of continuous algebraic functions on C, passing through the real roots of elements of H, and, if we did not use the equational constraint projection to eliminate (Subscript[x, k + 1]), slices of C×ℝ between two subsequent graphs. The algebraic function passing through a root a can be written as Root[(Subscript[g, i]),  p], where a is the p-th real root of (Subscript[g, i])(c, (Subscript[x, k + 1])), multiplicities counted. (If (Subscript[g, i]) is a reduct of an element of F we may use the element of  F instead, in order to get a smaller representation of the output.) The sample points in the cells extending C, which are graphs of algebraic functions, are points (c, a) where a is one of the roots of elements of H, annotated with all annotations of  c and all annotations of elements of F which vanish at  (c, a). The sample points in the cells extending C,  which are slices of  C×ℝ between graphs of algebraic functions are points (c, b), where b is a rational number between two subsequent  roots of elements of H, annotated with all annotations of c and all annotations of elements of F which vanish identically on C×ℝ.

For all projections other than the equational constraint projection, F contains all factors of polynomials of the input system that depend only on ((Subscript[x, 1]), …, (Subscript[x, k + 1])). For the equational constraint projection the procedure described above needs to be altered somewhat, to make sure that if a polynomial f of the input system vanishes on a cell C' extending C, then the sample point of C' is annotated with InputPolynomial(f). Let P be the set of factors of polynomials of the input system that depend on (Subscript[x, k + 1]), and do not depend on (Subscript[x, k + 2]), …, (Subscript[x, n]). Similarly as for F, we remove the leading terms of elements of P that we know are zero on C, and replace ((Subscript[x, 1]), …, (Subscript[x, k])) with ((Subscript[c, 1]), …, (Subscript[c, k]))  in those resulting polynomials that are not identically zero, obtaining a set Q of polynomials in (Subscript[x, k + 1]). We compute the complex roots of elements of Q. If we cannot compute the roots due to loss of precision, we revert to exact computation. Next, for every pair of an element of H and an element of Q that have a common real root (within the precision bounds), we check whether the number of common complex roots is implied by the principal subresultant coefficient annotations of c. If no, we revert to exact computation. If there are no ``unexpected'' common roots, we know that if a real root a of an element of H is equal to a real root of an element of Q within the precision bounds, then a is indeed a root of the element of Q, and so we can add annotation InputPolynomial(f) corresponding to the element of Q to the annotations of the sample point (c, a). Moreover, sample points of all cells extending C get InputPolynomial(f) annotations corresponding to elements of P that vanish identically on C.

To conclude the lifting of C, we remove those cells C' extending C whose sample points make the input system false. To this end we need to evaluate signs of each polynomial f of the input system which depends on (Subscript[x, k + 1]) and does not depend on (Subscript[x, k + 2]), …, (Subscript[x, n]) at the sample points. Let  c' = ((Subscript[c, 1]), …, (Subscript[c, k]), 
          (Subscript[c, k + 1])) be the sample point of C'. If the annotations of c' contain  InputPolynomial(f)  we know that the sign is zero. Otherwise we replace with (Subscript[x, 1]), …, (Subscript[x, k + 1]) in f with the coordinates of c'. If the we get an inexact zero, we revert to exact computations, else the sign of f on C' is the (nonzero) sign of the result.

Reverting to exact computations

Let us now describe how we ``revert to exact computations'' for a cell C. We need to compute the exact algebraic number values of ((Subscript[c, 1]), …, (Subscript[c, k])). We know the exact value of (Subscript[c, 1]), since in the first step of the lifting phase we used exact real root isolation. Suppose we have computed exact values ((Subscript[c, 1])', …, (Subscript[c, l])') of ((Subscript[c, 1]), …, (Subscript[c, l])). If (Subscript[c, l + 1]) is a rational number then we set (Subscript[c, l + 1])' = (Subscript[c, l + 1]). Otherwise we have a representation of (Subscript[c, l + 1])' as a p-th root of a polynomial with coefficients depending polynomially on((Subscript[c, 1])', …, (Subscript[c, l])'), so we get  (Subscript[c, l + 1]) by computing the real roots of this polynomial. We remember the remaining roots in case we need them later on. When we have found the exact values of coordinates of the sample point in C we continue lifting C using the classical version of CAD.

Remarks

When all coordinates of c are rational numbers all elements of H have exact rational coefficients, and so instead of using the method described above we can obtain their roots using exact real root isolation methods, and then set the algebraic number roots to a finite precision. Experiments suggest that this is most of the time
a better solution. The examples where computing with approximate numbers for rational sample points was faster involved low degree random polynomial systems with large integer coefficients.

Annotations telling which polynomial is zero on a given cell are also used by our implementation of CAD using exact algebraic number coordinates of sample points. They require no additional computations,just keeping
track of information, and allow to reduce the number of necessary sign computations.In some examples introducing the annotations reduced the CAD computation time by up to 25%.

Empirical Results

We compare our implementation of the version of CAD algorithm described above (APCAD, for arbitrary-precision floating point sample point coordinates used) with a version of CAD using exact algebraic sample point coordinates (ECAD). Both algorithms were implemented in the C kernel of Mathematica, and share the projection phase and the general outline of the lifting phase. Also whenever APCAD ``reverts to exact computation'', it uses the ECAD code.

The experiments have been conducted on a 700 MHz Pentium III computer, with 256 MB of RAM. All timings given in Time columns in tables are in seconds, the column Proj contains the times used by theprojection phase, the columns APCAD and ECAD contain times of the lifting phases of APCAD and ECAD. The Cells columns give the total number of ``fully lifted'' cells constructed during the lifting phase. A cell is ``fully lifted'' if it is a cell in (ℝ^n), where n is the total number of variables in the system, or if it is a lower dimensional cell that gets rejected, because its sample point makes the system false. For APCAD the Cells column gives the total number of ``fully lifted'' cells/the number of ``fully lifted'' cells that were computed with exact algebraic number sample points. The degree columns contain the maximal degrees of minimal polynomial of algebraic numbers that appeared in the computation.

Examples from applications

Stability of the Dormand-Prince fifth order embedded seven stage method

This is Example 4.4 from H. Hong, R. Liska, S. Steinberg, ``Testing Stability by Quantifier Elimination'', J. Symbolic Comp., 24 (1997), 161-188.

R[z_] := 
      1 + z + ((z^2)/2) + ((z^3)/6) + ((z^4)/24) + ((z^5)/120) + ((z^6)/600)
f = (Subscript[∀ , {x, y}, (x | y) ∈ Reals  &&  x < 0  &&  (x^2) + (y^2) < (99438/100000)])Expand[
            R[x + I y] R[x - I y]] < 1

Reduce[f]

({{Time

          Proj, Time

          APCAD, Time

          ECAD, Cells

          APCAD, Cells

          ECAD, Degree

          APCAD, Degree

          ECAD}, {1.52, 0.92, 233, 235/56, 235, 72, 72}})

Stability of a six-point upwind-biased second-order accurate scheme for approximating a two-dimensional advection equation

This is Example 5.4 from H. Hong, R. Liska, S. Steinberg, ``Testing Stability by Quantifier Elimination'', J. Symbolic Comp., 24 (1997), 161-188.

A = 
        (c^4) (a - b + 1) (a - b - 1) (((a - b))^2);
B = 
        2 (c^4) b (3 (a^2) b - 2 (a^2) - 2 a (b^2) + a + (b^3) - b) + 
          4 (c^3) a b ((a^2) - a + (b^2) - b) + 
          2 (c^2) a ((a^3) - 2 (a^2) b + 3 a (b^2) - a - 2 (b^2) + 
                b);
CC = 
        (c^4) (b^2) ((b^2) - 1) + 4 (c^3) a (b^2) (b - 1) + 
          2 (c^2) a b (3 a b - 2 a - 2 b + 1) + 
          4 c (a^2) b (a - 1) + (a^2) ((a^2) - 1);
DD = (c^2) R + 2 c S + T;
R = 
        8 (a^2) (b^2) - 12 (a^2) b + 5 (a^2) - 8 a (b^3) + 8 a (b^2) + 
          2 a b - 4 a + 4 (b^4) - 4 (b^3) - 3 (b^2) + 4 b;
S = 
        4 (a^3) b - 2 (a^3) - 4 (a^2) (b^2) - 2 (a^2) b + (a^2) + 
          4 a (b^3) - 2 a (b^2) + 2 a b - 2 (b^3) + (b^2);
T = 
        4 (a^4) - 8 (a^3) b - 4 (a^3) + 8 (a^2) (b^2) + 8 (a^2) b - 
          3 (a^2) - 12 a (b^2) + 2 a b + 4 a + 5 (b^2) - 
          4 b;

Reduced form of (5.14)

f = (Subscript[∀ , {a, b, c}, (a | b | c) ∈ Reals  &&  a ≥ 0  &&  b ≥ 0  &&  4 ((a^2) + (b^2)) < 1])(DD ≤ 0 || B ≤ 0)

Reduce[f]

Full form of (5.14)

g = (Subscript[∀ , {a, b, c}, (a | b | c) ∈ Reals  &&  0 ≤ a ≤ 1  &&  0 ≤ b ≤ 1])(A ≤ 0 && 
            CC ≤ 
              0 && (DD ≤ 0 || B ≤ 0))

Reduce[g]

({{Example

          , Time

          Proj, Time

          APCAD, Time

          ECAD, Cells

          APCAD, Cells

          ECAD, Degree

          APCAD, Degree

          ECAD}, {Reduce[f], 1.31, 10223, >36000, 5975/1111, ?, 480, ?}, {Reduce[g], 1.28, 18235, >36000, 6187/1056, ?, 480, ?}})

Robust multi-objective feedback design

This is Example 4.1 from P. Dorato, W. Yang, C. Abdallah, ``Robust Multi-Objective Feedback Design by Quantifier Elimination'', J. Symbolic Comp., 24 (1997), 153-160.

f = (Subscript[∀ , {p1, p2, w1, w2}, (p1 | p2 | w1 | w2 | q1) ∈ Reals  &&  16 ≤ 20 p1 ≤ 25  &&  16 ≤ 20 p2 ≤ 25  &&  0 ≤ w1 ≤ 2])
(p2 (1 + p1 q1) < 
              0 && -24 (w1^2) + (p2^2) ((((1 + p1 q1))^2) - 25) > 
              0 && (400 - (q1^2)) (w2^2) + 
                (p2^2) (400 (((1 + p1 q1))^2) - (q1^2)) > 
              0); 
Reduce[f, q1]

This is the second part of Example 4.2 from P. Dorato, W. Yang, C. Abdallah, ``Robust Multi-Objective Feedback Design by Quantifier Elimination'', J. Symbolic Comp., 24 (1997), 153-160. The first part requires showing that the formula is true for n/d=41/10, which both APCAD and ECAD do in 0.15 second.

f = (Subscript[∃ , {q1, q2}, q1 > 1 && q2 > 0])((Subscript[∀ , w, w ∈ Reals])(((n/d) - 
                            (q1^2)) (w^4) + (((n ((((1 + q1))^2) - 2 q2))/d) - ((q2^2) + (q1^2))) (w^2) + ((n/d) - 
                            1) (q2^2) ≥ 
                    0 && ((n/d) - 
                            (q1^2)) (w^4) + (((n ((((-1 + q1))^2) - 2 q2))/d) - ((q2^2) + (q1^2))) (w^2) + ((n/d) - 
                            1) (q2^2) ≥ 0));
Reduce[f /. (n/d) -> 4]

({{Example

          , Time

          Proj, Time

          APCAD, Time

          ECAD, Cells

          APCAD, Cells

          ECAD, Degree

          APCAD, Degree

          ECAD}, {4.1, 0.27, 28.8, 46.6, 66871/1877, 66254, 4, 4}, {4.2, 0.08, 1.52, 14.0, 282/27, 282, 24, 48}})

A geometric inequality

This is an inequality from D. Mitrinovic, J. E. Pecaric, V. Volenec, ``Recent Advances in Geometric Inequalities'', Kluwer Academic Publishers, 1989.

f = (Subscript[∀ , {a, b, c},  a ≥ 0  &&  b ≥ 0  &&  c ≥ 0  &&  a + c ≥ b  &&  b + c ≥ a  &&  a + b ≥ c])
(a^3) + (b^3) + (c^3) ≥ 
            3 a b c + 
              2 (Sqrt[9 + 6 (Sqrt[3])]) (a - b) (b - c) (-a + 
                    c);
Reduce[f]

Find constants that can replace 2 (Sqrt[9 + 6 (Sqrt[3])]) in this inequality.

g = (Subscript[∀ , {a, b, c}, a ≥ 0 && b ≥ 0 && c ≥ 0 && a + c ≥ b && b + c ≥ a && a + b ≥ c && k ∈ Reals])
(a^3) + 
            (b^3) + (c^3) ≥ 
          3 a b c + 
            k (a - b) (b - c) (-a + c); 
Reduce[g,  k]

ToRadicals[%] // ExpandAll

({{Example

          , Time

          Proj, Time

          APCAD, Time

          ECAD, Cells

          APCAD, Cells

          ECAD, Degree

          APCAD, Degree

          ECAD}, {Reduce[f], 0.35, 0.32, 16.4, 271/17, 271, 12, 48}, {Reduce[g], 0.35, 25.6, 324, 34035/505, 34035, 12, 108}})

Randomly generated systems

We have computed CAD of systems randomly generated in the following way. Each system is a conjunction of randomly generated polynomial equations and inequalities. Each equation and inequality contains 5 terms: one randomly generated 10-bit constant and 4 randomly generated terms with a prescribed number of variables, a prescribed maximal total degree, and randomly generated 10-bit coefficients. S(n,d,e,s,w) denotes a randomly generated system with n variables, maximal total degree d, e equations, s strict inequalities, and w  weak inequalities. Results given in the table are averages from 10 randomly generated examples of each type. Both algorithms were run with the same examples.

({{Example

          , Time

          Proj, Time

          APCAD, Time

          ECAD, Cells

          APCAD, Cells

          ECAD, Degree

          APCAD, Degree

          ECAD}, {S (2, 3, 1, 1, 1), 0.02, 0.075, 0.27, 19/0, 19, 7.1, 9.2}, {S (2, 3, 0, 2, 2), 0.04, 0.71, 2.13, 524/0, 524, 8.7, 23.9}, {S (2, 5, 1, 1, 1), 0.03, 0.16, 2.44, 24/0, 24, 16.5, 31.8}, {S (2, 5, 0, 2, 2), 0.1, 1.2, 14.2, 615/0, 615, 19.2, 70.1}, {S (2, 7, 1, 1, 1), 0.11, 0.49, 77.0, 36/0.6, 36, 35.3, 103}, {S (2, 7, 0, 2, 2), 0.34, 1.62, 377, 604/6, 604, 40.0, 229}, {S (2, 10, 1, 1, 1), 0.52, 2.21, >3600, 48/0.5, ?, 65.0, ?}, {S (2, 10, 0, 2, 2), 1.6, 23.2, >3600, 937/5, ?, 111, ?}, {S (3, 2, 2, 0, 0), 0.02, 0.19, 3.32, 26/1.2, 26, 4.0, 10.9}, {S (3, 2, 1, 1, 0), 0.02, 0.39, 4.99, 82/0.4, 82, 3.4, 15.7}, {S (3, 2, 0, 1, 1), 0.03, 0.74, 15.5, 472/2, 472, 4.0, 18.3}, {S (3, 3, 2, 0, 0), 0.08, 3.31, >3600, 123/3, ?, 15.0, ?}, {S (3, 3, 1, 1, 0), 0.24, 1.23, >3600, 201/3, ?, 17.0, ?}, {S (3, 3, 0, 1, 1), 0.24, 54.9, >3600, 1572/2, ?, 30.1, ?}})