BeginPackage["Groebner`"]               (* Author: David S. Ng   5/90 *)

GTrace::usage =
"Type \"Groebner`GTrace[s__]:=Print[s]\" to trace this package."

Spoly::usage =
"Spoly[f1_, f2_] computes the S-polynomial of the polynomials f1 and f2."

NormalForm::usage =
"NormalForm[g_, F_, V_] computes the normal form of the polynomial g
modulo the set F of polynomials with variables in V."

ReduceForm::usage =
"ReduceForm[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."

Basis::usage =
"Basis[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 < ... ."

Begin["`private`"]

singleterm[f_] := Not[TrueQ[Head[f] == Plus]];
(* Return True if the polynomial f has only one term, False otherwise *)
(* If Head[f] == Plus, then there are more than one term *)

numterm[f_] := If [singleterm[f], 1, Length[f]];
(* Compute the number of terms of the polynomial f *)

nthterm[f_, n_] := If [singleterm[f], f, f[[n]]];
(* Compute the n-th leading term of the polynomial f *)

lterm[f_] := If [singleterm[f], f, f[[Length[f]]]];
(* Compute the leading term of the polynomial f *)

SetAttributes[lterm, Listable];

lcoeff[f_] := lterm[f] /. Map[Function[#->1],Variables[f]];
(* Compute the leading coefficient of the polynomial f *)

lpp[f_] := If [f === 0, 0, lterm[f] / lcoeff[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. *)

Spoly[f1_, f2_] :=
(* Compute the S-polynomial of polynomials f1 and f2 *)
Block[ {s1 = lterm[f1], s2 = lterm[f2], g},
    g = Expand[(s2 f1 - s1 f2) / (PolynomialGCD[s1,s2] lcoeff[f2])];
    Return[g]
]

NormalForm[g_, F_, V_] := 
(* Compute a normal form of the polynomial g modulo F in variable set V *)
Block[ {h = g, i, j, leadterm = lterm[F]},
  (* Start with the larger terms of h; have larger reduction. *)
    For [i = numterm[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}, nthterm[h, i]/leadterm[[j]]]; *)
(* Because of the bug in PolynomialQ, we use our own routine instead *)
            If [polynomialq[nthterm[h, i]/leadterm[[j]], V],
            (* If leading term divides a term of h, then we have reduction *)
                h = Expand[h - nthterm[h, i]/leadterm[[j]] F[[j]]];
                i = numterm[h] + 1; (* add 1 to offset decrement *)
                (* After 1 reduction, start over again. *)
                Break[];
            ];
        ];
    ];
    Return[h]
]

ReduceForm[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/lcoeff[h]]]],
                                              (* normalize to monic *)
        {i, Length[F]}
    ];
    Return[G]
]

Basis[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];
            GTrace[StringForm["Spoly[g``,g``] = ``", i, j, h]];
            If [h =!= 0, (* i.e. h not identically 0 *)
                h = Expand[h/lcoeff[h]]; (* normalize to monic *)
                AppendTo[G, h];
            ];
        ];
    ];
    GTrace["Almost done...Performing Final Reduction."];
    Return[ReduceForm[G, V]]
]

End[] (* `private` *)
EndPackage[] (* Groebner` *)
Print["Package Groebner` loaded"]