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

)abbrev category PSETCAT PolynomialSetCategory
++ Author: Marc Moreno Maza
++ 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
++ References:
++ Description: A category for finite subsets of a polynomial ring.
++ Such a set is only regarded as a set of polynomials and not 
++ identified to the ideal it generates. So two distinct sets may 
++ generate the same the ideal. Furthermore, for \spad{R} being an 
++ integral domain, a set of polynomials may be viewed as a representation
++ of the ideal it generates in the polynomial ring \spad{(R)^(-1) P}, 
++ or the set of its zeros (described for instance by the radical of the 
++ previous ideal, or a split of the associated affine variety) and so on. 
++ So this category provides operations about those different notions.
++ Version: 2 

PolynomialSetCategory(R:Ring, E:OrderedAbelianMonoidSup,_
  VarSet:OrderedSet, P:RecursivePolynomialCategory(R,E,VarSet)): Category == 
    Join(SetCategory,Collection(P),CoercibleTo(List(P))) with
     finiteAggregate
     retractIfCan : List(P) -> Union($,"failed")
         ++ \axiom{retractIfCan(lp)} returns an element of the domain whose elements
         ++ are the members of \axiom{lp} if such an element exists, otherwise
         ++ \axiom{"failed"} is returned.
     retract : List(P) -> $
         ++ \axiom{retract(lp)} returns an element of the domain whose elements
         ++ are the members of \axiom{lp} if such an element exists, otherwise
         ++ an error is produced.
     mvar : $ -> VarSet
         ++  \axiom{mvar(ps)} returns the main variable of the non constant polynomial
         ++  with the greatest main variable, if any, else an error is returned.
     variables : $ -> List VarSet
         ++  \axiom{variables(ps)} returns the decreasingly sorted list of the
         ++  variables which are variables of some polynomial in \axiom{ps}.
     mainVariables : $  -> List VarSet
         ++ \axiom{mainVariables(ps)} returns the decreasingly sorted list of the
         ++ variables which are main variables of some polynomial in \axiom{ps}.
     mainVariable? : (VarSet,$) -> Boolean
         ++ \axiom{mainVariable?(v,ps)} returns true iff \axiom{v} is the main variable 
         ++ of some polynomial in \axiom{ps}.
     collectUnder : ($,VarSet) -> $
         ++ \axiom{collectUnder(ps,v)} returns the set consisting of the 
         ++ polynomials of \axiom{ps} with main variable less than \axiom{v}.
     collect : ($,VarSet) -> $
         ++ \axiom{collect(ps,v)}  returns the set consisting of the 
         ++ polynomials of \axiom{ps} with \axiom{v} as main variable.
     collectUpper : ($,VarSet) -> $
         ++ \axiom{collectUpper(ps,v)} returns the set consisting of the 
         ++ polynomials of \axiom{ps} with main variable greater than \axiom{v}.
     sort : ($,VarSet) -> Record(under:$,floor:$,upper:$)
         ++ \axiom{sort(v,ps)} returns \axiom{us,vs,ws} such that \axiom{us}
         ++ is \axiom{collectUnder(ps,v)}, \axiom{vs} is \axiom{collect(ps,v)}
         ++ and \axiom{ws} is \axiom{collectUpper(ps,v)}. 
     trivialIdeal?: $ -> Boolean
         ++ \axiom{trivialIdeal?(ps)} returns true iff \axiom{ps} does
         ++ not contain non-zero elements.
     if R has IntegralDomain
     then
       roughBase? : $ -> Boolean
           ++ \axiom{roughBase?(ps)} returns true iff for every pair \axiom{{p,q}}
           ++ of polynomials in \axiom{ps} their leading monomials are 
           ++ relatively prime.
       roughSubIdeal?  : ($,$) -> Boolean
           ++  \axiom{roughSubIdeal?(ps1,ps2)} returns true iff it can proved 
           ++ that all polynomials in  \axiom{ps1} lie in the ideal generated by 
           ++ \axiom{ps2} in  \axiom{\axiom{(R)^(-1) P}} without computing Groebner bases.
       roughEqualIdeals? : ($,$) -> Boolean
           ++ \axiom{roughEqualIdeals?(ps1,ps2)} returns true iff it can
           ++ proved that \axiom{ps1} and \axiom{ps2} generate the same ideal
           ++ in \axiom{(R)^(-1) P} without computing Groebner bases.
       roughUnitIdeal? : $ -> Boolean
           ++ \axiom{roughUnitIdeal?(ps)} returns true iff \axiom{ps} contains some
           ++ non null element lying in the base ring \axiom{R}.
       headRemainder : (P,$) -> Record(num:P,den:R)
           ++ \axiom{headRemainder(a,ps)} returns \axiom{[b,r]} such that the leading
           ++ monomial of \axiom{b} is reduced in the sense of Groebner bases w.r.t.
           ++ \axiom{ps} and \axiom{r*a - b} lies in the ideal generated by \axiom{ps}.
       remainder : (P,$) -> Record(rnum:R,polnum:P,den:R)
           ++ \axiom{remainder(a,ps)} returns \axiom{[c,b,r]} such that \axiom{b} is fully 
           ++ reduced in the sense of Groebner bases w.r.t. \axiom{ps}, 
           ++ \axiom{r*a - c*b} lies in the ideal generated by \axiom{ps}.
           ++ Furthermore, if \axiom{R} is a gcd-domain, \axiom{b} is primitive.
       rewriteIdealWithHeadRemainder : (List(P),$) -> List(P)
           ++ \axiom{rewriteIdealWithHeadRemainder(lp,cs)} returns \axiom{lr} such that
           ++ the leading monomial of every polynomial in \axiom{lr} is reduced
           ++ in the sense of Groebner bases w.r.t. \axiom{cs} and \axiom{(lp,cs)}
           ++ and \axiom{(lr,cs)} generate the same ideal in \axiom{(R)^(-1) P}.
       rewriteIdealWithRemainder : (List(P),$) -> List(P)
           ++ \axiom{rewriteIdealWithRemainder(lp,cs)} returns \axiom{lr} such that
           ++ every polynomial in \axiom{lr} is fully reduced in the sense 
           ++ of Groebner bases w.r.t. \axiom{cs} and \axiom{(lp,cs)} and 
           ++ \axiom{(lr,cs)} generate the same ideal in \axiom{(R)^(-1) P}.
       triangular? : $ -> Boolean
           ++ \axiom{triangular?(ps)} returns true iff \axiom{ps} is a triangular set,
           ++ i.e. two distinct polynomials have distinct main variables
           ++ and no constant lies in \axiom{ps}.

  add

     NNI ==> NonNegativeInteger
     B ==> Boolean

     elements: $ -> List(P)

     elements(ps:$):List(P) ==
       lp : List(P) := members(ps)$$

     variables1(lp:List(P)):(List VarSet) ==
       lvars : List(List(VarSet)) := [variables(p)$P for p in lp]
       sort(#1 > #2, removeDuplicates(concat(lvars)$List(VarSet)))

     variables2(lp:List(P)):(List VarSet) ==
       lvars : List(VarSet) := [mvar(p)$P for p in lp]
       sort(#1 > #2, removeDuplicates(lvars)$List(VarSet))

     variables (ps:$) ==
       variables1(elements(ps))

     mainVariables (ps:$) ==
       variables2(remove(ground?,elements(ps)))

     mainVariable? (v,ps) ==
       lp : List(P) := remove(ground?,elements(ps))
       while (not empty? lp) and (not (mvar(first(lp)) = v)) repeat
         lp := rest lp
       (not empty? lp)

     collectUnder (ps,v) ==
       lp : List P := elements(ps)
       lq : List P := []
       while (not empty? lp) repeat
         p := first lp
         lp := rest lp
         if (ground?(p)) or (mvar(p) < v)
           then
             lq := cons(p,lq)
       construct(lq)$$

     collectUpper (ps,v) ==
       lp : List P := elements(ps)
       lq : List P := []
       while (not empty? lp) repeat
         p := first lp
         lp := rest lp
         if (not ground?(p)) and (mvar(p) > v)
           then
             lq := cons(p,lq)
       construct(lq)$$

     collect (ps,v) ==
       lp : List P := elements(ps)
       lq : List P := []
       while (not empty? lp) repeat
         p := first lp
         lp := rest lp
         if (not ground?(p)) and (mvar(p) = v)
           then
             lq := cons(p,lq)
       construct(lq)$$

     sort (ps,v) ==
       lp : List P := elements(ps)
       us : List P := []
       vs : List P := []
       ws : List P := []
       while (not empty? lp) repeat
         p := first lp
         lp := rest lp
         if (ground?(p)) or (mvar(p) < v)
           then
             us := cons(p,us)
           else
             if (mvar(p) = v)
               then
                 vs := cons(p,vs)
               else
                 ws := cons(p,ws)
       [construct(us)$$,construct(vs)$$,construct(ws)$$]$Record(under:$,floor:$,upper:$)

     ps1 = ps2 ==
       {p for p in elements(ps1)} =$(Set P) {p for p in elements(ps2)}

     localInf? (p:P,q:P):B ==
       degree(p) <$E degree(q)

     localTriangular? (lp:List(P)):B ==
       lp := remove(zero?, lp)
       empty? lp => true
       any? (ground?, lp) => false
       lp := sort(mvar(#1)$P > mvar(#2)$P, lp)
       p,q : P
       p := first lp
       lp := rest lp
       while (not empty? lp) and (mvar(p) > mvar((q := first(lp)))) repeat
         p := q
         lp := rest lp
       empty? lp

     triangular? ps ==
       localTriangular? elements ps

     trivialIdeal? ps ==
       empty?(remove(zero?,elements(ps))$(List(P)))$(List(P))

     if R has IntegralDomain
     then

       roughUnitIdeal? ps ==
         any?(ground?,remove(zero?,elements(ps))$(List(P)))$(List P)

       relativelyPrimeLeadingMonomials? (p:P,q:P):B ==
         dp : E := degree(p)
         dq : E := degree(q)
         (sup(dp,dq)$E =$E dp +$E dq)@B

       roughBase? ps ==
         lp := remove(zero?,elements(ps))$(List(P))
         empty? lp => true
         rB? : B := true
         while (not empty? lp) and rB? repeat
           p := first lp
           lp := rest lp
           copylp := lp
           while (not empty? copylp) and rB? repeat
             rB? := relativelyPrimeLeadingMonomials?(p,first(copylp))
             copylp := rest copylp
         rB?

       roughSubIdeal?(ps1,ps2) ==
         lp: List(P) := rewriteIdealWithRemainder(elements(ps1),ps2)
         empty? (remove(zero?,lp))

       roughEqualIdeals? (ps1,ps2) ==
         ps1 =$$ ps2 => true
         roughSubIdeal?(ps1,ps2) and roughSubIdeal?(ps2,ps1)

     if (R has GcdDomain) and (VarSet has ConvertibleTo (Symbol))
     then

       LPR ==> List Polynomial R
       LS ==> List Symbol

       exactQuo(r:R,s:R):R ==
	 if R has EuclideanDomain then r quo$R s
         else (r exquo$R s)::R

       headRemainder (a,ps) ==
         lp1 : List(P) := remove(zero?, elements(ps))$(List(P))
         empty? lp1 => [a,1$R]
         any?(ground?,lp1) => [reductum(a),1$R]
         r : R := 1$R
         lp1 := sort(localInf?, reverse elements(ps))
         lp2 := lp1
         e : Union(E, "failed")
         while (not zero? a) and (not empty? lp2) repeat
           p := first lp2
           if ((e:= subtractIfCan(degree(a),degree(p))) case E)
             then
               g := gcd((lca := leadingCoefficient(a)),(lcp := leadingCoefficient(p)))$R
               (lca,lcp) := (exactQuo(lca,g),exactQuo(lcp,g))
               a := lcp * reductum(a) - monomial(lca, e::E)$P * reductum(p)
               r := r * lcp
               lp2 := lp1
             else
               lp2 := rest lp2
         [a,r]

       makeIrreducible! (frac:Record(num:P,den:R)):Record(num:P,den:R) ==
         g := gcd(frac.den,frac.num)$P
         one? g => frac
         frac.num := exactQuotient!(frac.num,g)
         frac.den := exactQuo(frac.den,g)
         frac

       remainder (a,ps) ==
         hRa := makeIrreducible! headRemainder (a,ps)
         a := hRa.num
         r : R := hRa.den
         zero? a => [1$R,a,r]
         b : P := monomial(1$R,degree(a))$P
         c : R := leadingCoefficient(a)
         while not zero?(a := reductum a) repeat
           hRa := makeIrreducible!  headRemainder (a,ps)
           a := hRa.num
           r := r * hRa.den
           g := gcd(c,(lca := leadingCoefficient(a)))$R
           b := ((hRa.den) * exactQuo(c,g)) * b + monomial(exactQuo(lca,g),degree(a))$P
           c := g
         [c,b,r]

       rewriteIdealWithHeadRemainder(ps,cs) ==
         trivialIdeal? cs => ps
         roughUnitIdeal? cs => [0$P]
         ps := remove(zero?,ps)
         empty? ps => ps
         any?(ground?,ps) => [1$P]
         rs : List P := []
         while not empty? ps repeat
           p := first ps
           ps := rest ps
           p := (headRemainder(p,cs)).num
           if not zero? p
             then 
               if ground? p
                 then
                   ps := []
                   rs := [1$P]
                 else
                   primitivePart! p
                   rs := cons(p,rs)
         removeDuplicates rs

       rewriteIdealWithRemainder(ps,cs) ==
         trivialIdeal? cs => ps
         roughUnitIdeal? cs => [0$P]
         ps := remove(zero?,ps)
         empty? ps => ps
         any?(ground?,ps) => [1$P]
         rs : List P := []
         while not empty? ps repeat
           p := first ps
           ps := rest ps
           p := (remainder(p,cs)).polnum
           if not zero? p
             then 
               if ground? p
                 then
                   ps := []
                   rs := [1$P]
                 else
                   rs := cons(unitCanonical(p),rs)
         removeDuplicates rs

)abbrev domain GPOLSET GeneralPolynomialSet
++ Author: Marc Moreno Maza
++ 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
++ References:
++ Description: A domain for polynomial sets.
++ Version: 1

GeneralPolynomialSet(R,E,VarSet,P) : Exports == Implementation where

  R:Ring
  VarSet:OrderedSet
  E:OrderedAbelianMonoidSup
  P:RecursivePolynomialCategory(R,E,VarSet)
  LP ==> List P
  PtoP ==> P -> P

  Exports ==  PolynomialSetCategory(R,E,VarSet,P)  with

     convert : LP -> $
       ++ \axiom{convert(lp)} returns the polynomial set whose members 
       ++ are the polynomials of \axiom{lp}.

     finiteAggregate
     shallowlyMutable

  Implementation == add

     Rep := List P

     construct lp ==
       (removeDuplicates(lp)$List(P))::$

     copy ps ==
       construct(copy(members(ps)$$)$LP)$$

     empty() ==
       []

     parts ps ==
       ps pretend LP

     map (f : PtoP, ps : $) : $ ==
       construct(map(f,members(ps))$LP)$$

     map! (f : PtoP, ps : $) : $  ==
       construct(map!(f,members(ps))$LP)$$

     member? (p,ps) ==
       member?(p,members(ps))$LP

     ps1 = ps2 ==
       {p for p in parts(ps1)} =$(Set P) {p for p in parts(ps2)}

     coerce(ps:$) : OutputForm ==
       lp : List(P) := sort(infRittWu?,members(ps))$(List P)
       brace([p::OutputForm for p in lp]$List(OutputForm))$OutputForm

     mvar ps ==
       empty? ps => error"Error from GPOLSET in mvar : #1 is empty"
       lv : List VarSet := variables(ps)
       empty? lv => error"Error from GPOLSET in mvar : every polynomial in #1 is constant"
       reduce(max,lv)$(List VarSet)

     retractIfCan(lp) ==
       (construct(lp))::Union($,"failed")

     coerce(ps:$) : (List P) ==
       ps pretend (List P)

     convert(lp:LP) : $ ==
       construct lp

