#==============================================================================
# IMPLEMENTATION IN GAUSS:  (procname_G)
# -----------------------------------------------------------------------------
# 1.) Implementation of Buchberger's algorithm for computing Grobner bases (GB)
# over a given domain S. GB(S,F) computes the reduced, mimimal Grobner basis
# for a list of polynomials F over the domain S. 
# Author : A. Pirklbauer July 1, 1992
# The implemented algorithm is :
#
#    procedure GB (F);	  			       #  F = {f1,...fn}
#      G := InterReduce(Q); k:= length(G);
#      B := {[i,j]: 1<=i<j<=k and Crit2(Gi,Gj)};
#      while B # {} do
#	 [i,j] := SelectPair (B,G);  B:= B - {[i,j]};
#	 if Crit1 ([i,j],B,G) then
#	   h := Reduce (SPoly(Gi,Gj),G);
#          if h # 0 then
#            G:= G + {h}; k:= k+1;
#            B:= B + {[i,k]: 1<=i<k and Crit2(Gi,Gk)};
#          endif;
#        endif;
#      endwhile;
#      redundant:= {h in G: R(h,G) # {h}};
#      return (InterReduce(G-redundant));
#    end GB;
#
#------------------------------------------------------------------------------

# infolevel[GB]:= 1:    # initialize global table infolevel. The user must 
                      # change this value to print more information. Example :
                      #   infolevel[GB] := 3;
                      #   GB(S,F);
                      # Used values are : 
                      #   1 => only result (final Grobner basis) is printed
                      #   2 => + statistics about reductions and time
                      #   3 => + each reduction step is printed
                      #   4 => + F after each reduction if new polynomials added
                      #   5 => + Now, also each S-Polynomial is printed
                      #   6 => + new pairs in B after each step (if added)

macro(ComputeReducers_G = `Gauss/GB/ComputeReducers`):
macro(SelectPoly_G = `Gauss/GB/SelectPoly`):
macro(Reduce_G = `Gauss/GB/Reduce_G`):
macro(ReduceAll_G = `Gauss/GB/ReduceAll_G`):
macro(NewBasis_G = `Gauss/GB/NewBasis_G`):
macro(InterReduce_G = `Gauss/GB/InterReduce_G`):
macro(Crit2_G = `Gauss/GB/Crit2_G`):
macro(Crit1_G = `Gauss/GB/Crit1_G`):
macro(TermLcm_G = `Gauss/GB/TermLcm_G`):
macro(SPoly_G = `Gauss/GB/SPoly_G`):
macro(SelectPair_G = `Gauss/GB/SelectPair_G`):
macro(InsertPair_G = `Gauss/GB/InsertPair_G`):

GB := proc (S,F)                 # S a domain, F a list/set of polynomials in S
  local k,B,memberB,HT,LCMHT,HC,G,G1,i,j,h,spoly,found,hc,ht,indB,lcmij,C,E,
        nec_reds, unnec_reds, total_time;

   if nargs <> 2 then ERROR(`bad number of arguments; call GB(S,F)`) fi;
   if not hasCategory(S,MultivariatePolynomial) then 
     ERROR(`Domain S must be MultivariatePolynomial`) 
   fi;
   C:= S[CoefficientRing]; E:= S[ExponentVector];
   if not hasCategory(C,GcdDomain) then
     ERROR(`CoefficientRing of S must be a GcdDomain`) 
   fi;
   G := [op(F)];                                              # convert to list
   if member(false, map(S[Type],G)) then 
     ERROR(`polynomials in F are not in domain S`)
   fi;

   nec_reds:=0; unnec_reds:=0; total_time:= time();           # for statistics
   B:= table(); memberB:= table(symmetric); LCMHT:= table();  # index set B
   userinfo(3,GB,print(`F start = `, map(S[Output],G)));
   userinfo(3,GB,print(`Interreducing initial basis ...`));
   G := InterReduce_G (S,G);            # reduce each element modulo the others
   userinfo(3,GB,print(`F after initial InterReduce = `, map(S[Output],G)));
   k := nops(G);                        # k:= number of polynomials in G;
   HC:=NULL; HT:=NULL;
   for i to k do                                  # HC[i] = head coeff of G[i]
     HC:=HC,S[Lcoeff](G[i],'ht');   HT:=HT,ht     # HT[i] = head term  of G[i]
   od;
   HC:=[HC]; HT:=[HT];
   
   # B := {[i,j]: 1<=i<j<=k and Crit2 (Gi,Gj)};
   for j from 2 to k do
     for i to (j-1) do
       if Crit2_G (E,HT,i,j) then InsertPair_G (E,i,j, B,memberB, LCMHT,HT) fi;
     od;
   od;
   userinfo(6,GB,print(`Initial pair set B = `, B));
   
   # while B <> {} do ... od;
   do
     indB:= map(op,[indices(B)]);                # indices: terms of Gi
     if indB = [] then break fi;                 # <=> while B <> {} do ... od;
     
     SelectPair_G (E,B,memberB,indB,'i','j');         # does also: B:=B-{[i,j]};
     lcmij:=LCMHT[i,j]; LCMHT[i,j]:= evaln(LCMHT[i,j]);

     if Crit1_G (E,HT,lcmij,i,j,memberB) then
       spoly:= SPoly_G(S,G[i],HC[i],HT[i],G[j],HC[j],HT[j],lcmij);   # S-Poly
       h := Reduce_G (S,spoly,G, HT,HC);                             # reduction
       userinfo(5,GB,print(`SPoly (`,i,j,`) = `, S[Output](spoly)));
       userinfo(3,GB,print(`Reduce(SPoly (`,i,j,`), G) = `, S[Output](h)));
       if (S[`<>`](h,S[0]) ) then                 # if (h <> 0) then ...
          
          # G:= G + {h}; k:= k+1;
          hc := S[Lcoeff](h,'ht');
          HC := [ op(HC), hc ];  HT := [ op(HT), ht ];
          G  := [ op(G), h ];  k:= k+1;
          userinfo(4,GB,print(`Basis F = `, map(S[Output],G)));
          
          # B:= B + {[i,k]: 1<=i<k and Crit2 (Gi,Gk)}; 
          for i to (k-1) do
            if Crit2_G(E,HT,i,k) then InsertPair_G(E,i,k,B,memberB, LCMHT,HT) fi
          od;
          userinfo(6,GB,print(`Pair set B = `, B));
         nec_reds := nec_reds + 1;
       else
         unnec_reds := unnec_reds + 1;
       fi;
     fi
   od;
   
   # redundant:= {h in G: R(h,G) # {h}};  G1:= G - redundant
   G1:=[];
   userinfo(3,GB,print(`Removing redundant elements ...`));
   for i to k do  found := false;
     for j to k while (not found) do
       found := (i<>j) and (E[Div](HT[i],HT[j])<>FAIL);   # if HT(j) | HT(i)
     od;
     if (not found) then G1:= [op(G1),G[i]] fi;
   od;

   userinfo(3,GB,print(`F before final InterReduce = `, map(S[Output],G1)));
   userinfo(3,GB,print(`Interreducing the basis ...`));

   G:= InterReduce_G (S,G1);             # reduce each element modulo the others
   G:= map(S[Primpart],G);
   total_time := time()-total_time;

   userinfo(3,GB,print(`F after final InterReduce = `, map(S[Output],G)));
   userinfo(2,GB,
     print(`Reductions :  necessary = `,nec_reds,` unnecessary = `,unnec_reds));
   userinfo(2,GB,print(`Total time : `,total_time));
   
   G;                                    # reduced, minimal Grobner basis
end:


#===========================================================================
#
# 2.) Full Reduction of a polynomial p modulo a set of polynomials Q
# Definition: the reducer set R(p,Q) := {q in Q such that HT(q) | HT(p) }
# The primitive part of the reduced polynomial is returned.
# Start with the whole polynomial. If no reducers exist, strip off leading
# monomial, otherwise, continue to reduce.
#
#    procedure Reduce(p,Q);                   # Q = {q1,q2,...,qn}
#      r:= p; s:= 0;
#      while r # 0 do
#        while R(r,Q) # {} do                 # while reducers (for r) exist :
#          q := SelectPoly(R(r,Q));           #   select a reducer from Q
#          r := r - (M(r)/M(q)) * q;          #   and reduce r
#        endwhile;
#        s:= s + M(r); r:= r - M(r);          # strip off leading monomial
#      endwhile;
#      return (s);
#    end Reduce;
#
#-------------------------------------------------------------------------------

ComputeReducers_G := proc (S,E,r,Q,htr,HT) local i,n,quo,reds;
  # precondition: HT[i] = headterm(Q[i]), htr := HT(r)
  if (S[`=`](r,S[0])) then RETURN ([]) fi;             # important ! (if r=0 ..)
  n := nops(Q); reds:= NULL;                           # reducer set R(r,Q)
  for i to n do
    quo := E[Div](htr,HT[i]);                          # exponent vector (term)
    if (quo <> FAIL) then                              # if HT(i) | HT(r)
      reds := reds, [i, quo] ;                         # quo = HT(r)/HT(i)
    fi;
  od;
  [reds];
end:

SelectPoly_G := proc(reducers,quo)      # precondition: reducers <> {}
  quo:= op(1, reducers)[2];
  op(1, reducers)[1];                   # take the first one (return index to Q)
end:

Reduce_G:= proc(S,poly,Q,HT,HC)    # precond.: HTi = hterm(Qi), HCi = hcoeff(Qi)
 local q,r,s,reds,quo,hcr,hcq,ht,htr,mr,HC1,HT1,i,k,E,C,scale,gcdqr,m1,m2;
  C:= S[CoefficientRing];  E:= S[ExponentVector];   # C: coefficients,  E: terms
  if (nargs = 3) then
    HC1:=NULL; HT1:=NULL; k:=nops(Q);
    for i to k do  HC1:=HC1,S[Lcoeff](Q[i],'ht'); HT1:=HT1,ht;  od;
    HC1:=[HC1]; HT1:=[HT1];
  else HC1:= HC; HT1:= HT; 
  fi;
  r:=poly; s:=S[0]; hcr:=S[Lcoeff](r,'htr');  # invariant H: hcr=HC(r),htr=HT(r)
  while (S[`<>`](r, S[0])) do  # H                    # while r <> 0
    reds := ComputeReducers_G (S,E,r,Q,htr,HT1);      # compute R(r,Q)
    scale := C[1];                                    # for rescaling s
    while reds <> [] do  # H
      q:= SelectPoly_G (reds,'quo'); hcq:= HC1[q];    # quo=HT(r)/HT(q), q=index
      
      #r:= r - (hcr/hcq)*quo*Q[q]);                   # reduce polynomial r :
      gcdqr := C[Gcd](hcq,hcr);                       # gcd=C[1] if C a field
      m1 := C[Div](hcq,gcdqr);                        # cofactor m1 = hcq/gcd 
      m2 := C[Div](hcr,gcdqr);                        # cofactor m2 = hcr/gcd 

      r := S[`.`](m1,r);                              # r := r * m1   (scale)
      scale := C[`*`](scale, m1);                     # update scalefactor
      k := S[`*`](S[List2Poly]([[m2,quo]]),Q[q]);     # k:= m2*(HT(r)/HT(q))*q
      r := S[`-`](r,k);                               # r:= r - k

      hcr := S[Lcoeff](r,'htr');                      # establish H: HC(r),HT(r)
      reds:= ComputeReducers_G (S,E,r,Q,htr,HT1);     # compute R(r,Q)
    od;
    if (S[`=`](r,S[0])) then break fi;                # IMPORTANT !  (if r=0 ..)
    mr := S[Lterm](r);                                # M(r) = HC(r)*HT(r)
    s := S[`.`](scale,s);                             # rescale :  s:=s*scale
    s := S[`+`](s,mr); r := S[`-`](r,mr);             # s:= s + mr;  r:= r - mr;
    hcr:= S[Lcoeff](r,'htr');                         # establish H: HC(r),HT(r)
  od;
  s:= S[Primpart](s);                                 # return primpart
  s;
end:


#===============================================================================
# 3.) InterReduce(F,X): Construction of a reduced ideal basis
#                  Ref: [1]: "Grobner Bases: An algorithmic method..."(Buchb85)
#                       [2]: Algorithm 10.3. in the Geddes textbook
#
#    procedure ReduceAll(F)          # remove any redundant elements
#      R:= F;  P:= {};
#      while R # {} do                     # P = {irreducible polynomials}
#        h:= an element of R;  R:= R-{h};  # R = {still reducible polynomials}
#        h:= Reduce (h,P);
#        if h # 0 then
#          P0 := {p in P: HT(h) | HT(p)};  # those which can be reduced by h
#          P := P - P0;                    # remove them from P = {irr. polys}
#          R := R + P0;                    # and add them to R = {red. polys}
#          P := P + {h};
#        fi;
#      end;
#      P;
#    end ReduceAll;
#  
#    procedure NewBasis (F)          # reduce each polynomial modulo the others
#      R := {};
#      for i from 1 to n do  f:= F[i];
#        h := Reduce (f, F - {f});         # reduce each F[i] modulo the others
#        if h <> 0 then R := R + {h} fi;
#      od;
#      R
#    end NewBasis;
#
#    procedure InterReduce (F)       # Construction of a reduced ideal basis
#      NewBasis (ReduceAll(F))
#    end InterReduce;
#
#-------------------------------------------------------------------------------

ReduceAll_G := proc(S,F)  local R,P,P0,h,i,HTh,HTp,E;
  R:= F;  P:= []; E := S[ExponentVector];
  while (R <> []) do
    h:= R[1]; R:= subsop(1=NULL,R);            # h:= an element of R; R:= R- {h}
    h:= Reduce_G (S,h,P);                      # h:= Reduce (h,P)
    if (S[`<>`](h,S[0]) ) then                 # if (h <> 0) then ...
      S[Lcoeff](h,'HTh');  P0:= NULL;
      for i to nops(P) do 
        S[Lcoeff](P[i],'HTp');
        if (E[Div](HTp,HTh) <> FAIL) then P0:= P0,i fi;     # if HT(h) | HT(p1)
      od;
      P0 := [P0];                              # P0 := {p in P: HT(h) | HT(p)};
      for i in P0 do R:=[op(R),P[i]] od;       # R := R + P0;
      for i in P0 do P:=subsop(i=NULL,P) od;   # P := P - P0;
      P := [op(P), h];                         # P := P + {h}
    fi;
  od;
  P;
end:

NewBasis_G := proc(S,F) local R,i,h,n,HT,HC,ht;
  n:= nops(F); HC:=NULL; HT:=NULL; R:=NULL;
  for i to n do HC:=HC,S[Lcoeff](F[i],'ht');  HT:= HT,ht  od;     # hcoeff/hterm
  HC:=[HC]; HT:=[HT];
  for i to n do
    h:=Reduce_G (S,F[i],subsop(i=NULL,F),subsop(i=NULL,HT),subsop(i=NULL,HC));
    if (S[`<>`](h,S[0]) ) then R:=R,h fi;    # if h<>0 then ...
  od;
  [R]
end:

InterReduce_G := proc (S,F) local R;
  R:= ReduceAll_G (S,F);
  NewBasis_G (S,R);
end:

#===============================================================================
# 4.) Auxiliary procedures for Buchbergers's Algorithm : 
#===============================================================================
   
Crit2_G := proc(E,HT,i,j)  local k,dim;      # returns (LCM(HTi,HTj) = HTi*HTj))
  dim := E[Dim];
  for k to dim do
   if (E[Degree](HT[i],k)<>0) and (E[Degree](HT[j],k)<>0) then RETURN (true) fi;
  od;
  false;           # i.e. returns (indets(HT[i]) intersect indets(HT[j]) <> {})
end:

Crit1_G := proc(E,HT,lcmij,i,j,memberB) local u,tmp;       # lcmij = LCMHT[i,j]
  # returns true if : ~E u: 1<=u<=k such that :
  #     {i<>u<>j:  [i,u],[u,j] not in B,  HT(u) | LCM ( HT(i),HT(j) )
  tmp:= true;
  for u from nops(HT) by -1 to 1 do
    if (u<>i and u<>j)        and                 #  i <> u <> j
       (memberB[i,u] <> true) and                 #  [i,u] not in B
       (memberB[u,j] <> true) and                 #  [u,j] not in B
       (E[Div](lcmij,HT[u])<>FAIL )               #  HT(u) | LCM ( HT(i),HT(j) )
    then
      tmp:= false;
      break;
    fi;
  od;
  tmp;
end:

TermLcm_G := proc(E,a,b) local k,dim,s;    # a,b are terms (exponent vectors)
  dim := E[Dim]; s:=NULL;
  for k to dim do
    s:= s, max(E[Degree](a,k),E[Degree](b,k));     # sequence of max's
  od;                                              # s a sequence => [s] a list
  E[List2Vect]([s]);                               # convert to exponent vector
end:

#-------------------------------------------------------------------------------
# SPoly, where HC's, HT's and the LCM(HT's) of the polynomials are known.
#-------------------------------------------------------------------------------

SPoly_G := proc (S,poly1,hc1,ht1, poly2,hc2,ht2, lcmij)
  local m1,m2,f1,f2,s1,s2,r,gcd12,C,E;
  # precondition: hc/ht(i) = HC/HT(poly(i));  lcmij = LCM (hc1,hc2)
  C:= S[CoefficientRing]; E:= S[ExponentVector];
  
  gcd12 := C[Gcd](hc1,hc2);                 # is always C[1] if S is a field
  m2 := C[Div](hc1,gcd12);                  # cofactor m2 :  m2 * gcd = hc1
  m1 := C[Div](hc2,gcd12);                  # cofactor m1 :  m1 * gcd = hc2
  
  f1 := E[Div](lcmij,ht1);                  # multiply HT(poly1) with f1  (term)
  f2 := E[Div](lcmij,ht2);                  # multiply HT(poly2) with f2  (term)
  
  s1 := S[`*`](S[List2Poly]([[m1,f1]]), poly1);      # s1 := f1*m1*poly1
  s2 := S[`*`](S[List2Poly]([[m2,f2]]), poly2);      # s2 := f2*m2*poly2
  r  := S[`-`](s1,s2);                               # r  := s1 - s2
  r;
end:


#------------------------------------------------------------------------------
#  Procedures operating on the index set B/memberB : SelectPair, InsertPair
#------------------------------------------------------------------------------

SelectPair_G := proc(E,B,memberB,indB,i,j) local i0,j0, minI,pair;  # B # {}

  # Selects a pair [i,j] from the set B such that
  #   LCM(HT(Gi),HT(Gj)) = min { LCM(HT(Gu),HT(Gv)): [u,v] in B  },
  # where the minimum is with respect to the ordering of terms (which is
  # defined in E) : "Normal Selection Strategy" (NSS).
  # Additionally, the pair is removed from B and memberB.
  
  minI := E[Min](op(indB));                # Indices of table B (indB) reflect
  pair := op(1,B[minI]);                   # ordering of the LCMHT's.
  i0:=pair[1]; j0:=pair[2]; i:=i0;j:=j0;   # Choosing the 'Min' implements NSS
  # remove the pair from B and memberB
  if nops(B[minI]) = 1 then B[minI]:= evaln(B[minI])  # last => remove entry
  else B[minI] := subsop(1=NULL,B[minI]);             # else take one element
  fi;
  memberB[i0,j0]:= evaln(memberB[i0,j0]);  
end:

InsertPair_G := proc(E,i,j, B,memberB, LCMHT,HT) local k;     # B:= B + [i,j]
 
  # Inserts a new pair [i,j into the sets B. The index of this table is the
  # LCM(HT(Gi), HT(Gj)) which lies in E. Each entry in B is a list of pairs
  # with the same LCM. This data structure is used in procedure 'SelectPair'
  # to implement the "Normal Selection Strategy" (see there).
  # Additionally the tables memberB and LCMHT are updated.
  
  memberB[i,j]:=true;
  LCMHT[i,j]:= TermLcm_G (E,HT[i],HT[j]);                # compute lcm(HTi,HTj)
  k:= LCMHT[i,j];                   # k:= index of table B (is exponent vector)
  if assigned(B[k]) then B[k]:= [op(B[k]), [i,j]]     # insert the pair [i,j]
  else B[k] := [ [i,j] ];
  fi;
end:

save `GB.m`;
done;