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

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra regset.spad}
\author{Marc Moreno Maza}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{category RSETCAT RegularTriangularSetCategory}
<<category RSETCAT RegularTriangularSetCategory>>=
)abbrev category RSETCAT RegularTriangularSetCategory
++ Author: Marc Moreno Maza
++ Date Created: 09/03/1998
++ Date Last Updated: 12/15/1998
++ Basic Functions:
++ Related Constructors:
++ Also See: essai Graphisme
++ AMS Classifications:
++ Keywords: polynomial, multivariate, ordered variables set
++ Description:
++ The category of regular triangular sets, introduced under
++ the name regular chains in [1] (and other papers).
++ In [3] it is proved that regular triangular sets and towers of simple
++ extensions of a field are equivalent notions.
++ In the following definitions, all polynomials and ideals
++ are taken from the polynomial ring \spad{k[x1,...,xn]} where \spad{k}
++ is the fraction field of \spad{R}.
++ The triangular set \spad{[t1,...,tm]} is regular 
++ iff for every \spad{i} the initial of \spad{ti+1} is invertible
++ in the tower of simple extensions associated with \spad{[t1,...,ti]}.
++ A family \spad{[T1,...,Ts]} of regular triangular sets 
++ is a split of Kalkbrener of a given ideal \spad{I}
++ iff the radical of \spad{I} is equal to the intersection
++ of the radical ideals generated by the saturated ideals 
++ of the \spad{[T1,...,Ti]}.
++ A family \spad{[T1,...,Ts]} of regular triangular sets 
++ is a split of Kalkbrener of a given triangular set \spad{T}
++ iff it is a split of Kalkbrener of the saturated ideal of \spad{T}.
++ Let \spad{K} be an algebraic closure of \spad{k}.
++ Assume that \spad{V} is finite with cardinality
++ \spad{n} and let \spad{A} be the affine space \spad{K^n}.
++ For a regular triangular set \spad{T} let denote by \spad{W(T)} the
++ set of regular zeros of \spad{T}.
++ A family \spad{[T1,...,Ts]} of regular triangular sets 
++ is a split of Lazard of a given subset \spad{S} of \spad{A}
++ iff the union of the \spad{W(Ti)} contains \spad{S} and
++ is contained in the closure of \spad{S} (w.r.t. Zariski topology).
++ A family \spad{[T1,...,Ts]} of regular triangular sets 
++ is a split of Lazard of a given triangular set \spad{T}
++ if it is a split of Lazard of \spad{W(T)}.
++ Note that if \spad{[T1,...,Ts]} is a split of Lazard of 
++ \spad{T} then it is also a split of Kalkbrener of \spad{T}.
++ The converse is false. 
++ This category provides operations related to both kinds of
++ splits, the former being related to ideals decomposition whereas 
++ the latter deals with varieties decomposition.
++ See the example illustrating the \spadtype{RegularTriangularSet} constructor
++ for more explanations about decompositions by means of regular triangular sets. \newline
++ References :
++  [1] M. KALKBRENER "Three contributions to elimination theory"
++      Phd Thesis, University of Linz, Austria, 1991.
++  [2] M. KALKBRENER "Algorithmic properties of polynomial rings"
++      Journal of Symbol. Comp. 1998
++  [3] P. AUBRY, D. LAZARD and M. MORENO MAZA "On the Theories
++      of Triangular Sets" Journal of Symbol. Comp. (to appear)
++  [4] M. MORENO MAZA "A new algorithm for computing triangular
++      decomposition of algebraic varieties" NAG Tech. Rep. 4/98.
++ Version: 2

RegularTriangularSetCategory(R:GcdDomain, E:OrderedAbelianMonoidSup,_
  V:OrderedSet,P:RecursivePolynomialCategory(R,E,V)):
         Category == 
   TriangularSetCategory(R,E,V,P) with 

     purelyAlgebraic?: (P,$) -> Boolean
         ++ \spad{purelyAlgebraic?(p,ts)} returns \spad{true} iff every
         ++ variable of \spad{p} is algebraic w.r.t. \spad{ts}.
     purelyTranscendental? : (P,$) -> Boolean
         ++ \spad{purelyTranscendental?(p,ts)} returns \spad{true} iff every
         ++ variable of \spad{p} is not algebraic w.r.t. \spad{ts}
     algebraicCoefficients?  : (P,$) -> Boolean
         ++ \spad{algebraicCoefficients?(p,ts)} returns \spad{true} iff every
         ++ variable of \spad{p} which is not the main one of \spad{p}
         ++ is algebraic w.r.t. \spad{ts}.
     purelyAlgebraic?: $ -> Boolean
         ++ \spad{purelyAlgebraic?(ts)} returns true iff for every algebraic 
         ++ variable \spad{v} of \spad{ts} we have 
         ++ \spad{algebraicCoefficients?(t_v,ts_v_-)} where \spad{ts_v}
         ++ is \axiomOpFrom{select}{TriangularSetCategory}(ts,v) and \spad{ts_v_-} is
         ++ \axiomOpFrom{collectUnder}{TriangularSetCategory}(ts,v).
     purelyAlgebraicLeadingMonomial?: (P, $) -> Boolean
         ++ \spad{purelyAlgebraicLeadingMonomial?(p,ts)} returns true iff
         ++ the main variable of any non-constant iterarted initial 
         ++ of \spad{p} is algebraic w.r.t. \spad{ts}. 
     invertibleElseSplit? : (P,$) -> Union(Boolean,List $)
         ++ \spad{invertibleElseSplit?(p,ts)} returns \spad{true} (resp.
         ++ \spad{false}) if \spad{p} is invertible in the tower 
         ++ associated with \spad{ts} or returns a split of Kalkbrener 
         ++ of \spad{ts}.
     invertible? : (P,$) -> List Record(val : Boolean, tower : $)
         ++ \spad{invertible?(p,ts)} returns \spad{lbwt} where \spad{lbwt.i}
         ++ is the result of \spad{invertibleElseSplit?(p,lbwt.i.tower)} and
         ++ the list of the \spad{(lqrwt.i).tower} is a split of Kalkbrener of \spad{ts}.
     invertible?: (P,$) -> Boolean
         ++ \spad{invertible?(p,ts)} returns true iff \spad{p} is invertible
         ++ in the tower associated with \spad{ts}.
     invertibleSet: (P,$) -> List $
         ++ \spad{invertibleSet(p,ts)} returns a split of Kalkbrener of the
         ++ quotient ideal of the ideal \axiom{I}  by \spad{p} where \spad{I} is 
         ++ the radical of saturated of \spad{ts}.
     lastSubResultantElseSplit: (P, P, $) -> Union(P,List $)
         ++ \spad{lastSubResultantElseSplit(p1,p2,ts)} returns either
         ++ \spad{g} a quasi-monic gcd of \spad{p1} and \spad{p2} w.r.t.
         ++ the \spad{ts} or a split of Kalkbrener of \spad{ts}.
         ++ This assumes that \spad{p1} and \spad{p2} have the same maim
         ++ variable and that this variable is greater that any variable
         ++ occurring in \spad{ts}. 
     lastSubResultant: (P, P, $) -> List Record(val : P, tower : $)
         ++ \spad{lastSubResultant(p1,p2,ts)} returns \spad{lpwt} such that
         ++ \spad{lpwt.i.val} is a quasi-monic gcd of \spad{p1} and \spad{p2} 
         ++ w.r.t. \spad{lpwt.i.tower}, for every \spad{i}, and such 
         ++ that the list of the \spad{lpwt.i.tower} is a split of Kalkbrener of 
         ++ \spad{ts}. Moreover, if \spad{p1} and \spad{p2} do not
         ++ have a non-trivial gcd w.r.t. \spad{lpwt.i.tower} then \spad{lpwt.i.val} 
         ++ is the resultant of these polynomials w.r.t. \spad{lpwt.i.tower}.
         ++ This assumes that \spad{p1} and \spad{p2} have the same maim
         ++ variable and that this variable is greater that any variable
         ++ occurring in \spad{ts}. 
     squareFreePart: (P,$) -> List Record(val : P, tower : $)
         ++ \spad{squareFreePart(p,ts)} returns \spad{lpwt} such that
         ++ \spad{lpwt.i.val} is a square-free polynomial
         ++ w.r.t. \spad{lpwt.i.tower}, this polynomial being associated with \spad{p}
         ++ modulo \spad{lpwt.i.tower}, for every \spad{i}. Moreover,
         ++ the list of the \spad{lpwt.i.tower} is a split
         ++ of Kalkbrener of \spad{ts}. 
         ++ WARNING: This assumes that \spad{p} is a non-constant polynomial such that
         ++ if \spad{p} is added to \spad{ts}, then the resulting set is a
         ++ regular triangular set. 
     intersect: (P,$) -> List $
         ++ \spad{intersect(p,ts)} returns the same as 
         ++ \spad{intersect([p],ts)}
     intersect: (List P, $) -> List $
         ++ \spad{intersect(lp,ts)} returns \spad{lts} a split of Lazard
         ++ of the intersection of the affine variety associated 
         ++ with \spad{lp} and the regular zero set of \spad{ts}.
     intersect: (List P, List $) -> List $
         ++ \spad{intersect(lp,lts)} returns the same as
         ++ \spad{concat([intersect(lp,ts) for ts in lts])|}
     intersect: (P, List $) -> List $
         ++ \spad{intersect(p,lts)} returns the same as
         ++ \spad{intersect([p],lts)}
     augment: (P,$) -> List $
         ++ \spad{augment(p,ts)} assumes that \spad{p} is a non-constant
         ++ polynomial whose main variable is greater than any variable
         ++ of \spad{ts}. This operation assumes also that if \spad{p} is 
         ++ added to \spad{ts} the resulting set, say \spad{ts+p}, is a
         ++ regular triangular set. Then it returns a split of Kalkbrener
         ++ of \spad{ts+p}. This may not be \spad{ts+p} itself, if for
         ++ instance \spad{ts+p} is required to be square-free.
     augment: (P,List $) -> List $
         ++ \spad{augment(p,lts)} returns the same as
         ++ \spad{concat([augment(p,ts) for ts in lts])}
     augment: (List P,$) -> List $   
         ++ \spad{augment(lp,ts)} returns \spad{ts} if \spad{empty? lp},
         ++ \spad{augment(p,ts)} if \spad{lp = [p]}, otherwise
         ++ \spad{augment(first lp, augment(rest lp, ts))}
     augment: (List P,List $) -> List $
         ++ \spad{augment(lp,lts)} returns the same as 
         ++ \spad{concat([augment(lp,ts) for ts in lts])}
     internalAugment: (P, $) -> $
         ++ \spad{internalAugment(p,ts)} assumes that \spad{augment(p,ts)}
         ++ returns a singleton and returns it.
     internalAugment: (List P, $) -> $
         ++ \spad{internalAugment(lp,ts)} returns \spad{ts} if \spad{lp}
         ++ is empty otherwise returns 
         ++ \spad{internalAugment(rest lp, internalAugment(first lp, ts))}
     extend: (P,$) -> List $
         ++ \spad{extend(p,ts)} assumes that \spad{p} is a non-constant
         ++ polynomial whose main variable is greater than any variable
         ++ of \spad{ts}. Then it returns a split of Kalkbrener
         ++ of \spad{ts+p}. This may not be \spad{ts+p} itself, if for
         ++ instance \spad{ts+p} is not a regular triangular set.
     extend: (P, List $)  -> List $
         ++ \spad{extend(p,lts)} returns the same as 
         ++ \spad{concat([extend(p,ts) for ts in lts])|}
     extend: (List P,$) -> List $   
         ++ \spad{extend(lp,ts)} returns \spad{ts} if \spad{empty? lp}
         ++ \spad{extend(p,ts)} if \spad{lp = [p]} else
         ++ \spad{extend(first lp, extend(rest lp, ts))}
     extend: (List P,List $) -> List $
         ++ \spad{extend(lp,lts)} returns the same as 
         ++ \spad{concat([extend(lp,ts) for ts in lts])|}
     zeroSetSplit: (List P, Boolean) -> List $
         ++ \spad{zeroSetSplit(lp,clos?)} returns \spad{lts} a split of Kalkbrener 
         ++ of the radical ideal associated with \spad{lp}.
         ++ If \spad{clos?} is false, it is also a decomposition of the
         ++ variety associated with \spad{lp} into the regular zero set of the \spad{ts} in \spad{lts}
         ++ (or, in other words, a split of Lazard of this variety).
         ++ See the example illustrating the \spadtype{RegularTriangularSet} constructor
         ++ for more explanations about decompositions by means of regular triangular sets. 

 add

     NNI ==> NonNegativeInteger
     INT ==> Integer
     LP ==> List P
     PWT ==> Record(val : P, tower : $)
     LpWT ==> Record(val : (List P), tower : $)
     Split ==> List $
     pack ==> PolynomialSetUtilitiesPackage(R,E,V,P)

     purelyAlgebraic?(p: P, ts: $): Boolean ==
       ground? p => true
       not algebraic?(mvar(p),ts) => false
       algebraicCoefficients?(p,ts)

     purelyTranscendental?(p:P,ts:$): Boolean  ==
       empty? ts => true
       lv : List V := variables(p)$P
       while (not empty? lv) and (not algebraic?(first(lv),ts)) repeat lv := rest lv
       empty? lv

     purelyAlgebraicLeadingMonomial?(p: P, ts: $): Boolean  ==
       ground? p => true
       algebraic?(mvar(p),ts) and purelyAlgebraicLeadingMonomial?(init(p), ts)

     algebraicCoefficients?(p:P,ts:$): Boolean  ==
       ground? p => true
       (not ground? init(p)) and not (algebraic?(mvar(init(p)),ts)) => false
       algebraicCoefficients?(init(p),ts) =>
         ground? tail(p) => true
         mvar(tail(p)) = mvar(p) => 
           algebraicCoefficients?(tail(p),ts)
         algebraic?(mvar(tail(p)),ts) => 
           algebraicCoefficients?(tail(p),ts)
         false
       false

     if V has Finite
     then
      purelyAlgebraic?(ts: $): Boolean ==
         empty? ts => true
         size()$V = #ts => true
         lp: LP := sort(infRittWu?,members(ts))
         i: NonNegativeInteger := size()$V
         for p in lp repeat
           v: V := mvar(p)
           (i = (lookup(v)$V)::NNI) => 
             i := subtractIfCan(i,1)::NNI
           univariate?(p)$pack => 
             i := subtractIfCan(i,1)::NNI
           not algebraicCoefficients?(p,collectUnder(ts,v)) =>
             return false
           i := subtractIfCan(i,1)::NNI
         true
           
     else

       purelyAlgebraic?(ts: $): Boolean ==
         empty? ts => true
         v: V := mvar(ts)
         p: P := select(ts,v)::P
         ts := collectUnder(ts,v)
         empty? ts => univariate?(p)$pack
         not purelyAlgebraic?(ts) => false
         algebraicCoefficients?(p,ts)

     augment(p:P,lts:List $) ==
       toSave: Split := []
       while not empty? lts repeat
         ts := first lts
         lts := rest lts
         toSave := concat(augment(p,ts),toSave)
       toSave

     augment(lp:LP,ts:$) ==
       toSave: Split := [ts]
       empty? lp => toSave
       lp := sort(infRittWu?,lp)
       while not empty? lp repeat
         p := first lp
         lp := rest lp
         toSave := augment(p,toSave)
       toSave

     augment(lp:LP,lts:List $) ==
       empty? lp => lts
       toSave: Split := []
       while not empty? lts repeat
         ts := first lts
         lts := rest lts
         toSave := concat(augment(lp,ts),toSave)
       toSave    

     extend(p:P,lts:List $) ==
       toSave : Split := []
       while not empty? lts repeat
         ts := first lts
         lts := rest lts
         toSave := concat(extend(p,ts),toSave)
       toSave

     extend(lp:LP,ts:$) ==
       toSave: Split := [ts]
       empty? lp => toSave
       lp := sort(infRittWu?,lp)
       while not empty? lp repeat
         p := first lp
         lp := rest lp
         toSave := extend(p,toSave)
       toSave

     extend(lp:LP,lts:List $) ==
       empty? lp => lts
       toSave: Split := []
       while not empty? lts repeat
         ts := first lts
         lts := rest lts
         toSave := concat(extend(lp,ts),toSave)
       toSave    

     intersect(lp:LP,lts:List $): List $  ==
       -- A VERY GENERAL default algorithm
       (empty? lp) or (empty? lts) => lts
       lp := [primitivePart(p) for p in lp]
       lp := removeDuplicates lp
       lp := remove(zero?,lp)
       any?(ground?,lp) => []
       toSee: List LpWT := [[lp,ts]$LpWT for ts in lts]
       toSave: List $ := []
       lp: LP
       p: P
       ts: $
       lus: List $
       while (not empty? toSee) repeat
         lpwt := first toSee; toSee := rest toSee
         lp := lpwt.val; ts := lpwt.tower
         empty? lp => toSave := cons(ts, toSave)
         p := first lp;  lp := rest lp
         lus := intersect(p,ts)
         toSee := concat([[lp,us]$LpWT for us in lus], toSee)
       toSave

     intersect(lp: LP,ts: $): List $  ==
       intersect(lp,[ts])

     intersect(p: P,lts: List $): List $  ==
       intersect([p],lts)

@
\section{package QCMPACK QuasiComponentPackage}
<<package QCMPACK QuasiComponentPackage>>=
)abbrev package QCMPACK QuasiComponentPackage
++ Author: Marc Moreno Maza
++   marc@nag.co.uk
++ Date Created: 08/30/1998
++ Date Last Updated: 12/16/1998
++ Basic Functions:
++ Related Constructors:
++ Also See: `tosedom.spad'
++ AMS Classifications:
++ Keywords:
++ Description: 
++ A package for removing redundant quasi-components and redundant
++ branches when decomposing a variety by means of quasi-components
++ of regular triangular sets. \newline
++ References :
++  [1] D. LAZARD "A new method for solving algebraic systems of 
++      positive dimension" Discr. App. Math. 33:147-160,1991
++  [2] M. MORENO MAZA "Calculs de pgcd au-dessus des tours
++      d'extensions simples et resolution des systemes d'equations
++      algebriques" These, Universite P.etM. Curie, Paris, 1997.
++  [3] M. MORENO MAZA "A new algorithm for computing triangular
++      decomposition of algebraic varieties" NAG Tech. Rep. 4/98.
++ Version: 3. 

QuasiComponentPackage(R,E,V,P,TS): Exports == Implementation where

  R : GcdDomain
  E : OrderedAbelianMonoidSup
  V : OrderedSet
  P : RecursivePolynomialCategory(R,E,V)
  TS : RegularTriangularSetCategory(R,E,V,P)
  N ==> NonNegativeInteger
  Z ==> Integer
  B ==> Boolean
  S ==> String
  LP ==> List P
  PtoP ==> P -> P
  PS ==> GeneralPolynomialSet(R,E,V,P)
  PWT ==> Record(val : P, tower : TS)
  BWT ==> Record(val : Boolean, tower : TS)
  LpWT ==> Record(val : (List P), tower : TS)
  Branch ==> Record(eq: List P, tower: TS, ineq: List P)
  UBF ==> Union(Branch,"failed")
  Split ==> List TS
  Key ==>  Record(left:TS, right:TS)
  Entry ==> Boolean
  H ==> TabulatedComputationPackage(Key, Entry)
  polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,P)

  Exports ==  with
     startTable!: (S,S,S) -> Void
         ++ \axiom{startTableGcd!(s1,s2,s3)} 
         ++ is an internal subroutine, exported only for developement.
     stopTable!: () -> Void
         ++ \axiom{stopTableGcd!()} 
         ++ is an internal subroutine, exported only for developement.
     supDimElseRittWu?: (TS,TS) -> Boolean
         ++ \axiom{supDimElseRittWu(ts,us)} returns true iff \axiom{ts}
         ++ has less elements than \axiom{us} otherwise if \axiom{ts}
         ++ has higher rank than \axiom{us} w.r.t. Riit and Wu ordering.
     algebraicSort: Split -> Split
         ++ \axiom{algebraicSort(lts)} sorts \axiom{lts} w.r.t 
         ++ \axiomOpFrom{supDimElseRittWu?}{QuasiComponentPackage}.
     moreAlgebraic?: (TS,TS) -> Boolean
         ++ \axiom{moreAlgebraic?(ts,us)} returns false iff \axiom{ts}
         ++ and \axiom{us} are both empty, or \axiom{ts}
         ++ has less elements than \axiom{us}, or some variable is
         ++ algebraic w.r.t. \axiom{us} and is not w.r.t. \axiom{ts}.
     subTriSet?: (TS,TS) -> Boolean
         ++ \axiom{subTriSet?(ts,us)} returns true iff \axiom{ts} is
         ++ a sub-set of \axiom{us}.
     subPolSet?: (LP, LP) -> Boolean
         ++ \axiom{subPolSet?(lp1,lp2)} returns true iff \axiom{lp1} is
         ++ a sub-set of \axiom{lp2}.
     internalSubPolSet?: (LP, LP) -> Boolean
         ++ \axiom{internalSubPolSet?(lp1,lp2)} returns true iff \axiom{lp1} is
         ++ a sub-set of \axiom{lp2} assuming that these lists are sorted
         ++ increasingly w.r.t. \axiomOpFrom{infRittWu?}{RecursivePolynomialCategory}.
     internalInfRittWu?: (LP, LP) -> Boolean
         ++ \axiom{internalInfRittWu?(lp1,lp2)}
         ++ is an internal subroutine, exported only for developement.
     infRittWu?: (LP, LP) -> Boolean
         ++ \axiom{infRittWu?(lp1,lp2)}
         ++ is an internal subroutine, exported only for developement.
     internalSubQuasiComponent?: (TS,TS) -> Union(Boolean,"failed")
         ++ \axiom{internalSubQuasiComponent?(ts,us)} returns a boolean \spad{b} value
         ++ if the fact that the regular zero set of \axiom{us} contains that of
         ++ \axiom{ts} can be decided (and in that case \axiom{b} gives this 
         ++ inclusion) otherwise returns \axiom{"failed"}.
     subQuasiComponent?: (TS,TS) -> Boolean
         ++ \axiom{subQuasiComponent?(ts,us)} returns true iff 
         ++ \axiomOpFrom{internalSubQuasiComponent?}{QuasiComponentPackage}
         ++ returs true.
     subQuasiComponent?: (TS,Split) -> Boolean
         ++ \axiom{subQuasiComponent?(ts,lus)} returns true iff
         ++ \axiom{subQuasiComponent?(ts,us)} holds for one \spad{us} in \spad{lus}.
     removeSuperfluousQuasiComponents: Split -> Split
         ++ \axiom{removeSuperfluousQuasiComponents(lts)} removes from \axiom{lts}
         ++ any \spad{ts} such that \axiom{subQuasiComponent?(ts,us)} holds for 
         ++ another \spad{us} in \axiom{lts}.
     subCase?: (LpWT,LpWT) -> Boolean
         ++ \axiom{subCase?(lpwt1,lpwt2)}
         ++ is an internal subroutine, exported only for developement.
     removeSuperfluousCases: List LpWT -> List LpWT
         ++ \axiom{removeSuperfluousCases(llpwt)}
         ++ is an internal subroutine, exported only for developement.
     prepareDecompose: (LP, List(TS),B,B) -> List Branch
         ++ \axiom{prepareDecompose(lp,lts,b1,b2)}
         ++ is an internal subroutine, exported only for developement.
     branchIfCan: (LP,TS,LP,B,B,B,B,B) -> Union(Branch,"failed")
         ++ \axiom{branchIfCan(leq,ts,lineq,b1,b2,b3,b4,b5)}
         ++ is an internal subroutine, exported only for developement.

  Implementation == add

     squareFreeFactors(lp: LP): LP == 
       lsflp: LP := []
       for p in lp repeat 
         lsfp := squareFreeFactors(p)$polsetpack
         lsflp := concat(lsfp,lsflp)
       sort(infRittWu?,removeDuplicates lsflp)

     startTable!(ok: S, ko: S, domainName: S): Void == 
       initTable!()$H
       if (not empty? ok) and (not empty? ko) then printInfo!(ok,ko)$H
       if (not empty? domainName) then startStats!(domainName)$H

     stopTable!(): Void ==   
       if makingStats?()$H then printStats!()$H
       clearTable!()$H

     supDimElseRittWu? (ts:TS,us:TS): Boolean ==
       #ts < #us => true
       #ts > #us => false
       lp1 :LP := members(ts)
       lp2 :LP := members(us)
       while (not empty? lp1) and (not infRittWu?(first(lp2),first(lp1))) repeat
         lp1 := rest lp1
         lp2 := rest lp2
       not empty? lp1

     algebraicSort (lts:Split): Split ==
       lts := removeDuplicates lts
       sort(supDimElseRittWu?,lts)

     moreAlgebraic?(ts:TS,us:TS): Boolean  ==
       empty? ts => empty? us 
       empty? us => true
       #ts < #us => false
       for p in (members us) repeat 
          not algebraic?(mvar(p),ts) => return false
       true

     subTriSet?(ts:TS,us:TS): Boolean  ==
       empty? ts => true
       empty? us => false
       mvar(ts) > mvar(us) => false
       mvar(ts) < mvar(us) => subTriSet?(ts,rest(us)::TS)
       first(ts)::P = first(us)::P => subTriSet?(rest(ts)::TS,rest(us)::TS)
       false

     internalSubPolSet?(lp1: LP, lp2: LP): Boolean  ==
       empty? lp1 => true
       empty? lp2 => false
       associates?(first lp1, first lp2) => 
         internalSubPolSet?(rest lp1, rest lp2)
       infRittWu?(first lp1, first lp2) => false
       internalSubPolSet?(lp1, rest lp2)

     subPolSet?(lp1: LP, lp2: LP): Boolean  ==
       lp1 := sort(infRittWu?, lp1)
       lp2 := sort(infRittWu?, lp2)
       internalSubPolSet?(lp1,lp2)

     infRittWu?(lp1: LP, lp2: LP): Boolean ==
       lp1 := sort(infRittWu?, lp1)
       lp2 := sort(infRittWu?, lp2)
       internalInfRittWu?(lp1,lp2)

     internalInfRittWu?(lp1: LP, lp2: LP): Boolean ==
       empty? lp1 => not empty? lp2
       empty? lp2 => false
       infRittWu?(first lp1, first lp2)$P => true
       infRittWu?(first lp2, first lp1)$P => false
       infRittWu?(rest lp1, rest lp2)$$

     subCase? (lpwt1:LpWT,lpwt2:LpWT): Boolean == 
       -- ASSUME lpwt.{1,2}.val is sorted w.r.t. infRittWu?
       not internalSubPolSet?(lpwt2.val, lpwt1.val) => false
       subQuasiComponent?(lpwt1.tower,lpwt2.tower)

     internalSubQuasiComponent?(ts:TS,us:TS): Union(Boolean,"failed") ==
       -- "failed" is false iff saturate(us) is radical
       subTriSet?(us,ts) => true
       not moreAlgebraic?(ts,us) => false::Union(Boolean,"failed")
       for p in (members us) repeat 
         mdeg(p) < mdeg(select(ts,mvar(p))::P) => 
           return("failed"::Union(Boolean,"failed"))
       for p in (members us) repeat 
         not zero? initiallyReduce(p,ts) =>
           return("failed"::Union(Boolean,"failed"))
       lsfp := squareFreeFactors(initials us)
       for p in lsfp repeat 
         not invertible?(p,ts)@B => 
           return(false::Union(Boolean,"failed"))
       true::Union(Boolean,"failed")

     subQuasiComponent?(ts:TS,us:TS): Boolean ==
       k: Key := [ts, us]
       e := extractIfCan(k)$H
       e case Entry => e::Entry
       ubf: Union(Boolean,"failed") := internalSubQuasiComponent?(ts,us)
       b: Boolean := (ubf case Boolean) and (ubf::Boolean)
       insert!(k,b)$H
       b

     subQuasiComponent?(ts:TS,lus:Split): Boolean ==
       for us in lus repeat
          subQuasiComponent?(ts,us)@B => return true
       false

     removeSuperfluousCases (cases:List LpWT) ==
       #cases < 2 => cases
       toSee := sort(supDimElseRittWu?(#1.tower,#2.tower),cases)
       lpwt1 : LpWT
       toSave,headmaxcases,maxcases,copymaxcases : List LpWT
       while not empty? toSee repeat
         lpwt1 := first toSee
         toSee := rest toSee
         toSave := []
         for lpwt2 in toSee repeat
            if subCase?(lpwt1,lpwt2) 
              then
                lpwt1 := lpwt2
              else
                if not subCase?(lpwt2,lpwt1) 
                  then
                    toSave := cons(lpwt2,toSave)
         if empty? maxcases
           then
             headmaxcases := [lpwt1]
             maxcases := headmaxcases
           else
             copymaxcases := maxcases
             while (not empty? copymaxcases) and _
               (not subCase?(lpwt1,first(copymaxcases))) repeat
                 copymaxcases := rest copymaxcases
             if empty? copymaxcases
               then
                 setrest!(headmaxcases,[lpwt1])
                 headmaxcases := rest headmaxcases
         toSee := reverse toSave
       maxcases

     removeSuperfluousQuasiComponents(lts: Split): Split ==
       lts := removeDuplicates lts
       #lts < 2 => lts
       toSee := algebraicSort lts
       toSave,headmaxlts,maxlts,copymaxlts : Split
       while not empty? toSee repeat
         ts := first toSee
         toSee := rest toSee
         toSave := []
         for us in toSee repeat
            if subQuasiComponent?(ts,us)@B
              then
                ts := us
              else
                if not subQuasiComponent?(us,ts)@B 
                  then
                    toSave := cons(us,toSave)
         if empty? maxlts
           then
             headmaxlts := [ts]
             maxlts := headmaxlts
           else
             copymaxlts := maxlts
             while (not empty? copymaxlts) and _
               (not subQuasiComponent?(ts,first(copymaxlts))@B) repeat
                 copymaxlts := rest copymaxlts
             if empty? copymaxlts
               then
                 setrest!(headmaxlts,[ts])
                 headmaxlts := rest headmaxlts
         toSee := reverse toSave
       algebraicSort maxlts

     removeAssociates (lp:LP):LP ==
       removeDuplicates [primitivePart(p) for p in lp]

     branchIfCan(leq: LP,ts: TS,lineq: LP, b1:B,b2:B,b3:B,b4:B,b5:B):UBF ==
        -- ASSUME pols in leq are squarefree and mainly primitive
        -- if b1 then CLEAN UP leq
        -- if b2 then CLEAN UP lineq
        -- if b3 then SEARCH for ZERO in lineq with leq
        -- if b4 then SEARCH for ZERO in lineq with ts
        -- if b5 then SEARCH for ONE in leq with lineq
        if b1 
          then 
            leq := removeAssociates(leq)
            leq := remove(zero?,leq)
            any?(ground?,leq) => 
              return("failed"::Union(Branch,"failed"))
        if b2
          then
            any?(zero?,lineq) =>
              return("failed"::Union(Branch,"failed"))
            lineq := removeRedundantFactors(lineq)$polsetpack
        if b3
          then
            ps: PS := construct(leq)$PS
            for q in lineq repeat
              zero? remainder(q,ps).polnum =>
                return("failed"::Union(Branch,"failed"))
        (empty? leq) or (empty? lineq) => ([leq, ts, lineq]$Branch)::UBF
        if b4
          then
            for q: local in lineq repeat
              zero? initiallyReduce(q,ts) => 
                return("failed"::Union(Branch,"failed"))
        if b5
          then
            newleq: LP := []
            for p in leq repeat
              for q: local in lineq repeat
                if mvar(p) = mvar(q)
                  then
                    g := gcd(p,q)
                    newp := (p exquo g)::P
                    ground? newp => 
                      return("failed"::Union(Branch,"failed"))
                    newleq := cons(newp,newleq)
                  else
                    newleq := cons(p,newleq)
            leq := newleq
        leq := sort(infRittWu?, removeDuplicates leq)
        ([leq, ts, lineq]$Branch)::UBF

     prepareDecompose(lp: LP, lts: List(TS), b1: B, b2: B): List Branch ==
       -- if b1 then REMOVE REDUNDANT COMPONENTS in lts
       -- if b2 then SPLIT the input system with squareFree
       lp := sort(infRittWu?, remove(zero?,removeAssociates(lp)))
       any?(ground?,lp) => []
       empty? lts => []
       if b1 then lts := removeSuperfluousQuasiComponents lts
       not b2 =>
         [[lp,ts,squareFreeFactors(initials ts)]$Branch for ts in lts]
       toSee: List Branch 
       lq: LP := []         
       toSee := [[lq,ts,squareFreeFactors(initials ts)]$Branch for ts in lts]
       empty? lp => toSee
       for p in lp repeat
         lsfp := squareFreeFactors(p)$polsetpack
         branches: List Branch := []
         lq := []
         for f in lsfp repeat
           for branch in toSee repeat
             leq : LP := branch.eq
             ts := branch.tower
             lineq : LP := branch.ineq
             ubf1: UBF := branchIfCan(leq,ts,lq,false,false,true,true,true)@UBF
             ubf1 case "failed" => "leave"
             ubf2: UBF := branchIfCan([f],ts,lineq,false,false,true,true,true)@UBF
             ubf2 case "failed" => "leave"
             leq := sort(infRittWu?,removeDuplicates concat(ubf1.eq,ubf2.eq))
             lineq := sort(infRittWu?,removeDuplicates concat(ubf1.ineq,ubf2.ineq))
             newBranch := branchIfCan(leq,ts,lineq,false,false,false,false,false)
             branches:= cons(newBranch::Branch,branches)
           lq := cons(f,lq)
         toSee := branches
       sort(supDimElseRittWu?(#1.tower,#2.tower),toSee)

@
\section{package RSETGCD RegularTriangularSetGcdPackage}
<<package RSETGCD RegularTriangularSetGcdPackage>>=
)abbrev package RSETGCD RegularTriangularSetGcdPackage
++ Author: Marc Moreno Maza (marc@nag.co.uk)
++ Date Created: 08/30/1998
++ Date Last Updated: 12/15/1998
++ Basic Functions:
++ Related Constructors:
++ Also See: 
++ AMS Classifications:
++ Keywords:
++ Description: 
++ An internal package for computing gcds and resultants of univariate
++ polynomials with coefficients in a tower of simple extensions of a field.\newline
++ References :
++  [1] M. MORENO MAZA and R. RIOBOO "Computations of gcd over
++      algebraic towers of simple extensions" In proceedings of AAECC11
++      Paris, 1995.
++  [2] M. MORENO MAZA "Calculs de pgcd au-dessus des tours
++      d'extensions simples et resolution des systemes d'equations
++      algebriques" These, Universite P.etM. Curie, Paris, 1997.
++  [3] M. MORENO MAZA "A new algorithm for computing triangular
++      decomposition of algebraic varieties" NAG Tech. Rep. 4/98.
++ Version: 4.

RegularTriangularSetGcdPackage(R,E,V,P,TS): Exports == Implementation where

  R : GcdDomain
  E : OrderedAbelianMonoidSup
  V : OrderedSet
  P : RecursivePolynomialCategory(R,E,V)
  TS : RegularTriangularSetCategory(R,E,V,P)
  N ==> NonNegativeInteger
  Z ==> Integer
  B ==> Boolean
  S ==> String
  LP ==> List P
  PtoP ==> P -> P
  PS ==> GeneralPolynomialSet(R,E,V,P)
  PWT ==> Record(val : P, tower : TS)
  BWT ==> Record(val : Boolean, tower : TS)
  LpWT ==> Record(val : (List P), tower : TS)
  Branch ==> Record(eq: List P, tower: TS, ineq: List P)
  UBF ==> Union(Branch,"failed")
  Split ==> List TS
  KeyGcd ==> Record(arg1: P, arg2: P, arg3: TS, arg4: B)
  EntryGcd ==> List PWT
  HGcd ==> TabulatedComputationPackage(KeyGcd, EntryGcd)
  KeyInvSet ==> Record(arg1: P, arg3: TS)
  EntryInvSet ==> List TS
  HInvSet ==> TabulatedComputationPackage(KeyInvSet, EntryInvSet)
  polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,P)
  quasicomppack ==> QuasiComponentPackage(R,E,V,P,TS)

  Exports ==  with
     startTableGcd!: (S,S,S) -> Void
         ++ \axiom{startTableGcd!(s1,s2,s3)} 
         ++ is an internal subroutine, exported only for developement.
     stopTableGcd!: () -> Void
         ++ \axiom{stopTableGcd!()} 
         ++ is an internal subroutine, exported only for developement.
     startTableInvSet!: (S,S,S) -> Void
         ++ \axiom{startTableInvSet!(s1,s2,s3)} 
         ++ is an internal subroutine, exported only for developement.
     stopTableInvSet!: () -> Void
         ++ \axiom{stopTableInvSet!()} is an internal subroutine,
         ++ exported only for developement.
     prepareSubResAlgo: (P,P,TS) -> List LpWT 
         ++ \axiom{prepareSubResAlgo(p1,p2,ts)} 
         ++ is an internal subroutine, exported only for developement.
     internalLastSubResultant: (P,P,TS,B,B) -> List PWT 
         ++ \axiom{internalLastSubResultant(p1,p2,ts,inv?,break?)}
         ++ is an internal subroutine, exported only for developement.
     internalLastSubResultant: (List LpWT,V,B) -> List PWT 
         ++ \axiom{internalLastSubResultant(lpwt,v,flag)} is an internal
         ++ subroutine, exported only for developement.
     integralLastSubResultant: (P,P,TS) -> List PWT 
         ++ \axiom{integralLastSubResultant(p1,p2,ts)} 
         ++ is an internal subroutine, exported only for developement.
     toseLastSubResultant: (P,P,TS) -> List PWT 
         ++ \axiom{toseLastSubResultant(p1,p2,ts)} has the same specifications as
         ++ \axiomOpFrom{lastSubResultant}{RegularTriangularSetCategory}.
     toseInvertible?: (P,TS) -> B
         ++ \axiom{toseInvertible?(p1,p2,ts)} has the same specifications as
         ++ \axiomOpFrom{invertible?}{RegularTriangularSetCategory}.
     toseInvertible?: (P,TS) -> List BWT
         ++ \axiom{toseInvertible?(p1,p2,ts)} has the same specifications as
         ++ \axiomOpFrom{invertible?}{RegularTriangularSetCategory}.
     toseInvertibleSet: (P,TS) -> Split
         ++ \axiom{toseInvertibleSet(p1,p2,ts)} has the same specifications as
         ++ \axiomOpFrom{invertibleSet}{RegularTriangularSetCategory}.
     toseSquareFreePart: (P,TS) -> List PWT 
         ++ \axiom{toseSquareFreePart(p,ts)} has the same specifications as
         ++ \axiomOpFrom{squareFreePart}{RegularTriangularSetCategory}.

  Implementation == add

     startTableGcd!(ok: S, ko: S, domainName: S): Void == 
       initTable!()$HGcd
       printInfo!(ok,ko)$HGcd
       startStats!(domainName)$HGcd

     stopTableGcd!(): Void ==   
       if makingStats?()$HGcd then printStats!()$HGcd
       clearTable!()$HGcd

     startTableInvSet!(ok: S, ko: S, domainName: S): Void == 
       initTable!()$HInvSet
       printInfo!(ok,ko)$HInvSet
       startStats!(domainName)$HInvSet

     stopTableInvSet!(): Void ==   
       if makingStats?()$HInvSet then printStats!()$HInvSet
       clearTable!()$HInvSet

     toseInvertible?(p:P,ts:TS): Boolean == 
       q := primitivePart initiallyReduce(p,ts)
       zero? q => false
       normalized?(q,ts) => true
       v := mvar(q)
       not algebraic?(v,ts) => 
         toCheck: List BWT := toseInvertible?(p,ts)@(List BWT)
         for bwt in toCheck repeat
           bwt.val = false => return false
         return true
       ts_v := select(ts,v)::P
       ts_v_- := collectUnder(ts,v)
       lgwt := internalLastSubResultant(ts_v,q,ts_v_-,false,true)
       for gwt in lgwt repeat
         g := gwt.val; 
         (not ground? g) and (mvar(g) = v) => 
           return false
       true
       
     toseInvertible?(p:P,ts:TS): List BWT ==
       q := primitivePart initiallyReduce(p,ts)
       zero? q => [[false,ts]$BWT]
       normalized?(q,ts) => [[true,ts]$BWT]
       v := mvar(q)
       not algebraic?(v,ts) => 
         lbwt: List BWT := []
         toCheck: List BWT := toseInvertible?(init(q),ts)@(List BWT)
         for bwt in toCheck repeat
           bwt.val => lbwt := cons(bwt,lbwt)
           newq := removeZero(q,bwt.tower)
           zero? newq => lbwt := cons(bwt,lbwt)
           lbwt := concat(toseInvertible?(newq,bwt.tower)@(List BWT), lbwt)
         return lbwt
       ts_v := select(ts,v)::P
       ts_v_- := collectUnder(ts,v)
       ts_v_+ := collectUpper(ts,v)
       lgwt := internalLastSubResultant(ts_v,q,ts_v_-,false,false)
       lbwt: List BWT := []
       for gwt in lgwt repeat
         g := gwt.val; ts := gwt.tower
         (ground? g) or (mvar(g) < v) => 
           ts := internalAugment(ts_v,ts)
           ts := internalAugment(members(ts_v_+),ts)
           lbwt := cons([true, ts]$BWT,lbwt)
         g := mainPrimitivePart g
         ts_g := internalAugment(g,ts)
         ts_g := internalAugment(members(ts_v_+),ts_g)
         -- USE internalAugment with parameters ??
         lbwt := cons([false, ts_g]$BWT,lbwt)
         h := lazyPquo(ts_v,g)
         (ground? h) or (mvar(h) < v) => "leave"
         h := mainPrimitivePart h
         ts_h := internalAugment(h,ts)
         ts_h := internalAugment(members(ts_v_+),ts_h)
         -- USE internalAugment with parameters ??
         -- CAN BE OPTIMIZED if the input tower is separable
         inv := toseInvertible?(q,ts_h)@(List BWT)
         lbwt := concat([bwt for bwt in inv | bwt.val],lbwt)
       sort(#1.val < #2.val,lbwt)

     toseInvertibleSet(p:P,ts:TS): Split ==
       k: KeyInvSet := [p,ts]
       e := extractIfCan(k)$HInvSet
       e case EntryInvSet => e::EntryInvSet
       q := primitivePart initiallyReduce(p,ts)
       zero? q => []
       normalized?(q,ts) => [ts]
       v := mvar(q)
       toSave: Split := []
       not algebraic?(v,ts) => 
         toCheck: List BWT := toseInvertible?(init(q),ts)@(List BWT)
         for bwt in toCheck repeat
           bwt.val => toSave := cons(bwt.tower,toSave)
           newq := removeZero(q,bwt.tower)
           zero? newq => "leave"
           toSave := concat(toseInvertibleSet(newq,bwt.tower), toSave)
         toSave := removeDuplicates toSave
         return algebraicSort(toSave)$quasicomppack
       ts_v := select(ts,v)::P
       ts_v_- := collectUnder(ts,v)
       ts_v_+ := collectUpper(ts,v)
       lgwt := internalLastSubResultant(ts_v,q,ts_v_-,false,false)
       for gwt in lgwt repeat
         g := gwt.val; ts := gwt.tower
         (ground? g) or (mvar(g) < v) => 
           ts := internalAugment(ts_v,ts)
           ts := internalAugment(members(ts_v_+),ts)
           toSave := cons(ts,toSave)
         g := mainPrimitivePart g
         h := lazyPquo(ts_v,g)
         h := mainPrimitivePart h
         (ground? h) or (mvar(h) < v) => "leave"
         ts_h := internalAugment(h,ts)
         ts_h := internalAugment(members(ts_v_+),ts_h)
         inv := toseInvertibleSet(q,ts_h)
         toSave := removeDuplicates concat(inv,toSave)
       toSave := algebraicSort(toSave)$quasicomppack
       insert!(k,toSave)$HInvSet
       toSave

     toseSquareFreePart_wip(p:P, ts: TS): List PWT ==
     -- ASSUME p is not constant and mvar(p) > mvar(ts)
     -- ASSUME init(p) is invertible w.r.t. ts
     -- ASSUME p is mainly primitive
       one? mdeg(p) => [[p,ts]$PWT]
       v := mvar(p)$P
       q: P := mainPrimitivePart D(p,v)
       lgwt: List PWT := internalLastSubResultant(p,q,ts,true,false)
       lpwt : List PWT := []
       sfp : P
       for gwt in lgwt repeat
         g := gwt.val; us := gwt.tower
         (ground? g) or (mvar(g) < v) =>
           lpwt := cons([p,us],lpwt)
         g := mainPrimitivePart g
         sfp := lazyPquo(p,g)
         sfp := mainPrimitivePart stronglyReduce(sfp,us)
         lpwt := cons([sfp,us],lpwt)
       lpwt

     toseSquareFreePart_base(p:P, ts: TS): List PWT == [[p,ts]$PWT]

     toseSquareFreePart(p:P, ts: TS): List PWT == toseSquareFreePart_wip(p,ts)

     prepareSubResAlgo(p1:P,p2:P,ts:TS): List LpWT ==
       -- ASSUME mvar(p1) = mvar(p2) > mvar(ts) and mdeg(p1) >= mdeg(p2)
       -- ASSUME init(p1) invertible modulo ts !!!
       toSee: List LpWT := [[[p1,p2],ts]$LpWT]
       toSave: List LpWT := []
       v := mvar(p1)
       while (not empty? toSee) repeat
         lpwt := first toSee; toSee := rest toSee
         p1 := lpwt.val.1; p2 := lpwt.val.2
         ts := lpwt.tower
         lbwt := toseInvertible?(leadingCoefficient(p2,v),ts)@(List BWT)
         for bwt in lbwt repeat
           (bwt.val = true) and positive? degree(p2,v) =>
             p3 := prem(p1, -p2)
             s: P := init(p2)**(mdeg(p1) - mdeg(p2))::N
             toSave := cons([[p2,p3,s],bwt.tower]$LpWT,toSave)
           -- p2 := initiallyReduce(p2,bwt.tower)
           newp2 := primitivePart initiallyReduce(p2,bwt.tower)
           (bwt.val = true) =>
             -- toSave := cons([[p2,0,1],bwt.tower]$LpWT,toSave)
             toSave := cons([[p2,0,1],bwt.tower]$LpWT,toSave)
           -- zero? p2 => 
           zero? newp2 => 
             toSave := cons([[p1,0,1],bwt.tower]$LpWT,toSave)
           -- toSee := cons([[p1,p2],ts]$LpWT,toSee)
           toSee := cons([[p1,newp2],bwt.tower]$LpWT,toSee)
       toSave

     integralLastSubResultant(p1:P,p2:P,ts:TS): List PWT ==
       -- ASSUME mvar(p1) = mvar(p2) > mvar(ts) and mdeg(p1) >= mdeg(p2)
       -- ASSUME p1 and p2 have no algebraic coefficients
       lsr := lastSubResultant(p1, p2)
       ground?(lsr) => [[lsr,ts]$PWT]
       mvar(lsr) < mvar(p1) => [[lsr,ts]$PWT]
       gi1i2 := gcd(init(p1),init(p2))
       ex: Union(P,"failed") := (gi1i2 * lsr) exquo$P init(lsr)
       ex case "failed" => [[lsr,ts]$PWT]
       [[ex::P,ts]$PWT]
            
     internalLastSubResultant(p1:P,p2:P,ts:TS,b1:B,b2:B): List PWT ==
       -- ASSUME mvar(p1) = mvar(p2) > mvar(ts) and mdeg(p1) >= mdeg(p2)
       -- if b1 ASSUME init(p2) invertible w.r.t. ts
       -- if b2 BREAK with the first non-trivial gcd 
       k: KeyGcd := [p1,p2,ts,b2]
       e := extractIfCan(k)$HGcd
       e case EntryGcd => e::EntryGcd
       toSave: List PWT 
       empty? ts => 
         toSave := integralLastSubResultant(p1,p2,ts)
         insert!(k,toSave)$HGcd
         return toSave
       toSee: List LpWT 
       if b1
         then
           p3 := prem(p1, -p2)
           s: P := init(p2)**(mdeg(p1) - mdeg(p2))::N
           toSee := [[[p2,p3,s],ts]$LpWT]
         else
           toSee := prepareSubResAlgo(p1,p2,ts)
       toSave := internalLastSubResultant(toSee,mvar(p1),b2)
       insert!(k,toSave)$HGcd
       toSave

     internalLastSubResultant(llpwt: List LpWT,v:V,b2:B): List PWT ==
       toReturn: List PWT := []; toSee: List LpWT; 
       while (not empty? llpwt) repeat
         toSee := llpwt; llpwt := []
         -- CONSIDER FIRST the vanishing current last subresultant
         for lpwt in toSee repeat 
           p1 := lpwt.val.1; p2 := lpwt.val.2; s := lpwt.val.3; ts := lpwt.tower
           lbwt := toseInvertible?(leadingCoefficient(p2,v),ts)@(List BWT)
           for bwt in lbwt repeat
             bwt.val = false => 
               toReturn := cons([p1,bwt.tower]$PWT, toReturn)
               b2 and positive?(degree(p1,v)) => return toReturn
             llpwt := cons([[p1,p2,s],bwt.tower]$LpWT, llpwt)
         empty? llpwt => "leave"
         -- CONSIDER NOW the branches where the computations continue
         toSee := llpwt; llpwt := []
         lpwt := first toSee; toSee := rest toSee
         p1 := lpwt.val.1; p2 := lpwt.val.2; s := lpwt.val.3
         delta: N := (mdeg(p1) - degree(p2,v))::N
         p3: P := LazardQuotient2(p2, leadingCoefficient(p2,v), s, delta)
         zero?(degree(p3,v)) =>
           toReturn := cons([p3,lpwt.tower]$PWT, toReturn)
           for lpwt: free in toSee repeat 
             toReturn := cons([p3,lpwt.tower]$PWT, toReturn)
         (p1, p2) := (p3, next_subResultant2(p1, p2, p3, s))
         s := leadingCoefficient(p1,v)
         llpwt := cons([[p1,p2,s],lpwt.tower]$LpWT, llpwt)
         for lpwt:local in toSee repeat 
           llpwt := cons([[p1,p2,s],lpwt.tower]$LpWT, llpwt)
       toReturn

     toseLastSubResultant(p1:P,p2:P,ts:TS): List PWT ==
       ground? p1 => 
         error"in toseLastSubResultantElseSplit$TOSEGCD  : bad #1"
       ground? p2 => 
         error"in toseLastSubResultantElseSplit$TOSEGCD : bad #2"
       not (mvar(p2) = mvar(p1)) => 
         error"in toseLastSubResultantElseSplit$TOSEGCD : bad #2"
       algebraic?(mvar(p1),ts) =>
         error"in toseLastSubResultantElseSplit$TOSEGCD : bad #1"
       not initiallyReduced?(p1,ts) => 
         error"in toseLastSubResultantElseSplit$TOSEGCD : bad #1"
       not initiallyReduced?(p2,ts) => 
         error"in toseLastSubResultantElseSplit$TOSEGCD : bad #2"
       purelyTranscendental?(p1,ts) and purelyTranscendental?(p2,ts) =>
         integralLastSubResultant(p1,p2,ts)
       if mdeg(p1) < mdeg(p2) then 
          (p1, p2) := (p2, p1)
          if odd?(mdeg(p1)) and odd?(mdeg(p2)) then p2 := - p2
       internalLastSubResultant(p1,p2,ts,false,false)

@
\section{package RSDCMPK RegularSetDecompositionPackage}
<<package RSDCMPK RegularSetDecompositionPackage>>=
)abbrev package RSDCMPK RegularSetDecompositionPackage
++ Author: Marc Moreno Maza
++ Date Created: 09/16/1998
++ Date Last Updated: 12/16/1998
++ Basic Functions:
++ Related Constructors:
++ Also See: 
++ AMS Classifications:
++ Keywords:
++ Description: 
++ A package providing a new algorithm for solving polynomial systems
++ by means of regular chains. Two ways of solving are proposed:
++ in the sense of Zariski closure (like in Kalkbrener's algorithm)
++ or in the sense of the regular zeros (like in Wu, Wang or Lazard
++ methods). This algorithm is valid for nay type
++ of regular set. It does not care about the way a polynomial is
++ added in an regular set, or how two quasi-components are compared
++ (by an inclusion-test), or how the invertibility test is made in
++ the tower of simple extensions associated with a regular set.
++ These operations are realized respectively by the domain \spad{TS}
++ and the packages \axiomType{QCMPACK}(R,E,V,P,TS) and \axiomType{RSETGCD}(R,E,V,P,TS).
++ The same way it does not care about the way univariate polynomial
++ gcd (with coefficients in the tower of simple extensions associated 
++ with a regular set) are computed. The only requirement is that these
++ gcd need to have invertible initials (normalized or not).
++ WARNING. There is no need for a user to call diectly any operation
++ of this package since they can be accessed by the domain \axiom{TS}.
++ Thus, the operations of this package are not documented.\newline
++ References :
++  [1] M. MORENO MAZA "A new algorithm for computing triangular
++      decomposition of algebraic varieties" NAG Tech. Rep. 4/98.
++ Version: 5. Same as 4 but Does NOT use any unproved criteria.

RegularSetDecompositionPackage(R,E,V,P,TS): Exports == Implementation where

  R : GcdDomain
  E : OrderedAbelianMonoidSup
  V : OrderedSet
  P : RecursivePolynomialCategory(R,E,V)
  TS : RegularTriangularSetCategory(R,E,V,P)
  N ==> NonNegativeInteger
  Z ==> Integer
  B ==> Boolean
  LP ==> List P
  PS ==> GeneralPolynomialSet(R,E,V,P)
  PWT ==> Record(val : P, tower : TS)
  BWT ==> Record(val : Boolean, tower : TS)
  LpWT ==> Record(val : (List P), tower : TS)
  Wip ==> Record(done: Split, todo: List LpWT)
  Branch ==> Record(eq: List P, tower: TS, ineq: List P)
  UBF ==> Union(Branch,"failed")
  Split ==> List TS
  iprintpack ==> InternalPrintPackage()
  polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,P)
  quasicomppack ==> QuasiComponentPackage(R,E,V,P,TS)
  regsetgcdpack ==> RegularTriangularSetGcdPackage(R,E,V,P,TS)

  Exports ==  with

     KrullNumber: (LP, Split) -> N
     numberOfVariables: (LP, Split) -> N
     algebraicDecompose: (P,TS,B) -> Record(done: Split, todo: List LpWT) 
     transcendentalDecompose: (P,TS,N) -> Record(done: Split, todo: List LpWT) 
     transcendentalDecompose: (P,TS) -> Record(done: Split, todo: List LpWT) 
     internalDecompose: (P,TS,N,B) -> Record(done: Split, todo: List LpWT)
     internalDecompose: (P,TS,N) -> Record(done: Split, todo: List LpWT)
     internalDecompose: (P,TS) -> Record(done: Split, todo: List LpWT)
     decompose: (LP, Split, B, B) -> Split
     decompose: (LP, Split, B, B, B, B, B) -> Split
     upDateBranches: (LP,Split,List LpWT,Wip,N) -> List LpWT
     convert: Record(val: List P,tower: TS) -> String
     printInfo: (List Record(val: List P,tower: TS), N) -> Void

  Implementation == add

     KrullNumber(lp: LP, lts: Split): N ==
       ln: List N := [#(ts) for ts in lts]
       n := #lp + reduce(max,ln)

     numberOfVariables(lp: LP, lts: Split): N ==
       lv: List V := variables([lp]$PS)
       for ts in lts repeat lv := concat(variables(ts), lv)
       # removeDuplicates(lv)

     algebraicDecompose(p: P, ts: TS, clos?: B): Record(done: Split, todo: List LpWT) ==
       ground? p =>
         error " in algebraicDecompose$REGSET: should never happen !"
       v := mvar(p); n := #ts
       ts_v_- := collectUnder(ts,v)
       ts_v_+ := collectUpper(ts,v)
       ts_v := select(ts,v)::P
       if mdeg(p) < mdeg(ts_v)
         then 
           lgwt := internalLastSubResultant(ts_v,p,ts_v_-,true,false)$regsetgcdpack
         else
           lgwt := internalLastSubResultant(p,ts_v,ts_v_-,true,false)$regsetgcdpack
       lts: Split := []
       llpwt: List LpWT := []
       for gwt in lgwt repeat
         g := gwt.val; us := gwt.tower
         zero? g => 
           error " in algebraicDecompose$REGSET: should never happen !!"
         ground? g => "leave"
         if mvar(g) = v then lts := concat(augment(members(ts_v_+),augment(g,us)),lts) 
         h := leadingCoefficient(g,v)
         b: Boolean := purelyAlgebraic?(us)
         lsfp := squareFreeFactors(h)$polsetpack
         lus := augment(members(ts_v_+),augment(ts_v,us)@Split)
         for f in lsfp repeat
           ground? f => "leave"
           b and purelyAlgebraic?(f,us) => "leave"
           for vs in lus repeat
              llpwt := cons([[f,p],vs]$LpWT, llpwt)
       [lts,llpwt]

     transcendentalDecompose(p: P, ts: TS,bound: N): Record(done: Split, todo: List LpWT) ==
       lts: Split
       if #ts < bound 
         then
           lts := augment(p,ts)
         else
           lts := []
       llpwt: List LpWT := []
       [lts,llpwt]

     transcendentalDecompose(p: P, ts: TS): Record(done: Split, todo: List LpWT) ==
       lts: Split:= augment(p,ts)
       llpwt: List LpWT := []
       [lts,llpwt]

     internalDecompose(p: P, ts: TS,bound: N,clos?:B): Record(done: Split, todo: List LpWT) ==
       clos? => internalDecompose(p,ts,bound)
       internalDecompose(p,ts)

     internalDecompose(p: P, ts: TS,bound: N): Record(done: Split, todo: List LpWT) ==
       -- ASSUME p not constant
       llpwt: List LpWT := []
       lts: Split := []
       -- EITHER mvar(p) is null
       if (not zero? tail(p)) and (not ground? (lmp := leastMonomial(p)))
         then
           llpwt := cons([[mvar(p)::P],ts]$LpWT,llpwt)
           p := (p exquo lmp)::P
       ip := squareFreePart init(p); tp := tail p
       p := mainPrimitivePart p
       -- OR init(p) is null or not
       lbwt := invertible?(ip,ts)@(List BWT)
       for bwt in lbwt repeat
         bwt.val =>
           if algebraic?(mvar(p),bwt.tower) 
             then 
               rsl := algebraicDecompose(p,bwt.tower,true)
             else
               rsl := transcendentalDecompose(p,bwt.tower,bound)
           lts := concat(rsl.done,lts)
           llpwt := concat(rsl.todo,llpwt)
           -- purelyAlgebraicLeadingMonomial?(ip,bwt.tower) => "leave"  -- UNPROVED CRITERIA
           purelyAlgebraic?(ip,bwt.tower) and purelyAlgebraic?(bwt.tower) => "leave" -- SAFE
           (not ground? ip) => 
             zero? tp => llpwt := cons([[ip],bwt.tower]$LpWT, llpwt)
             (not ground? tp) => llpwt := cons([[ip,tp],bwt.tower]$LpWT, llpwt)
         riv := removeZero(ip,bwt.tower)
         (zero? riv) =>
           zero? tp => lts := cons(bwt.tower,lts)
           (not ground? tp) => llpwt := cons([[tp],bwt.tower]$LpWT, llpwt)
         llpwt := cons([[riv * mainMonomial(p) + tp],bwt.tower]$LpWT, llpwt)
       [lts,llpwt]

     internalDecompose(p: P, ts: TS): Record(done: Split, todo: List LpWT) ==
       -- ASSUME p not constant
       llpwt: List LpWT := []
       lts: Split := []
       -- EITHER mvar(p) is null
       if (not zero? tail(p)) and (not ground? (lmp := leastMonomial(p)))
         then
           llpwt := cons([[mvar(p)::P],ts]$LpWT,llpwt)
           p := (p exquo lmp)::P
       ip := squareFreePart init(p); tp := tail p
       p := mainPrimitivePart p
       -- OR init(p) is null or not
       lbwt := invertible?(ip,ts)@(List BWT)
       for bwt in lbwt repeat
         bwt.val =>
           if algebraic?(mvar(p),bwt.tower) 
             then 
               rsl := algebraicDecompose(p,bwt.tower,false)
             else
               rsl := transcendentalDecompose(p,bwt.tower)
           lts := concat(rsl.done,lts)
           llpwt :=  concat(rsl.todo,llpwt)
           purelyAlgebraic?(ip,bwt.tower) and purelyAlgebraic?(bwt.tower) => "leave"
           (not ground? ip) => 
             zero? tp => llpwt := cons([[ip],bwt.tower]$LpWT, llpwt)
             (not ground? tp) => llpwt := cons([[ip,tp],bwt.tower]$LpWT, llpwt)
         riv := removeZero(ip,bwt.tower)
         (zero? riv) =>
           zero? tp => lts := cons(bwt.tower,lts)
           (not ground? tp) => llpwt := cons([[tp],bwt.tower]$LpWT, llpwt)
         llpwt := cons([[riv * mainMonomial(p) + tp],bwt.tower]$LpWT, llpwt)
       [lts,llpwt]

     decompose(lp: LP, lts: Split, clos?: B, info?: B): Split ==
       decompose(lp,lts,false,false,clos?,true,info?)

     convert(lpwt: LpWT): String ==
       ls: List String := ["<", string((#(lpwt.val))::Z), ",", string((#(lpwt.tower))::Z), ">" ]
       concat ls

     printInfo(toSee: List LpWT, n: N): Void ==
       lpwt := first toSee
       s: String := concat ["[", string((#toSee)::Z), " ", convert(lpwt)@String]
       m: N := #(lpwt.val)
       toSee := rest toSee
       for lpwt: local in toSee repeat
         m := m + #(lpwt.val)
         s := concat [s, ",", convert(lpwt)@String]
       s := concat [s, " -> |", string(m::Z), "|; {", string(n::Z),"}]"]
       iprint(s)$iprintpack

     decompose(lp: LP, lts: Split, cleanW?: B, sqfr?: B, clos?: B, rem?: B, info?: B): Split ==
       -- if cleanW? then REMOVE REDUNDANT COMPONENTS in lts
       -- if sqfr? then SPLIT the system with SQUARE-FREE FACTORIZATION
       -- if clos? then SOLVE in the closure sense 
       -- if rem? then REDUCE the current p by using remainder
       -- if info? then PRINT info
       empty? lp => lts
       branches: List Branch := prepareDecompose(lp,lts,cleanW?,sqfr?)$quasicomppack
       empty? branches => []
       toSee: List LpWT := [[br.eq,br.tower]$LpWT for br in branches]
       toSave: Split := []
       if clos? then bound := KrullNumber(lp,lts) else bound := numberOfVariables(lp,lts)
       while (not empty? toSee) repeat
         if info? then printInfo(toSee,#toSave)
         lpwt := first toSee; toSee := rest toSee
         lp := lpwt.val; ts := lpwt.tower
         empty? lp => 
           toSave := cons(ts, toSave)
         p := first lp;  lp := rest lp
         if rem? and (not ground? p) and (not empty? ts)  
            then 
              p := remainder(p,ts).polnum
         p := removeZero(p,ts)
         zero? p => toSee := cons([lp,ts]$LpWT, toSee)
         ground? p => "leave"
         rsl := internalDecompose(p,ts,bound,clos?)
         toSee := upDateBranches(lp,toSave,toSee,rsl,bound)
       removeSuperfluousQuasiComponents(toSave)$quasicomppack

     upDateBranches(leq:LP,lts:Split,current:List LpWT,wip: Wip,n:N): List LpWT ==
       newBranches: List LpWT := wip.todo
       newComponents: Split := wip.done
       branches1, branches2:  List LpWT 
       branches1 := []; branches2  := []
       for branch in newBranches repeat
         us := branch.tower
         #us > n => "leave"
         newleq := sort(infRittWu?,concat(leq,branch.val))
         --foo := rewriteSetWithReduction(newleq,us,initiallyReduce,initiallyReduced?)
         --any?(ground?,foo)  => "leave"
         branches1 := cons([newleq,us]$LpWT, branches1)
       for us in newComponents repeat
         #us > n => "leave"
         subQuasiComponent?(us,lts)$quasicomppack => "leave"
         --newleq := leq
         --foo := rewriteSetWithReduction(newleq,us,initiallyReduce,initiallyReduced?)
         --any?(ground?,foo)  => "leave"
         branches2 := cons([leq,us]$LpWT, branches2)
       empty? branches1 => 
         empty? branches2 => current
         concat(branches2, current)
       branches := concat [branches2, branches1, current]
       -- branches := concat(branches,current)
       removeSuperfluousCases(branches)$quasicomppack

@
\section{domain REGSET RegularTriangularSet}
Several domain constructors implement regular triangular sets (or regular
chains). Among them {\bf RegularTriangularSet} and 
{\bf SquareFreeRegularTriangularSet}. They also implement an algorithm
by Marc Moreno Maza for computing triangular decompositions of polynomial
systems. This method is refined in the package {\bf LazardSetSolvingPackage}
in order to produce decompositions by means of Lazard triangular sets.
<<domain REGSET RegularTriangularSet>>=
)abbrev domain REGSET RegularTriangularSet
++ Author: Marc Moreno Maza
++ Date Created: 08/25/1998
++ Date Last Updated: 16/12/1998
++ Basic Functions:
++ Related Constructors:
++ Also See: 
++ AMS Classifications:
++ Keywords:
++ Description: 
++ This domain provides an implementation of regular chains.
++ Moreover, the operation \axiomOpFrom{zeroSetSplit}{RegularTriangularSetCategory}
++ is an implementation of a new algorithm for solving polynomial systems by
++ means of regular chains.\newline
++ References :
++  [1] M. MORENO MAZA "A new algorithm for computing triangular
++      decomposition of algebraic varieties" NAG Tech. Rep. 4/98.
++ Version: Version 11. 

RegularTriangularSet(R,E,V,P) : Exports == Implementation where

  R : GcdDomain
  E : OrderedAbelianMonoidSup
  V : OrderedSet
  P : RecursivePolynomialCategory(R,E,V)
  N ==> NonNegativeInteger
  Z ==> Integer
  B ==> Boolean
  LP ==> List P
  PtoP ==> P -> P
  PS ==> GeneralPolynomialSet(R,E,V,P)
  PWT ==> Record(val : P, tower : $)
  BWT ==> Record(val : Boolean, tower : $)
  LpWT ==> Record(val : (List P), tower : $)
  Split ==> List $
  iprintpack ==> InternalPrintPackage()
  polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,P)
  quasicomppack ==> QuasiComponentPackage(R,E,V,P,$)
  regsetgcdpack ==> RegularTriangularSetGcdPackage(R,E,V,P,$)
  regsetdecomppack ==> RegularSetDecompositionPackage(R,E,V,P,$)

  Exports ==  RegularTriangularSetCategory(R,E,V,P) with

     internalAugment: (P,$,B,B,B,B,B) -> List $
       ++ \axiom{internalAugment(p,ts,b1,b2,b3,b4,b5)}
       ++ is an internal subroutine, exported only for developement.
     zeroSetSplit: (LP, B, B) -> Split
       ++ \axiom{zeroSetSplit(lp,clos?,info?)} has the same specifications as
       ++ \axiomOpFrom{zeroSetSplit}{RegularTriangularSetCategory}.
       ++ Moreover, if \axiom{clos?} then solves in the sense of the Zariski closure
       ++ else solves in the sense of the regular zeros. If \axiom{info?} then
       ++ do print messages during the computations.
     zeroSetSplit: (LP, B, B, B, B) -> Split
       ++ \axiom{zeroSetSplit(lp,b1,b2.b3,b4)} 
       ++ is an internal subroutine, exported only for developement.
     internalZeroSetSplit: (LP, B, B, B) -> Split
       ++ \axiom{internalZeroSetSplit(lp,b1,b2,b3)}
       ++ is an internal subroutine, exported only for developement.
     pre_process: (LP, B, B) -> Record(val: LP, towers: Split)
       ++ \axiom{pre_process(lp,b1,b2)} 
       ++ is an internal subroutine, exported only for developement.

  Implementation == add

     Rep == LP

     copy ts ==
       per(copy(rep(ts))$LP)
     empty() ==
       per([])
     empty?(ts:$) ==
       empty?(rep(ts))
     members ts ==
       rep(ts)
     map (f : PtoP, ts : $) : $ ==
       construct(map(f,rep(ts))$LP)$$
     map! (f : PtoP, ts : $) : $  ==
       construct(map!(f,rep(ts))$LP)$$
     member? (p,ts) ==
       member?(p,rep(ts))$LP
     roughUnitIdeal? ts ==
       false
     coerce(ts:$) : OutputForm ==
       lp : List(P) := reverse(rep(ts))
       brace([p::OutputForm for p in lp]$List(OutputForm))$OutputForm
     mvar ts ==
       empty? ts => error "mvar$REGSET: #1 is empty"
       mvar(first(rep(ts)))$P
     first ts ==
       empty? ts => "failed"::Union(P,"failed")
       first(rep(ts))::Union(P,"failed")
     last ts ==
       empty? ts => "failed"::Union(P,"failed")
       last(rep(ts))::Union(P,"failed")
     rest ts ==
       empty? ts => "failed"::Union($,"failed")
       per(rest(rep(ts)))::Union($,"failed")
     coerce(ts:$) : (List P) ==
       rep(ts)

     collectUpper (ts,v) ==
       empty? ts => ts
       lp := rep(ts)
       newlp : Rep := []
       while (not empty? lp) and (mvar(first(lp)) > v) repeat
         newlp := cons(first(lp),newlp)
         lp := rest lp
       per(reverse(newlp))

     collectUnder (ts,v) ==
       empty? ts => ts
       lp := rep(ts)
       while (not empty? lp) and (mvar(first(lp)) >= v) repeat
         lp := rest lp
       per(lp)

     construct(lp:List(P)) ==
       ts : $ := per([])
       empty? lp => ts
       lp := sort(infRittWu?,lp)
       while not empty? lp repeat
         eif := extendIfCan(ts,first(lp))
         not (eif case $) =>
           error"in construct : List P -> $  from REGSET : bad #1"
         ts := eif::$
         lp := rest lp
       ts

     extendIfCan(ts:$,p:P) ==
       ground? p => "failed"::Union($,"failed")       
       empty? ts => 
         p := primitivePart p
         (per([p]))::Union($,"failed")
       not (mvar(ts) < mvar(p)) => "failed"::Union($,"failed")
       invertible?(init(p),ts)@Boolean => 
         (per(cons(p,rep(ts))))::Union($,"failed")
       "failed"::Union($,"failed")

     removeZero(p:P, ts:$): P ==
       (ground? p) or (empty? ts) => p
       v := mvar(p)
       ts_v_- := collectUnder(ts,v)
       if algebraic?(v,ts) 
         then
           q := lazyPrem(p,select(ts,v)::P)
           zero? q => return q
           zero? removeZero(q,ts_v_-) => return 0
       empty? ts_v_- => p
       q: P := 0
       while positive? degree(p,v) repeat
          q := removeZero(init(p),ts_v_-) * mainMonomial(p) + q
          p := tail(p)
       q + removeZero(p,ts_v_-)

     internalAugment(p:P,ts:$): $ ==
       -- ASSUME that adding p to ts DOES NOT require any split
       ground? p => error "in internalAugment$REGSET: ground? #1"
       first(internalAugment(p,ts,false,false,false,false,false))

     internalAugment(lp:List(P),ts:$): $ ==
       -- ASSUME that adding p to ts DOES NOT require any split
       empty? lp => ts
       internalAugment(rest lp, internalAugment(first lp, ts))

     internalAugment(p:P,ts:$,rem?:B,red?:B,prim?:B,sqfr?:B,extend?:B): Split ==
       -- ASSUME p is not a constant
       -- ASSUME mvar(p) is not algebraic w.r.t. ts
       -- ASSUME init(p) invertible modulo ts
       -- if rem? then REDUCE p by remainder
       -- if prim? then REPLACE p by its main primitive part
       -- if sqfr? then FACTORIZE SQUARE FREE p over R
       -- if extend? DO NOT ASSUME every pol in ts_v_+ is invertible modulo ts
       v := mvar(p)
       ts_v_- := collectUnder(ts,v)
       ts_v_+ := collectUpper(ts,v)
       if rem? then p := remainder(p,ts_v_-).polnum
       -- if rem? then p := reduceByQuasiMonic(p,ts_v_-)
       if red? then p := removeZero(p,ts_v_-)
       if prim? then p := mainPrimitivePart p
       if sqfr?
         then
           lsfp := squareFreeFactors(p)$polsetpack
           lts: Split := [per(cons(f,rep(ts_v_-))) for f in lsfp]
         else
           lts: Split := [per(cons(p,rep(ts_v_-)))]
       extend? => extend(members(ts_v_+),lts)
       [per(concat(rep(ts_v_+),rep(us))) for us in lts]

     augment(p:P,ts:$): List $ ==
       ground? p => error "in augment$REGSET: ground? #1"
       algebraic?(mvar(p),ts) => error "in augment$REGSET: bad #1"
       -- ASSUME init(p) invertible modulo ts
       -- DOES NOT ASSUME anything else.
       -- THUS reduction, mainPrimitivePart and squareFree are NEEDED
       internalAugment(p,ts,true,true,true,true,true)

     extend(p:P,ts:$): List $ ==
       ground? p => error "in extend$REGSET: ground? #1"
       v := mvar(p)
       not (mvar(ts) < mvar(p)) => error "in extend$REGSET: bad #1"
       lts: List($) := []
       split: List($) := invertibleSet(init(p),ts)
       for us in split repeat
         lts := concat(augment(p,us),lts)
       lts

     invertible?(p:P,ts:$): Boolean == 
       toseInvertible?(p,ts)$regsetgcdpack
       
     invertible?(p:P,ts:$): List BWT ==
       toseInvertible?(p,ts)$regsetgcdpack

     invertibleSet(p:P,ts:$): Split ==
       toseInvertibleSet(p,ts)$regsetgcdpack

     lastSubResultant(p1:P,p2:P,ts:$): List PWT ==
       toseLastSubResultant(p1,p2,ts)$regsetgcdpack

     squareFreePart(p:P, ts: $): List PWT ==
       toseSquareFreePart(p,ts)$regsetgcdpack

     intersect(p:P, ts: $): List($) == decompose([p], [ts], false, false)$regsetdecomppack

     intersect(lp: LP, lts: List($)): List($) == decompose(lp, lts, false, false)$regsetdecomppack
        -- SOLVE in the regular zero sense 
        -- and DO NOT PRINT info

     zeroSetSplit(lp:List(P)) == zeroSetSplit(lp,true,false)
        -- by default SOLVE in the closure sense 
        -- and DO NOT PRINT info

     zeroSetSplit(lp:List(P), clos?: B) == zeroSetSplit(lp,clos?, false)
        -- DO NOT PRINT info

     zeroSetSplit(lp:List(P), clos?: B, info?: B) ==
       -- if clos? then SOLVE in the closure sense 
       -- if info? then PRINT info
       -- by default USE hash-tables
       -- and PREPROCESS the input system
       zeroSetSplit(lp,true,clos?,info?,true)

     zeroSetSplit(lp:List(P),hash?:B,clos?:B,info?:B,prep?:B) == 
       -- if hash? then USE hash-tables
       -- if info? then PRINT information
       -- if clos? then SOLVE in the closure sense
       -- if prep? then PREPROCESS the input system
       if hash? 
         then
           s1, s2, s3, dom1, dom2, dom3: String
           e: String := empty()$String
           if info? then (s1,s2,s3) := ("w","g","i") else (s1,s2,s3) := (e,e,e)
           if info? 
             then 
               (dom1, dom2, dom3) := ("QCMPACK", "REGSETGCD: Gcd", "REGSETGCD: Inv Set")
             else
               (dom1, dom2, dom3) := (e,e,e)
           startTable!(s1,"W",dom1)$quasicomppack
           startTableGcd!(s2,"G",dom2)$regsetgcdpack
           startTableInvSet!(s3,"I",dom3)$regsetgcdpack
       lts := internalZeroSetSplit(lp,clos?,info?,prep?)
       if hash? 
         then
           stopTable!()$quasicomppack
           stopTableGcd!()$regsetgcdpack
           stopTableInvSet!()$regsetgcdpack
       lts

     internalZeroSetSplit(lp:LP,clos?:B,info?:B,prep?:B) ==
       -- if info? then PRINT information
       -- if clos? then SOLVE in the closure sense
       -- if prep? then PREPROCESS the input system
       if prep?
         then
           pp := pre_process(lp,clos?,info?)
           lp := pp.val
           lts := pp.towers
         else
           ts: $ := [[]]
           lts := [ts]
       lp := remove(zero?, lp)
       any?(ground?, lp) => []
       empty? lp => lts
       empty? lts => lts
       lp := sort(infRittWu?,lp)
       clos? => decompose(lp,lts, clos?, info?)$regsetdecomppack
       -- IN DIM > 0 with clos? the following is false ...
       for p in lp repeat
         lts := decompose([p],lts, clos?, info?)$regsetdecomppack
       lts

     largeSystem?(lp:LP): Boolean == 
       -- Gonnet and Gerdt and not Wu-Wang.2
       #lp > 16 => true
       #lp < 13 => false
       lts: List($) := []
       (#lp :: Z - numberOfVariables(lp,lts)$regsetdecomppack :: Z) > 3

     smallSystem?(lp:LP): Boolean == 
       -- neural, Vermeer, Liu, and not f-633 and not Hairer-2
       #lp < 5

     mediumSystem?(lp:LP): Boolean == 
       -- f-633 and not Hairer-2
       lts: List($) := []
       (numberOfVariables(lp,lts)$regsetdecomppack :: Z - #lp :: Z) < 2

     lin?(p:P):Boolean == ground?(init(p)) and one?(mdeg(p))

     pre_process(lp:LP,clos?:B,info?:B): Record(val: LP, towers: Split) ==
       -- if info? then PRINT information
       -- if clos? then SOLVE in the closure sense
       ts: $ := [[]]; 
       lts: Split := [ts]
       empty? lp => [lp,lts]
       lp1: List P := []
       lp2: List P := []
       for p in lp repeat 
          ground? (tail p) => lp1 := cons(p, lp1)
          lp2 := cons(p, lp2)
       lts: Split := decompose(lp1,[ts],clos?,info?)$regsetdecomppack
       probablyZeroDim?(lp)$polsetpack =>
          largeSystem?(lp) => return [lp2,lts]
          if #lp > 7
            then 
              -- Butcher (8,8) + Wu-Wang.2 (13,16) 
              lp2 := crushedSet(lp2)$polsetpack
              lp2 := remove(zero?,lp2)
              any?(ground?,lp2) => return [lp2, lts]
              lp3 := [p for p in lp2 | lin?(p)]
              lp4 := [p for p in lp2 | not lin?(p)]
              if clos?
                then 
                  lts := decompose(lp4,lts, clos?, info?)$regsetdecomppack
                else
                  lp4 := sort(infRittWu?,lp4)
                  for p in lp4 repeat
                    lts := decompose([p],lts, clos?, info?)$regsetdecomppack
              lp2 := lp3
            else
              lp2 := crushedSet(lp2)$polsetpack
              lp2 := remove(zero?,lp2)
              any?(ground?,lp2) => return [lp2, lts]
          if clos?
            then
              lts := decompose(lp2,lts, clos?, info?)$regsetdecomppack
            else
              lp2 := sort(infRittWu?,lp2)
              for p in lp2 repeat
                lts := decompose([p],lts, clos?, info?)$regsetdecomppack
          lp2 := []
          return [lp2,lts]
       smallSystem?(lp) => [lp2,lts]
       mediumSystem?(lp) => [crushedSet(lp2)$polsetpack,lts]
       lp3 := [p for p in lp2 | lin?(p)]
       lp4 := [p for p in lp2 | not lin?(p)]
       if clos?
         then 
           lts := decompose(lp4,lts, clos?, info?)$regsetdecomppack
         else
           lp4 := sort(infRittWu?,lp4)
           for p in lp4 repeat
             lts := decompose([p],lts, clos?, info?)$regsetdecomppack
       if clos?
         then 
           lts := decompose(lp3,lts, clos?, info?)$regsetdecomppack
         else
           lp3 := sort(infRittWu?,lp3)
           for p in lp3 repeat
             lts := decompose([p],lts, clos?, info?)$regsetdecomppack
       lp2 := []
       return [lp2,lts]

@
\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>>

<<category RSETCAT RegularTriangularSetCategory>>
<<package QCMPACK QuasiComponentPackage>>
<<package RSETGCD RegularTriangularSetGcdPackage>>
<<package RSDCMPK RegularSetDecompositionPackage>>
<<domain REGSET RegularTriangularSet>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}