CAD Adjacency Computation Using Validated Numerics

Adam Strzeboński, Wolfram Research, Inc .

Abstract

We will present an algorithm for computation of cell adjacencies for well-based cylindrical algebraic decomposition. The algorithm determines cell adjacency information using validated numerical methods similar to those used in CAD construction, thus computing CAD with adjacency information in time comparable to that of computing CAD without adjacency information. Our implementation in Mathematica uses cell adjacency information to compute topological operations e.g. closure, boundary, and connected components. Other applications include computing topological properties e.g. homology groups, visualization and path planning.

Code

CAFView

cview[cad_,vars_,rul_,dvars_,ind_]:=
Module[{k=Length[ind],cc,ee,a,b,ans=False,ha,hb,x,na,rr},
ha=Hue[k/Length[vars],1/3];
hb=Hue[(k+1)/Length[vars],1/3];
If[rul=!={},
a=cad/.rul;
If[a===True||a===False,Return[Row[{a},Background->ha,Frame->True]]]];
x=dvars[[k+1]];
cc=If[Head[cad]===Or,List@@cad,{cad}];
If[k==Length[vars]-1,
Return[Or@@(Row[{Tooltip[#,N[#/.rul]]},Background->ha,Frame->True]&/@cc)]];
ee=cc[[1]];
If[Head[ee]===And,ee=ee[[1]]];
If[!FreeQ[ee,Alternatives@@Drop[vars,k+1]],
Return[Row[{True},Background->ha,Frame->True] &&Row[{Opener[Dynamic[Evaluate[x]]],PaneSelector[{False->Row[{Subscript["b",Sequence@@ind,1]},Frame->True,Background->hb],True->Row[{cview[cad,vars,rul,dvars,Append[ind,1]]},Frame->True]},Dynamic[Evaluate[x]],ImageSize->Automatic]}]]];
Do[ee=cc[[i]];
{a,b}=If[Head[ee]===And, {ee[[1]],Drop[ee,1]},{ee,True}];
na=N[a/.rul];
rr=If[na===True,rul,{}];
If[na=!=False,ans=ans||Row[{Tooltip[a,na]},Background->ha,Frame->True] &&Row[{Opener[Dynamic[Evaluate[x]]],PaneSelector[{False->Row[{Subscript["b",Sequence@@ind,i]},Frame->True,Background->hb],True->Row[{cview[b,vars,rr,dvars,Append[ind,i]]},Frame->True]},Dynamic[Evaluate[x]],ImageSize->Automatic]}]],
{i,Length[cc]}];
ans]

CAFView[cad_,vars_,rul_:{}]:=With[{dvars=Table[Unique["dynamicx"],{Length[vars]}]},cview[cad,vars,rul,dvars,{}]]

CAD picture

plusform[f_List]:=If[Length[f]==1,f[[1]],HoldForm[Plus[##]]&[Sequence@@f]]
timesform[{f_,g_}]:=If[NumericQ[f]||NumericQ[g],f g,HoldForm[f g]]
polyform[f_,vars_]:=plusform[(timesform[{polyform[#[[2]],Rest[vars]],vars[[1]]^#[[1,1]]}]&/@CoefficientRules[f,vars[[1]]])]
polyform[f_,{}]:=f

rdisp[Root[f_,i_]]:=Row[{(Subscript[Root, y,i]),polyform[f[y],{y,x}]}]

fig1[pq_]:=Module[{ineqs,pl,cad,bds,pts,ff,pp,n,cbs},
ff={-5 (((y-3))^3)+7 (x^2) (y-3)+1,(x^4)+2 (((y-4))^4)-100};
ineqs=And@@Thread[ff<0];
pl=RegionPlot[ineqs,{x,-4.5,4.5},{y,-0.5,7.75},AspectRatio->Automatic,PlotPoints->100];
cad=CylindricalDecomposition[ineqs,{x,y}];
bds=Select[cad/.{Or->List,And->List,HoldPattern[Inequality[a_,__,b_]]:>{a,b}},ListQ[#[[1]]]&];
pts=Sort[Union@@(First/@bds),Less];
pp=If[pq,N[Point[{#,y}/.Solve[Or@@Thread[ff==0]/.x->#,y,Reals]]&/@pts],{}];
n=Length[bds];
cbs=Table[{Graphics[{c=Hue[2(i-1)/n],Thickness[0.005],Line[{{pts[[i]],0},{pts[[i+1]],0}}]}],Plot[Tooltip[#,rdisp[#]],{x,bds[[i,1,1]],bds[[i,1,2]]},ColorFunction->(c&),PlotStyle->Thickness[0.005],PlotRange->{{-4.5,4.5},{-0.5,7.75}}]&/@Flatten[bds[[i,2]]]},{i,n}];
Show[{pl,cbs,Graphics[{{Thickness[0.002],Gray,Line[{{#,-0.5},{#,7.75}}]&/@pts},{PointSize[Large],Red,pp}}]},Axes->True,Frame->False]]

CADAdjacencyPlot code

$offset=0.00001;

linepts[Line[ll_]]:=If[ListQ[ll[[1,1]]],Join@@ll,ll]

curve[f_,x_,a_,b_,offset_,opts_]:=Union@@(linepts/@Cases[Plot@@{f,{x,a+(b-a)*offset,b+(a-b)*offset},PlotRange->All,Exclusions->None,Sequence@@opts},_Line,Infinity])

cellplot[True,{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]:=Polygon[{{x1,y1},{x2,y1},{x2,y2},{x1,y2}}]

cellplot[cond_,{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]/;FreeQ[cond,x]:=
cellplot[Inequality[x1,Less,x,Less,x2]&&cond,{{x,x1,x2},{y,y1,y2}},opts]

cellplot[x_>a_,{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]:=
Polygon[{{a,y1},{x2,y1},{x2,y2},{a,y2}}]

cellplot[x_<b_,{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]:=
Polygon[{{x1,y1},{b,y1},{b,y2},{x1,y2}}]

cellplot[HoldPattern[Inequality[a_,Less,x_,Less,b_]],{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]:=
Polygon[{{a,y1},{b,y1},{b,y2},{a,y2}}]

cellplot[x_>a_&&cond_,{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]:=
cellplot[Inequality[a,Less,x,Less,x2]&&cond,{{x,x1,x2},{y,y1,y2}},opts]

cellplot[x_<b_&&cond_,{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]:=
cellplot[Inequality[x1,Less,x,Less,b]&&cond,{{x,x1,x2},{y,y1,y2}},opts]

cellplot[HoldPattern[Inequality[a_,Less,x_,Less,b_]]&&y_>r_,{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]:=
Module[{up,dn},
dn=curve[r,x,a,b,$offset,opts];
up={{a,y2},{b,y2}};
Polygon[Join[dn,Reverse[up]]]]

cellplot[HoldPattern[Inequality[a_,Less,x_,Less,b_]]&&y_<s_,{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]:=
Module[{up,dn},
dn={{a,y1},{b,y1}};
up=curve[s,x,a,b,$offset,opts];
Polygon[Join[dn,Reverse[up]]]]

cellplot[HoldPattern[Inequality[a_,Less,x_,Less,b_]]&&HoldPattern[Inequality[r_,Less,y_,Less,s_]],{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]:=
Module[{up,dn},
dn=curve[r,x,a,b,$offset,opts];
up=curve[s,x,a,b,$offset,opts];
Polygon[Join[dn,Reverse[up]]]]

cellplot[x_==a_,{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]:=
Line[{{a,y1},{a,y2}}]

cellplot[x_==a_&&y_>r_,{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]:=
Line[{{a,N[r/.x->a]},{a,y2}}]

cellplot[x_==a_&&y_<s_,{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]:=
Line[{{a,y1},{a,N[s/.x->a]}}]

cellplot[x_==a_&&HoldPattern[Inequality[r_,Less,y_,Less,s_]],{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]:=
Line[{{a,N[r/.x->a]},{a,N[s/.x->a]}}]

cellplot[HoldPattern[Inequality[a_,Less,x_,Less,b_]]&&y_==r_,{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]:=
Line[curve[r,x,a,b,$offset,opts]]

cellplot[x_==a_&&y_==r_,{{x_,x1_,x2_},{y_, y1_,y2_}},opts_]:=
Point[{a,N[r/.x->a]}]

ptcoords[x_==a_&&y_==r_,{x_,y_}]:={a,N[r/.x->a]}
ptcoords[other__]:=None

cellformat[{c_Point,ind_,form_,inQ_}]:=Tooltip[{PointSize[Large],If[inQ,Hue[Random[]],GrayLevel[0.5]],c},{ind,form}]

cellformat[{c_Line,ind_,form_,inQ_}]:=Tooltip[{Thickness[Large],If[inQ,Hue[Random[]],GrayLevel[0.7]],c},{ind,form}]

cellformat[{c_Polygon,ind_,form_,inQ_}]:=Tooltip[{If[inQ,Hue[Random[]],GrayLevel[0.9]],c},{ind,form}]

CADAdjacencyPlot[ineqs_,{x_, y_},op___]:=
Module[{opts={op},x1,x2,y1,y2,d,cells,lo,hi,inQ,n,cc,ee,vv},
{cells,lo,hi,inQ}=Reduce`CADAdjacencies[ineqs,{x,y}];
cc=Select[ptcoords[#,{x,y}]&/@cells,#=!=None&];
If[cc==={},
{x1,x2,y1,y2}={-1,1,-1,1},
cc=Transpose[cc];
{x1,x2,y1,y2}={Min[cc[[1]]],Max[cc[[1]]],Min[cc[[2]]],Max[cc[[2]]]};
d=If[x1<x2,0.25(x2-x1),1];
x1-=d;
x2+=d;
d=If[y1<y2,0.25(y2-y1),1];
y1-=d;
y2+=d];
n=Length[cells];
cc=Table[{cellplot[cells[[i]],{{x,x1,x2},{y,y1,y2}},opts],i,cells[[i]],inQ[[i]]},{i,n}];
cc=Join[Select[cc,Head[#[[1]]]===Polygon&],Select[cc,Head[#[[1]]]===Line&],Select[cc,Head[#[[1]]]===Point&]];
SeedRandom[1234];
cc=cellformat/@cc;
vv=Select[Range[n],inQ[[#]]&];
ee=Join@@Table[If[inQ[[i]],DirectedEdge[i,#]&/@Intersection[lo[[i]],vv],{}],{i,n}];
GraphicsRow[{Graphics[cc],Graph[vv,ee,VertexLabels->"Name"]}]]

IntervalRoots picture

SeedRandom[1234];
cc=RandomReal[{-10,10},{7,2}];
rr=RandomReal[{1,2},7];
dd=MapThread[Disk[#1,#2]&,{cc,rr}];
pp=Table[RandomPoint[Disk[cc[[i]],0.7rr[[i]]],RandomInteger[{2,5}]],{i,7}];
qq=RandomPoint[Annulus[{0,0},{21,30}],5];
fig2=Graphics[{{Opacity[0.3],Blue,dd},{PointSize[Large],{Red,Point[qq]},{Red,Point[Rest[#]],Green,Point[First[#]]}&/@pp},{Red,Circle[{0,0},20]}}];

CADAdjacency picture

SeedRandom[1234];
{xl,xu}={12,21};
n0=5;
{x0,x1}={16,19};
cc=Table[-50+100/(n0+2)i,{i,n0+1}];
cc=Drop[cc,{(n0+1)/2}];
cc0=RandomReal[{#-2,#+2}]&/@cc;
rr0=RandomReal[{3,5},n0];
rec0=MapThread[Rectangle[{xl,#1-#2},{xu,#1+#2}]&,{cc0,rr0}];
y0=Table[RandomReal[{cc0[[i]]-0.5rr0[[i]],cc0[[i]]+0.5rr0[[i]]}],{i,n0}];
p0=Point[{x0,#}&/@y0];
ran[{a_,b_}]:=If[Random[]>0.5,{Max[#],Min[#]},{RankedMin[#,3]}]&[RandomReal[{a,b},5]];
cc1=Join[Join@@Table[ran[{cc0[[i]]-0.5rr0[[i]],cc0[[i]]+0.5rr0[[i]]}],{i,n0}],{-52.7,54.5,57.9}];
n1=Length[cc1];
rr1=RandomReal[{0.5,0.7},n1];
rec1=MapThread[Rectangle[{x1-0.7,#1-#2},{x1+1,#1+#2}]&,{cc1,rr1}];
y1=Table[RandomReal[{cc1[[i]]-0.75rr1[[i]],cc1[[i]]+0.75rr1[[i]]}],{i,n1}];
p1=Point[{x1,#}&/@y1];
ll={Line[{{xl,-50},{xu,-50}}],Line[{{xl,50},{xu,50}}]};
fig3=Graphics[{{Opacity[0.3],Blue,rec0},{Opacity[0.9],Yellow,rec1},{PointSize[Medium],{Red,p1},{Green,p0}},{Red,ll},{Opacity[0.5],Red,Thick,Line[{{x0,0},{40,0}}]},{Thick,Blue,Line[{{xl,0},{xu,0}}],Yellow,Line[{{x1-0.7,0},{x1+1,0}}]},{PointSize[Medium],{Red,Point[{x1,0}]},{Green,Point[{x0,0}]}},{Text[C',{x0,-3}],Text[C,{30,-3}]}},PlotRange->{{-60,60},{-60,60}},Axes->True];

Semialgebraic sets

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

S((Subscript[x, 1]),...,(Subscript[x, n]))== (Underoverscript[∨, (i=1), k]) (Underoverscript[∧, (j=1), l]) (Subscript[f, i,j])((Subscript[x, 1]),...,(Subscript[x, n])) (Subscript[ρ, 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.

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+1)) 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 semialgebraic function, and (Subscript[r, 1]) and (Subscript[r, 2]) are continuous semialgebraic functions, -∞, or ∞, and (Subscript[r, 1])<(Subscript[r, 2]) on (Subscript[C, k]).

fig1[False]

• By a semialgebraic function we mean a function

(Subscript[Root, (Subscript[x, k+1]),p])f:(Subscript[C, k])∋((Subscript[x, 1]),...,(Subscript[x, k]))⟶(Subscript[Root, (Subscript[x, k+1]),p])f((Subscript[x, 1]),...,(Subscript[x, k]))∈ℝ

where

f==(Subscript[c, 0])(Subscript[x, k+1]^m)+(Subscript[c, 1]) (Subscript[x, k+1]^(m-1))+…+(Subscript[c, m])∈ℝ[(Subscript[x, 1]),…,(Subscript[x, k]),(Subscript[x, k+1])]

is a polynomial and

(Subscript[Root, (Subscript[x, k+1]),p])f((Subscript[x, 1]),...,(Subscript[x, k]))

is the p-th real root (multiplicities counted) of f treated as a univariate polynomial in (Subscript[x, k+1]).

• Polynomials whose roots bound cells extending C are delineable on C.
A set {(Subscript[f, 1]),...,(Subscript[f, l])}⊆ℝ[(Subscript[x, 1]),…,(Subscript[x, k])][(Subscript[x, k+1])] of polynomials is delineable on C⊆(ℝ^k) if all (Subscript[f, i]) are degree-invariant on C, the real roots of (Subscript[f, i ]) are continuous functions on C, have constant multiplicities on C, and any two of the root functions are either always equal or always different on C.

Cylindrical Algebraic Decomposition

A cylindrical algebraic decomposition (CAD) of (ℝ^n) is a finite collection D of pairwise disjoint cylindrically arranged cells in  (ℝ^n) such that (⋃_(C∈D))C=(ℝ^n).

A CAD D of (ℝ^n) is consistent with a semialgebraic set A if A is a union of cells of D.

The CAD algorithm, introduced by Collins (1975), allows to compute a CAD of (ℝ^n) consistent with a semialgebraic set given by a system of polynomial equations and inequalities.

ineqs=5 (((y-3))^3)-7 (x^2) (y-3)>1&&(x^4)+2 (((y-4))^4)≤100;
CylindricalDecomposition[ineqs,{x,y}]

cad=CylindricalDecomposition[ineqs,{x,y},"Function"]

Reduce`CADPlot[cad,{x,-4,4},{y,0.5,7}]

comp=CylindricalDecomposition[cad,{x,y},"ComponentsFunction"]

GraphicsRow[(Reduce`CADPlot[#,{x,-4,4},{y,0.5,7}]&/@comp)]

CELL Adjacencies

Let (Subscript[C, 1]),(Subscript[C, 2])∈D. (Subscript[C, 1]) and (Subscript[C, 2]) are adjacent if (Subscript[C, 1])≠(Subscript[C, 2]) and (Subscript[C, 1])(Subscript[∪C, 2]) is connected.

CADAdjacencyPlot[5 (((y-3))^3)-7 (x^2) (y-3)>1&&(x^4)+2 (((y-4))^4)≤100,{x,y}]

Well-based CAD

Let D be a CAD of (ℝ^n). Let (Subscript[Π, k])(D) denote the CAD of (ℝ^k) obtained by projecting the cells of D on the space of the first k coordinates.

A CAD D of (ℝ^n) is well-based (Schwartz & Sharir, 1983) if each polynomial whose root bounds a cell in (Subscript[Π, k+1])(D) has a constant sign on all cells in (Subscript[Π, k+1])(D) and does not vanish identically on any cell in (Subscript[Π, k])(D), for 1≤k<n.

Each semialgebraic set has a well-based CAD after a linear change of coordinates.

Excluded case 1

(-1<x<0∧y==(Subscript[Root, y,2])(y^2)+(x^2)-1)∨x==0 is a CAD representation of (x<0∧y>0∧(x^2)+(y^2)==1)∨x==0

Show[{Plot[Sqrt[1-x^2],{x,-1.2,0},Axes->False,PlotStyle->Thick],Graphics[{Green,Thick,Line[{{0,0},{0,1.2}}]}]},PlotRange->All,AspectRatio->Automatic]

This CAD can be refined to a well based CAD by subdividing the line x==0 with the roots of (y^2)+(x^2)-1 at x==0.

Excluded case 2

(x==0∧y==0)∨(x>0∧y>0∧z==(Subscript[Root, z,1])y z-x) is a CAD representation of x≥0∧z>0∧y z==x

Show[{Plot3D[x/y,{x,0,1},{y,0,1},Axes->False,PlotPoints->100],Graphics3D[{Green,Thickness[0.01],Line[{{0,0,0},{0,0,10}}],Red,Line[{{0,0,0},{0,1,0}}]}]}]

Properties of well-based CAD

Let D be a well-based CAD of (ℝ^n). Then:

• (Subscript[Π, k])(D) is a well-based CAD of (ℝ^k).

• The closure of a cell of D is a union of cells of D.

• Lower-dimensional cells adjacent to a section must be sections defined by the same polynomial.

More precisely, if C and C' are adjacent cells in (Subscript[Π, k])(D), dim(C')<dim(C), and

S = {(a,(Subscript[a, k+1])):a∈C∧(Subscript[a, k+1])==(Subscript[Root, f,p])(a)}

is a cell in (Subscript[Π, k+1])(D) then either

(OverBar[S])⋂(C'×ℝ) = {(b,(Subscript[b, k+1])):b∈C'∧(Subscript[b, k+1])==(Subscript[Root, f,q])(b)}

or ∀b∈C' (Underscript[lim, a∈C,a->b])(Subscript[Root, f,p])(a)=-∞, or ∀b∈C' (Underscript[lim, a∈C,a->b])(Subscript[Root, f,p])(a)=∞.

• All cell adjacencies can be determined by adjacencies between cells whose dimensions differ by one.

More precisely, if (Subscript[C, k]) and (Subscript[C, l]) are adjacent cells with k=dim((Subscript[C, k]))≤ dim((Subscript[C, l]))=l then there exist cells (Subscript[C, k+1]),...,(Subscript[C, l-1]) such that dim((Subscript[C, j]))=j  and (Subscript[C, j])⊂(OverBar[(Subscript[C, j+1])]) for k≤j<l.

Finding cell adjacencies

Algorithm CADAdjacency (general idea)
Input: A well-based CAD D of (ℝ^n).
Output: The set of all pairs (C,C') of adjacent cells of D such that dim(C)-dim(C')=1.

• Let (Subscript[A, 1]) be the set of all pairs of adjacent intervals and points in (Subscript[Π, 1])(D).

• For k=1 to n-1

• Set (Subscript[A, k+1])=∅.

• For each section

S = {(a,(Subscript[a, k+1])):a∈C∧(Subscript[a, k+1])==(Subscript[Root, f,p])(a)}

in (Subscript[Π, k+1])(D) and each cell C'  such that (C,C')∈(Subscript[A, k]) check which of 
	  the three possibilities is true:

(OverBar[S])⋂(C'×ℝ) = S' = {(b,(Subscript[b, k+1])):b∈C'∧(Subscript[b, k+1])==(Subscript[Root, f,q])(b)},

or ∀b∈C' (Underscript[lim, a∈C,a->b])(Subscript[Root, f,p])(a)=-∞,

or ∀b∈C' (Underscript[lim, a∈C,a->b])(Subscript[Root, f,p])(a)=∞.

If the first possibility is true, add (S,S')  to (Subscript[A, k+1]).

• For each sector S in (Subscript[Π, k+1])(D) bounded by (possibly infinite) sections (Subscript[S, 1]) and (Subscript[S, 2])

• Add (S,(Subscript[S, 1]))  and (S,(Subscript[S, 2]))  to (Subscript[A, k+1]) (if finite).

• For each cell C' such that (C,C')∈(Subscript[A, k]) let (Subscript[S, 1]^(')) and (Subscript[S, 2]^(')) be 
		  the (possibly infinite) sections over C' adjacent to (Subscript[S, 1]) and (Subscript[S, 2]). 
		  For all sectors S' between (Subscript[S, 1]^(')) and (Subscript[S, 2]^(')) add (S,S')  to (Subscript[A, k+1]).

• Return (Subscript[A, n]).

The main computational step

Given a pair (C,C') of adjacent cells with dim(C)-dim(C')=1 and a section

S = {(a,(Subscript[a, k+1])):a∈C∧(Subscript[a, k+1])==(Subscript[Root, f,p])(a)}

find the (possibly infinite) section S' over C' adjacent to S. That is decide whether

(OverBar[S])⋂(C'×ℝ) = S' = {(b,(Subscript[b, k+1])):b∈C'∧(Subscript[b, k+1])==(Subscript[Root, f,q])(b)},

or ∀b∈C' (Underscript[lim, a∈C,a->b])(Subscript[Root, f,p])(a)=-∞ (then S'=-∞),

or ∀b∈C' (Underscript[lim, a∈C,a->b])(Subscript[Root, f,p])(a)=∞ (then S'=∞).

IntervalRoots

Algorithm IntervalRoots (specification)
Input: Discs (Subscript[Z, 0]),..., (Subscript[Z, N])⊂ℂ.
Output: Discs (Subscript[D, 1]),..., (Subscript[D, m])⊂ℂ, positive integers (Subscript[k, 1]),..., (Subscript[k, m]), and a positive radius R such that:

• (Subscript[D, i])⋂(Subscript[D, j])=∅ for i≠j,

• (Subscript[D, i])⊂Δ(0,R),

• for any (Subscript[a, 0])∈(Subscript[Z, 0]), ..., (Subscript[a, N])∈(Subscript[Z, N]) the polynomial f=(Subscript[a, N])(z^N)+...+(Subscript[a, 0]) has exactly (Subscript[k, i]) roots (counted with multiplicity) in (Subscript[D, i]) and f has no roots in Δ(0,R)\(⋃_(i=1)^n)(Subscript[D, i]).

The algorithm may return Failed.

Property:
For any g=(Subscript[b, N])(z^N)+...+(Subscript[b, 0])∈ℂ[z]\{0} and any M>0 there exists ϵ>0 such that if (Subscript[b, i])∈(Subscript[Z, i]) and ρ((Subscript[Z, i]))<ϵ for all 0≤i≤N then IntervalRoots outputs  (Subscript[D, 1]),..., (Subscript[D, m]), (Subscript[k, 1]),..., (Subscript[k, m]), and R such that

• ρ((Subscript[D, i]))<1/M,

• g has exactly one root in (Subscript[D, i]) (necessarily of multiplicity (Subscript[k, i])),

• g has no roots in ℂ\(⋃_(i=1)^n)(Subscript[D, i]),

• R>M.

fig2

The main computational step (general idea)

Given a pair (C,C') of adjacent cells with dim(C)-dim(C')=1 and a section

S = {(a,(Subscript[a, k+1])):a∈C∧(Subscript[a, k+1])==(Subscript[Root, f,p])(a)}

find the (possibly infinite) section S' over C' adjacent to S. That is decide whether

(OverBar[S])⋂(C'×ℝ) = S' = {(b,(Subscript[b, k+1])):b∈C'∧(Subscript[b, k+1])==(Subscript[Root, f,q])(b)},

or ∀b∈C' (Underscript[lim, a∈C,a->b])(Subscript[Root, f,p])(a)=-∞ (then S'=-∞),

or ∀b∈C' (Underscript[lim, a∈C,a->b])(Subscript[Root, f,p])(a)=∞ (then S'=∞).

• Find a sample point in each cell of the CAD D of (ℝ^n) such that each coordinate is given by a isolating disc satisfying the IntervalRoots property for all projection polynomials.

• For k=1 to n and each pair (C,C') of adjacent cells with dim(C)-dim(C')=1 find an additional sample point in C with coordinates given by isolating discs contained in the isolating discs of the main sample point in C'.

fig3

Empirical results

An algorithm computing CAD cell adjacencies has been implemented in C, as a part of the kernel of Mathematica. The implementation takes a quantified system of polynomial equations and inequalities S and uses Mathematica’s multi-algorithm implementation of CAD computation to find a CAD D of (ℝ^n) consistent with the solution set A of S. If D is well-based then the implementation finds the cell adjacencies. The implementation is geared towards solving a specified topological problem, e.g. finding the boundary or the connected components of A, hence it avoids computing cell adjacencies for cells that are known not to belong to the closure of A. The current implementation also works for non-well-based problems in (ℝ^3) using ideas from D. S. Arnon, G. E. Collins, and S. McCallum, Adjacency Algorithm for Cylindrical Algebraic Decomposition of Three-dimensional Space, J. Symbolic Comp. 5 (1988), 163–187.

The experiments have been conducted on a Linux laptop with a 4-core 2.7 GHz Intel Core i7 processor and 16 GB of RAM. The reported CPU time is a total from all cores used. For each example we give three timings. (Subscript[t, CAD]) is the computation time of constructing a CAD consistent with the solution set the input system. (Subscript[t, SP]) is the time of refining the CAD to a well-based CAD D of (ℝ^n) and of constructing the “main” sample points in the CAD cells. Our implementation refines the CAD while constructing sample points, which is why we cannot give separate timings. The third timing, (Subscript[t, ADJ]) is the time of computing cell adjacency information. We also report the dimension dim of the embedding space, the number (Subscript[N, CELL]) of cells in the CAD of A, the number (Subscript[N, ADJ]) of computed pairs of adjacent cells whose dimensions differ by one, and the number (Subscript[N, CC]) of connected components of A.

Union of two unit balls in  (ℝ^n)

Find cell adjacencies for a CAD of the union of two unit balls in  (ℝ^n): B({0,...,0}, 1)⋃B({1,...,1}, 1).

Note that for n≤3 the balls have full-dimensional intersection, for n=4 they touch at one point, and for n>4 they are disjoint.

Import["C:\\Users\\adams\\Desktop\\Table1.jpg"]

Wilson’s benchmark problems

Here we used modified versions of examples from Wilson’s benchmark set. Of the 77 examples we selected 63 that involved at least 3 variables and we used quantifier-free versions of the examples. In 21 of the examples the system was not well-based and involved more than 3 variables, hence our algorithm did not apply. 7 examples did not finish in 600 seconds. Of the 35 examples for which our implementation succeeded, 29 were well-based and 6 were not well-based and in (ℝ^3). On average, (Subscript[t, CAD]) took 55% of the total time, (Subscript[t, SP]) took 34%, and (Subscript[t, ADJ]) took 11%. Five examples with the largest number of cells are given in Table 2.  All but the third example are well-based.

Import["C:\\Users\\adams\\Desktop\\Table2.jpg"]

Intersections of solids with 9(((x+y+z))^2)>(z^2)+1

Here we took the 32 3D solids that appear in Mathematica SolidData and intersected each of them with the solution set of 9 (((x+y+z))^2)>(z^2)+1. All 32 examples were well-based and in all our implementation succeeded. On average, (Subscript[t, CAD]) took 37% of the total time, (Subscript[t, SP]) took 47%, and (Subscript[t, ADJ]) took 16%. The five solids which resulted in the largest number of cells are:

Steinmetz 6-solid: 2 (x^2)+(((y-z))^2)≤2∧2 (x^2)+(((y+z))^2)≤2∧2 (y^2)+(((x-z))^2)≤2∧2 (y^2)+(((x+z))^2)≤2∧(((x-y))^2)+2 (z^2)≤2∧(((x+y))^2)+2 (z^2)≤2

Sphericon: ((x^2)+(y^2)≤(((-1+|z|))^2)∧x≥0∧-1≤z≤1)∨((x^2)+(z^2)≤(((-1+|y|))^2)∧x≤0∧-1≤y≤1)

Steinmetz 4-solid: (x^2)+(y^2)≤1∧9 (x^2)+(y^2)+8 (z^2)≤9+((164 y z)/29)∧3 (x^2)+((284 x y)/41)+7 (y^2)+((82 y z)/29)+8 (z^2)≤9+((49 x z)/10)∧3 (x^2)+7 (y^2)+((49 x z)/10)+((82 y z)/29)+8 (z^2)≤9+((284 x y)/41)

Solid capsule: (x^2)+(y^2)+(((-(1/2)+z))^2)≤1∨(x^2)+(y^2)+((((1/2)+z))^2)≤1∨(-(1/2)≤z≤(1/2)∧(x^2)+(y^2)≤1)

Reuleaux tetrahedron: (x^2)+(y^2)+((((19/31)+z))^2)≤1∧(((-(15/26)+x))^2)+(y^2)+(((-(9/44)+z))^2)≤1∧((((11/38)+x))^2)+(((-(1/2)+y))^2)+(((-(9/44)+z))^2)≤1∧((((11/38)+x))^2)+((((1/2)+y))^2)+(((-(9/44)+z))^2)≤1

Import["C:\\Users\\adams\\Desktop\\Table3.jpg"]

Thank you for listening!