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

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra triset.spad}
\author{Marc Moreno Maza}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject

\section{category TSETCAT TriangularSetCategory}

<<category TSETCAT TriangularSetCategory>>=
import Boolean
import NonNegativeInteger
import OutputForm
import List
)abbrev category TSETCAT TriangularSetCategory
++ Author: Marc Moreno Maza (marc@nag.co.uk)
++ Date Created: 04/26/1994
++ Date Last Updated: 12/15/1998
++ 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 \axiom{R} be an integral domain and \axiom{V} a finite ordered set of 
++ variables, say \axiom{X1 < X2 < ... < Xn}.  
++ A set \axiom{S} of polynomials in \axiom{R[X1,X2,...,Xn]} is triangular
++ if no elements of \axiom{S} lies in \axiom{R}, and if two distinct 
++ elements of \axiom{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 \axiom{P} is reduced w.r.t a non-constant polynomial 
++ \axiom{Q} if the degree of \axiom{P} in the main variable of \axiom{Q} 
++ is less than the main degree of \axiom{Q}.
++ A polynomial \axiom{P} is reduced w.r.t a triangular set \axiom{T}
++ if it is reduced w.r.t. every polynomial of \axiom{T}. \newline
++ References :
++  [1] P. AUBRY, D. LAZARD and M. MORENO MAZA "On the Theories
++      of Triangular Sets" Journal of Symbol. Comp. (to appear)
++ Version: 4.

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

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

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

@

\section{domain GTSET GeneralTriangularSet}

<<domain GTSET GeneralTriangularSet>>=
)abbrev domain GTSET GeneralTriangularSet
++ Author: Marc Moreno Maza (marc@nag.co.uk)
++ Date Created: 10/06/1995
++ Date Last Updated: 06/12/1996
++ Basic Functions:
++ Related Constructors:
++ Also See: 
++ AMS Classifications:
++ Keywords:
++ Description: 
++ A domain constructor of the category \axiomType{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. (to appear)
++ 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

     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 variavles 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")

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

     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 => "leave"
         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       
       rs: LP
       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 (one? first(rs)))
         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)@B
           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 EuclideanDomain) and (R has CharacteristicZero)
     then
       irreducibleFactors lp ==
         newlp : LP := []
         lrrz : List RRZ
         rrz : RRZ
         fp : FP
         while not empty? lp repeat
           p := first lp
           lp := rest lp
           fp := factor(p)$pf
           lrrz := factors(fp)$FP
           lf := remove(ground?,[rrz.factor for rrz in lrrz])
           newlp := concat(lf,newlp)
         removeDuplicates newlp

       lazyIrreducibleFactors lp ==
         lp := removeRedundantFactors(lp)
         newlp : LP := []
         lrrz : List RRZ
         rrz : RRZ
         fp : FP
         while not empty? lp repeat
           p := first lp
           lp := rest lp
           fp := factor(p)$pf
           lrrz := factors(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 factors(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)@B
                   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)@B
           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)@B
           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
       #lp = 1$N => 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   

@
\section{domain WUTSET WuWenTsunTriangularSet}
<<domain WUTSET WuWenTsunTriangularSet>>=
)abbrev domain WUTSET WuWenTsunTriangularSet
++ Author: Marc Moreno Maza (marc@nag.co.uk)
++ Date Created: 11/18/1995
++ Date Last Updated: 12/15/1998
++ Basic Functions:
++ Related Constructors:
++ Also See: 
++ AMS Classifications:
++ Keywords:
++ Description: A domain constructor of the category \axiomType{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 \axiomOpFrom{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
  GPS ==> GeneralPolynomialSet(R,E,V,P)
  RBT ==> Record(bas:$,top:LP)
  RUL ==> Record(chs:Union($,"failed"),rfs: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")
        ++ \axiom{medialSet(ps,redOp?,redOp)} returns \axiom{bs} a basic set
        ++ (in Wu Wen Tsun sense w.r.t the reduction-test \axiom{redOp?})
        ++ of some set generating the same ideal as \axiom{ps} (with 
        ++ rank not higher than  any basic set of \axiom{ps}), if no non-zero 
        ++ constant polynomials appear during the computatioms, else 
        ++ \axiom{"failed"} is returned. In the former case, \axiom{bs} has to be 
        ++ understood as a candidate for being a characteristic set of \axiom{ps}.
        ++ In the original algorithm, \axiom{bs} is simply a basic set of \axiom{ps}.
     medialSet: LP -> Union($,"failed")
        ++ \axiom{medial(ps)} returns the same as 
        ++ \axiom{medialSet(ps,initiallyReduced?,initiallyReduce)}.
     characteristicSet : (LP,((P,P)->B),((P,P)->P)) -> Union($,"failed")
        ++ \axiom{characteristicSet(ps,redOp?,redOp)} returns a non-contradictory
        ++ characteristic set of \axiom{ps} in Wu Wen Tsun sense w.r.t the 
        ++ reduction-test \axiom{redOp?} (using \axiom{redOp} to reduce 
        ++ polynomials w.r.t a \axiom{redOp?} basic set), if no
        ++ non-zero constant polynomial appear during those reductions,
        ++ else \axiom{"failed"} is returned.
        ++ The operations \axiom{redOp} and \axiom{redOp?} must satisfy 
        ++ the following conditions: \axiom{redOp?(redOp(p,q),q)} holds
        ++ for every polynomials \axiom{p,q} and there exists an integer
        ++ \axiom{e} and a polynomial \axiom{f} such that we have
        ++ \axiom{init(q)^e*p = f*q + redOp(p,q)}.
     characteristicSet: LP -> Union($,"failed")
        ++ \axiom{characteristicSet(ps)} returns the same as 
        ++ \axiom{characteristicSet(ps,initiallyReduced?,initiallyReduce)}.
     characteristicSerie  : (LP,((P,P)->B),((P,P)->P)) -> List $
        ++ \axiom{characteristicSerie(ps,redOp?,redOp)} returns a list \axiom{lts}
        ++ of triangular sets such that the zero set of \axiom{ps} is the
        ++ union of the regular zero sets of the members of \axiom{lts}.
        ++ This is made by the Ritt and Wu Wen Tsun process applying
        ++ the operation \axiom{characteristicSet(ps,redOp?,redOp)}
        ++ to compute characteristic sets in Wu Wen Tsun sense.
     characteristicSerie: LP ->  List $
        ++ \axiom{characteristicSerie(ps)} returns the same as 
        ++ \axiom{characteristicSerie(ps,initiallyReduced?,initiallyReduce)}.

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

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

     Rep == LP

     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")@B
         if not contradiction
           then
             bs := (rec::RBT).bas
             rs := (rec::RBT).top
             rs :=  rewriteIdealWithRemainder(rs,bs)
             contradiction := ((not empty? rs) and (one? first(rs)))
             if (not empty? rs) and (not contradiction)
               then
                 rs := rewriteSetWithReduction(rs,bs,redOp,redOp?)
                 contradiction := ((not empty? rs) and (one? first(rs)))
         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 (one? first(rs)))
       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" :: Union($,"failed")
       ((foo::RBT).bas) :: Union($,"failed")

     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")@B
         if not contradiction
           then
             ms := (rec::RBT).bas
             qs := (rec::RBT).top
             qs := rewriteIdealWithRemainder(qs,ms)
             contradiction := ((not empty? qs) and (one? first(qs))) 
             if not contradiction
               then
                 rs :=  rewriteSetWithReduction(qs,ms,lazyPrem,reduced?)
                 contradiction := ((not empty? rs) and (one? first(rs)))
             if  (not contradiction) and (not empty? rs)
               then
                 qs := removeDuplicates(concat(rs,concat(members(ms),qs)))
       contradiction => "failed"::Union($,"failed")
       ms::Union($,"failed")

     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,[nil()$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 := [[nil()$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"::Union($,"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)

@
\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 TSETCAT TriangularSetCategory>>
<<domain GTSET GeneralTriangularSet>>
<<package PSETPK PolynomialSetUtilitiesPackage>>
<<domain WUTSET WuWenTsunTriangularSet>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}