

last edited 2 years ago by test1 
1 2 3 4  
Editor: test1
Time: 2016/09/01 17:48:53 GMT+0 

Note: 
removed: CylindricalAlgebraicDecompositionUtilities \begin{spad} )ab pack CADU CylindricalAlgebraicDecompositionUtilities  CylindricalAlgebraicDecompositionUtilities(R,P) : PUB == PRIV where   Tese are some standard tools which are needed to compute with univariate  polynomials.   A gcd basis for a set of polynomials is a set of pairwise relatively prime  polynomials which all divide the original set and whose product is the  same than the product of the original set.   A square free basis for a list of polynomials is just a list  of square free polynomials which are pairwise relatively primes and have  the same roots than the original polynomials.   R : GcdDomain  P : UnivariatePolynomialCategory(R)   PUB == with  squareFreeBasis : List(P) > List(P)  ++  gcdBasis : List(P) > List(P)  ++ decompose a list of polynomials into pairwise relatively  ++ prime polynomials  gcdBasisAdd : (P,List(P)) > List(P)  ++ add one polynomial to list of pairwise relatively prime polynomials   PRIV == add   squareFreeBasis(lpols) ==  lpols = [] => []  pol := first(lpols)  sqpol := unitCanonical(squareFreePart(pol))  gcdBasis(cons(sqpol,squareFreeBasis(rest(lpols))))   gcdBasisAdd(p,lpols) ==  (degree(p) = 0) => lpols  null lpols => [unitCanonical p]  p1 := first(lpols)  g := gcd(p,p1)  (degree(g) = 0) => cons(p1,gcdBasisAdd(p,rest lpols))  p := (p exquo g)::P  p1 := (p1 exquo g)::P  basis := gcdBasisAdd(p,rest(lpols))  if degree(p1) > 0 then basis := cons(p1,basis)  gcdBasisAdd(g,basis)    gcdBasis(lpols) ==  (#lpols <= 1) => lpols  basis := gcdBasis(rest lpols)  gcdBasisAdd(first(lpols),basis) \end{spad}  SimpleCell \begin{spad} )ab domain SCELL SimpleCell   This domain is made to work with onedimensional semialgebraic cells  ie either an algebraic point, either an interval. The point is specified as  specification of an algebraic value.  SimpleCell(TheField,ThePols) : PUB == PRIV where  TheField : RealClosedField  ThePols : UnivariatePolynomialCategory(TheField)  O ==> OutputForm  B ==> Boolean  Z ==> Integer  N ==> NonNegativeInteger    VARS ==> VariationsPackage(TheField,ThePols)  VARS ==> RealPolynomialUtilitiesPackage(TheField,ThePols)  LF ==> List(TheField)   PUB == CoercibleTo(O) with   allSimpleCells : (ThePols,Symbol) > List %  allSimpleCells : (List(ThePols),Symbol) > List %  hasDimension? : % > B  samplePoint : % > TheField  stablePol : % > ThePols  variableOf : % > Symbol  separe : (LF,TheField,TheField) > LF  pointToCell : (TheField,B,Symbol) > %   PRIV == add   Rep := Record(samplePoint:TheField,  hasDim:B,  varOf:Symbol)   samplePoint(c) == c.samplePoint   stablePol(c) == error "Prout"   hasDimension?(c) == c.hasDim   variableOf(c) == c.varOf   coerce(c:%):O ==  o : O := ((c.varOf)::O) = ((c.samplePoint)::O)  brace [o,(c.hasDim)::O]   separe(liste,gauche,droite) ==  milieu : TheField := (gauche + droite) / (2::TheField)  liste = [] => [milieu]  #liste = 1 => [gauche,first(liste),droite]  nbe := first(liste)  lg :List(TheField) := []  ld :List(TheField) := rest(liste)  sg := sign(milieunbe)  while sg > 0 repeat  lg := cons(nbe,lg)  ld = [] => return(separe(reverse(lg),gauche,milieu))  nbe := first(ld)  sg := sign(milieunbe)  ld := rest(ld)  sg < 0 =>  append(separe(reverse(lg),gauche,milieu),  rest(separe(cons(nbe,ld),milieu,droite)))  newDroite := (gauche+milieu)/(2::TheField)  null lg =>  newGauche := (milieu+droite)/(2::TheField)  while newGauche >= first(ld) repeat  newGauche := (milieu+newGauche)/(2::TheField)  append([gauche,milieu],separe(ld,newGauche,droite))  while newDroite <= first(lg) repeat  newDroite := (newDroite+milieu)/(2::TheField)  newGauche := (milieu+droite)/(2::TheField)  null ld => append(separe(reverse(lg),gauche,newDroite),[milieu,droite])  while newGauche >= first(ld) repeat  newGauche := (milieu+newGauche)/(2::TheField)  append(separe(reverse(lg),gauche,newDroite),  cons(milieu,separe(ld,newGauche,droite)))    pointToCell(sp,hasDim?,varName) ==  [sp,hasDim?,varName]$Rep   allSimpleCells(p:ThePols,var:Symbol) ==  allSimpleCells([p],var)   PACK ==> CylindricalAlgebraicDecompositionUtilities(TheField,ThePols)  allSimpleCells(lp:List(ThePols),var:Symbol) ==  lp1 := gcdBasis(lp)$PACK  null(lp1) => [pointToCell(0,true,var)]  b := ("max" / [ boundOfCauchy(p)$VARS for p in lp1 ])::TheField  l := "append" / [allRootsOf(makeSUP(unitCanonical(p))) for p in lp1]  l := sort(l)  l1 := separe(l,b,b)  res : List(%) := [pointToCell(first(l1),true,var)]  l1 := rest(l1)  while not(null(l1)) repeat  res := cons(pointToCell(first(l1),false,var),res)  l1 := rest(l1)  l1 = [] => return(error "Liste vide")  res := cons(pointToCell(first(l1),true,var),res)  l1 := rest(l1)  reverse! res \end{spad}  Cell \begin{spad} )ab domain CELL Cell  Cell(TheField) : PUB == PRIV where  TheField : RealClosedField   ThePols ==> Polynomial(TheField)   O ==> OutputForm  B ==> Boolean  Z ==> Integer  N ==> NonNegativeInteger  BUP ==> SparseUnivariatePolynomial(TheField)  SCELL ==> SimpleCell(TheField,BUP)    PUB == CoercibleTo(O) with   samplePoint : % > List(TheField)  dimension : % > N  hasDimension? : (%,Symbol) > B  makeCell : List(SCELL) > %  makeCell : (SCELL,%) > %  mainVariableOf : % > Symbol  variablesOf : % > List Symbol  projection : % > Union(%,"failed")   PRIV == add   Rep := List(SCELL)   coerce(c:%):O ==  paren [sc::O for sc in c]   projection(cell) ==  null cell => error "projection: should not appear"  r := rest(cell)  null r => "failed"  r   makeCell(l:List(SCELL)) == l   makeCell(scell,toAdd) == cons(scell,toAdd)   mainVariableOf(cell) ==  null(cell) =>  error "Should not appear"  variableOf(first(cell))   variablesOf(cell) ==  null(cell) => []  cons(mainVariableOf(cell),variablesOf(rest(cell)::%))   dimension(cell) ==  null(cell) => 0  hasDimension?(first(cell)) => 1+dimension(rest(cell))  dimension(rest(cell))   hasDimension?(cell,var) ==  null(cell) =>  error "Should not appear"  sc : SCELL := first(cell)  v := variableOf(sc)  v = var => hasDimension?(sc)  v < var => false  v > var => true  error "Caca Prout"   samplePoint(cell) ==  null(cell) => []  cons(samplePoint(first(cell)),samplePoint(rest(cell))) \end{spad}  CylindricalAlgebraicDecompositionPackage \begin{spad} )ab pack CAD CylindricalAlgebraicDecompositionPackage  CylindricalAlgebraicDecompositionPackage(TheField) : PUB == PRIV where   TheField : RealClosedField   ThePols ==> Polynomial(TheField)  P ==> ThePols  BUP ==> SparseUnivariatePolynomial(TheField)  RUP ==> SparseUnivariatePolynomial(ThePols)   Z ==> Integer  N ==> NonNegativeInteger   CELL ==> Cell(TheField)  SCELL ==> SimpleCell(TheField,BUP)   PUB == with   cylindricalDecomposition: List P > List CELL  cylindricalDecomposition: (List(P),List(Symbol)) > List CELL  projectionSet: (List RUP) > List P  coefficientSet: RUP > List P  discriminantSet : List RUP > List(P)  resultantSet : List RUP > List P  principalSubResultantSet : (RUP,RUP) > List P  specialise : (List(ThePols),CELL) > List(BUP)   PRIV == add   cylindricalDecomposition(lpols) ==  lv : List(Symbol) := []  for pol in lpols repeat  ground?(pol) => "next pol"  lv := removeDuplicates(append(variables(pol),lv))  lv := reverse(sort(lv))  cylindricalDecomposition(lpols,lv)   cylindricalDecomposition(lpols,lvars) ==  lvars = [] => error("CAD: cylindricalDecomposition: empty list of vars")  mv := first(lvars)  lv := rest(lvars)  lv = [] =>  lp1 := [ univariate(pol) for pol in lpols ]  scells := allSimpleCells(lp1,mv)$SCELL  [ makeCell([scell]) for scell in scells ]  lpols1 := projectionSet [univariate(pol,mv) for pol in lpols]  previousCad := cylindricalDecomposition(lpols1,lv)  res : List(CELL) := []  for cell in previousCad repeat  lspec := specialise(lpols,cell)  scells := allSimpleCells(lspec,mv)  res := append(res,[makeCell(scell,cell) for scell in scells])  res    PACK1 ==> CylindricalAlgebraicDecompositionUtilities(ThePols,RUP)  PACK2 ==> CylindricalAlgebraicDecompositionUtilities(TheField,BUP)   specialise(lpols,cell) ==  lpols = [] => error("CAD: specialise: empty list of pols")  sp := samplePoint(cell)  vl := variablesOf(cell)  res : List(BUP) := []  for pol in lpols repeat  p1 := univariate(eval(pol,vl,sp))  zero?(p1) => return(error "Bad sepcialise")  degree(p1) = 0 => "next pol"  res := cons(p1,res)  res    coefficientSet(pol) ==  res : List(ThePols) := []  for c in coefficients(pol) repeat  ground?(c) => return(res)  res := cons(c,res)  res   SUBRES ==> SubResultantPackage(ThePols,RUP)  discriminantSet(lpols) ==  res : List(ThePols) := []  for p in lpols repeat  v := subresultantVector(p,differentiate(p))$SUBRES  not(zero?(degree(v.0))) => return(error "Bad discriminant")  d : ThePols := leadingCoefficient(v.0)  d := discriminant p  zero?(d) => return(error "Non Square Free polynomial")  if not(ground? d) then res := cons(d,res)  res   principalSubResultantSet(p,q) ==  if degree(p) < degree(q)  then (p,q) := (q,p)  if degree(p) = degree(q)  then (p,q) := (q,pseudoRemainder(p, q))  v := subresultantVector(p,q)$SUBRES  [coefficient(v.i,i) for i in 0..(((#v)2)::N)]   resultantSet(lpols) ==  res : List(ThePols) := []  laux := lpols  for p in lpols repeat  laux := rest(laux)  for q in laux repeat  r : ThePols := first(principalSubResultantSet(p,q))  r := resultant(p,q)  zero?(r) => return(error "Non relatively prime polynomials")  if not(ground? r) then res := cons(r,res)  res   projectionSet(lpols) ==  res : List(ThePols) := []  for p in lpols repeat  c := content(p)  ground?(c) => "next p"  res := cons(c,res)  lp1 := [primitivePart p for p in lpols]  f : ((RUP,RUP) > Boolean) := (x,y)+>(degree(x) <= degree(y))  lp1 := sort(f,lp1)  lsqfrb := squareFreeBasis(lp1)$PACK1  lsqfrb := sort(f,lsqfrb)  for p in lp1 repeat  res := append(res,coefficientSet(p))  res := append(res,discriminantSet(lsqfrb))  append(res,resultantSet(lsqfrb)) \end{spad} removed: )lib SCELL CELL CAD
Note: This package is currently included in FriCAS.
On Wednesday, October 12, 2005 5:29 PM Renaud Rioboo wrote:
Dear Axiom fans,
In the last few years there have been a regain of interest in some Computer Algebra users circles about Cylindrical Algebraic Decomposition (aka CAD). As some of you may know my thesis was about CAD and I am regularly receiving queries about CAD and the Axiom implementation I wrote.
I already privately gave sources or ran particular problems for some people but I have never officially released my sources because I do not think they are at production level and many work has still to be done on them.
Recently Martin Rubey thought it could be a good idea to release them and I have of course no objection on that point nor on sharing these sources with the community. If there is some interest on it you may include it into the Axiom distribution with the same permissions than the rest of the software.
While I wait for my lab to give me the necessary permissions to be able to use Axiom and export this software you may find it's sources at the unlisted url
I would like to add just a few words about that work which is a straightforward implementation of the standard papers by Arnon, Collins and McCallum? which are more than 20 years old now.
My work on Cad is 15 years old and was developped for the defense of my thesis where CAD was an sample application for real algebraic numbers manipulations. By the time I wrote it I also needed some subresultant calculations which were not in Axiom (Scratchpad by that time) and the required machinery to work with real algebraic numbers.
Along the times, real algebraic numbers were included into Axiom and subresultant calculations which are used in several places of the Axiom Library were modified to all use Lionel Ducos's package which enables a performance gain. While making these modifications to Axiom I tried to keep my CAD package uptodate with the Axiom Library and have used it as a test for other packages. It thus should compile under recent Axiom versions and should compute something that resembles a CAD.
Don't expect this package to be the absolute CAD program, I never found time to further work on it in order to include many enhancements which are present in other CAD packages. While Axiom versions evolved I had to remove what I thought were enhancements and remove some support for generalization. As of today it only includes one projection method which is the McCallum? projection and no reconstruction nor adjacency's are taken into account.
However it did reasonabily compare with more recent techniques such as Rational Univariate Representation (aka RUR) for simple cases though it could not produce results for some more complicated cases. By the time this comparison was made Axiom had some severe memory restrictions which seem to have disappear now. I thus think that many objections to using CAD could now be revised.
For hard problems this package should thus not be worse than any other CAD package except that you cannot fall into an optimzation drawback :) There is no optimization !
I will of course try to maintain the package to the best of my ability. One thing I should do is include comments describing the different functions, but even that requires going into the code and the algorithms which are quite old for me now.
Let me know if there is some interest on it.
Best regards,
Renaud
Test Input
Ran := RECLOS(FRAC INT)
(1) 
Cad := CAD(Ran)
(2) 
p1 : POLY(Ran) := x^2+y^21
(3) 
p2 : POLY(Ran) := yx^2
(4) 
lpols : List(POLY(Ran)) := [p1,p2]
(5) 
cad := cylindricalDecomposition(lpols)$Cad
(6) 
precision(30)
(7) 
ls := [ samplePoint cell for cell in cad]
(8) 
lf := [ [relativeApprox(coord,2^(precision()))::Float for coord in point] for point in ls]
(9) 
lp := [ point(term::List(DoubleFloat))$Point(DoubleFloat) for term in lf ]
(10) 
[ sign(eval(p1,variables(p1), samplePoint(cell))::Ran) for cell in cad ]
(11) 