BeginPackage["Groebner`"]               (* Author: David S. Ng   	5/90 *)   
  										(* modified by: K.-P. Ndf	6/96 *)	   
  gbTrace::usage =   
  	"Type \"Groebner`gbTrace[s__]:=Print[s]\" to trace this package."   
  	(* sollte besser als Option eingebaut werden *)   
    
  sPoly::usage =   
  	"sPoly[f1_, f2_] computes the S-polynomial of the polynomials f1 and f2."   
    
  monicForm::usage =  
  	"monicForm[F_,V_] divides any polynomial in F by it's leading coefficient."  
     
  normalForm::usage =   
  	"normalForm[g_, F_, V_] computes the normal form of the polynomial g   
  modulo the set F of polynomials with variables in V."   
    
  reducedForm::usage =   
  	"reducedForm[F_, V_] computes the reduced form of the set F of   
  polynomials with variables in V, i.e. every f in F is   
  in normal form modulo F - {f} and is normalized to monic."   
      
  stairForm::usage =   
  	" (* T = *) stairForm[F_, V_] computes the recursively reduced form of the set F of   
  polynomials with variables in V, i.e. every f in T is finaly  
  in normal form modulo T - {f} and is normalized to monic."   
     
  reducedGBasis::usage =   
  	"reducedGBasis[F_, V_] computes the Reduced Groebner Basis of    
  the set F of polynomials with variables in V using Mathematica's    
  built-in lexicographic ordering, i.e.   
  1 < x < x^2 < ... < y < x y < x^2 y < ... < y^2 < x y^2 < ... < z < ... ."   
     
  monomialQ::usage =    
  	"monomialQ[f] gives True is f is a monomial    
  and False otherwise"   
     
  leadingTerm::usage=   
  	"leadingTerm[f] gives the leading term of the polynomial f"   
     
  Begin["`Private`"]   
     
  singleMonomialQ[f_] := Not[TrueQ[Head[f] == Plus]];   
  (* Return True if the polynomial f has only one monomial,   
  False otherwise *)   
  (* If Head[f] == Plus, then there are more than one term *)   
     
  nbOfMonomials[f_] := If [singleMonomialQ[f], 1, Length[f]];   
  (* Compute the number of terms of the polynomial f *)   
     
  nthMonomial[f_, n_] := If [singleMonomialQ[f], f, f[[n]]];   
  (* Compute the n-th monomial of the polynomial f *)   
     
  leadingMonomial[f_] := If [singleMonomialQ[f], f, f[[Length[f]]]];   
  (* Compute the leading monomial of the polynomial f *)   
     
  SetAttributes[leadingMonomial, Listable];   
     
  leadingCoefficient[f_] :=    
  	leadingMonomial[f] /. Map[Function[#->1],Variables[f]];   
  (* Compute the leading coefficient of the polynomial f *)   
     
  leadingTerm[f_] := If [f === 0, 0,    
  	leadingMonomial[f] / leadingCoefficient[f]];   
  (* Compute the leading power product of the polynomial f *)   
     
  polynomialQ[f_, V_] :=   
  Intersection[Variables[Denominator[ExpandDenominator[Together[f]]]], V] == {};   
  (* Determines if the given rational function f is actually   
     a polynomial with variables in V. *)   
    
  monicForm[F_,V_] := Expand[#/leadingCoefficient[#]]& /@ F;   
  (* normalize to monic *)   
          
  sPoly[f1_, f2_] :=   
  (* Compute the S-polynomial of polynomials f1 and f2 *)   
  Block[ {s1 = leadingMonomial[f1], s2 = leadingMonomial[f2], g},   
      g = Expand[(s2 f1 - s1 f2) /    
      	(PolynomialGCD[s1,s2] leadingCoefficient[f2])];   
      Return[g]   
  ];   
  (* (1) sPoly is here not symmetric *)   
  (* (2) Trace the evaluation and test if sPoly can be improved    
  		according rational number operations *)   
     
  normalForm[g_, F_, V_] :=    
  (* Compute a normal form of the polynomial g modulo F in variable set V *)   
  Block[ {h = g, i, j, headMonomial = leadingMonomial[F]},   
    (* Start with the larger terms of h; have larger reduction. *)   
      For [i = nbOfMonomials[h], i > 0, i--,   
          If [h === 0, Break[]];   
          For [j = 1, j <= Length[F], j++,   
          (* Try every polynomial f until reduction is found *)   
  (* Print["     ", {i,j}, nthMonomial[h, i]/headMonomial[[j]]]; *)   
  (* Because of the bug in PolynomialQ, we use our own routine instead *)   
   (* which is also much quicker *)   
              If [polynomialQ[nthMonomial[h, i]/headMonomial[[j]], V],   
              (* If leading term divides a term of h, then we have reduction *)   
                  h = Expand[h - nthMonomial[h, i]/headMonomial[[j]] F[[j]]];   
                  i = nbOfMonomials[h] + 1; (* add 1 to offset decrement *)   
                  (* After 1 reduction, start over again. *)   
  				(* "Could this not be improved?" *)   
                  Break[];   
              ];   
          ];   
      ];   
      Return[h]   
  ];   
     
  reducedForm[F_, V_] :=   
  (* Compute the reduced form for F with variable set V: i.e.   
     Every f in F is in normal form modulo F - {f} and is normalized to monic. *)   
  Block[ {G = {}, h, i},   
      Do [   
          h = normalForm[F[[i]], Drop[F, {i, i}], V];   
          If [h =!= 0, AppendTo[G, Expand[h/leadingCoefficient[h]]]],   
                                                (* normalize to monic *)   
          {i, Length[F]}   
      ];   
      Return[G]   
  ];   
     
  stairForm[F_, V_] :=   
  (* Replaces recursively the reduced form for F with variable set V: i.e.   
     Every f in F is in normal form modulo F - {f} and is normalized to monic *)   
  Block[ {H = {}, G = Expand[F], h, i},   
      While[ H =!= G,   
  		G = Union[G];   
  		If[G[[1]] === 0, G = Drop[G, 1]];  
  			Print[G];  
  		(* F = G in the next row did not work! *)  
  		G = monicForm[G,V];  
  		H = G;  
          Do[G[[i]] = normalForm[H[[i]], Drop[H, {i, i}], V];,   
          {i, Length[H]}  
  		];    
      ];   
      Return[G]   
  ];   
    
  reducedGBasis[F_, V_] :=   
  (* Compute the Reduced Groebner Basis of the set of polynomials F with   
     variables in V using Mathematica's built-in lexicographic ordering, i.e.   
     1 < x < x^2 < ... < y < x y < x^2 y < ... < y^2 < x y^2 < ... < z < ... .   
     Note that no parameter is allowed in the coefficients of the polynomials. *)   
  Block[ {G, B, h, i, j},   
      G = Expand[F];   
      For [i = 2, i <= Length[G], i++,   
          For [j = 1, j < i, j++,   
              h = normalForm[sPoly[G[[i]], G[[j]]], G, V];   
              gbTrace[StringForm["Spoly[g``,g``] = ``", i, j, h]];   
              If [h =!= 0, (* i.e. h not identically 0 *)   
                  h = Expand[h/leadingCoefficient[h]]; (* normalize to monic *)   
                  AppendTo[G, h];   
              ];   
          ];   
      ];   
      gbTrace["Almost done...Performing Final Reduction."];   
      Return[ReducedForm[G, V]]   
  ]   
     
  End[] (* `Private` *)   
  EndPackage[] (* Groebner` *)   
  Print["Package Groebner` loaded"]