triset.spad

)abbrev category TSETCAT TriangularSetCategory
++ Author: Marc Moreno Maza (marc@nag.co.uk)
++ Date Created: 04/26/1994
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords: polynomial, multivariate, ordered variables set
++ Description:
++ The category of triangular sets of multivariate polynomials
++ with coefficients in an integral domain.
++ Let \spad{R} be an integral domain and \spad{V} a finite ordered set of
++ variables, say \spad{X1 < X2 < ... < Xn}.
++ A set \spad{S} of polynomials in \spad{R[X1, X2, ..., Xn]} is triangular
++ if no elements of \spad{S} lies in \spad{R}, and if two distinct
++ elements of \spad{S} have distinct main variables.
++ Note that the empty set is a triangular set. A triangular set is not
++ necessarily a (lexicographical) Groebner basis and the notion of
++ reduction related to triangular sets is based on the recursive view
++ of polynomials. We recall this notion here and refer to [1] for more details.
++ A polynomial \spad{P} is reduced w.r.t a non-constant polynomial
++ \spad{Q} if the degree of \spad{P} in the main variable of \spad{Q}
++ is less than the main degree of \spad{Q}.
++ A polynomial \spad{P} is reduced w.r.t a triangular set \spad{T}
++ if it is reduced w.r.t. every polynomial of \spad{T}. \newline
++ References:
++  [1] P. AUBRY, D. LAZARD and M. MORENO MAZA "On the Theories
++      of Triangular Sets" Journal of Symbol. Comp. 28:105-124, 1999.
++ Version: 4.

TriangularSetCategory(R : IntegralDomain, E : OrderedAbelianMonoidSup, _
 V : OrderedSet, P : RecursivePolynomialCategory(R, E, V)):
         Category ==
   Join(PolynomialSetCategory(R, E, V, P), finiteAggregate, shallowlyMutable) with

     infRittWu? : (%, %) -> Boolean
         ++ \spad{infRittWu?(ts1, ts2)} returns true iff \spad{ts2} has higher rank
         ++ than \spad{ts1} in Wu Wen Tsun sense.
     basicSet : (List P,((P,P)->Boolean)) -> Union(Record(bas:%,top:List P),"failed")
         ++ \spad{basicSet(ps, redOp?)} returns \spad{[bs, ts]} where
         ++ \spad{concat(bs, ts)} is \spad{ps} and \spad{bs}
         ++ is a basic set in Wu Wen Tsun sense of \spad{ps} w.r.t
         ++ the reduction-test \spad{redOp?}, if no non-zero constant
         ++ polynomial lie in \spad{ps}, otherwise \spad{"failed"} is returned.
     basicSet : (List P,(P->Boolean),((P,P)->Boolean)) -> Union(Record(bas:%,top:List P),"failed")
         ++ \spad{basicSet(ps, pred?, redOp?)} returns the same as \spad{basicSet(qs, redOp?)}
         ++ where \spad{qs} consists of the polynomials of \spad{ps}
         ++ satisfying property \spad{pred?}.
     initials : % -> List P
         ++ \spad{initials(ts)} returns the list of the non-constant initials
         ++ of the members of \spad{ts}.
     degree : % -> NonNegativeInteger
         ++ \spad{degree(ts)} returns the product of main degrees of the
         ++ members of \spad{ts}.
     quasiComponent : % -> Record(close : List P, open : List P)
         ++ \spad{quasiComponent(ts)} returns \spad{[lp, lq]} where \spad{lp} is the list
         ++ of the members of \spad{ts} and \spad{lq}is \spad{initials(ts)}.
     normalized? : (P, %) -> Boolean
         ++ \spad{normalized?(p, ts)} returns true iff \spad{p} and all its iterated initials
         ++ have degree zero w.r.t. the main variables of the polynomials of \spad{ts}
     normalized? : %  -> Boolean
         ++ \spad{normalized?(ts)} returns true iff for every \spad{p} in \spad{ts} we have
         ++ \spad{normalized?(p, us)} where \spad{us} is \spad{collectUnder(ts, mvar(p))}.
     reduced? : (P, %, ((P, P) -> Boolean)) -> Boolean
         ++ \spad{reduced?(p, ts, redOp?)} returns true iff \spad{p} is reduced w.r.t.
         ++ in the sense of the operation \spad{redOp?}, that is if for every \spad{t} in
         ++ \spad{ts} \spad{redOp?(p, t)} holds.
     stronglyReduced? : (P, %) -> Boolean
         ++ \spad{stronglyReduced?(p, ts)} returns true iff \spad{p}
         ++ is reduced w.r.t. \spad{ts}.
     headReduced? : (P, %) -> Boolean
         ++ \spad{headReduced?(p, ts)} returns true iff the head of \spad{p} is
         ++ reduced w.r.t. \spad{ts}.
     initiallyReduced? : (P, %) -> Boolean
         ++ \spad{initiallyReduced?(p, ts)} returns true iff \spad{p} and all its iterated initials
         ++ are reduced w.r.t. to the elements of \spad{ts} with the same main variable.
     autoReduced? : (%, ((P, List(P)) -> Boolean)) -> Boolean
         ++ \spad{autoReduced?(ts, redOp?)} returns true iff every element of \spad{ts} is
         ++ reduced w.r.t to every other in the sense of \spad{redOp?}
     stronglyReduced? : % -> Boolean
         ++ \spad{stronglyReduced?(ts)} returns true iff every element of \spad{ts} is
         ++ reduced w.r.t to any other element of \spad{ts}.
     headReduced? : % -> Boolean
         ++ headReduced?(ts) returns true iff the head of every element of \spad{ts} is
         ++ reduced w.r.t to any other element of \spad{ts}.
     initiallyReduced? : % -> Boolean
         ++ initiallyReduced?(ts) returns true iff for every element \spad{p} of \spad{ts}
         ++ \spad{p} and all its iterated initials are reduced w.r.t. to the other elements
         ++ of \spad{ts} with the same main variable.
     reduce : (P, %, ((P, P) -> P), ((P, P) -> Boolean) ) -> P
         ++ \spad{reduce(p, ts, redOp, redOp?)} returns a polynomial \spad{r}  such that
         ++ \spad{redOp?(r, p)} holds for every \spad{p} of \spad{ts}
         ++ and there exists some product \spad{h} of the initials of the members
         ++ of \spad{ts} such that \spad{h*p - r} lies in the ideal generated by \spad{ts}.
         ++ The operation \spad{redOp} must satisfy the following conditions.
         ++ For every \spad{p} and \spad{q} we have  \spad{redOp?(redOp(p, q), q)}
         ++ and there exists an integer \spad{e} and a polynomial \spad{f} such that
         ++ \spad{init(q)^e*p = f*q + redOp(p, q)}.
     rewriteSetWithReduction : (List P, %, ((P, P) -> P), ((P, P) -> Boolean) ) -> List P
         ++ \spad{rewriteSetWithReduction(lp, ts, redOp, redOp?)} returns a list \spad{lq} of
         ++ polynomials such that \spad{[reduce(p, ts, redOp, redOp?) for p in lp]} and \spad{lp}
         ++ have the same zeros inside the regular zero set of \spad{ts}. Moreover, for every
         ++ polynomial \spad{q} in \spad{lq} and every polynomial \spad{t} in \spad{ts}
         ++ \spad{redOp?(q, t)} holds and there exists a polynomial \spad{p}
         ++ in the ideal generated by \spad{lp} and a product \spad{h} of \spad{initials(ts)}
         ++ such that \spad{h*p - r} lies in the ideal generated by \spad{ts}.
         ++ The operation \spad{redOp} must satisfy the following conditions.
         ++ For every \spad{p} and \spad{q} we have \spad{redOp?(redOp(p, q), q)}
         ++ and there exists an integer \spad{e} and a polynomial \spad{f}
         ++ such that \spad{init(q)^e*p = f*q + redOp(p, q)}.
     stronglyReduce : (P, %) -> P
         ++ \spad{stronglyReduce(p, ts)} returns a polynomial \spad{r}  such that
         ++ \spad{stronglyReduced?(r, ts)} holds and there exists some product
         ++ \spad{h} of \spad{initials(ts)}
         ++ such that \spad{h*p - r} lies in the ideal generated by \spad{ts}.
     headReduce : (P, %) -> P
         ++ \spad{headReduce(p, ts)} returns a polynomial \spad{r}  such that \spad{headReduce?(r, ts)}
         ++ holds and there exists some product \spad{h} of \spad{initials(ts)}
         ++ such that \spad{h*p - r} lies in the ideal generated by \spad{ts}.
     initiallyReduce : (P, %) -> P
         ++ \spad{initiallyReduce(p, ts)} returns a polynomial \spad{r}
         ++ such that  \spad{initiallyReduced?(r, ts)}
         ++ holds and there exists some product \spad{h} of \spad{initials(ts)}
         ++ such that \spad{h*p - r} lies in the ideal generated by \spad{ts}.
     removeZero : (P, %) -> P
         ++ \spad{removeZero(p, ts)} returns \spad{0} if \spad{p} reduces
         ++ to \spad{0} by pseudo-division w.r.t \spad{ts} otherwise
         ++ returns a polynomial \spad{q} computed from \spad{p}
         ++ by removing any coefficient in \spad{p} reducing to \spad{0}.
     collectQuasiMonic : % -> %
         ++ \spad{collectQuasiMonic(ts)} returns the subset of \spad{ts}
         ++ consisting of the polynomials with initial in \spad{R}.
     reduceByQuasiMonic : (P, %) -> P
         ++ \spad{reduceByQuasiMonic(p, ts)} returns the same as
         ++ \spad{remainder(p, collectQuasiMonic(ts)).polnum}.
     zeroSetSplit : List P -> List %
         ++ \spad{zeroSetSplit(lp)} returns a list \spad{lts} of triangular sets such that
         ++ the zero set of \spad{lp} is the union of the closures of the regular zero sets
         ++ of the members of \spad{lts}.
     zeroSetSplitIntoTriangularSystems : List P -> List Record(close : %, open : List P)
         ++ \spad{zeroSetSplitIntoTriangularSystems(lp)} returns a list of triangular
         ++ systems \spad{[[ts1, qs1], ..., [tsn, qsn]]} such that the zero set of \spad{lp}
         ++ is the union of the closures of the \spad{W_i} where \spad{W_i} consists
         ++ of the zeros of \spad{ts} which do not cancel any polynomial in \spad{qsi}.
     first : % -> Union(P,"failed")
         ++ \spad{first(ts)} returns the polynomial of \spad{ts} with greatest main variable
         ++ if \spad{ts} is not empty, otherwise returns \spad{"failed"}.
     last : % -> Union(P,"failed")
         ++ \spad{last(ts)} returns the polynomial of \spad{ts} with smallest main variable
         ++ if \spad{ts} is not empty, otherwise returns \spad{"failed"}.
     rest : % -> Union(%,"failed")
         ++ \spad{rest(ts)} returns the polynomials of \spad{ts} with smaller main variable
         ++ than \spad{mvar(ts)} if \spad{ts} is not empty, otherwise returns "failed"
     algebraicVariables : % -> List(V)
         ++ \spad{algebraicVariables(ts)} returns the decreasingly sorted list of the main
         ++ variables of the polynomials of \spad{ts}.
     algebraic? : (V, %) -> Boolean
         ++ \spad{algebraic?(v, ts)} returns true iff \spad{v} is the main variable of some
         ++ polynomial in \spad{ts}.
     select : (%,V) -> Union(P,"failed")
         ++ \spad{select(ts, v)} returns the polynomial of \spad{ts} with \spad{v} as
         ++ main variable, if any.
     extendIfCan : (%,P) -> Union(%,"failed")
         ++ \spad{extendIfCan(ts, p)} returns a triangular set which encodes the simple
         ++ extension by \spad{p} of the extension of the base field defined by \spad{ts},
         ++ according to the properties of triangular sets of the current domain.
         ++ If the required properties do not hold then "failed" is returned.
         ++ This operation encodes in some sense the properties of the
         ++ triangular sets of the current category. Is is used to implement
         ++ the \spad{construct} operation to guarantee that every triangular
         ++ set build from a list of polynomials has the required properties.
     extend : (%, P) -> %
         ++ \spad{extend(ts, p)} returns a triangular set which encodes the simple
         ++ extension by \spad{p} of the extension of the base field defined by \spad{ts},
         ++ according to the properties of triangular sets of the current category
         ++ If the required properties do not hold an error is returned.
     if V has Finite
     then
       coHeight : % -> NonNegativeInteger
           ++ \spad{coHeight(ts)} returns \spad{size()\$V} minus \spad{\#ts}.
  add

     B ==> Boolean
     RBT ==> Record(bas : %, top : List P)

     import from Integer
     import from List(V)
     import from List(P)

     ts : % = us : % ==
       empty?(ts)$% => empty?(us)$%
       empty?(us)$% => false
       first(ts)::P =$P first(us)::P => rest(ts)::% =$% rest(us)::%
       false

     infRittWu?(ts, us) ==
       empty?(us)$% => not empty?(ts)$%
       empty?(ts)$% => false
       p : P := (last(ts))::P
       q : P := (last(us))::P
       infRittWu?(p, q)$P => true
       supRittWu?(p, q)$P => false
       v : V := mvar(p)
       infRittWu?(collectUpper(ts, v), collectUpper(us, v))$%

     reduced?(p, ts, redOp?) ==
       lp : List P := members(ts)
       while (not empty? lp) and (redOp?(p, first(lp))) repeat
         lp := rest lp
       empty? lp

     basicSet(ps, redOp?) ==
       ps := remove(zero?, ps)
       any?(ground?,ps) => "failed"::Union(RBT,"failed")
       ps := sort(infRittWu?, ps)
       p, b : P
       bs := empty()$%
       ts : List P := []
       while not empty? ps repeat
         b := first(ps)
         bs := extend(bs, b)$%
         ps := rest ps
         while (not empty? ps) and (not reduced?((p := first(ps)), bs, redOp?)) repeat
           ts := cons(p, ts)
           ps := rest ps
       ([bs,ts]$RBT)::Union(RBT,"failed")

     basicSet(ps, pred?, redOp?) ==
       ps := remove(zero?, ps)
       any?(ground?,ps) => "failed"::Union(RBT,"failed")
       gps : List P := []
       bps : List P := []
       while not empty? ps repeat
         p := first ps
         ps := rest ps
         if pred?(p)
           then
             gps := cons(p, gps)
           else
             bps := cons(p, bps)
       gps := sort(infRittWu?, gps)
       p, b : P
       bs := empty()$%
       ts : List P := []
       while not empty? gps repeat
         b := first(gps)
         bs := extend(bs, b)$%
         gps := rest gps
         while (not empty? gps) and (not reduced?((p := first(gps)), bs, redOp?)) repeat
           ts := cons(p, ts)
           gps := rest gps
       ts := sort(infRittWu?, concat(ts, bps))
       ([bs,ts]$RBT)::Union(RBT,"failed")

     initials ts ==
       lip : List P := []
       empty? ts => lip
       lp := members(ts)
       while not empty? lp repeat
          p := first(lp)
          if not ground?((ip := init(p)))
            then
              lip := cons(primPartElseUnitCanonical(ip), lip)
          lp := rest lp
       removeDuplicates lip

     degree ts ==
       empty? ts => 0$NonNegativeInteger
       lp := members ts
       d : NonNegativeInteger := mdeg(first lp)
       while not empty? (lp := rest lp) repeat
         d := d * mdeg(first lp)
       d

     quasiComponent ts ==
       [members(ts), initials(ts)]

     normalized?(p, ts) ==
       normalized?(p, members(ts))$P

     stronglyReduced? (p, ts) ==
       reduced?(p, members(ts))$P

     headReduced? (p, ts) ==
       stronglyReduced?(head(p), ts)

     initiallyReduced? (p, ts) ==
       lp : List (P) := members(ts)
       red : Boolean := true
       while (not empty? lp) and (not ground?(p)$P) and red repeat
         while (not empty? lp) and (mvar(first(lp)) > mvar(p)) repeat
           lp := rest lp
         if (not empty? lp)
           then
             if  (mvar(first(lp)) = mvar(p))
               then
                 if reduced?(p, first(lp))
                   then
                     lp := rest lp
                     p := init(p)
                   else
                     red := false
               else
                 p := init(p)
       red

     reduce(p : P, ts, redOp, redOp?) ==
       (empty? ts) or (ground? p) => p
       ts0 := ts
       while (not empty? ts) and (not ground? p) repeat
          reductor := (first ts)::P
          ts := (rest ts)::%
          if not redOp?(p, reductor)
            then
              p := redOp(p, reductor)
              ts := ts0
       p

     rewriteSetWithReduction(lp, ts, redOp, redOp?) ==
       trivialIdeal? ts => lp
       lp := remove(zero?, lp)
       empty? lp => lp
       any?(ground?, lp) => [1$P]
       rs : List P := []
       while not empty? lp repeat
         p := first lp
         lp := rest lp
         p := primPartElseUnitCanonical reduce(p, ts, redOp, redOp?)
         if not zero? p
           then
             if ground? p
               then
                 lp := []
                 rs := [1$P]
               else
                 rs := cons(p, rs)
       removeDuplicates rs

     stronglyReduce(p, ts) ==
       reduce (p, ts, lazyPrem, reduced?)

     headReduce(p, ts) ==
       reduce (p, ts, headReduce, headReduced?)

     initiallyReduce(p, ts) ==
       reduce (p, ts, initiallyReduce, initiallyReduced?)

     removeZero(p, ts) ==
       (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_-)

     reduceByQuasiMonic(p, ts) ==
       (ground? p) or (empty? ts) => p
       remainder(p, collectQuasiMonic(ts)).polnum

     autoReduced?(ts : %, redOp? : ((P, List(P)) -> Boolean)) ==
       empty? ts => true
       lp : List (P) := members(ts)
       p : P := first(lp)
       lp := rest lp
       while (not empty? lp) and redOp?(p, lp) repeat
          p := first lp
          lp := rest lp
       empty? lp

     stronglyReduced? ts ==
       autoReduced? (ts, reduced?)

     normalized? ts ==
       autoReduced? (ts, normalized?)

     headReduced? ts ==
       autoReduced? (ts, headReduced?)

     initiallyReduced?  ts ==
       autoReduced? (ts, initiallyReduced?)

     mvar ts ==
       empty? ts => error"Error from TSETCAT in mvar : #1 is empty"
       mvar((first(ts))::P)$P

     first ts ==
       empty? ts => "failed"::Union(P,"failed")
       lp : List(P) := sort(supRittWu?, members(ts))$(List P)
       first(lp)::Union(P,"failed")

     last ts ==
       empty? ts => "failed"::Union(P,"failed")
       lp : List(P) := sort(infRittWu?, members(ts))$(List P)
       first(lp)::Union(P,"failed")

     rest ts ==
       empty? ts => "failed"::Union(%,"failed")
       lp : List(P) := sort(supRittWu?, members(ts))$(List P)
       construct(rest(lp))::Union(%,"failed")

     coerce (ts : %) : List(P) ==
       sort(supRittWu?, members(ts))$(List P)

     algebraicVariables ts ==
       [mvar(p) for p in members(ts)]

     algebraic? (v, ts) ==
       member?(v, algebraicVariables(ts))

     select(ts : %, v : V) : Union(P, "failed") ==
       lp : List (P) := sort(supRittWu?, members(ts))$(List P)
       while (not empty? lp) and (not (v = mvar(first lp))) repeat
         lp := rest lp
       empty? lp => "failed"::Union(P,"failed")
       (first lp)::Union(P,"failed")

     collectQuasiMonic ts ==
       lp : List(P) := members(ts)
       newlp : List(P) := []
       while (not empty? lp) repeat
         if ground? init(first(lp)) then newlp := cons(first(lp), newlp)
         lp := rest lp
       construct(newlp)

     collectUnder (ts, v) ==
       lp : List (P) := sort(supRittWu?, members(ts))$(List P)
       while (not empty? lp) and (not (v > mvar(first lp))) repeat
         lp := rest lp
       construct(lp)

     collectUpper  (ts, v) ==
       lp1 : List(P) := sort(supRittWu?, members(ts))$(List P)
       lp2 : List(P) := []
       while (not empty? lp1) and  (mvar(first lp1) > v) repeat
         lp2 := cons(first(lp1), lp2)
         lp1 := rest lp1
       construct(reverse lp2)

     construct(lp : List(P)) ==
       rif := retractIfCan(lp)@Union(%,"failed")
       not (rif case %) => error"in construct : LP -> % from TSETCAT : bad arg"
       rif

     retractIfCan(lp : List(P)) ==
       empty? lp => (empty()$%)::Union(%,"failed")
       lp := sort(supRittWu?, lp)
       rif := retractIfCan(rest(lp))@Union(%,"failed")
       not (rif case %) => error"in retractIfCan : LP -> ... from TSETCAT : bad arg"
       extendIfCan(rif@%, first(lp))@Union(%, "failed")

     extend(ts : %, p : P) : % ==
       eif := extendIfCan(ts,p)@Union(%,"failed")
       not (eif case %) =>
           error "in extend : (%, P) -> % from TSETCAT : bad args"
       eif

     if V has Finite
     then

       coHeight ts ==
         n := size()$V
         m := #(members ts)
         subtractIfCan(n, m)$NonNegativeInteger::NonNegativeInteger


)abbrev domain GTSET GeneralTriangularSet
++ Author: Marc Moreno Maza (marc@nag.co.uk)
++ Date Created: 10/06/1995
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Description:
++ A domain constructor of the category \spadtype{TriangularSetCategory}.
++ The only requirement for a list of polynomials to be a member of such
++ a domain is the following: no polynomial is constant and two distinct
++ polynomials have distinct main variables. Such a triangular set may
++ not be auto-reduced or consistent. Triangular sets are stored
++ as sorted lists w.r.t. the main variables of their members but they
++ are displayed in reverse order.\newline
++ References:
++  [1] P. AUBRY, D. LAZARD and M. MORENO MAZA "On the Theories
++      of Triangular Sets" Journal of Symbol. Comp. 28:105-124, 1999.
++ Version: 1

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

  R : IntegralDomain
  E : OrderedAbelianMonoidSup
  V : OrderedSet
  P : RecursivePolynomialCategory(R, E, V)
  N ==> NonNegativeInteger
  Z ==> Integer
  B ==> Boolean
  LP ==> List P
  PtoP ==> P -> P

  Exports ==  TriangularSetCategory(R, E, V, P)

  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

     -- the following assume that rep(ts) is decreasingly sorted
     -- w.r.t. the main variables of the polynomials in rep(ts)
     coerce(ts : %) : OutputForm ==
       lp : List(P) := reverse(rep(ts))
       brace([p::OutputForm for p in lp]$List(OutputForm))$OutputForm
     mvar ts ==
       empty? ts => error"failed in mvar : % -> V from GTSET"
       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)

     -- for another domain of TSETCAT build on this domain GTSET
     -- the following operations must be redefined
     extendIfCan(ts : %, p : P) ==
       ground? p => "failed"::Union(%,"failed")
       empty? ts => (per([unitCanonical(p)]$LP))::Union(%,"failed")
       not (mvar(ts) < mvar(p)) => "failed"::Union(%,"failed")
       (per(cons(p,rep(ts))))::Union(%,"failed")

)abbrev package PSETPK PolynomialSetUtilitiesPackage
++ Author: Marc Moreno Maza (marc@nag.co.uk)
++ Date Created: 12/01/1995
++ SPARC Version
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Examples:
++ References:
++ Description:
++ This package provides modest routines for polynomial system solving.
++ The aim of many of the operations of this package is to remove certain
++ factors in some polynomials in order to avoid unnecessary computations
++ in algorithms involving splitting techniques by partial factorization.
++ Version: 3

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

  R : IntegralDomain
  E : OrderedAbelianMonoidSup
  V : OrderedSet
  P : RecursivePolynomialCategory(R, E, V)
  N ==> NonNegativeInteger
  Z ==> Integer
  B ==> Boolean
  LP ==> List P
  FP ==> Factored P
  T ==> GeneralTriangularSet(R, E, V, P)
  RBT ==> Record(bas : T, top : LP)
  RUL ==> Record(chs:Union(T,"failed"),rfs:LP)
  GPS ==> GeneralPolynomialSet(R, E, V, P)
  pf ==> MultivariateFactorize(V, E, R, P)

  Exports ==  with

     removeRedundantFactors : LP -> LP
        ++ \spad{removeRedundantFactors(lp)} returns \spad{lq} such that if
        ++ \spad{lp = [p1, ..., pn]} and \spad{lq = [q1, ..., qm]}
        ++ then the product \spad{p1*p2*...*pn} vanishes iff the product \spad{q1*q2*...*qm} vanishes,
        ++ and the product of degrees of the \spad{qi} is not greater than
        ++ the one of the \spad{pj}, and no polynomial in \spad{lq}
        ++ divides another polynomial in \spad{lq}. In particular,
        ++ polynomials lying in the base ring \spad{R} are removed.
        ++ Moreover, \spad{lq} is sorted w.r.t \spad{infRittWu?}.
        ++ Furthermore, if R is gcd-domain, the polynomials in \spad{lq} are
        ++ pairwise without common non trivial factor.
     removeRedundantFactors : (P, P) -> LP
        ++ \spad{removeRedundantFactors(p, q)} returns the same as
        ++ \spad{removeRedundantFactors([p, q])}
     removeSquaresIfCan : LP -> LP
        ++ \spad{removeSquaresIfCan(lp)} returns
        ++ \spad{removeDuplicates [squareFreePart(p)$P for p in lp]}
        ++ if \spad{R} is gcd-domain else returns \spad{lp}.
     unprotectedRemoveRedundantFactors : (P, P) -> LP
        ++ \spad{unprotectedRemoveRedundantFactors(p, q)} returns the same as
        ++ \spad{removeRedundantFactors(p, q)} but does assume that neither
        ++ \spad{p} nor \spad{q} lie in the base ring \spad{R} and assumes that
        ++ \spad{infRittWu?(p, q)} holds. Moreover, if \spad{R} is gcd-domain,
        ++ then \spad{p} and \spad{q} are assumed to be square free.
     removeRedundantFactors : (LP, P) -> LP
        ++ \spad{removeRedundantFactors(lp, q)} returns the same as
        ++ \spad{removeRedundantFactors(cons(q, lp))} assuming
        ++ that \spad{removeRedundantFactors(lp)} returns \spad{lp}
        ++ up to replacing some polynomial \spad{pj} in \spad{lp}
        ++ by some some polynomial \spad{qj} associated to \spad{pj}.
     removeRedundantFactors : (LP, LP) -> LP
        ++ \spad{removeRedundantFactors(lp, lq)} returns the same as
        ++ \spad{removeRedundantFactors(concat(lp, lq))} assuming
        ++ that \spad{removeRedundantFactors(lp)} returns \spad{lp}
        ++ up to replacing some polynomial \spad{pj} in \spad{lp}
        ++ by some polynomial \spad{qj} associated to \spad{pj}.
     removeRedundantFactors : (LP, LP, (LP -> LP)) -> LP
        ++ \spad{removeRedundantFactors(lp, lq, remOp)} returns the same as
        ++ \spad{concat(remOp(removeRoughlyRedundantFactorsInPols(lp, lq)), lq)}
        ++ assuming that \spad{remOp(lq)} returns \spad{lq} up to similarity.
     certainlySubVariety? : (LP, LP) -> B
        ++ \spad{certainlySubVariety?(newlp, lp)} returns true iff for every \spad{p}
        ++ in \spad{lp} the remainder of \spad{p} by \spad{newlp} using the division algorithm
        ++ of Groebner techniques is zero.
     possiblyNewVariety? : (LP, List LP) -> B
        ++ \spad{possiblyNewVariety?(newlp, llp)} returns true iff for every \spad{lp}
        ++ in \spad{llp} certainlySubVariety?(newlp, lp) does not hold.
     probablyZeroDim? : LP -> B
        ++ \spad{probablyZeroDim?(lp)} returns true iff the number of polynomials
        ++ in \spad{lp} is not smaller than the number of variables occurring
        ++ in these polynomials.
     selectPolynomials : ((P -> B), LP) -> Record(goodPols : LP, badPols : LP)
        ++ \spad{selectPolynomials(pred?, ps)} returns \spad{gps, bps} where
        ++ \spad{gps} is a list of the polynomial \spad{p} in \spad{ps}
        ++ such that \spad{pred?(p)} holds and \spad{bps} are the other ones.
     selectOrPolynomials : (List (P -> B), LP) -> Record(goodPols : LP, badPols : LP)
        ++ \spad{selectOrPolynomials(lpred?, ps)} returns \spad{gps, bps} where
        ++ \spad{gps} is a list of the polynomial \spad{p} in \spad{ps}
        ++ such that \spad{pred?(p)} holds for some \spad{pred?} in \spad{lpred?}
        ++ and \spad{bps} are the other ones.
     selectAndPolynomials : (List (P -> B), LP) -> Record(goodPols : LP, badPols : LP)
        ++ \spad{selectAndPolynomials(lpred?, ps)} returns \spad{gps, bps} where
        ++ \spad{gps} is a list of the polynomial \spad{p} in \spad{ps}
        ++ such that \spad{pred?(p)} holds for every \spad{pred?} in \spad{lpred?}
        ++ and \spad{bps} are the other ones.
     quasiMonicPolynomials : LP -> Record(goodPols : LP, badPols : LP)
        ++ \spad{quasiMonicPolynomials(lp)} returns \spad{qmps, nqmps} where
        ++ \spad{qmps} is a list of the quasi-monic polynomials in \spad{lp}
        ++ and \spad{nqmps} are the other ones.
     univariate? : P -> B
        ++ \spad{univariate?(p)} returns true iff \spad{p} involves one and
        ++ only one variable.
     univariatePolynomials : LP -> Record(goodPols : LP, badPols : LP)
        ++ \spad{univariatePolynomials(lp)} returns \spad{ups, nups} where
        ++ \spad{ups} is a list of the univariate polynomials,
        ++ and \spad{nups} are the other ones.
     linear? : P -> B
        ++ \spad{linear?(p)} returns true iff \spad{p} does not lie
        ++ in the base ring \spad{R} and has main degree \spad{1}.
     linearPolynomials : LP -> Record(goodPols : LP, badPols : LP)
        ++ \spad{linearPolynomials(lp)} returns \spad{lps, nlps} where
        ++ \spad{lps} is a list of the linear polynomials in lp,
        ++ and  \spad{nlps} are the other ones.
     bivariate? : P -> B
        ++ \spad{bivariate?(p)} returns true iff \spad{p} involves two and
        ++ only two variables.
     bivariatePolynomials : LP -> Record(goodPols : LP, badPols : LP)
        ++ \spad{bivariatePolynomials(lp)} returns \spad{bps, nbps} where
        ++ \spad{bps} is a list of the bivariate polynomials,
        ++ and \spad{nbps} are the other ones.
     removeRoughlyRedundantFactorsInPols : (LP, LP) -> LP
        ++ \spad{removeRoughlyRedundantFactorsInPols(lp, lf)} returns
        ++ \spad{newlp}where \spad{newlp} is obtained from \spad{lp}
        ++ by removing in every polynomial \spad{p} of \spad{lp}
        ++ any occurrence of a polynomial \spad{f} in \spad{lf}.
        ++ This may involve a lot of exact-quotients computations.
     removeRoughlyRedundantFactorsInPols : (LP, LP, B) -> LP
        ++ \spad{removeRoughlyRedundantFactorsInPols(lp, lf, opt)} returns
        ++ the same as \spad{removeRoughlyRedundantFactorsInPols(lp, lf)}
        ++ if \spad{opt} is \spad{false} and if the previous operation
        ++ does not return any non null and constant polynomial,
        ++ else return \spad{[1]}.
     removeRoughlyRedundantFactorsInPol : (P, LP) -> P
        ++ \spad{removeRoughlyRedundantFactorsInPol(p, lf)} returns the same as
        ++ removeRoughlyRedundantFactorsInPols([p], lf, true)
     interReduce : LP -> LP
        ++ \spad{interReduce(lp)} returns \spad{lq} such that \spad{lp}
        ++ and \spad{lq} generate the same ideal and no polynomial
        ++ in \spad{lq} is reducible by the others in the sense
        ++ of Groebner bases. Since no assumptions are required
        ++ the result may depend on the ordering the reductions are
        ++ performed.
     roughBasicSet : LP -> Union(Record(bas:T,top:LP),"failed")
        ++ \spad{roughBasicSet(lp)} returns the smallest (with Ritt-Wu
        ++ ordering) triangular set contained in \spad{lp}.
     crushedSet : LP -> LP
        ++ \spad{crushedSet(lp)} returns \spad{lq} such that \spad{lp} and
        ++ and \spad{lq} generate the same ideal and no rough basic
        ++ sets reduce (in the sense of Groebner bases) the other
        ++ polynomials in \spad{lq}.
     rewriteSetByReducingWithParticularGenerators : (LP, (P->B), ((P, P)->B), ((P, P)->P)) -> LP
        ++ \spad{rewriteSetByReducingWithParticularGenerators(lp, pred?, redOp?, redOp)}
        ++ returns \spad{lq} where \spad{lq} is computed by the following
        ++ algorithm. Chose a basic set w.r.t. the reduction-test \spad{redOp?}
        ++ among the polynomials satisfying property \spad{pred?},
        ++ if it is empty then leave, else reduce the other polynomials by
        ++ this basic set w.r.t. the reduction-operation \spad{redOp}.
        ++ Repeat while another basic set with smaller rank can be computed.
        ++ See code. If \spad{pred?} is \spad{quasiMonic?} the ideal is unchanged.
     rewriteIdealWithQuasiMonicGenerators : (LP, ((P, P)->B), ((P, P)->P)) -> LP
        ++ \spad{rewriteIdealWithQuasiMonicGenerators(lp, redOp?, redOp)} returns
        ++ \spad{lq} where \spad{lq} and \spad{lp} generate
        ++ the same ideal in \spad{R^(-1) P} and \spad{lq}
        ++ has rank not higher than the one of \spad{lp}.
        ++ Moreover, \spad{lq} is computed by reducing \spad{lp}
        ++ w.r.t. some basic set of the ideal generated by
        ++ the quasi-monic polynomials in \spad{lp}.
     if R has GcdDomain
     then
       squareFreeFactors : P -> LP
          ++ \spad{squareFreeFactors(p)} returns the square-free factors of \spad{p}
          ++ over \spad{R}
       univariatePolynomialsGcds : LP -> LP
          ++ \spad{univariatePolynomialsGcds(lp)} returns \spad{lg} where
          ++ \spad{lg} is a list of the gcds of every pair in \spad{lp}
          ++ of univariate polynomials in the same main variable.
       univariatePolynomialsGcds : (LP, B) -> LP
          ++ \spad{univariatePolynomialsGcds(lp, opt)} returns the same as
          ++ \spad{univariatePolynomialsGcds(lp)} if \spad{opt} is
          ++ \spad{false} and if the previous operation does not return
          ++ any non null and constant polynomial, else return \spad{[1]}.
       removeRoughlyRedundantFactorsInContents : (LP, LP) -> LP
          ++ \spad{removeRoughlyRedundantFactorsInContents(lp, lf)} returns
          ++ \spad{newlp}where \spad{newlp} is obtained from \spad{lp}
          ++ by removing in the content of every polynomial of \spad{lp}
          ++ any occurence of a polynomial \spad{f} in \spad{lf}. Moreover,
          ++ squares over \spad{R} are first removed in the content
          ++ of every polynomial of \spad{lp}.
       removeRedundantFactorsInContents : (LP, LP) -> LP
          ++ \spad{removeRedundantFactorsInContents(lp, lf)} returns \spad{newlp}
          ++ where \spad{newlp} is obtained from \spad{lp} by removing
          ++ in the content of every polynomial of \spad{lp} any non trivial
          ++ factor of any polynomial \spad{f} in \spad{lf}. Moreover,
          ++ squares over \spad{R} are first removed in the content
          ++ of every polynomial of \spad{lp}.
       removeRedundantFactorsInPols : (LP, LP) -> LP
          ++ \spad{removeRedundantFactorsInPols(lp, lf)} returns \spad{newlp}
          ++ where \spad{newlp} is obtained from \spad{lp} by removing
          ++ in every polynomial \spad{p} of \spad{lp} any non trivial
          ++ factor of any polynomial \spad{f} in \spad{lf}. Moreover,
          ++ squares over \spad{R} are first removed in every
          ++ polynomial \spad{lp}.
     if (R has PolynomialFactorizationExplicit) and (R has CharacteristicZero)
     then
       irreducibleFactors : LP -> LP
          ++ \spad{irreducibleFactors(lp)} returns \spad{lf} such that if
          ++ \spad{lp = [p1, ..., pn]} and \spad{lf = [f1, ..., fm]} then
          ++ \spad{p1*p2*...*pn=0} means \spad{f1*f2*...*fm=0}, and the \spad{fi}
          ++ are irreducible over \spad{R} and are pairwise distinct.
       lazyIrreducibleFactors : LP -> LP
          ++ \spad{lazyIrreducibleFactors(lp)} returns \spad{lf} such that if
          ++ \spad{lp = [p1, ..., pn]} and \spad{lf = [f1, ..., fm]} then
          ++ \spad{p1*p2*...*pn=0} means \spad{f1*f2*...*fm=0}, and the \spad{fi}
          ++ are irreducible over \spad{R} and are pairwise distinct.
          ++ The algorithm tries to avoid factorization into irreducible
          ++ factors as far as possible and makes previously use of gcd
          ++ techniques over \spad{R}.
       removeIrreducibleRedundantFactors : (LP, LP) -> LP
          ++ \spad{removeIrreducibleRedundantFactors(lp, lq)} returns the same
          ++ as \spad{irreducibleFactors(concat(lp, lq))} assuming
          ++ that \spad{irreducibleFactors(lp)} returns \spad{lp}
          ++ up to replacing some polynomial \spad{pj} in \spad{lp}
          ++ by some polynomial \spad{qj} associated to \spad{pj}.

  Implementation ==  add

     autoRemainder : T -> List(P)

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

     selectPolynomials  (pred?, ps) ==
       gps : LP := []
       bps : LP := []
       while not empty? ps repeat
         p := first ps
         ps := rest ps
         if pred?(p)
           then
             gps := cons(p, gps)
           else
             bps := cons(p, bps)
       gps := sort(infRittWu?, gps)
       bps := sort(infRittWu?, bps)
       [gps, bps]

     selectOrPolynomials (lpred?, ps) ==
       gps : LP := []
       bps : LP := []
       while not empty? ps repeat
         p := first ps
         ps := rest ps
         clpred? :=  lpred?
         while (not empty? clpred?) and (not (first clpred?)(p)) repeat
           clpred? :=  rest clpred?
         if not empty?(clpred?)
           then
             gps := cons(p, gps)
           else
             bps := cons(p, bps)
       gps := sort(infRittWu?, gps)
       bps := sort(infRittWu?, bps)
       [gps, bps]

     selectAndPolynomials (lpred?, ps) ==
       gps : LP := []
       bps : LP := []
       while not empty? ps repeat
         p := first ps
         ps := rest ps
         clpred? :=  lpred?
         while (not empty? clpred?) and ((first clpred?)(p)) repeat
           clpred? :=  rest clpred?
         if empty?(clpred?)
           then
             gps := cons(p, gps)
           else
             bps := cons(p, bps)
       gps := sort(infRittWu?, gps)
       bps := sort(infRittWu?, bps)
       [gps, bps]

     linear? p ==
       ground? p => false
       (mdeg(p) = 1)

     linearPolynomials  ps ==
       selectPolynomials(linear?, ps)

     univariate? p ==
       ground? p => false
       not(ground?(init(p))) => false
       tp := tail(p)
       ground?(tp) => true
       not (mvar(p) = mvar(tp)) => false
       univariate?(tp)

     univariatePolynomials ps ==
       selectPolynomials(univariate?, ps)

     bivariate? p ==
       ground? p => false
       ground? tail(p) => univariate?(init(p))
       vp := mvar(p)
       vtp := mvar(tail(p))
       ((ground? init(p)) and (vp = vtp)) => bivariate? tail(p)
       ((ground? init(p)) and (vp > vtp)) => univariate? tail(p)
       not univariate?(init(p)) => false
       vip := mvar(init(p))
       vip > vtp => false
       vip = vtp => univariate? tail(p)
       vtp < vp => false
       zero? degree(tail(p), vip) => univariate? tail(p)
       bivariate? tail(p)

     bivariatePolynomials ps ==
       selectPolynomials(bivariate?, ps)

     quasiMonicPolynomials ps ==
       selectPolynomials(quasiMonic?, ps)

     removeRoughlyRedundantFactorsInPols (lp, lf, opt) ==
       empty? lp => lp
       newlp : LP := []
       stop : B := false
       lp := remove(zero?, lp)
       lf := sort(infRittWu?, lf)
       test : Union(P,"failed")
       while (not empty? lp) and (not stop) repeat
         p := first lp
         lp := rest lp
         copylf := lf
         while (not empty? copylf) and (not ground? p) and (not (mvar(p) < mvar(first copylf))) repeat
           f := first copylf
           copylf := rest copylf
           while (((test := p exquo$P f)) case P) repeat
             p := test::P
         stop := opt and ground?(p)
         newlp := cons(unitCanonical(p), newlp)
       stop => [1$P]
       newlp

     removeRoughlyRedundantFactorsInPol(p, lf) ==
       zero? p => p
       lp : LP := [p]
       first removeRoughlyRedundantFactorsInPols (lp, lf, true()$B)

     removeRoughlyRedundantFactorsInPols (lp, lf) ==
       removeRoughlyRedundantFactorsInPols (lp, lf, false()$B)

     possiblyNewVariety?(newlp, llp) ==
       while (not empty? llp) and _
        (not certainlySubVariety?(newlp, first(llp))) repeat
         llp := rest llp
       empty? llp

     certainlySubVariety?(lp, lq) ==
       gs := construct(lp)$GPS
       while (not empty? lq) and _
        (zero? (remainder(first(lq), gs)$GPS).polnum) repeat
         lq := rest lq
       empty? lq

     probablyZeroDim?(lp : List P) : Boolean ==
       m := #lp
       lv : List V := variables(first lp)
       while not empty? (lp := rest lp) repeat
         lv := concat(variables(first lp), lv)
       n := #(removeDuplicates lv)
       not (n > m)

     interReduce(lp : LP) : LP ==
       ps := lp
       rs : List(P) := []
       repeat
         empty? ps => return rs
         ps := sort(supRittWu?, ps)
         p := first ps
         ps := rest ps
         r := remainder(p, [ps]$GPS).polnum
         zero? r => iterate
         ground? r => return []
         associates?(r, p) => rs := cons(r, rs)
         ps := concat(ps, cons(r, rs))
         rs := []

     roughRed?(p : P, q : P) : B ==
       ground? p => false
       ground? q => true
       mvar(p) > mvar(q)

     roughBasicSet(lp) == basicSet(lp, roughRed?)$T

     autoRemainder(ts : T) : List(P) ==
       empty? ts => members(ts)
       lp := sort(infRittWu?, reverse members(ts))
       newlp : List(P) := [primPartElseUnitCanonical first(lp)]
       lp := rest(lp)
       while not empty? lp repeat
         p := (remainder(first(lp), construct(newlp)$GPS)$GPS).polnum
         if not zero? p
           then
             if ground? p
               then
                 newlp := [1$P]
                 lp := []
               else
                 newlp := cons(p, newlp)
                 lp := rest(lp)
           else
             lp := rest(lp)
       newlp

     crushedSet(lp) ==
       rec := roughBasicSet(lp)
       contradiction := (rec case "failed")@B
       finished : B := false
       while (not finished) and (not contradiction) repeat
         bs := (rec::RBT).bas
         rs := (rec::RBT).top
         rs :=  rewriteIdealWithRemainder(rs, bs)$T
         contradiction := ((not empty? rs) and (first(rs) = 1))
         if not contradiction
           then
             rs := concat(rs, autoRemainder(bs))
             rec := roughBasicSet(rs)
             contradiction := (rec case "failed")@B
             not contradiction => finished := not infRittWu?((rec::RBT).bas, bs)
       contradiction => [1$P]
       rs

     rewriteSetByReducingWithParticularGenerators (ps, pred?, redOp?, redOp) ==
       rs : LP := remove(zero?, ps)
       any?(ground?, rs) => [1$P]
       contradiction : B := false
       bs1 : T := empty()$T
       rec : Union(RBT,"failed")
       ar : Union(T, List(P))
       stop : B := false
       while (not contradiction) and (not stop) repeat
         rec := basicSet(rs, pred?, redOp?)$T
         bs2 : T := (rec::RBT).bas
         rs := (rec::RBT).top
         -- ar := autoReduce(bs2, lazyPrem, reduced?)@Union(T, List(P))
         ar := bs2::Union(T, List(P))
         if (ar case T) then
             bs2 := ar@T
             if infRittWu?(bs2, bs1)
               then
                 rs := rewriteSetWithReduction(rs, bs2, redOp, redOp?)$T
                 bs1 := bs2
               else
                 stop := true
             rs := concat(members(bs2), rs)
         else
             rs := concat(ar::LP, rs)
         if any?(ground?, rs)
           then
             contradiction := true
             rs := [1$P]
       rs

     removeRedundantFactors (lp : LP, lq : LP, remOp : (LP -> LP)) ==
       -- ASSUME remOp(lp) returns lp up to similarity
       lq := removeRoughlyRedundantFactorsInPols(lq, lp, false)
       lq := remOp lq
       sort(infRittWu?, concat(lp, lq))

     removeRedundantFactors (lp : LP, lq : LP) ==
       lq := removeRoughlyRedundantFactorsInPols(lq, lp, false)
       lq := removeRedundantFactors lq
       sort(infRittWu?, concat(lp, lq))

     if (R has PolynomialFactorizationExplicit) and (R has CharacteristicZero)
     then
       irreducibleFactors lp ==
         newlp : LP := []
         fp : FP
         while not empty? lp repeat
           p := first lp
           lp := rest lp
           fp := factor(p)$pf
           lrrz := factorList(fp)$FP
           lf := remove(ground?, [rrz.factor for rrz in lrrz])
           newlp := concat(lf, newlp)
         removeDuplicates newlp

       lazyIrreducibleFactors lp ==
         lp := removeRedundantFactors(lp)
         newlp : LP := []
         fp : FP
         while not empty? lp repeat
           p := first lp
           lp := rest lp
           fp := factor(p)$pf
           lrrz := factorList(fp)$FP
           lf := remove(ground?, [rrz.factor for rrz in lrrz])
           newlp := concat(lf, newlp)
         newlp

       removeIrreducibleRedundantFactors (lp : LP, lq : LP) ==
         -- ASSUME lp only contains irreducible factors over R
         lq := removeRoughlyRedundantFactorsInPols(lq, lp, false)
         lq := irreducibleFactors lq
         sort(infRittWu?, concat(lp, lq))

     if R has GcdDomain
     then

       squareFreeFactors(p : P) ==
         sfp : Factored P := squareFree(p)$P
         lsf : List P := [foo.factor for foo in factorList(sfp)]
         lsf

       univariatePolynomialsGcds (ps, opt) ==
         lg : LP := []
         pInV : LP
         stop : B := false
         ps := sort(infRittWu?, ps)
         p, g : P
         v : V
         while (not empty? ps) and (not stop) repeat
           while (not empty? ps) and (not univariate?((p := first(ps)))) repeat
             ps := rest ps
           if not empty? ps
             then
               v := mvar(p)$P
               pInV := [p]
               while (not empty? ps) and (mvar((p := first(ps))) = v) repeat
                 if (univariate?(p))
                   then
                     pInV := cons(p, pInV)
                 ps := rest ps
               g := gcd(pInV)$P
               stop := opt and (ground? g)
               lg := cons(g, lg)
         stop => [1$P]
         lg

       univariatePolynomialsGcds ps ==
         univariatePolynomialsGcds (ps, false)

       removeSquaresIfCan lp ==
         empty? lp => lp
         removeDuplicates [squareFreePart(p)$P for p in lp]

       rewriteIdealWithQuasiMonicGenerators (ps, redOp?, redOp) ==
         ups := removeSquaresIfCan(univariatePolynomialsGcds(ps, true))
         ps := removeDuplicates concat(ups, ps)
         rewriteSetByReducingWithParticularGenerators(ps, quasiMonic?, redOp?, redOp)

       removeRoughlyRedundantFactorsInContents (ps, lf) ==
         empty? ps => ps
         newps : LP := []
         p, newp, cp, newcp, f, g : P
         test : Union(P,"failed")
         copylf : LP
         while not empty? ps repeat
           p := first ps
           ps := rest ps
           cp := mainContent(p)$P
           newcp := squareFreePart(cp)$P
           newp := (p exquo$P cp)::P
           if not ground? newcp
             then
               copylf := [f for f in lf | mvar(f) <= mvar(newcp)]
               while (not empty? copylf) and (not ground? newcp) repeat
                 f := first copylf
                 copylf := rest copylf
                 test := (newcp exquo$P f)
                 if (test case P) then
                     newcp := test@P
           if ground? newcp
             then
               newp := unitCanonical(newp)
             else
               newp := unitCanonical(newp * newcp)
           newps := cons(newp, newps)
         newps

       removeRedundantFactorsInContents (ps, lf) ==
         empty? ps => ps
         newps : LP := []
         p, newp, cp, newcp, f, g : P
         while not empty? ps repeat
           p := first ps
           ps := rest ps
           cp := mainContent(p)$P
           newcp := squareFreePart(cp)$P
           newp := (p exquo$P cp)::P
           if not ground? newcp
             then
               copylf := lf
               while (not empty? copylf) and (not ground? newcp) repeat
                 f := first copylf
                 copylf := rest copylf
                 g := gcd(newcp, f)$P
                 if not ground? g
                   then
                     newcp := (newcp exquo$P g)::P
           if ground? newcp
             then
               newp := unitCanonical(newp)
             else
               newp := unitCanonical(newp * newcp)
           newps := cons(newp, newps)
         newps

       removeRedundantFactorsInPols (ps, lf) ==
         empty? ps => ps
         newps : LP := []
         p, newp, cp, newcp, f, g : P
         while not empty? ps repeat
           p := first ps
           ps := rest ps
           cp := mainContent(p)$P
           newcp := squareFreePart(cp)$P
           newp := (p exquo$P cp)::P
           newp := squareFreePart(newp)$P
           copylf := lf
           while not empty? copylf repeat
             f := first copylf
             copylf := rest copylf
             if not ground? newcp
               then
                 g := gcd(newcp, f)$P
                 if not ground? g
                   then
                     newcp := (newcp exquo$P g)::P
             if not ground? newp
               then
                 g := gcd(newp, f)$P
                 if not ground? g
                   then
                     newp := (newp exquo$P g)::P
           if ground? newcp
             then
               newp := unitCanonical(newp)
             else
               newp := unitCanonical(newp * newcp)
           newps := cons(newp, newps)
         newps

       removeRedundantFactors (a : P, b : P) : LP ==
         a := primPartElseUnitCanonical(squareFreePart(a))
         b := primPartElseUnitCanonical(squareFreePart(b))
         if not infRittWu?(a, b)
           then
            (a, b) := (b, a)
         if ground? a
           then
             if ground? b
               then
                 return([])
               else
                 return([b])
           else
             if ground? b
               then
                 return([a])
               else
                 return(unprotectedRemoveRedundantFactors(a, b))

       unprotectedRemoveRedundantFactors (a, b) ==
         c := b exquo$P a
         if (c case P) then
             d : P := c@P
             if ground? d then
                 return([a])
             else
                 return([a, d])
         else
             g : P := gcd(a, b)$P
             if ground? g then
                 return([a, b])
             else
                 return([g, (a exquo$P g)::P, (b exquo$P g)::P])

     else

       removeSquaresIfCan lp ==
         lp

       rewriteIdealWithQuasiMonicGenerators (ps, redOp?, redOp) ==
         rewriteSetByReducingWithParticularGenerators(ps, quasiMonic?, redOp?, redOp)

       removeRedundantFactors (a : P, b : P) ==
         a := primPartElseUnitCanonical(a)
         b := primPartElseUnitCanonical(b)
         if not infRittWu?(a, b)
           then
            (a, b) := (b, a)
         if ground? a
           then
             if ground? b
               then
                 return([])
               else
                 return([b])
           else
             if ground? b
               then
                 return([a])
               else
                 return(unprotectedRemoveRedundantFactors(a, b))

       unprotectedRemoveRedundantFactors (a, b) ==
         c := b exquo$P a
         if (c case P) then
             d : P := c@P
             if ground? d then
                 return([a])
             else
                 if infRittWu?(d, a) then (a, d) := (d, a)
                 return(unprotectedRemoveRedundantFactors(a, d))
         else
              return([a, b])

     removeRedundantFactors (lp : LP) ==
       lp := remove(ground?, lp)
       lp := removeDuplicates [primPartElseUnitCanonical(p) for p in lp]
       lp := removeSquaresIfCan lp
       lp := removeDuplicates [unitCanonical(p) for p in lp]
       empty? lp => lp
       size?(lp, 1$N)$(List P) => lp
       lp := sort(infRittWu?, lp)
       p : P := first lp
       lp := rest lp
       base : LP := unprotectedRemoveRedundantFactors(p, first lp)
       top : LP := rest lp
       while not empty? top repeat
         p := first top
         base := removeRedundantFactors(base, p)
         top := rest top
       base

     removeRedundantFactors (lp : LP, a : P) ==
       lp := remove(ground?, lp)
       lp := sort(infRittWu?, lp)
       ground? a => lp
       empty? lp => [a]
       toSee : LP := lp
       toSave : LP := []
       while not empty? toSee repeat
         b := first toSee
         toSee := rest toSee
         if not infRittWu?(b, a)
           then
             (c, d) := (a, b)
           else
             (c, d) := (b, a)
         rrf := unprotectedRemoveRedundantFactors(c, d)
         empty? rrf => error"in removeRedundantFactors : (LP,P) -> LP from PSETPK"
         c := first rrf
         rrf := rest rrf
         if empty? rrf
           then
             if associates?(c, b)
               then
                 toSave := concat(toSave, toSee)
                 a := b
                 toSee := []
               else
                 a := c
                 toSee := concat(toSave, toSee)
                 toSave := []
           else
             d := first rrf
             rrf := rest rrf
             if empty? rrf
               then
                 if associates?(c, b)
                   then
                     toSave := concat(toSave, [b])
                     a := d
                   else
                     if associates?(d, b)
                       then
                         toSave := concat(toSave, [b])
                         a := c
                       else
                         toSave := removeRedundantFactors(toSave, c)
                         a := d
               else
                 e := first rrf
                 not empty? rest(rrf) => error"in removeRedundantFactors:(LP,P)->LP from PSETPK"
                 -- ASSUME that neither c, nor d, nor e may be associated to b
                 toSave := removeRedundantFactors(toSave, c)
                 toSave := removeRedundantFactors(toSave, d)
                 a := e
         if empty? toSee
           then
             toSave := sort(infRittWu?, cons(a, toSave))
       toSave

)abbrev domain WUTSET WuWenTsunTriangularSet
++ Author: Marc Moreno Maza (marc@nag.co.uk)
++ Date Created: 11/18/1995
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Description: A domain constructor of the category \spadtype{GeneralTriangularSet}.
++ The only requirement for a list of polynomials to be a member of such
++ a domain is the following: no polynomial is constant and two distinct
++ polynomials have distinct main variables. Such a triangular set may
++ not be auto-reduced or consistent. The \spadopFrom{construct}{WuWenTsunTriangularSet} operation
++ does not check the previous requirement. Triangular sets are stored
++ as sorted lists w.r.t. the main variables of their members.
++ Furthermore, this domain exports operations dealing with the
++ characteristic set method of Wu Wen Tsun and some optimizations
++ mainly proposed by Dong Ming Wang.\newline
++ References:
++  [1] W. T. WU "A Zero Structure Theorem for polynomial equations solving"
++       MM Research Preprints, 1987.
++  [2] D. M. WANG "An implementation of the characteristic set method in Maple"
++       Proc. DISCO'92. Bath, England.
++ Version: 3

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

  R : IntegralDomain
  E : OrderedAbelianMonoidSup
  V : OrderedSet
  P : RecursivePolynomialCategory(R, E, V)
  N ==> NonNegativeInteger
  Z ==> Integer
  B ==> Boolean
  LP ==> List P
  A ==> FiniteEdge P
  H ==> FiniteSimpleHypergraph P
  RBT ==> Record(bas : %, top : LP)
  pa ==> PolynomialSetUtilitiesPackage(R, E, V, P)
  NLpT ==> SplittingNode(LP, %)
  ALpT ==> SplittingTree(LP, %)
  O ==> OutputForm
  OP ==> OutputPackage

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

     medialSet : (LP,((P,P)->B),((P,P)->P)) -> Union(%,"failed")
        ++ \spad{medialSet(ps, redOp?, redOp)} returns \spad{bs} a basic set
        ++ (in Wu Wen Tsun sense w.r.t the reduction-test \spad{redOp?})
        ++ of some set generating the same ideal as \spad{ps} (with
        ++ rank not higher than  any basic set of \spad{ps}), if no non-zero
        ++ constant polynomials appear during the computations, else
        ++ \spad{"failed"} is returned. In the former case, \spad{bs} has to be
        ++ understood as a candidate for being a characteristic set of \spad{ps}.
        ++ In the original algorithm, \spad{bs} is simply a basic set of \spad{ps}.
     medialSet : LP -> Union(%,"failed")
        ++ \spad{medial(ps)} returns the same as
        ++ \spad{medialSet(ps, initiallyReduced?, initiallyReduce)}.
     characteristicSet : (LP,((P,P)->B),((P,P)->P)) -> Union(%,"failed")
        ++ \spad{characteristicSet(ps, redOp?, redOp)} returns a non-contradictory
        ++ characteristic set of \spad{ps} in Wu Wen Tsun sense w.r.t the
        ++ reduction-test \spad{redOp?} (using \spad{redOp} to reduce
        ++ polynomials w.r.t a \spad{redOp?} basic set), if no
        ++ non-zero constant polynomial appear during those reductions,
        ++ else \spad{"failed"} is returned.
        ++ The operations \spad{redOp} and \spad{redOp?} must satisfy
        ++ the following conditions: \spad{redOp?(redOp(p, q), q)} holds
        ++ for every polynomials \spad{p, q} and there exists an integer
        ++ \spad{e} and a polynomial \spad{f} such that we have
        ++ \spad{init(q)^e*p = f*q + redOp(p, q)}.
     characteristicSet : LP -> Union(%,"failed")
        ++ \spad{characteristicSet(ps)} returns the same as
        ++ \spad{characteristicSet(ps, initiallyReduced?, initiallyReduce)}.
     characteristicSerie  : (LP, ((P, P)->B), ((P, P)->P)) -> List %
        ++ \spad{characteristicSerie(ps, redOp?, redOp)} returns a list \spad{lts}
        ++ of triangular sets such that the zero set of \spad{ps} is the
        ++ union of the regular zero sets of the members of \spad{lts}.
        ++ This is made by the Ritt and Wu Wen Tsun process applying
        ++ the operation \spad{characteristicSet(ps, redOp?, redOp)}
        ++ to compute characteristic sets in Wu Wen Tsun sense.
     characteristicSerie : LP ->  List %
        ++ \spad{characteristicSerie(ps)} returns the same as
        ++ \spad{characteristicSerie(ps, initiallyReduced?, initiallyReduce)}.

  Implementation == GeneralTriangularSet(R, E, V, P) add

     removeSquares : % -> Union(%,"failed")

     Rep ==> LP

     rep(s : %) : Rep == s pretend Rep
     per(l : Rep) : % == l pretend %

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

     medialSetWithTrace (ps:LP,redOp?:((P,P)->B),redOp:((P,P)->P)):Union(RBT,"failed") ==
       qs := rewriteIdealWithQuasiMonicGenerators(ps, redOp?, redOp)$pa
       contradiction : B := any?(ground?, ps)
       contradiction => "failed"::Union(RBT,"failed")
       rs : LP := qs
       bs : %
       while (not empty? rs) and (not contradiction) repeat
         rec := basicSet(rs, redOp?)
         contradiction := (rec case "failed")
         if not contradiction then
             bs := (rec::RBT).bas
             rs := (rec::RBT).top
             rs :=  rewriteIdealWithRemainder(rs, bs)
             contradiction := ((not empty? rs) and (first(rs) = 1))
             if (not empty? rs) and (not contradiction)
               then
                 rs := rewriteSetWithReduction(rs, bs, redOp, redOp?)
                 contradiction := ((not empty? rs) and (first(rs) = 1))
         if (not empty? rs) and (not contradiction) then
             rs := removeDuplicates concat(rs, members(bs))
             rs := rewriteIdealWithQuasiMonicGenerators(rs, redOp?, redOp)$pa
             contradiction := ((not empty? rs) and (first(rs) = 1))
       contradiction => "failed"::Union(RBT,"failed")
       ([bs,qs]$RBT)::Union(RBT,"failed")

     medialSet(ps : LP, redOp? : ((P, P)->B), redOp : ((P, P)->P)) ==
       foo : Union(RBT,"failed") := medialSetWithTrace(ps,redOp?,redOp)
       (foo case "failed") => "failed"
       (foo@RBT).bas

     medialSet(ps : LP) == medialSet(ps, initiallyReduced?, initiallyReduce)

     characteristicSetUsingTrace(ps:LP,redOp?:((P,P)->B),redOp:((P,P)->P)):Union(%,"failed") ==
       ps := removeAssociates ps
       ps := remove(zero?, ps)
       contradiction : B := any?(ground?, ps)
       contradiction => "failed"::Union(%,"failed")
       rs : LP := ps
       qs : LP := ps
       ms : %
       while (not empty? rs) and (not contradiction) repeat
         rec := medialSetWithTrace (qs, redOp?, redOp)
         contradiction := (rec case "failed")
         if not contradiction then
             ms := (rec::RBT).bas
             qs := (rec::RBT).top
             qs := rewriteIdealWithRemainder(qs, ms)
             contradiction := ((not empty? qs) and (first(qs) = 1))
             if not contradiction then
                 rs :=  rewriteSetWithReduction(qs, ms, lazyPrem, reduced?)
                 contradiction := ((not empty? rs) and (first(rs) = 1))
             if  (not contradiction) and (not empty? rs) then
                 qs := removeDuplicates(concat(rs, concat(members(ms), qs)))
       contradiction => "failed"
       ms

     characteristicSet(ps : LP, redOp? : ((P, P)->B), redOp : ((P, P)->P)) ==
       characteristicSetUsingTrace(ps, redOp?, redOp)

     characteristicSet(ps : LP) == characteristicSet(ps, initiallyReduced?, initiallyReduce)

     characteristicSerie(ps : LP, redOp? : ((P, P)->B), redOp : ((P, P)->P)) ==
       a := [[ps, empty()$%]$NLpT]$ALpT
       while ((esl := extractSplittingLeaf(a)) case ALpT) repeat
          ps := value(value(esl::ALpT)$ALpT)$NLpT
          charSet? := characteristicSetUsingTrace(ps, redOp?, redOp)
          if not (charSet? case %) then
                setvalue!(esl::ALpT, [[]$LP, empty()$%, true]$NLpT)
                updateStatus!(a)
          else
                cs := (charSet?)::%
                lics := initials(cs)
                lics := removeRedundantFactors(lics)$pa
                lics := sort(infRittWu?, lics)
                if empty? lics then
                      setvalue!(esl::ALpT, [ps, cs, true]$NLpT)
                      updateStatus!(a)
                else
                      ln : List NLpT := [[[]$LP, cs, true]$NLpT]
                      while not empty? lics repeat
                         newps := cons(first(lics), concat(cs::LP, ps))
                         lics := rest lics
                         newps := removeDuplicates newps
                         newps := sort(infRittWu?, newps)
                         ln := cons([newps, empty()$%, false]$NLpT, ln)
                      splitNodeOf!(esl::ALpT, a, ln)
       remove(empty()$%, conditions(a))

     characteristicSerie(ps : LP) ==  characteristicSerie (ps, initiallyReduced?, initiallyReduce)

     if R has GcdDomain
     then

       removeSquares (ts:%):Union(%,"failed") ==
         empty?(ts)$% => ts::Union(%,"failed")
         p := (first ts)::P
         rsts : Union(%,"failed")
         rsts := removeSquares((rest ts)::%)
         not(rsts case %) => "failed"
         newts := rsts@%
         empty? newts =>
           p := squareFreePart(p)
           (per([primitivePart(p)]$LP))::Union(%,"failed")
         zero? initiallyReduce(init(p),newts) => "failed"::Union(%,"failed")
         p := primitivePart(removeZero(p, newts))
         ground? p => "failed"::Union(%,"failed")
         not (mvar(newts) < mvar(p)) => "failed"::Union(%,"failed")
         p := squareFreePart(p)
         (per(cons(unitCanonical(p),rep(newts))))::Union(%,"failed")

       zeroSetSplit lp ==
         lts : List % := characteristicSerie(lp, initiallyReduced?, initiallyReduce)
         lts := removeDuplicates(lts)$(List %)
         newlts : List % := []
         while not empty? lts repeat
           ts := first lts
           lts := rest lts
           iic := removeSquares(ts)
           if iic case % then
               newlts := cons(iic@%, newlts)
         newlts := removeDuplicates(newlts)$(List %)
         sort(infRittWu?, newlts)

     else

       zeroSetSplit lp ==
         lts : List % := characteristicSerie(lp, initiallyReduced?, initiallyReduce)
         sort(infRittWu?, removeDuplicates lts)

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