path: root/src/algebra/qalgset.spad.pamphlet

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra qalgset.spad}
\author{William Sit}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{domain QALGSET QuasiAlgebraicSet}
<<domain QALGSET QuasiAlgebraicSet>>=
)abbrev domain QALGSET QuasiAlgebraicSet
++ Author:  William Sit
++ Date Created: March 13, 1992
++ Date Last Updated: June 12, 1992 
++ Basic Operations:
++ Related Constructors:GroebnerPackage
++ See Also: QuasiAlgebraicSet2
++ AMS Classifications:
++ Keywords: Zariski closed sets, quasi-algebraic sets
++ 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) : C == T
 where
   R         :  GcdDomain
   Expon     :  OrderedAbelianMonoidSup
   Var       :  OrderedSet
   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
 
   C == 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.
   T  == 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)
       not zero? n => "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 | not (overset?(s,lset))]
 
     overset?(p,qlist) ==
       empty? qlist => false
       or/[part?(brace(q)$Set(Dpoly), brace(p)$Set(Dpoly))$Set(Dpoly) 
             for q in qlist]

@
\section{package QALGSET2 QuasiAlgebraicSet2}
<<package QALGSET2 QuasiAlgebraicSet2>>=
)abbrev package QALGSET2 QuasiAlgebraicSet2
++ Author:  William Sit
++ Date Created: March 13, 1992
++ Date Last Updated: June 12, 1992
++ Basic Operations:
++ Related Constructors:GroebnerPackage, IdealDecompositionPackage,
++                      PolynomialIdeals
++ See Also: QuasiAlgebraicSet
++ AMS Classifications:
++ Keywords: Zariski closed sets, quasi-algebraic sets
++ References:William Sit, "An Algorithm for Parametric Linear Systems"
++            J. Sym. Comp., April, 1992
++ Description:
++ \spadtype{QuasiAlgebraicSet2} adds a function \spadfun{radicalSimplify}
++ which uses \spadtype{IdealDecompositionPackage} to simplify
++ the representation of a quasi-algebraic set.  A quasi-algebraic set
++ 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).
++ Quasi-algebraic sets are implemented in the domain
++ \spadtype{QuasiAlgebraicSet}, where two simplification routines are
++ provided:
++ \spadfun{idealSimplify} and \spadfun{simplify}.
++ The function
++ \spadfun{radicalSimplify} is added
++ for comparison study only.  Because the domain
++ \spadtype{IdealDecompositionPackage} provides facilities for
++ computing with radical ideals, it is necessary to restrict
++ the ground ring to the domain \spadtype{Fraction Integer},
++ and the polynomial ring to be of type
++ \spadtype{DistributedMultivariatePolynomial}.
++ The routine \spadfun{radicalSimplify} uses these to compute groebner
++ basis of radical ideals and
++ is inefficient and restricted when compared to the
++ two in \spadtype{QuasiAlgebraicSet}.
QuasiAlgebraicSet2(vl,nv) : C == T where
   vl       :   List Symbol
   nv       :   NonNegativeInteger
   R        ==> Integer
   F        ==> Fraction R
   Var      ==> OrderedVariableList vl
   NNI      ==> NonNegativeInteger
   Expon    ==> DirectProduct(nv,NNI)
   Dpoly    ==> DistributedMultivariatePolynomial(vl,F)
   QALG     ==> QuasiAlgebraicSet(F, Var, Expon, Dpoly)
   newExpon ==> DirectProduct(#newvl, NNI)
   newPoly  ==> DistributedMultivariatePolynomial(newvl,F)
   newVar   ==> OrderedVariableList newvl
   Status   ==> Union(Boolean,"failed") -- empty or not, or don't know
 
   C == with
     radicalSimplify:QALG -> QALG
       ++ radicalSimplify(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
       ++ using groebner basis of radical ideals
   T  == add
                ----  Local Functions  ----
     ts:=new()$Symbol
     newvl:=concat(ts, vl)
     tv:newVar:=(variable ts)::newVar
     npoly         :     Dpoly            ->  newPoly
     oldpoly       :     newPoly          ->  Union(Dpoly,"failed")
     f             :     Var              ->  newPoly
     g             :     newVar           ->  Dpoly
 
     import PolynomialIdeals(F,newExpon,newVar,newPoly)
     import GroebnerPackage(F,Expon,Var,Dpoly)
     import GroebnerPackage(F,newExpon,newVar,newPoly)
     import IdealDecompositionPackage(newvl,#newvl)
     import QuasiAlgebraicSet(F, Var, Expon, Dpoly)
     import PolynomialCategoryLifting(Expon,Var,F,Dpoly,newPoly)
     import PolynomialCategoryLifting(newExpon,newVar,F,newPoly,Dpoly)
     f(v:Var):newPoly ==
       variable((convert v)@Symbol)@Union(newVar,"failed")::newVar
         ::newPoly
     g(v:newVar):Dpoly ==
       v = tv => 0
       variable((convert v)@Symbol)@Union(Var,"failed")::Var::Dpoly
 
     npoly(p:Dpoly) : newPoly ==  map(f #1, #1::newPoly, p)
 
     oldpoly(q:newPoly) : Union(Dpoly,"failed") ==
       (x:=mainVariable q) case "failed" => (leadingCoefficient q)::Dpoly
       (x::newVar = tv) => "failed"
       map(g #1,#1::Dpoly, q)
 
     radicalSimplify x ==
       status(x)$QALG = true => x     -- x is empty
       z0:=definingEquations x
       n0:=definingInequation x
       t:newPoly:= coerce(tv)$newPoly
       tp:newPoly:= t * (npoly n0) - 1$newPoly
       gen:List newPoly:= concat(tp, [npoly g for g in z0])
       id:=ideal gen
       ngb:=generators radical(id)
       member? (1$newPoly, ngb) => empty()$QALG
       gb:List Dpoly:=nil
       while not empty? ngb repeat
         if ((k:=oldpoly ngb.first) case Dpoly) then gb:=concat(k, gb)
         ngb:=ngb.rest
       y:=quasiAlgebraicSet(gb, primitivePart normalForm(n0, gb))
       setStatus(y,false::Status)

@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>

<<domain QALGSET QuasiAlgebraicSet>>
<<package QALGSET2 QuasiAlgebraicSet2>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}