polset.spad

)abbrev category PSETCAT PolynomialSetCategory
++ Author: Marc Moreno Maza
++ Date Created: 04/26/1994
++ 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)),
         RetractableFrom(List(P)), finiteAggregate) with
     mvar : % -> VarSet
         ++  \spad{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
         ++  \spad{variables(ps)} returns the decreasingly sorted list of the
         ++  variables which are variables of some polynomial in \spad{ps}.
     mainVariables : %  -> List VarSet
         ++ \spad{mainVariables(ps)} returns the decreasingly sorted list of the
         ++ variables which are main variables of some polynomial in \spad{ps}.
     mainVariable? : (VarSet, %) -> Boolean
         ++ \spad{mainVariable?(v, ps)} returns true iff \spad{v} is the main variable
         ++ of some polynomial in \spad{ps}.
     collectUnder : (%, VarSet) -> %
         ++ \spad{collectUnder(ps, v)} returns the set consisting of the
         ++ polynomials of \spad{ps} with main variable less than \spad{v}.
     collect : (%, VarSet) -> %
         ++ \spad{collect(ps, v)}  returns the set consisting of the
         ++ polynomials of \spad{ps} with \spad{v} as main variable.
     collectUpper : (%, VarSet) -> %
         ++ \spad{collectUpper(ps, v)} returns the set consisting of the
         ++ polynomials of \spad{ps} with main variable greater than \spad{v}.
     sort : (%, VarSet) -> Record(under : %, floor : %, upper : %)
         ++ \spad{sort(v, ps)} returns \spad{us, vs, ws} such that \spad{us}
         ++ is \spad{collectUnder(ps, v)}, \spad{vs} is \spad{collect(ps, v)}
         ++ and \spad{ws} is \spad{collectUpper(ps, v)}.
     trivialIdeal? : % -> Boolean
         ++ \spad{trivialIdeal?(ps)} returns true iff \spad{ps} does
         ++ not contain non-zero elements.
     if R has IntegralDomain
     then
       roughBase? : % -> Boolean
           ++ \spad{roughBase?(ps)} returns true iff for every pair \spad{{p, q}}
           ++ of polynomials in \spad{ps} their leading monomials are
           ++ relatively prime.
       roughSubIdeal?  : (%, %) -> Boolean
           ++  \spad{roughSubIdeal?(ps1, ps2)} returns true iff it can proved
           ++ that all polynomials in  \spad{ps1} lie in the ideal generated by
           ++ \spad{ps2} in \spad{(R)^(-1) P} without computing Groebner bases.
       roughEqualIdeals? : (%, %) -> Boolean
           ++ \spad{roughEqualIdeals?(ps1, ps2)} returns true iff it can
           ++ proved that \spad{ps1} and \spad{ps2} generate the same ideal
           ++ in \spad{(R)^(-1) P} without computing Groebner bases.
       roughUnitIdeal? : % -> Boolean
           ++ \spad{roughUnitIdeal?(ps)} returns true iff \spad{ps} contains some
           ++ non null element lying in the base ring \spad{R}.
       headRemainder : (P, %) -> Record(num : P, den : R)
           ++ \spad{headRemainder(a, ps)} returns \spad{[b, r]} such that the leading
           ++ monomial of \spad{b} is reduced in the sense of Groebner bases w.r.t.
           ++ \spad{ps} and \spad{r*a - b} lies in the ideal generated by \spad{ps}.
       remainder : (P, %) -> Record(rnum : R, polnum : P, den : R)
           ++ \spad{remainder(a, ps)} returns \spad{[c, b, r]} such that \spad{b} is fully
           ++ reduced in the sense of Groebner bases w.r.t. \spad{ps},
           ++ \spad{r*a - c*b} lies in the ideal generated by \spad{ps}.
           ++ Furthermore, if \spad{R} is a gcd-domain, \spad{b} is primitive.
       rewriteIdealWithHeadRemainder : (List(P), %) -> List(P)
           ++ \spad{rewriteIdealWithHeadRemainder(lp, cs)} returns \spad{lr} such that
           ++ the leading monomial of every polynomial in \spad{lr} is reduced
           ++ in the sense of Groebner bases w.r.t. \spad{cs} and \spad{(lp, cs)}
           ++ and \spad{(lr, cs)} generate the same ideal in \spad{(R)^(-1) P}.
       rewriteIdealWithRemainder : (List(P), %) -> List(P)
           ++ \spad{rewriteIdealWithRemainder(lp, cs)} returns \spad{lr} such that
           ++ every polynomial in \spad{lr} is fully reduced in the sense
           ++ of Groebner bases w.r.t. \spad{cs} and \spad{(lp, cs)} and
           ++ \spad{(lr, cs)} generate the same ideal in \spad{(R)^(-1) P}.
       triangular? : % -> Boolean
           ++ \spad{triangular?(ps)} returns true iff \spad{ps} is a triangular set,
           ++ i.e. two distinct polynomials have distinct main variables
           ++ and no constant lies in \spad{ps}.
       iexactQuo : (R, R) -> R
           ++ iexactQuo(x, y) should be local but conditional

  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((z1 : VarSet, z2 : VarSet) : Boolean +-> z1 > z2, removeDuplicates(concat(lvars)$List(VarSet)))

     variables2(lp : List(P)) : (List VarSet) ==
       lvars : List(VarSet) := [mvar(p)$P for p in lp]
       sort((z1 : VarSet, z2 : VarSet) : Boolean +-> z1 > z2, 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 ==
       set([p for p in elements(ps1)]) =$(Set P)
          set([p for p in elements(ps2)])

     exactQuo : (R, R) -> R

     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((z1 : P, z2 : P) : Boolean +-> mvar(z1)$P > mvar(z2)$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

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

       exactQuo(r : R, s : R) : R == iexactQuo(r, s)

       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
         (g = 1) => 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
++ 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 ==  Join(PolynomialSetCategory(R, E, VarSet, P), finiteAggregate, _
       shallowlyMutable)  with

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

  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 ==
       set([p for p in parts(ps1)]) =$(Set P)
         set([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

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