axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / QALGSET.spad

)abbrev domain QALGSET QuasiAlgebraicSet
++ Author:  William Sit
++ Date Created: March 13, 1992
++ Date Last Updated: June 12, 1992 
++ References:William Sit, "An Algorithm for Parametric Linear Systems"
++            J. Sym. Comp., April, 1992
++ Description:
++ \spadtype{QuasiAlgebraicSet} constructs a domain representing
++ quasi-algebraic sets, which is the intersection of a Zariski
++ closed set, defined as the common zeros of a given list of
++ polynomials (the defining polynomials for equations), and a principal
++ Zariski open set, defined as the complement of the common
++ zeros of a polynomial f (the defining polynomial for the inequation).
++ This domain provides simplification of a user-given representation
++ using groebner basis computations.
++ There are two simplification routines: the first function
++ \spadfun{idealSimplify}  uses groebner
++ basis of ideals alone, while the second, \spadfun{simplify} uses both
++ groebner basis and factorization.  The resulting defining equations L
++ always form a groebner basis, and the resulting defining
++ inequation f is always reduced.  The function \spadfun{simplify} may
++ be applied several times if desired.   A third simplification
++ routine \spadfun{radicalSimplify} is provided in
++ \spadtype{QuasiAlgebraicSet2}  for comparison study only,
++ as it is inefficient compared to the other two, as well as is
++ restricted to only certain coefficient domains.  For detail analysis
++ and a comparison of the three methods, please consult the reference
++ cited.
++
++ A polynomial function q defined on the quasi-algebraic set
++ is equivalent to its reduced form with respect to L.  While
++ this may be obtained using the usual normal form
++ algorithm, there is no canonical form for q.
++
++ The ordering in groebner basis computation is determined by
++ the data type of the input polynomials.  If it is possible
++ we suggest to use refinements of total degree orderings.

QuasiAlgebraicSet(R, Var,Expon,Dpoly) : SIG == CODE where
  R : GcdDomain
  Var : OrderedSet
  Expon : OrderedAbelianMonoidSup
  Dpoly : PolynomialCategory(R,Expon,Var)

  NNI      ==> NonNegativeInteger
  newExpon ==> Product(NNI,Expon)
  newPoly  ==> PolynomialRing(R,newExpon)
  Ex       ==> OutputForm
  mrf      ==> MultivariateFactorize(Var,Expon,R,Dpoly)
  Status   ==> Union(Boolean,"failed") -- empty or not, or don't know
 
  SIG ==> Join(SetCategory, CoercibleTo OutputForm) with
     --- should be Object instead of SetCategory, bug in LIST Object ---
     --- equality is not implemented ---

     empty : () -> $
       ++ empty() returns the empty quasi-algebraic set

     quasiAlgebraicSet : (List Dpoly, Dpoly) -> $
       ++ quasiAlgebraicSet(pl,q) returns the quasi-algebraic set
       ++ with defining equations p = 0 for p belonging to the list pl, and
       ++ defining inequation q ^= 0.

     status : $ -> Status
       ++ status(s) returns true if the quasi-algebraic set is empty,
       ++ false if it is not, and "failed" if not yet known

     setStatus : ($, Status) -> $
       ++ setStatus(s,t) returns the same representation for s, but
       ++ asserts the following: if t is true, then s is empty,
       ++ if t is false, then s is non-empty, and if t = "failed",
       ++ then no assertion is made (that is, "don't know").
       ++ Note: for internal use only, with care.

     empty? : $ -> Boolean
       ++ empty?(s) returns
       ++ true if the quasialgebraic set s has no points,
       ++ and false otherwise.

     definingEquations : $ -> List Dpoly
       ++ definingEquations(s) returns a list of defining polynomials
       ++ for equations, that is, for the Zariski closed part of s.

     definingInequation : $ -> Dpoly
       ++ definingInequation(s) returns a single defining polynomial for the
       ++ inequation, that is, the Zariski open part of s.

     idealSimplify : $ -> $
       ++ idealSimplify(s) returns a different and presumably simpler
       ++ representation of s with the defining polynomials for the
       ++ equations
       ++ forming a groebner basis, and the defining polynomial for the
       ++ inequation reduced with respect to the basis, using Buchberger's
       ++ algorithm.

     if (R has EuclideanDomain) and (R has CharacteristicZero) then

       simplify : $ -> $
         ++ simplify(s) returns a different and presumably simpler
         ++ representation of s with the defining polynomials for the
         ++ equations
         ++ forming a groebner basis, and the defining polynomial for the
         ++ inequation reduced with respect to the basis, using a heuristic
         ++ algorithm based on factoring.

  CODE ==> add

     Rep := Record(status:Status,zero:List Dpoly, nzero:Dpoly)
     x:$
 
     import GroebnerPackage(R,Expon,Var,Dpoly)
     import GroebnerPackage(R,newExpon,Var,newPoly)
     import GroebnerInternalPackage(R,Expon,Var,Dpoly)
 
                       ----  Local Functions  ----
 
     minset   : List List Dpoly -> List List Dpoly
     overset? : (List Dpoly, List List Dpoly) -> Boolean
     npoly    : Dpoly            ->  newPoly
     oldpoly  : newPoly          ->  Union(Dpoly,"failed")
 
 
     if (R has EuclideanDomain) and (R has CharacteristicZero) then
       factorset (y:Dpoly):List Dpoly ==
         ground? y => []
         [j.factor for j in factors factor$mrf  y]
 
       simplify x ==
         if x.status case "failed" then
           x:=quasiAlgebraicSet(zro:=groebner x.zero, redPol(x.nzero,zro))
         (pnzero:=x.nzero)=0 => empty()
         nzro:=factorset pnzero
         mset:=minset [factorset p for p in x.zero]
         mset:=[setDifference(s,nzro) for s in mset]
         zro:=groebner [*/s for s in mset]
         member? (1$Dpoly, zro) => empty()
         [x.status, zro, primitivePart redPol( */nzro, zro)]
 
     npoly(f:Dpoly) : newPoly ==
       zero? f => 0
       monomial(leadingCoefficient f,makeprod(0,degree f))$newPoly +
             npoly(reductum f)
 
     oldpoly(q:newPoly) : Union(Dpoly,"failed") ==
       q=0$newPoly => 0$Dpoly
       dq:newExpon:=degree q
       n:NNI:=selectfirst (dq)
       n^=0 => "failed"
       ((g:=oldpoly reductum q) case "failed") => "failed"
       monomial(leadingCoefficient q,selectsecond dq)$Dpoly + (g::Dpoly)
 
     coerce x ==
       x.status = true => "Empty"::Ex
       bracket [[hconcat(f::Ex, " = 0"::Ex) for f in x.zero ]::Ex,
                 hconcat( x.nzero::Ex, " != 0"::Ex)]
 
     empty? x ==
       if x.status case "failed" then x:=idealSimplify x
       x.status :: Boolean
 
     empty() == [true::Status, [1$Dpoly], 0$Dpoly]

     status x == x.status

     setStatus(x,t) == [t,x.zero,x.nzero]

     definingEquations x == x.zero

     definingInequation x == x.nzero

     quasiAlgebraicSet(z0,n0) == ["failed", z0, n0]
 
     idealSimplify x ==
       x.status case Boolean => x
       z0:= x.zero
       n0:= x.nzero
       empty? z0 => [false, z0, n0]
       member? (1$Dpoly, z0) => empty()
       tp:newPoly:=(monomial(1,makeprod(1,0$Expon))$newPoly * npoly n0)-1
       ngb:=groebner concat(tp, [npoly g for g in z0])
       member? (1$newPoly, ngb) => empty()
       gb:List Dpoly:=nil
       while not empty? ngb repeat
         if ((f:=oldpoly ngb.first) case Dpoly) then gb:=concat(f, gb)
         ngb:=ngb.rest
       [false::Status, gb, primitivePart redPol(n0, gb)]
 
 
     minset lset ==
       empty? lset => lset
       [s for s  in lset | ^(overset?(s,lset))]
 
     overset?(p,qlist) ==
       empty? qlist => false
       or/[(brace$(Set Dpoly) q) <$(Set Dpoly) (brace$(Set Dpoly) p) _
           for q in qlist]