[Axiom-developer] CAD package from Renaud Rioboo

[daly]

Waldek,

Here is Renaud's CAD package for your new release.

Tim

=======================================================================
)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)

)ab domain SCELL SimpleCell

-- This domain is made to work with one-dimensional semi-algebraic 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(milieu-nbe)
       while sg > 0 repeat
         lg := cons(nbe,lg)
         ld = [] => return(separe(reverse(lg),gauche,milieu))
         nbe := first(ld)
         sg := sign(milieu-nbe)
         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

)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)))

)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))
         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)
         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))
           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) := (degree(#1) <= degree(#2))
       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))