Initial population.

1 files changed, 1606 insertions, 0 deletions

diff --git a/src/algebra/sregset.spad.pamphlet b/src/algebra/sregset.spad.pamphlet
new file mode 100644
index 00000000..19886881
--- /dev/null
+++ b/src/algebra/sregset.spad.pamphlet
@@ -0,0 +1,1606 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra sregset.spad}
+\author{Marc Moreno Maza}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{category SFRTCAT SquareFreeRegularTriangularSetCategory}
+<<category SFRTCAT SquareFreeRegularTriangularSetCategory>>=
+)abbrev category SFRTCAT SquareFreeRegularTriangularSetCategory
+++ Author: Marc Moreno Maza
+++ Date Created: 09/03/1996
+++ Date Last Updated: 09/10/1998  
+++ Basic Functions:
+++ Related Constructors:
+++ Also See: essai Graphisme
+++ AMS Classifications:
+++ Keywords: polynomial, multivariate, ordered variables set
+++ Description:
+++ The category of square-free regular triangular sets.
+++ A regular triangular set \spad{ts} is square-free if
+++ the gcd of any polynomial \spad{p} in \spad{ts} and 
+++ \spad{differentiate(p,mvar(p))} w.r.t. 
+++ \axiomOpFrom{collectUnder}{TriangularSetCategory}(ts,\axiomOpFrom{mvar}{RecursivePolynomialCategory}(p))
+++ has degree zero w.r.t. \spad{mvar(p)}. Thus any square-free regular
+++ set defines a tower of square-free simple extensions.\newline
+++ References :
+++  [1] D. LAZARD "A new method for solving algebraic systems of 
+++      positive dimension" Discr. App. Math. 33:147-160,1991
+++  [2] M. KALKBRENER "Algorithmic properties of polynomial rings"
+++      Habilitation Thesis, ETZH, Zurich, 1995.
+++  [3] M. MORENO MAZA "A new algorithm for computing triangular
+++      decomposition of algebraic varieties" NAG Tech. Rep. 4/98.
+
+
+
+SquareFreeRegularTriangularSetCategory(R:GcdDomain,E:OrderedAbelianMonoidSup,_
+  V:OrderedSet,P:RecursivePolynomialCategory(R,E,V)):
+         Category == 
+   RegularTriangularSetCategory(R,E,V,P) 
+
+@
+\section{package SFQCMPK SquareFreeQuasiComponentPackage}
+<<package SFQCMPK SquareFreeQuasiComponentPackage>>=
+)abbrev package SFQCMPK SquareFreeQuasiComponentPackage
+++ Author: Marc Moreno Maza
+++ Date Created: 09/23/1998
+++ Date Last Updated: 12/16/1998
+++ Basic Functions:
+++ Related Constructors:
+++ Also See: 
+++ AMS Classifications:
+++ Keywords:
+++ Description: 
+++   A internal 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
+++      Tech. Report (PoSSo project)
+++  [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: 1. 
+
+SquareFreeQuasiComponentPackage(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)
+  SQUAREFREE ==> SquareFreeRegularTriangularSetCategory(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 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?(ts,us)}{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
+       void()
+
+     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)
+
+     if TS has SquareFreeRegularTriangularSetCategory(R,E,V,P)
+     then
+
+       internalSubQuasiComponent?(ts:TS,us:TS): Union(Boolean,"failed") ==
+         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 
+           b: B := invertible?(p,ts)$TS
+           not b => 
+             return(false::Union(Boolean,"failed"))
+         true::Union(Boolean,"failed")
+
+     else
+
+       internalSubQuasiComponent?(ts:TS,us:TS): Union(Boolean,"failed") ==
+         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? reduceByQuasiMonic(p,ts) =>
+             return("failed"::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,lpwt2 : 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 in lineq repeat
+              zero? initiallyReduce(q,ts) => 
+                return("failed"::Union(Branch,"failed"))
+        if b5
+          then
+            newleq: LP := []
+            for p in leq repeat
+              for q 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 SFRGCD SquareFreeRegularTriangularSetGcdPackage}
+<<package SFRGCD SquareFreeRegularTriangularSetGcdPackage>>=
+)abbrev package SFRGCD SquareFreeRegularTriangularSetGcdPackage
+++ Author: Marc Moreno Maza
+++ Date Created: 09/23/1998
+++ Date Last Updated: 10/01/1998
+++ Basic Functions:
+++ Related Constructors:
+++ Also See: 
+++ AMS Classifications:
+++ Keywords:
+++ Description: 
+++ A internal package for computing gcds and resultants of univariate polynomials
+++ with coefficients in a tower of simple extensions of a field.
+++ There is no need to use directly this package since its main operations are
+++ available from \spad{TS}. \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: 1. 
+
+SquareFreeRegularTriangularSetGcdPackage(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)
+  iprintpack ==> InternalPrintPackage()
+  polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,P)
+  quasicomppack ==> SquareFreeQuasiComponentPackage(R,E,V,P,TS)
+
+  SQUAREFREE ==> SquareFreeRegularTriangularSetCategory(R,E,V,P)
+
+  Exports ==  with
+     startTableGcd!: (S,S,S) -> Void
+     stopTableGcd!: () -> Void
+     startTableInvSet!: (S,S,S) -> Void
+     stopTableInvSet!: () -> Void
+     stosePrepareSubResAlgo: (P,P,TS) -> List LpWT 
+     stoseInternalLastSubResultant: (P,P,TS,B,B) -> List PWT 
+     stoseInternalLastSubResultant: (List LpWT,V,B) -> List PWT 
+     stoseIntegralLastSubResultant: (P,P,TS) -> List PWT 
+     stoseLastSubResultant: (P,P,TS) -> List PWT 
+     stoseInvertible?: (P,TS) -> B
+     stoseInvertible?_sqfreg: (P,TS) -> List BWT
+     stoseInvertibleSet_sqfreg: (P,TS) -> Split
+     stoseInvertible?_reg: (P,TS) -> List BWT
+     stoseInvertibleSet_reg: (P,TS) -> Split
+     stoseInvertible?: (P,TS) -> List BWT
+     stoseInvertibleSet: (P,TS) -> Split
+     stoseSquareFreePart: (P,TS) -> List PWT 
+
+  Implementation == add
+
+     startTableGcd!(ok: S, ko: S, domainName: S): Void == 
+       initTable!()$HGcd
+       printInfo!(ok,ko)$HGcd
+       startStats!(domainName)$HGcd
+       void()
+
+     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
+       void()
+
+     stopTableInvSet!(): Void ==   
+       if makingStats?()$HInvSet then printStats!()$HInvSet
+       clearTable!()$HInvSet
+
+     stoseInvertible?(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 := stoseInvertible?(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 := stoseInternalLastSubResultant(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
+
+     stosePrepareSubResAlgo(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 := stoseInvertible?(leadingCoefficient(p2,v),ts)@(List BWT)
+         for bwt in lbwt repeat
+           (bwt.val = true) and (degree(p2,v) > 0) =>
+             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],bwt.tower]$LpWT,toSee)
+           toSee := cons([[p1,newp2],bwt.tower]$LpWT,toSee)
+       toSave
+
+     stoseIntegralLastSubResultant(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]
+            
+     stoseInternalLastSubResultant(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 := stoseIntegralLastSubResultant(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 := stosePrepareSubResAlgo(p1,p2,ts)
+       toSave := stoseInternalLastSubResultant(toSee,mvar(p1),b2)
+       insert!(k,toSave)$HGcd
+       toSave
+
+     stoseInternalLastSubResultant(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 := stoseInvertible?(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 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 in toSee repeat 
+           llpwt := cons([[p1,p2,s],lpwt.tower]$LpWT, llpwt)
+       toReturn
+
+     stoseLastSubResultant(p1:P,p2:P,ts:TS): List PWT ==
+       ground? p1 => 
+         error"in stoseLastSubResultantElseSplit$SFRGCD  : bad #1"
+       ground? p2 => 
+         error"in stoseLastSubResultantElseSplit$SFRGCD : bad #2"
+       not (mvar(p2) = mvar(p1)) => 
+         error"in stoseLastSubResultantElseSplit$SFRGCD : bad #2"
+       algebraic?(mvar(p1),ts) =>
+         error"in stoseLastSubResultantElseSplit$SFRGCD : bad #1"
+       not initiallyReduced?(p1,ts) => 
+         error"in stoseLastSubResultantElseSplit$SFRGCD : bad #1"
+       not initiallyReduced?(p2,ts) => 
+         error"in stoseLastSubResultantElseSplit$SFRGCD : bad #2"
+       purelyTranscendental?(p1,ts) and purelyTranscendental?(p2,ts) =>
+         stoseIntegralLastSubResultant(p1,p2,ts)
+       if mdeg(p1) < mdeg(p2) then 
+          (p1, p2) := (p2, p1)
+          if odd?(mdeg(p1)) and odd?(mdeg(p2)) then p2 := - p2
+       stoseInternalLastSubResultant(p1,p2,ts,false,false)
+
+     stoseSquareFreePart_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]
+       mdeg(p) = 1 => [[p,ts]$PWT]
+       v := mvar(p)$P
+       q: P := mainPrimitivePart D(p,v)
+       lgwt: List PWT := stoseInternalLastSubResultant(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
+
+     stoseSquareFreePart_base(p:P, ts: TS): List PWT == [[p,ts]$PWT]
+
+     stoseSquareFreePart(p:P, ts: TS): List PWT == stoseSquareFreePart_wip(p,ts)
+       
+     stoseInvertible?_sqfreg(p:P,ts:TS): List BWT ==
+       --iprint("+")$iprintpack
+       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 := stoseInvertible?_sqfreg(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(stoseInvertible?_sqfreg(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 := stoseInternalLastSubResultant(ts_v,q,ts_v_-,false,false)
+       lbwt: List BWT := []
+       lts, lts_g, lts_h: Split
+       for gwt in lgwt repeat
+         g := gwt.val; ts := gwt.tower
+         (ground? g) or (mvar(g) < v) => 
+           lts := augment(ts_v,ts)$TS
+           lts := augment(members(ts_v_+),lts)$TS
+           for ts in lts repeat
+             lbwt := cons([true, ts]$BWT,lbwt)
+         g := mainPrimitivePart g
+         lts_g := augment(g,ts)$TS
+         lts_g := augment(members(ts_v_+),lts_g)$TS
+         -- USE stoseInternalAugment with parameters ??
+         for ts_g in lts_g repeat
+           lbwt := cons([false, ts_g]$BWT,lbwt)
+         h := lazyPquo(ts_v,g)
+         (ground? h) or (mvar(h) < v) => "leave"
+         h := mainPrimitivePart h
+         lts_h := augment(h,ts)$TS
+         lts_h := augment(members(ts_v_+),lts_h)$TS
+         -- USE stoseInternalAugment with parameters ??
+         for ts_h in lts_h repeat
+           lbwt := cons([true, ts_h]$BWT,lbwt)
+       sort(#1.val < #2.val,lbwt)
+
+     stoseInvertibleSet_sqfreg(p:P,ts:TS): Split ==
+       --iprint("*")$iprintpack
+       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 := stoseInvertible?_sqfreg(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(stoseInvertibleSet_sqfreg(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 := stoseInternalLastSubResultant(ts_v,q,ts_v_-,false,false)
+       lts, lts_h: Split
+       for gwt in lgwt repeat
+         g := gwt.val; ts := gwt.tower
+         (ground? g) or (mvar(g) < v) => 
+           lts := augment(ts_v,ts)$TS
+           lts := augment(members(ts_v_+),lts)$TS
+           toSave := concat(lts,toSave)
+         g := mainPrimitivePart g
+         h := lazyPquo(ts_v,g)
+         h := mainPrimitivePart h
+         (ground? h) or (mvar(h) < v) => "leave"
+         lts_h := augment(h,ts)$TS
+         lts_h := augment(members(ts_v_+),lts_h)$TS
+         toSave := concat(lts_h,toSave)
+       toSave := algebraicSort(toSave)$quasicomppack
+       insert!(k,toSave)$HInvSet
+       toSave
+       
+     stoseInvertible?_reg(p:P,ts:TS): List BWT ==
+       --iprint("-")$iprintpack
+       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 := stoseInvertible?_reg(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(stoseInvertible?_reg(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 := stoseInternalLastSubResultant(ts_v,q,ts_v_-,false,false)
+       lbwt: List BWT := []
+       lts, lts_g, lts_h: Split
+       for gwt in lgwt repeat
+         g := gwt.val; ts := gwt.tower
+         (ground? g) or (mvar(g) < v) => 
+           lts := augment(ts_v,ts)$TS
+           lts := augment(members(ts_v_+),lts)$TS
+           for ts in lts repeat
+             lbwt := cons([true, ts]$BWT,lbwt)
+         g := mainPrimitivePart g
+         lts_g := augment(g,ts)$TS
+         lts_g := augment(members(ts_v_+),lts_g)$TS
+         -- USE internalAugment with parameters ??
+         for ts_g in lts_g repeat
+           lbwt := cons([false, ts_g]$BWT,lbwt)
+         h := lazyPquo(ts_v,g)
+         (ground? h) or (mvar(h) < v) => "leave"
+         h := mainPrimitivePart h
+         lts_h := augment(h,ts)$TS
+         lts_h := augment(members(ts_v_+),lts_h)$TS
+         -- USE internalAugment with parameters ??
+         for ts_h in lts_h repeat
+           inv := stoseInvertible?_reg(q,ts_h)@(List BWT)
+           lbwt := concat([bwt for bwt in inv | bwt.val],lbwt)
+       sort(#1.val < #2.val,lbwt)
+
+     stoseInvertibleSet_reg(p:P,ts:TS): Split ==
+       --iprint("/")$iprintpack
+       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 := stoseInvertible?_reg(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(stoseInvertibleSet_reg(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 := stoseInternalLastSubResultant(ts_v,q,ts_v_-,false,false)
+       lts, lts_h: Split
+       for gwt in lgwt repeat
+         g := gwt.val; ts := gwt.tower
+         (ground? g) or (mvar(g) < v) => 
+           lts := augment(ts_v,ts)$TS
+           lts := augment(members(ts_v_+),lts)$TS
+           toSave := concat(lts,toSave)
+         g := mainPrimitivePart g
+         h := lazyPquo(ts_v,g)
+         h := mainPrimitivePart h
+         (ground? h) or (mvar(h) < v) => "leave"
+         lts_h := augment(h,ts)$TS
+         lts_h := augment(members(ts_v_+),lts_h)$TS
+         for ts_h in lts_h repeat
+           inv := stoseInvertibleSet_reg(q,ts_h)
+           toSave := removeDuplicates concat(inv,toSave)
+       toSave := algebraicSort(toSave)$quasicomppack
+       insert!(k,toSave)$HInvSet
+       toSave
+
+     if TS has SquareFreeRegularTriangularSetCategory(R,E,V,P)
+     then
+
+       stoseInvertible?(p:P,ts:TS): List BWT == stoseInvertible?_sqfreg(p,ts)
+
+       stoseInvertibleSet(p:P,ts:TS): Split == stoseInvertibleSet_sqfreg(p,ts)
+
+     else
+       
+       stoseInvertible?(p:P,ts:TS): List BWT == stoseInvertible?_reg(p,ts)
+ 
+       stoseInvertibleSet(p:P,ts:TS): Split == stoseInvertibleSet_reg(p,ts)
+
+@
+\section{package SRDCMPK SquareFreeRegularSetDecompositionPackage}
+<<package SRDCMPK SquareFreeRegularSetDecompositionPackage>>=
+)abbrev package SRDCMPK SquareFreeRegularSetDecompositionPackage
+++ Author: Marc Moreno Maza
+++ Date Created: 09/23/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 provided:
+++ 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-
+++ Moreno 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 \spad{QCMPPK(R,E,V,P,TS)} and \spad{RSETGCD(R,E,V,P,TS)}.
+++ The same way it does not care about the way univariate polynomial
+++ gcds (with coefficients in the tower of simple extensions associated 
+++ with a regular set) are computed. The only requirement is that these
+++ gcds 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 \axiomType{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: 2. Does not use any unproved criteria.
+
+SquareFreeRegularSetDecompositionPackage(R,E,V,P,TS): Exports == Implementation where
+
+  R : GcdDomain
+  E : OrderedAbelianMonoidSup
+  V : OrderedSet
+  P : RecursivePolynomialCategory(R,E,V)
+  TS : SquareFreeRegularTriangularSetCategory(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 ==> SquareFreeQuasiComponentPackage(R,E,V,P,TS)
+  regsetgcdpack ==> SquareFreeRegularTriangularSetGcdPackage(R,E,V,P,TS)
+
+  Exports ==  with
+
+     KrullNumber: (LP, Split) -> N
+     numberOfVariables: (LP, Split) -> N
+     algebraicDecompose: (P,TS) -> 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): 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
+       lgwt: List PWT
+       if mdeg(p) < mdeg(ts_v)
+         then 
+           lgwt := stoseInternalLastSubResultant(ts_v,p,ts_v_-,true,false)$regsetgcdpack
+         else
+           lgwt := stoseInternalLastSubResultant(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"
+         h := leadingCoefficient(g,v)
+         lus := augment(members(ts_v_+),augment(ts_v,us)$TS)$TS
+         lsfp := squareFreeFactors(h)$polsetpack
+         for f in lsfp repeat
+           ground? f => "leave"
+           for vs in lus repeat
+             llpwt := cons([[f,p],vs]$LpWT, llpwt)
+         n < #us => 
+           error " in algebraicDecompose$REGSET: should never happen !!!"
+         mvar(g) = v => 
+           lts := concat(augment(members(ts_v_+),augment(g,us)$TS)$TS,lts)         
+       [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)$TS
+         else
+           lts := []
+       llpwt: List LpWT := []
+       [lts,llpwt]
+
+     transcendentalDecompose(p: P, ts: TS): Record(done: Split, todo: List LpWT) ==
+       lts: Split:= augment(p,ts)$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: List BWT := stoseInvertible?_sqfreg(ip,ts)$regsetgcdpack
+       for bwt in lbwt repeat
+         bwt.val =>
+           if algebraic?(mvar(p),bwt.tower) 
+             then 
+               rsl := algebraicDecompose(p,bwt.tower)
+             else
+               rsl := transcendentalDecompose(p,bwt.tower,bound)
+           lts := concat(rsl.done,lts)
+           llpwt :=  concat(rsl.todo,llpwt)
+           (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: List BWT := stoseInvertible?_sqfreg(ip,ts)$regsetgcdpack
+       for bwt in lbwt repeat
+         bwt.val =>
+           if algebraic?(mvar(p),bwt.tower) 
+             then 
+               rsl := algebraicDecompose(p,bwt.tower)
+             else
+               rsl := transcendentalDecompose(p,bwt.tower)
+           lts := concat(rsl.done,lts)
+           llpwt :=  concat(rsl.todo,llpwt)
+           (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 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
+       void()
+
+     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 SREGSET SquareFreeRegularTriangularSet}
+<<domain SREGSET SquareFreeRegularTriangularSet>>=
+)abbrev domain SREGSET SquareFreeRegularTriangularSet
+++ 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 square-free regular chains.
+++ Moreover, the operation \axiomOpFrom{zeroSetSplit}{SquareFreeRegularTriangularSetCategory}
+++ 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: 2
+
+SquareFreeRegularTriangularSet(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 ==> SquareFreeQuasiComponentPackage(R,E,V,P,$)
+  regsetgcdpack ==> SquareFreeRegularTriangularSetGcdPackage(R,E,V,P,$)
+  regsetdecomppack ==> SquareFreeRegularSetDecompositionPackage(R,E,V,P,$)
+
+  Exports ==  SquareFreeRegularTriangularSetCategory(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} 
+       ++ from \spadtype{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
+
+     rep(s:$):Rep == s pretend Rep
+     per(l:Rep):$ == l pretend $
+
+     copy ts ==
+       per(copy(rep(ts))$LP)
+     empty() ==
+       per([])
+     empty?(ts:$) ==
+       empty?(rep(ts))
+     parts ts ==
+       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
+     unitIdealIfCan() ==
+       "failed"::Union($,"failed")
+     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$SREGSET: #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 SREGSET : bad #1"
+         ts := eif::$
+         lp := rest lp
+       ts
+
+     extendIfCan(ts:$,p:P) ==
+       ground? p => "failed"::Union($,"failed")       
+       empty? ts => 
+         p := squareFreePart primitivePart p
+         (per([p]))::Union($,"failed")
+       not (mvar(ts) < mvar(p)) => "failed"::Union($,"failed")
+       invertible?(init(p),ts)@Boolean => 
+         lts: Split := augment(p,ts)
+         #lts ~= 1 => "failed"::Union($,"failed")
+         (first lts)::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$SREGSET: 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
+       lts: Split
+       if sqfr?
+         then
+           lts: Split := []
+           lsfp := squareFreeFactors(p)$polsetpack
+           for f in lsfp repeat
+             (ground? f) or (mvar(f) < v) => "leave"
+             lpwt := squareFreePart(f,ts_v_-)
+             for pwt in lpwt repeat 
+               sfp := pwt.val; us := pwt.tower
+               lts := cons( per(cons(pwt.val, rep(pwt.tower))), lts)
+         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$SREGSET: ground? #1"
+       algebraic?(mvar(p),ts) => error "in augment$SREGSET: 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$SREGSET: ground? #1"
+       v := mvar(p)
+       not (mvar(ts) < mvar(p)) => error "in extend$SREGSET: bad #1"
+       split: List($) := invertibleSet(init(p),ts)
+       lts: List($) := []
+       for us in split repeat
+         lts := concat(augment(p,us),lts)
+       lts
+
+     invertible?(p:P,ts:$): Boolean == 
+       stoseInvertible?(p,ts)$regsetgcdpack
+       
+     invertible?(p:P,ts:$): List BWT ==
+       stoseInvertible?_sqfreg(p,ts)$regsetgcdpack
+
+     invertibleSet(p:P,ts:$): Split ==
+       stoseInvertibleSet_sqfreg(p,ts)$regsetgcdpack
+
+     lastSubResultant(p1:P,p2:P,ts:$): List PWT ==
+       stoseLastSubResultant(p1,p2,ts)$regsetgcdpack
+
+     squareFreePart(p:P, ts: $): List PWT ==
+       stoseSquareFreePart(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
+
+     decompose(p:P, ts: $): List($) == decompose([p], [ts], true, false)$regsetdecomppack
+
+     decompose(lp: LP, lts: List($)): List($) == decompose(lp, lts, true, false)$regsetdecomppack
+        -- SOLVE in the closure 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)
+
+     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 not 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))
+     lin?(p:P):Boolean == ground?(init(p)) and (mdeg(p) = 1)
+
+     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 SFRTCAT SquareFreeRegularTriangularSetCategory>>
+<<package SFQCMPK SquareFreeQuasiComponentPackage>>
+<<package SFRGCD SquareFreeRegularTriangularSetGcdPackage>>
+<<package SRDCMPK SquareFreeRegularSetDecompositionPackage>>
+<<domain SREGSET SquareFreeRegularTriangularSet>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}