axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / PSETPK.spad

)abbrev package PSETPK PolynomialSetUtilitiesPackage
++ Author: Marc Moreno Maza (marc@nag.co.uk)
++ Date Created: 12/01/1995
++ Date Last Updated: 12/15/1998
++ 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.

PolynomialSetUtilitiesPackage (R,E,V,P) : SIG == CODE 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)

  SIG ==> 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}.
       
  CODE ==> 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 => "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       
       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)@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
       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