--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 package GBEUCLID EuclideanGroebnerBasisPackage
++ Authors: Gebauer, Moeller
++ Date Created: 12-1-86
++ Date Last Updated: 2-28-91
++ Basic Functions:
++ Related Constructors: Ideal, IdealDecompositionPackage, GroebnerPackage
++ Also See:
++ AMS Classifications:
++ Keywords: groebner basis, polynomial ideal, euclidean domain
++ References:
++ Description: \spadtype{EuclideanGroebnerBasisPackage} computes groebner
++ bases for polynomial ideals over euclidean domains.
++ The basic computation provides
++ a distinguished set of generators for these ideals.
++ This basis allows an easy test for membership: the operation
++ \spadfun{euclideanNormalForm} returns zero on ideal members. The string 
++ "info" and "redcrit" can be given as additional args to provide 
++ incremental information during the computation. If "info" is given,
++  a computational summary is given for each s-polynomial. If "redcrit" 
++ is given, the reduced critical pairs are printed. The term ordering
++ is determined by the polynomial type used. Suggested types include
++ \spadtype{DistributedMultivariatePolynomial},
++ \spadtype{HomogeneousDistributedMultivariatePolynomial},
++ \spadtype{GeneralDistributedMultivariatePolynomial}.
 
EuclideanGroebnerBasisPackage(Dom, Expon, VarSet, Dpol): T == C where
 
 Dom: EuclideanDomain
 Expon: OrderedAbelianMonoidSup
 VarSet: OrderedSet
 Dpol: PolynomialCategory(Dom, Expon, VarSet)
 
 T== with
 
     euclideanNormalForm: (Dpol, List(Dpol) )  ->  Dpol
       ++ euclideanNormalForm(poly,gb) reduces the polynomial poly modulo the
       ++ precomputed groebner basis gb giving a canonical representative
       ++ of the residue class.
     euclideanGroebner: List(Dpol) -> List(Dpol)
       ++ euclideanGroebner(lp) computes a groebner basis for a polynomial ideal
       ++ over a euclidean domain generated by the list of polynomials lp.
     euclideanGroebner: (List(Dpol), String) -> List(Dpol)
       ++ euclideanGroebner(lp, infoflag) computes a groebner basis 
       ++ for a polynomial ideal over a euclidean domain
       ++ generated by the list of polynomials lp.
       ++ During computation, additional information is printed out
       ++ if infoflag is given as 
       ++ either "info" (for summary information) or
       ++ "redcrit" (for reduced critical pairs)
     euclideanGroebner: (List(Dpol), String, String ) -> List(Dpol)
       ++ euclideanGroebner(lp, "info", "redcrit") computes a groebner basis
       ++ for a polynomial ideal generated by the list of polynomials lp.
       ++ If the second argument is "info", a summary is given of the critical pairs.
       ++ If the third argument is "redcrit", critical pairs are printed.
 C== add
   Ex ==> OutputForm
   lc ==> leadingCoefficient
   red ==> reductum

   import OutputForm
 
   ------  Definition list of critPair
   ------  lcmfij is now lcm of headterm of poli and polj
   ------  lcmcij is now lcm of of lc poli and lc polj
 
   critPair ==>Record(lcmfij: Expon, lcmcij: Dom, poli:Dpol, polj: Dpol )
   Prinp    ==> Record( ci:Dpol,tci:Integer,cj:Dpol,tcj:Integer,c:Dpol,
                tc:Integer,rc:Dpol,trc:Integer,tH:Integer,tD:Integer)
 
   ------  Definition of intermediate functions
 
   strongGbasis: (List(Dpol), Integer, Integer) -> List(Dpol)
   eminGbasis: List(Dpol) -> List(Dpol)
   ecritT: (critPair ) -> Boolean
   ecritM: (Expon, Dom, Expon, Dom) -> Boolean
   ecritB: (Expon, Dom, Expon, Dom, Expon, Dom) -> Boolean
   ecrithinH: (Dpol, List(Dpol)) -> Boolean
   ecritBonD: (Dpol, List(critPair)) -> List(critPair)
   ecritMTondd1:(List(critPair)) -> List(critPair)
   ecritMondd1:(Expon, Dom, List(critPair)) -> List(critPair)
   crithdelH: (Dpol, List(Dpol)) -> List(Dpol)
   eupdatF: (Dpol, List(Dpol) ) -> List(Dpol)
   updatH: (Dpol, List(Dpol), List(Dpol), List(Dpol) ) -> List(Dpol)
   sortin: (Dpol, List(Dpol) ) -> List(Dpol)
   eRed: (Dpol, List(Dpol), List(Dpol) )  ->  Dpol
   ecredPol: (Dpol, List(Dpol) ) -> Dpol
   esPol: (critPair) -> Dpol
   updatD: (List(critPair), List(critPair)) -> List(critPair)
   lepol: Dpol -> Integer
   prinshINFO : Dpol -> Void
   prindINFO: (critPair, Dpol, Dpol,Integer,Integer,Integer) -> Integer
   prinpolINFO: List(Dpol) -> Void
   prinb: Integer -> Void
 
   ------    MAIN ALGORITHM GROEBNER ------------------------
   euclideanGroebner( Pol: List(Dpol) ) ==
     eminGbasis(strongGbasis(Pol,0,0))
 
   euclideanGroebner( Pol: List(Dpol), xx1: String) ==
     xx1 = "redcrit" =>
       eminGbasis(strongGbasis(Pol,1,0))
     xx1 = "info" =>
       eminGbasis(strongGbasis(Pol,2,1))
     print("   "::Ex)
     print("WARNING: options are - redcrit and/or info - "::Ex)
     print("         you didn't type them correct"::Ex)
     print("         please try again"::Ex)
     print("   "::Ex)
     []
 
   euclideanGroebner( Pol: List(Dpol), xx1: String, xx2: String) ==
     (xx1 = "redcrit" and xx2 = "info") or
      (xx1 = "info" and xx2 = "redcrit")   =>
       eminGbasis(strongGbasis(Pol,1,1))
     xx1 = "redcrit" and xx2 = "redcrit" =>
       eminGbasis(strongGbasis(Pol,1,0))
     xx1 = "info" and xx2 = "info" =>
       eminGbasis(strongGbasis(Pol,2,1))
     print("   "::Ex)
     print("WARNING:  options are - redcrit and/or info - "::Ex)
     print("          you didn't type them correct"::Ex)
     print("          please try again "::Ex)
     print("   "::Ex)
     []
 
   ------    calculate basis
 
   strongGbasis(Pol: List(Dpol),xx1: Integer, xx2: Integer ) ==
     dd1, D : List(critPair)
 
     ---------   create D and Pol
 
     Pol1:= sort( (degree #1 > degree #2) or
                    ((degree #1 = degree #2 ) and
                       sizeLess?(leadingCoefficient #2,leadingCoefficient #1)),
                 Pol)
     Pol:= [first(Pol1)]
     H:= Pol
     Pol1:= rest(Pol1)
     D:= nil
     while not null Pol1 repeat
        h:= first(Pol1)
        Pol1:= rest(Pol1)
        en:= degree(h)
        lch:= lc h
        dd1:= [[sup(degree(x), en), lcm(leadingCoefficient x, lch), x, h]$critPair
            for x in Pol]
        D:= updatD(ecritMTondd1(sort((#1.lcmfij < #2.lcmfij) or
                                      (( #1.lcmfij = #2.lcmfij ) and
                                        ( sizeLess?(#1.lcmcij,#2.lcmcij)) ),
                                     dd1)), ecritBonD(h,D))
        Pol:= cons(h, eupdatF(h, Pol))
        ((en = degree(first(H))) and (leadingCoefficient(h) = leadingCoefficient(first(H)) ) ) =>
              " go to top of while "
        H:= updatH(h,H,crithdelH(h,H),[h])
        H:= sort((degree #1 > degree #2) or
                ((degree #1 = degree #2 ) and
                  sizeLess?(leadingCoefficient #2,leadingCoefficient #1)), H)
     D:= sort((#1.lcmfij < #2.lcmfij) or
             (( #1.lcmfij = #2.lcmfij ) and
               ( sizeLess?(#1.lcmcij,#2.lcmcij)) ) ,D)
     xx:= xx2
 
     --------  loop
 
     while not null D repeat
         D0:= first D
         ep:=esPol(D0)
         D:= rest(D)
         eh:= ecredPol(eRed(ep,H,H),H)
         if xx1 = 1 then
               prinshINFO(eh)
         eh = 0 =>
              if xx2 = 1 then
                  ala:= prindINFO(D0,ep,eh,#H, #D, xx)
                  xx:= 2
              " go to top of while "
         eh := unitCanonical eh
         e:= degree(eh)
         leh:= lc eh
         dd1:= [[sup(degree(x), e), lcm(leadingCoefficient x, leh), x, eh]$critPair
            for x in Pol]
         D:= updatD(ecritMTondd1(sort( (#1.lcmfij <
              #2.lcmfij) or (( #1.lcmfij = #2.lcmfij ) and
               ( sizeLess?(#1.lcmcij,#2.lcmcij)) ), dd1)), ecritBonD(eh,D))
         Pol:= cons(eh,eupdatF(eh,Pol))
         not ecrithinH(eh,H) or
           ((e = degree(first(H))) and (leadingCoefficient(eh) = leadingCoefficient(first(H)) ) ) =>
              if xx2 = 1 then
                  ala:= prindINFO(D0,ep,eh,#H, #D, xx)
                  xx:= 2
              " go to top of while "
         H:= updatH(eh,H,crithdelH(eh,H),[eh])
         H:= sort( (degree #1 > degree #2) or
             ((degree #1 = degree #2 ) and
                 sizeLess?(leadingCoefficient #2,leadingCoefficient #1)), H)
         if xx2 = 1 then
           ala:= prindINFO(D0,ep,eh,#H, #D, xx)
           xx:= 2
           " go to top of while "
     if xx2 = 1 then
       prinpolINFO(Pol)
       print("    THE GROEBNER BASIS over EUCLIDEAN DOMAIN"::Ex)
     if xx1 = 1 and not one? xx2 then
       print("    THE GROEBNER BASIS over EUCLIDEAN DOMAIN"::Ex)
     H
 
             --------------------------------------
 
             --- erase multiple of e in D2 using crit M
 
   ecritMondd1(e: Expon, c: Dom, D2: List(critPair))==
      null D2 => nil
      x:= first(D2)
      ecritM(e,c, x.lcmfij, lcm(leadingCoefficient(x.poli), leadingCoefficient(x.polj)))
         => ecritMondd1(e, c, rest(D2))
      cons(x, ecritMondd1(e, c, rest(D2)))
 
            -------------------------------
 
   ecredPol(h: Dpol, F: List(Dpol) ) ==
        h0:Dpol:= 0
        null F => h
        while h ~= 0 repeat
           h0:= h0 + monomial(leadingCoefficient(h),degree(h))
           h:= eRed(red(h), F, F)
        h0
             ----------------------------
 
             --- reduce dd1 using crit T and crit M
 
   ecritMTondd1(dd1: List(critPair))==
           null dd1 => nil
           f1:= first(dd1)
           s1:= #(dd1)
           cT1:= ecritT(f1)
           s1= 1 and cT1 => nil
           s1= 1 => dd1
           e1:= f1.lcmfij
           r1:= rest(dd1)
           f2:= first(r1)
           e1 = f2.lcmfij and f1.lcmcij = f2.lcmcij =>
              cT1 =>   ecritMTondd1(cons(f1, rest(r1)))
              ecritMTondd1(r1)
           dd1 := ecritMondd1(e1, f1.lcmcij, r1)
           cT1 => ecritMTondd1(dd1)
           cons(f1, ecritMTondd1(dd1))
 
             -----------------------------
 
             --- erase elements in D fullfilling crit B
 
   ecritBonD(h:Dpol, D: List(critPair))==
         null D => nil
         x:= first(D)
         x1:= x.poli
         x2:= x.polj
         ecritB(degree(h), leadingCoefficient(h), degree(x1),leadingCoefficient(x1),degree(x2),leadingCoefficient(x2)) =>
           ecritBonD(h, rest(D))
         cons(x, ecritBonD(h, rest(D)))
 
             -----------------------------
 
             --- concat F and h and erase multiples of h in F
 
   eupdatF(h: Dpol, F: List(Dpol)) ==
       null F => nil
       f1:= first(F)
       ecritM(degree h, leadingCoefficient(h), degree f1, leadingCoefficient(f1))
           => eupdatF(h, rest(F))
       cons(f1, eupdatF(h, rest(F)))
 
             -----------------------------
             --- concat H and h and erase multiples of h in H
 
   updatH(h: Dpol, H: List(Dpol), Hh: List(Dpol), Hhh: List(Dpol)) ==
       null H => append(Hh,Hhh)
       h1:= first(H)
       hlcm:= sup(degree(h1), degree(h))
       plc:= extendedEuclidean(leadingCoefficient(h), leadingCoefficient(h1))
       hp:= monomial(plc.coef1,subtractIfCan(hlcm, degree(h))::Expon)*h +
            monomial(plc.coef2,subtractIfCan(hlcm, degree(h1))::Expon)*h1
       (ecrithinH(hp, Hh) and ecrithinH(hp, Hhh)) =>
         hpp:= append(rest(H),Hh)
         hp:= ecredPol(eRed(hp,hpp,hpp),hpp)
         updatH(h, rest(H), crithdelH(hp,Hh),cons(hp,crithdelH(hp,Hhh)))
       updatH(h, rest(H), Hh,Hhh)
 
             --------------------------------------------------
             ---- delete elements in cons(h,H)
 
   crithdelH(h: Dpol, H: List(Dpol))==
        null H => nil
        h1:= first(H)
        dh1:= degree h1
        dh:= degree h
        ecritM(dh, lc h, dh1, lc h1) => crithdelH(h, rest(H))
        dh1 = sup(dh,dh1) =>
           plc:= extendedEuclidean( lc h1, lc h)
           cons(plc.coef1*h1 + monomial(plc.coef2,subtractIfCan(dh1,dh)::Expon)*h,
               crithdelH(h,rest(H)))
        cons(h1, crithdelH(h,rest(H)))
 
   eminGbasis(F: List(Dpol)) ==
        null F => nil
        newbas := eminGbasis rest F
        cons(ecredPol( first(F), newbas),newbas)
 
             ------------------------------------------------
             --- does h belong to H
 
   ecrithinH(h: Dpol, H: List(Dpol))==
        null H  => true
        h1:= first(H)
        ecritM(degree h1, lc h1, degree h, lc h) => false
        ecrithinH(h, rest(H))
 
            -----------------------------
            --- calculate  euclidean S-polynomial of a critical pair
 
   esPol(p:critPair)==
      Tij := p.lcmfij
      fi := p.poli
      fj := p.polj
      lij:= lcm(leadingCoefficient(fi), leadingCoefficient(fj))
      red(fi)*monomial((lij exquo leadingCoefficient(fi))::Dom,
                        subtractIfCan(Tij, degree fi)::Expon) -
        red(fj)*monomial((lij exquo leadingCoefficient(fj))::Dom,
                         subtractIfCan(Tij, degree fj)::Expon)
 
            ----------------------------
 
            --- euclidean reduction mod F
 
   eRed(s: Dpol, H: List(Dpol), Hh: List(Dpol)) ==
     ( s = 0 or null H ) => s
     f1:= first(H)
     ds:= degree s
     lf1:= leadingCoefficient(f1)
     ls:= leadingCoefficient(s)
     e: Union(Expon, "failed")
     (((e:= subtractIfCan(ds, degree f1))  case "failed" ) or sizeLess?(ls, lf1) ) =>
        eRed(s, rest(H), Hh)
     sdf1:= divide(ls, lf1)
     q1:= sdf1.quotient
     sdf1.remainder = 0 =>
        eRed(red(s) - monomial(q1,e)*reductum(f1), Hh, Hh)
     eRed(s -(monomial(q1, e)*f1), rest(H), Hh)
 
            ----------------------------
 
            --- crit T  true, if e1 and e2 are disjoint
 
   ecritT(p: critPair) ==
          pi:= p.poli
          pj:= p.polj
          ci:= lc pi
          cj:= lc pj
          (p.lcmfij = degree pi + degree pj) and  (p.lcmcij = ci*cj)
 
            ----------------------------
 
            --- crit M - true, if lcm#2 multiple of lcm#1
 
   ecritM(e1: Expon, c1: Dom, e2: Expon, c2: Dom) ==
     en: Union(Expon, "failed")
     ((en:=subtractIfCan(e2, e1)) case "failed") or
       ((c2 exquo c1) case "failed") => false
     true
            ----------------------------
 
            --- crit B - true, if eik is a multiple of eh and eik not equal
            ---          lcm(eh,ei) and eik not equal lcm(eh,ek)
 
   ecritB(eh:Expon, ch: Dom, ei:Expon, ci: Dom, ek:Expon, ck: Dom) ==
       eik:= sup(ei, ek)
       cik:= lcm(ci, ck)
       ecritM(eh, ch, eik, cik) and
             not ecritM(eik, cik, sup(ei, eh), lcm(ci, ch)) and
                not ecritM(eik, cik, sup(ek, eh), lcm(ck, ch))
 
            -------------------------------
 
            --- reduce p1 mod lp
 
   euclideanNormalForm(p1: Dpol, lp: List(Dpol))==
       eRed(p1, lp, lp)
 
            ---------------------------------
 
            ---  insert element in sorted list
 
   sortin(p1: Dpol, lp: List(Dpol))==
      null lp => [p1]
      f1:= first(lp)
      elf1:= degree(f1)
      ep1:= degree(p1)
      ((elf1 < ep1) or ((elf1 = ep1) and
        sizeLess?(leadingCoefficient(f1),leadingCoefficient(p1)))) =>
         cons(f1,sortin(p1, rest(lp)))
      cons(p1,lp)
 
   updatD(D1: List(critPair), D2: List(critPair)) ==
      null D1 => D2
      null D2 => D1
      dl1:= first(D1)
      dl2:= first(D2)
      (dl1.lcmfij  <  dl2.lcmfij) => cons(dl1, updatD(D1.rest, D2))
      cons(dl2, updatD(D1, D2.rest))
 
            ----  calculate number of terms of polynomial
 
   lepol(p1:Dpol)==
      n: Integer
      n:= 0
      while p1 ~= 0 repeat
         n:= n + 1
         p1:= red(p1)
      n
 
            ----  print blanc lines
 
   prinb(n: Integer)==
        for i 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 =>
         print("  There is 1  Groebner Basis Polynomial "::Ex)
         prinb(2)
       print("  There are "::Ex)
       prinb(1)
       print(n::Ex)
       prinb(1)
       print("  Groebner Basis Polynomials. "::Ex)
       prinb(2)
 

