gbintern.spad

)abbrev package GBINTERN GroebnerInternalPackage
++ Author:
++ Keywords:
++ References:
++   R. Gebauer, M. Moeller,  On an Installation of Buchberger's Algorithm,
++   JSC (1988) 6, 275-286.
++ Description
++ This package provides low level tools for Groebner basis computations
GroebnerInternalPackage(Dom, Expon, Dpol) : T == C where
 Dom :   GcdDomain
 Expon : OrderedAbelianMonoidSup
 Dpol : Join(FiniteAbelianMonoidRing(Dom, Expon),
             VariablesCommuteWithCoefficients)
 NNI    ==> NonNegativeInteger
   ------  Definition of Record critPair and Prinp

 critPair ==> Record(lcmfij : Expon, totdeg : NonNegativeInteger,
                      poli : Dpol, polj : Dpol)
 sugarPol ==> Record(totdeg : NonNegativeInteger, pol : Dpol)
 Prinp    ==> Record(ci : Dpol, tci : Integer, cj : Dpol, tcj : Integer,
                     c : Dpol, tc : Integer, rc : Dpol, trc : Integer,
                     tF : Integer, tD : Integer)
 Prinpp   ==> Record(ci : Dpol, tci : Integer, cj : Dpol, tcj : Integer,
                     c : Dpol, tc : Integer, rc : Dpol, trc : Integer,
                     tF : Integer, tDD : Integer, tDF : Integer)
 T== with

     credPol :  (Dpol, List(Dpol))  -> Dpol
        ++ credPol \undocumented
     redPol :   (Dpol, List(Dpol))  -> Dpol
        ++ redPol \undocumented
     gbasis :  (List(Dpol), Integer, Integer) -> List(Dpol)
        ++ gbasis(pols, f1, f2) computes a Groebner basis of pols.
     gbasisExtend: (List(Dpol), List(Dpol), Integer, Integer) -> List(Dpol)
        ++ gbasisExtend(pols, gb, f1, f2) computes a Groebner basis of
        ++ concat(pols, gb) where gb is assumed to be already a Groebner
        ++ basis, i.e. no S-polynomials among elements of gb will ever be
        ++ considered.
     critT :  critPair   -> Boolean
        ++ critT \undocumented
     critM :  (Expon, Expon) -> Boolean
        ++ critM \undocumented
     critB :  (Expon, Expon, Expon, Expon) -> Boolean
        ++ critB \undocumented
     critBonD :  (Dpol, List(critPair)) -> List(critPair)
        ++ critBonD \undocumented
     critMTonD1 : (List(critPair)) -> List(critPair)
        ++ critMTonD1 \undocumented
     critMonD1 :  (Expon, List(critPair)) -> List(critPair)
        ++ critMonD1 \undocumented
     redPo : (Dpol, List(Dpol) )  ->  Record(poly : Dpol, mult : Dom)
        ++ redPo \undocumented
     hMonic :  Dpol  -> Dpol
        ++ hMonic \undocumented
     updatF : (Dpol, NNI, List(sugarPol) ) -> List(sugarPol)
        ++ updatF \undocumented
     sPol :  critPair  -> Dpol
        ++ sPol \undocumented
     updatD : (List(critPair), List(critPair)) -> List(critPair)
        ++ updatD \undocumented
     minGbasis : List(Dpol) -> List(Dpol)
        ++ minGbasis \undocumented
     lepol : Dpol -> Integer
        ++ lepol \undocumented
     prinshINFO : Dpol -> Void
        ++ prinshINFO \undocumented
     prindINFO : (critPair, Dpol, Dpol, Integer, Integer, Integer) -> Integer
        ++ prindINFO \undocumented
     fprindINFO : (critPair, Dpol, Dpol, Integer, Integer, Integer
                 , Integer) ->  Integer
        ++ fprindINFO \undocumented
     prinpolINFO : List(Dpol) -> Void
        ++ prinpolINFO \undocumented
     prinb : Integer-> Void
        ++ prinb \undocumented
     critpOrder : (critPair, critPair) -> Boolean
        ++ critpOrder \undocumented
     makeCrit : (sugarPol, Dpol, NonNegativeInteger) -> critPair
        ++ makeCrit \undocumented
     virtualDegree : Dpol -> NonNegativeInteger
        ++ virtualDegree \undocumented

 C== add
   Ex ==> OutputForm
   import from OutputForm

   ------  Definition of intermediate functions
   if Dpol has totalDegree : Dpol -> NonNegativeInteger then
     virtualDegree p == totalDegree p
   else
     virtualDegree p == 0

   ------  ordering of critpairs

   critpOrder(cp1, cp2) ==
     cp1.totdeg < cp2.totdeg => true
     cp2.totdeg < cp1.totdeg => false
     cp1.lcmfij < cp2.lcmfij

   ------    creating a critical pair

   makeCrit(sp1, p2, totdeg2) ==
     p1 := sp1.pol
     deg := sup(degree(p1), degree(p2))
     e1 := subtractIfCan(deg, degree(p1))::Expon
     e2 := subtractIfCan(deg, degree(p2))::Expon
     tdeg := max(sp1.totdeg + virtualDegree(monomial(1, e1)),
                 totdeg2 + virtualDegree(monomial(1, e2)))
     [deg, tdeg, p1, p2]$critPair

   ------    calculate basis

   gbasis(Pol: List(Dpol), xx1: Integer, xx2: Integer) ==
       gbasisExtend(Pol, [], xx1, xx2)

   gbasisExtend(Pol: List(Dpol), gb: List(Dpol), xx1: Integer, xx2: Integer) ==
     D, D1 : List(critPair)
     ---------   create D and Pol

     sgb: List(sugarPol) := [[virtualDegree x, hMonic x] for x in gb]
     Pol1 := sort((z1, z2) +-> degree z1 > degree z2, Pol)
     basPols := updatF(hMonic(first Pol1), virtualDegree(first Pol1), sgb)
     Pol1 := rest(Pol1)
     D := []
     while not(empty?(Pol1)) repeat
        h := hMonic(first(Pol1))
        Pol1 := rest(Pol1)
        toth := virtualDegree h
        D1 := [makeCrit(x, h, toth) for x in basPols]
        D := updatD(critMTonD1(sort(critpOrder, D1)),
                   critBonD(h, D))
        basPols := updatF(h, toth, basPols)
     D := sort(critpOrder, D)
     xx := xx2
     --------  loop

     redPols := [x.pol for x in basPols]
     while not(empty?(D)) repeat
         D0 := first D
         s := hMonic(sPol(D0))
         D := rest(D)
         h := hMonic(redPol(s, redPols))
         if xx1 = 1  then
              prinshINFO(h)
         h = 0  =>
          if xx2 = 1 then
           prindINFO(D0, s, h, # basPols, # D, xx)
           xx := 2
          " go to top of while "
         degree(h) = 0 =>
           D := []
           if xx2 = 1 then
            prindINFO(D0, s, h, # basPols, # D, xx)
            xx := 2
           basPols := updatF(h, 0, [])
           break
         D1 := [makeCrit(x, h, D0.totdeg) for x in basPols]
         D := updatD(critMTonD1(sort(critpOrder, D1)),
                   critBonD(h, D))
         basPols := updatF(h, D0.totdeg, basPols)
         redPols := concat(redPols, h)
         if xx2 = 1 then
            prindINFO(D0, s, h, # basPols, # D, xx)
            xx := 2
     Pol := [x.pol for x in basPols]
     if xx2 = 1 then
       prinpolINFO(Pol)
       messagePrint("    THE GROEBNER BASIS POLYNOMIALS")
     if xx1 = 1 and xx2 ~= 1 then
       messagePrint("    THE GROEBNER BASIS POLYNOMIALS")
     Pol

             --------------------------------------

             --- erase multiple of e in D2 using crit M

   critMonD1(e : Expon, D2 : List(critPair))==
      [x for x in D2 | not critM(e, x.lcmfij)]

             ----------------------------

             --- reduce D1 using crit T and crit M

   critMTonD1(D1 : List(critPair))==
           empty?(D1) => []
           f1 := first(D1)
           s1 := #(D1)
           cT1 := critT(f1)
           s1= 1 and cT1 => []
           s1= 1 => D1
           e1 := f1.lcmfij
           r1 := rest(D1)
           e1 = (first r1).lcmfij  =>
              cT1 =>   critMTonD1(cons(f1, rest(r1)))
              critMTonD1(r1)
           D1 := critMonD1(e1, r1)
           cT1 => critMTonD1(D1)
           cons(f1, critMTonD1(D1))


             -----------------------------

             --- erase elements in D fullfilling crit B

   critBonD(h : Dpol, D : List(critPair))==
       d := degree(h)
       [x for x in D | not critB(d, x.lcmfij, degree(x.poli), degree(x.polj))]

             -----------------------------

             --- concat F and h and erase multiples of h in F

   updatF(h : Dpol, deg : NNI, F : List(sugarPol)) ==
       d := degree h
       concat!([f for f in F | not critM(d, degree(f.pol))], [deg, h]$sugarPol)

             -----------------------------

             --- merge ordered critical pair lists D1 and D2

   updatD(D1 : List(critPair), D2 : List(critPair)) ==
      empty?(D1) => D2
      empty?(D2) => D1
      res: List(critPair) := []
      while not empty? D1 and not empty? D2 repeat
         dl1 := first D1
         dl2 := first D2
         if critpOrder(dl1, dl2) then
            res := cons(dl1, res)
            D1 := rest D1
         else
            res := cons(dl2, res)
            D2 := rest D2
      for e in D1 repeat res := cons(e, res)
      for e in D2 repeat res := cons(e, res)
      reverse! res

            -----------------------------

            --- remove gcd from pair of coefficients

   gcdCo(c1 : Dom, c2 : Dom) : Record(co1 : Dom, co2 : Dom) ==
      d := gcd(c1, c2)
      [(c1 exquo d)::Dom, (c2 exquo d)::Dom]

            --- calculate S-polynomial of a critical pair

   sPol(p : critPair)==
      Tij := p.lcmfij
      fi := p.poli
      fj := p.polj
      cc := gcdCo(leadingCoefficient fi, leadingCoefficient fj)
      reductum(fi)*monomial(cc.co2, subtractIfCan(Tij, degree fi)::Expon) -
        reductum(fj)*monomial(cc.co1, subtractIfCan(Tij, degree fj)::Expon)

            ----------------------------

            --- reduce critpair polynomial mod F
            --- iterative version

   redPo(s : Dpol, F : List(Dpol)) ==
      m : Dom := 1
      Fh := F
      while not(s = 0 or empty?(F)) repeat
        f1 := first(F)
        s1 := degree(s)
        e : Union(Expon, "failed")
        (e := subtractIfCan(s1, degree(f1))) case Expon  =>
           cc := gcdCo(leadingCoefficient f1, leadingCoefficient s)
           s := cc.co1*reductum(s) - monomial(cc.co2, e)*reductum(f1)
           m := m*cc.co1
           F := Fh
        F := rest F
      [s, m]

   redPol(s : Dpol, F : List(Dpol)) ==  credPol(redPo(s, F).poly, F)

            ----------------------------

            --- crit T  true, if e1 and e2 are disjoint

   critT(p : critPair) == p.lcmfij =  (degree(p.poli) + degree(p.polj))

            ----------------------------

            --- crit M - true, if lcm#2 multiple of lcm#1

   critM(e1 : Expon, e2 : Expon) ==
         en : Union(Expon, "failed")
         (en := subtractIfCan(e2, e1)) case Expon

            ----------------------------

            --- crit B - true, if eik is a multiple of eh and eik not equal
            ---          lcm(eh, ei) and eik not equal lcm(eh, ek)

   critB(eh : Expon, eik : Expon, ei : Expon, ek : Expon) ==
       critM(eh, eik) and (eik ~= sup(eh, ei)) and (eik ~= sup(eh, ek))

            ----------------------------

            ---  make polynomial monic case Domain a Field

   hMonic(p : Dpol) ==
        p= 0 => p
        -- inv(leadingCoefficient(p))*p
        primitivePart p

            -----------------------------

            ---  reduce all terms of h mod F  (iterative version )

   credPol(h : Dpol, F : List(Dpol) ) ==
        empty?(F) => h
        h0 : Dpol := monomial(leadingCoefficient h, degree h)
        while (h := reductum h) ~= 0 repeat
           hred := redPo(h, F)
           h := hred.poly
           h0 := (hred.mult)*h0 + monomial(leadingCoefficient(h), degree h)
        h0

            -------------------------------

            ----  calculate minimal basis for ordered F

   minGbasis(F : List(Dpol)) ==
        empty?(F) => []
        newbas := minGbasis rest F
        cons(hMonic credPol( first(F), newbas), newbas)

            -------------------------------

            ----  calculate number of terms of polynomial

   lepol(p1 : Dpol)==
      n : Integer
      n := 0
      while p1 ~= 0 repeat
         n := n + 1
         p1 := reductum(p1)
      n

            ----  print blanc lines

   prinb(n : Integer)==
      for x in 1..n repeat
         messagePrint("    ")

            ----  print reduced critpair polynom

   prinshINFO(h : Dpol)==
           prinb(2)
           messagePrint(" reduced Critpair - Polynom :")
           prinb(2)
           print(h::Ex)
           prinb(2)

            -------------------------------

            ----  print info string

   prindINFO(cp : critPair, ps : Dpol, ph : Dpol, i1 : Integer,
             i2 : Integer, n : Integer) ==
       ll : List Prinp
       a : Dom
       cpi := cp.poli
       cpj := cp.polj
       if n = 1 then
        prinb(1)
        messagePrint("you choose option  -info-  ")
        messagePrint("abbrev. for the following information strings are")
        messagePrint("  ci  =>  Leading monomial  for critpair calculation")
        messagePrint("  tci =>  Number of terms of polynomial i")
        messagePrint("  cj  =>  Leading monomial  for critpair calculation")
        messagePrint("  tcj =>  Number of terms of polynomial j")
        messagePrint("  c   =>  Leading monomial of critpair polynomial")
        messagePrint("  tc  =>  Number of terms of critpair polynomial")
        messagePrint("  rc  =>  Leading monomial of redcritpair polynomial")
        messagePrint("  trc =>  Number of terms of redcritpair polynomial")
        messagePrint("  tF  =>  Number of polynomials in reduction list F")
        messagePrint("  tD  =>  Number of critpairs still to do")
        prinb(4)
        n := 2
       prinb(1)
       a := 1
       ph = 0  =>
          ps = 0 =>
            ll := [[monomial(a, degree(cpi)), lepol(cpi),
                  monomial(a, degree(cpj)),
                   lepol(cpj), ps, 0, ph, 0, i1, i2]$Prinp]
            print(ll::Ex)
            prinb(1)
            n
          ll := [[monomial(a, degree(cpi)), lepol(cpi),
             monomial(a, degree(cpj)), lepol(cpj), monomial(a, degree(ps)),
              lepol(ps), ph, 0, i1, i2]$Prinp]
          print(ll::Ex)
          prinb(1)
          n
       ll := [[monomial(a, degree(cpi)), lepol(cpi),
            monomial(a, degree(cpj)), lepol(cpj), monomial(a, degree(ps)),
             lepol(ps), monomial(a, degree(ph)), lepol(ph), i1, i2]$Prinp]
       print(ll::Ex)
       prinb(1)
       n

            -------------------------------

            ----  print the groebner basis polynomials

   prinpolINFO(pl : List(Dpol))==
       n : Integer
       n := # pl
       prinb(1)
       n = 1 =>
         messagePrint("  There is 1  Groebner Basis Polynomial ")
         prinb(2)
       messagePrint("  There are ")
       prinb(1)
       print(n::Ex)
       prinb(1)
       messagePrint("  Groebner Basis Polynomials. ")
       prinb(2)

   fprindINFO(cp : critPair, ps : Dpol, ph : Dpol, i1 : Integer,
             i2 : Integer, i3 : Integer, n : Integer) ==
       ll : List Prinpp
       a : Dom
       cpi := cp.poli
       cpj := cp.polj
       if n = 1 then
        prinb(1)
        messagePrint("you choose option  -info-  ")
        messagePrint("abbrev. for the following information strings are")
        messagePrint("  ci  =>  Leading monomial  for critpair calculation")
        messagePrint("  tci =>  Number of terms of polynomial i")
        messagePrint("  cj  =>  Leading monomial  for critpair calculation")
        messagePrint("  tcj =>  Number of terms of polynomial j")
        messagePrint("  c   =>  Leading monomial of critpair polynomial")
        messagePrint("  tc  =>  Number of terms of critpair polynomial")
        messagePrint("  rc  =>  Leading monomial of redcritpair polynomial")
        messagePrint("  trc =>  Number of terms of redcritpair polynomial")
        messagePrint("  tF  =>  Number of polynomials in reduction list F")
        messagePrint("  tD  =>  Number of critpairs still to do")
        messagePrint("  tDF =>  Number of subproblems still to do")
        prinb(4)
        n := 2
       prinb(1)
       a := 1
       ph = 0  =>
          ps = 0 =>
            ll := [[monomial(a, degree(cpi)), lepol(cpi),
              monomial(a, degree(cpj)),
               lepol(cpj), ps, 0, ph, 0, i1, i2, i3]$Prinpp]
            print(ll::Ex)
            prinb(1)
            n
          ll := [[monomial(a, degree(cpi)), lepol(cpi),
            monomial(a, degree(cpj)), lepol(cpj), monomial(a, degree(ps)),
             lepol(ps), ph, 0, i1, i2, i3]$Prinpp]
          print(ll::Ex)
          prinb(1)
          n
       ll := [[monomial(a, degree(cpi)), lepol(cpi),
            monomial(a, degree(cpj)), lepol(cpj), monomial(a, degree(ps)),
             lepol(ps), monomial(a, degree(ph)), lepol(ph), i1, i2, i3]$Prinpp]
       print(ll::Ex)
       prinb(1)
       n

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