Practical computations with Gröbner bases

Daniel Lichtblau
Copyright ©  2000

Abstract

Since their invention by Buchberger in 1965, Gröbner bases have become a pervasive tool in computational mathematics. In this paper we show several sorts of Gröbner basis computations in Mathematica. We discuss computations of approximate Gröbner bases, with applications to solving systems and implicitization. We also show how to use the built-in GroebnerBasis function to compute polynomial greatest common divisors over various fields, series inverses, syzygy modules, and more. Some of this is by now classical; the emphasis is on effective use of off-the-shelf technology, although possibly some of the methods are new. We will point out areas where existing functionality needs enhancement, and note some directions for future work.

Introduction

Invented by Buchberger (in his 1965 thesis) to solve the ideal membership problem algorithmically, Gröbner bases have come to be used throughout computational mathematics, especially in the realm of symbolic algebra. For a large set of examples, one might look at the survey article [Buchberger 1985]. Indeed, that paper provides an excellent cross-section of the wealth of Gröbner basis techniques already in use at that time, as well as providing abundant references to the early history of the technology. A related reference is [Buchberger 1988] which covers nicely the applications of Gröbner bases to computational geometry that were known at that time.  In this paper we will consider several applications of interest, old and new, that are motivated by quite practical concerns. We will, for example, see some ways in which use of Gröbner bases may enhance speed. Another general consideration will be ease of programming in situations where the methods presented might not be the fastest available.

Specific topics discussed below include:
(i) Numeric Gröbner bases.
(ii) Basis conversion from one term order to another.
(iii) Implicitization of parametrized varieties.
(iv) Working with primitive elements.
(v) Finding partial reciprocals of power series.
(vi) Computing polynomial greatest common divisors.
(vii) Computing Gröbner bases over modules.
(vii) Computing conversion matrices.
(ix) Computating syzygy modules.
(x) Computating univariate polynomial matrix Hermite normal forms.
(xi) Finding dimensions of certain vector spaces and ideals.
Several related topics are also presented.

We will assume familiarity with the basic notions and notations of Gröbner bases. The reader may find relevant background material in [Buchberger 1985; Cox, Little, and O'Shea 1992; Becker, Weispfenning, and Kredel 1993; Adams and Loustaunau 1994]. Emphasis below will be placed on informal explanations for the various algorithms and ideas presented, rather than formal proofs, as by and large the theory is well understood but the applications are not universally known.

1.	Numerical Gröbner bases and some methods for basis conversion

An example

We begin with the following set of equations, taken from ISSAC '97 challenge problem #6 [Trott, 1997]. All timings in this paper were taken from tests run in version 4 of Mathematica [Wolfram 1999]. The machine has a 300Mhz Pentium chip with 128 Mb of RAM, and runs under the Linux operating system.

polys={8 (w^2)+5 w x-4 w y+2 w z+3 w+5 (x^2)+2 x y-7 x z-7 x+7 (y^2)-8 y z-7 y+7 (z^2)-8 z+8,3 (w^2)-5 w x-3 w y-6 w z+9 w+4 (x^2)+2 x y-2 x z+7 x+9 (y^2)+6 y z+5 y+7 (z^2)+7 z+5,-2 (w^2)+9 w x+9 w y-7 w z-4 w+8 (x^2)+9 x y-3 x z+8 x+6 (y^2)-7 y z+4 y-6 (z^2)+8 z+2,7 (w^2)+5 w x+3 w y-5 w z-5 w+2 (x^2)+9 x y-7 x z+4 x-4 (y^2)-5 y z+6 y-4 (z^2)-9 z+2};vars={x,y,z,w};

We will find a lexicographic Gröbner basis. This is relatively slow. Soon we will see some short cuts.

Timing[gblex = GroebnerBasis[polys, vars];]

Let us find the shape of the Gröbner basis. The routine below is a crude indicator. Better would be something that returns all maximal power products. This is not too hard to code, but is more than our present purposes require.

shape[gb_, vars_] := Map[Exponent[#,vars]&, gb]

shape[gblex, vars]

The shape is as follows. We have a univariate polynomial of degree 16 in w, a linear polynomial in z with constant term a degree 15 polynomial in w, and similarly for y and x.

Using a lexicographic Gröbner basis to solve a system of polynomials.

Let us use this Gröbner basis to find the roots of the original set of polynomials. The idea is simple (see e.g. [Cox, Little, and O'Shea 1992, chapter 2, section 8]). We already saw that the basis is in effect triangulated, that is, we have a polynomial in the last variable, a polynomial in the last two variables, etc. (similar to a triangulated matrix). We will find the roots of the univariate polynomial, then back-substitute these into the three linear polynomials. Were the structure of these polynomials not linear in the variables of interest, it might not be so simple, though the basic idea would be the same. Note that Solve will handle the root extraction just fine. We reverse the order of the variables to show that Solve will automatically choose a sensible reordering for this example (we know this because it does not spend much time recomputing a new Gröbner basis).

Timing[soln = Solve[gblex==0, Reverse[vars]];]

Let us check the solution. It turns out that we require Mathematica's extended precision (often called "significance arithmetic") because  machine arithmetic will not suffice. It is also the case that numerical evaluation near zero is slow, and for a very good reason. (Question: How can you tell when you have a result of the desired precision? Answer: You cannot, because zero has no precision. Internally Mathematica must eventually switch from precision to accuracy, a compromise that is necessary to keep programs running in finite time). We gain a bit of speed by the trick of adding 1 to each polynomial, numerically evaluating, then subtracting it again. We see that every component, after chopping, is indeed zero. So our solution set is correct.

Off[$MaxExtraPrecision::meprec]

Timing[Max[Abs[N[polys + 1 /. soln, 100] - 1 // Chop]]]

Approximate Gröbner bases

In Mathematica we can find an approximate Gröbner basis. This is often much faster because the coefficients, unlike exact integers, do not grow in the number of bits needed to represent them. Below we will find a Gröbner basis using a starting precision of 250.

Timing[gbapprox =
	GroebnerBasis[polys, vars, CoefficientDomain->InexactNumbers[250]];]

So this is much faster. Now let us use this instead of the exact gblex to get approximate roots of polys.

Timing[numsoln = Solve[gbapprox==0, vars];]
Timing[Max[Abs[N[polys + 1 /. numsoln, 100] - 1 // Chop]]]

This took a bit longer in Solve for reasons related to the use of bignums and significance (rather than machine) arithmetic, but overall time including the basis computation is clearly much less. Actually by using machine arithmetic we can get solutions an order of magnitude faster, albeit with considerably higher residual errors. Note that these residuals of course can be made arbitrarily large or small simply by rescaling the input polynomials. As the coefficients are all between one and ten in absolute value the residuals below are within a factor of ten of those we would have by normalizing so that the largest coefficients are 1.

Timing[numsoln = Solve[N[gbapprox]==0, vars];]
Timing[Max[Abs[N[polys + 1 /. numsoln, 100] - 1 // Chop]]]

The fact that these often can be computed relatively quickly a hint that approximate Gröbner bases might be quite useful in technology, especially in finding approximate roots to zero-dimensional varieties. Actually one might prefer to compute using a term order that is less prone to intermediate swell and huge results. Typically the total degree reverse-lexicographic order is suitable. One can compute an approximate Gröbner basis in this ordering, then find endomorphism maps in the appropriate finite dimensional vector space, compute eigenvalues, and from there find roots to the original polynomial system. This method, spelled out in much more detail e.g. in [Cox 1998, section 6; Lichtblau 2000a] and references therein, is the core of the numeric algebraic equation solver  in Mathematica. There is some discussion of the use of significance arithmetic in computation of Gröbner bases in [Lichtblau 2000a]. A recent discussion of significance arithmetic in general is in [Sofroniou and Spaletta 2000].

Simple basis conversion

We start this section with the statement that lexicographic (lex) Gröbner bases are often of interest but hard to compute. It is often considerably faster to compute a basis with respect to a different term order, for example degree-reverse-lexicographic (degrevlex). For many ideals, in particular the one above, there is a simple method to convert to a lexicographic basis. After computing a basis in degrevlex order we can use linear algebra. It is in essence a simple form of the method in [Faugére,  Gianni,  Lazard, and Mora 1993], one that works whenever the ideal satisfies the hypotheses of the shape lemma [Becker, Marinari, Mora, and Traverso 1994]. Specifically, we assume the lex basis has the form

{univariate polynomial of known degree d in last variable,
 polynomial linear in next-to-last variable with constant term is of degree at most (d-1)  in last variable,
 polynomial linear in third-from-last variable with constant term of degree at most (d-1) degree in last variable, ...}

To begin the process, we form a polynomial in the last variable, of degree d, having new variables to represent undetermined coefficients.  Specifically, we will write (Subscript[a, 0])+(Subscript[a, 1])*w + ... + (Subscript[a, deg])*(w^d). We want this polynomial to reduce to zero, and this will give us linear equations in the coefficients. To get these equations we first call PolynomialReduce to reduce this with respect to the basis we have. The basis has only numeric coefficients, so each term in the result is a power product in {x,y,z,w} with a coefficient that is still linear in the new variables. Setting these coefficients to zero gives a linear system in the new variables, we solve it to get the values for these coefficients that give us a univariate polynomial in w. This goes into our lex basis under construction. Similar method for finding the remaining polynomials in the lexicographic basis.

One sees from gbapprox that we have a basis whose shape is amenable to this simple conversion method. What happens when this is not so? More on this later. For now note that we need to get the degree (easy) and to decide which is our last variable (easier). And we still need a degrevlex basis.

var = Last[vars];
deg = Exponent[First[gbapprox], var];
len = Length[vars];
Timing[gbdrlex = GroebnerBasis[polys, vars,
	MonomialOrder->DegreeReverseLexicographic];]

It was easy to get a basis in degrevlex order. We use the definitions below and the basis conversion is now significantly faster than a direct lex basis computation. One sees that the code is a bit more complicated than that of most other examples herein. The reason is that wee do some extra work in order to take optimal advantage of LinearSolve.

Off[General::spell];
Off[General::spell1];

adjustCoefficients[poly_Plus] := Module[
	{lis=Apply[List,poly] /. a_?NumberQ*b_ :> a, gcd},
	lis = Select[lis, NumberQ];
	gcd = Apply[GCD, lis];
	Expand[poly/gcd]
	]

basisConvert[vars_,var_,basis_,deg_]:=Module[{newvars,nv,newpoly1,np1,fpoly,rhsterms,result},newvars=Array[nv,deg,0];
newpoly1=Table[var^j,{j,0,deg-1}];
np1=Map[Last,PolynomialReduce[Append[newpoly1,var^deg],basis,vars,MonomialOrder->DegreeReverseLexicographic]];
fpoly=Last[np1];
np1=newvars.Drop[np1,-1];
newpoly1=newvars.newpoly1;
rhsterms=Reverse[Drop[vars,-1]];
result=findPolynomials[newpoly1,Prepend[rhsterms,var^deg],np1,Prepend[rhsterms,fpoly],vars,basis,newvars];
Map[adjustCoefficients,result]]

findPolynomials[origpoly_,origleads_,poly_,leads_,vars_,gb_,lvars_]:=Module[{jac,mults,rhs,gcd,denoms,coeffs,newpoly,dvars,extra,dd},dvars=Array[dd,Length[leads]];
extra=dvars.leads;
monoms=Map[Last,First[Internal`DistributedTermsList[poly+extra,vars,MonomialOrder->DegreeReverseLexicographic]]];
jac=Outer[D,monoms,lvars];
rhs=-Outer[D,monoms,dvars];
denoms=Denominator[jac];
mults=Map[Apply[LCM,#]&,denoms];
jac*=mults;
rhs*=mults;
gcd=Apply[GCD,Flatten[jac]];
jac/=gcd;
rhs/=gcd;
jac=Transpose[jac];
mults=Map[Apply[GCD,#]&,jac];
jac/=mults;
jac=Transpose[jac];
coeffs=LinearSolve[jac,rhs];
coeffs/=mults;
(origpoly/.Thread[lvars->coeffs])+origleads]

Timing[lexgb1 = basisConvert[vars, var, gbdrlex, deg];]

Clearly this is a significant improvement in speed over the direct computation. We check that the result is in agreement with lexgb. First they must have the same number of elements.

Length[lexgb1]==Length[gblex]

Now we divide corresponding elements, call Together on the result, and check that they are the same up to multiplicative units. Since we normalize so that the polynomial coefficients do not have a common denominator, the units must in fact be ±1.

Map[Together, lexgb1/gblex]

The nice thing is that often we can guess the correct shape from the degrevlex basis. We assume it has the shape above, and the degree deg is simply the size of the normal set (this is the set of power products that is not reducible by the Gröbner basis). This size is invariant for the ideal, that is, independent of the specific basis. This is because the vector space ℚ[variables]/ideal is independent of any Gröbner basis for the ideal, and this normal set is a basis for that vector space, so has size independent of the Gröbner basis. How does one find this size given a non-lexicographic basis? That is left as an exercise for the reader. Of importance is that any ideal that is zero dimensional (as an ideal), radical and in general position with respect to the last variable will have this sort of shape, hence will be amenable to this simple variant of FGLM basis conversion.

Suppose our basis did not have the shape needed for the code above to obtain a lexicographic basis? In thet case one might  try the general version of the FGLM method, or else another conversion method such as the Gröbner walk [Collart, Kalkbrenner, and Mall 1997]. Moreover, provided the ideal is radical and zero-dimensional, a generic linear change of coordinates suffices to put it into the appropriate shape [Gianni, Trager, and Zacharias 1988].

Another way to obtain an exact Gröbner basis

We start with an inexact basis computed to high precision, then rationalize it to obtain an exact one.  By trial-and-error one may learn that for the example above, an initial precision of 1000 digits is sufficient whereas 800 is not. Including the check for correctness, this process takes under 15 seconds. To make this into a general method, one would increment the precision until the Rationalize result stabilized, then try the check. This would of course add substantially to the run time.

The code below uses the first and last elements in our internal ordering, because at least one of these is not the highest in the term order (this one being normalized to have a coefficient of inexact 1). We typically get a good approximation of the common denominator, up to a small factor, by looking at any other coefficient, then we use it to clear denominators up to the small remaining factor. What remains will be much faster to rationalize.

fastRationalize[polys_, vars_] := Module[
	{lmonoms, coeffs, denoms},
	lmonoms = Map[{First[#],Last[#]}&, polys];
	coeffs = lmonoms /. Thread[vars->1];
	denoms = Denominator[Rationalize[coeffs]];
	denoms = Map[Apply[LCM,#]&, denoms];
	Rationalize[Expand[denoms*polys]]
	]

Timing[gbapprox1000 =
	GroebnerBasis[polys, vars,
	CoefficientDomain->InexactNumbers[1000]];]
Timing[gb1000 = fastRationalize[gbapprox1000, vars];]
Timing[lexgb2 = Map[adjustCoefficients, gb1000];]

We now check that this is correct. We do so by verifying that each polynomial is in the ideal. For this purpose it suffices to reduce by a known Gröbner basis, and so we use gbdrlex. Actually we should also reduce gbdrlex by this putative basis in order to get full confirmation, but provided we have independent knowledge that the size of the normal set is 16 we can actually dispense with this step.

Timing[Map[(PolynomialReduce[#, gbdrlex, vars, (* want all 0's *)
	MonomialOrder->DegreeReverseLexicographic] [[2]])&, lexgb2]]
Map[Together, lexgb2/lexgb1]

Minimal polynomials via Gröbner bases

Often one can use Gröbner technology to simplify nontrivial algebraic numbers. For example, the value below is the surface area of the 11-gonal antiprism with height = side length = 1 [Weisstein 2000].

s=(11/2)*(Cot[Pi/11]+Sqrt[3]);

Now suppose we wish to obtain the minimal polynomial for s. Converting to exponentials, we can express Cot[Pi/11] as a rational function in terms of a primitive 22th root of -1. From basic cyclotomic polynomial theory we know this is of degree 10 over the rationals [Dummit and Foote  1991]. As Sqrt[3] is of degree 2 over the rationals. Hence s must satisfy a polynomial of degree no larger than 20. We will use a technique from [Buchberger 1985, section 6] to find that polynomial. In essence this is a precursor to the FGLM basis conversion algorithm. Other nice descriptions of the method may be found in [Becker, Weispfenning, and Kredel 1993, chapter 7; Adams and Loustaunau 1994, chapter 2, section 6]. In brief, we replace radicals by variables, set up polynomial identities for those variables, and then rewrite s taking into account all algebraic dependencies. We set this up as below by (i) explicitly rewriting in terms of radicals,  (ii) replacing the radicals with new variables, (iii) forming polynomials to relate those variables to rational values, and (iv) forming a polynomial expression for our result. We use the fact that Sqrt[-1] is the 11th power of one of our radicals (otherwise we will have an algebraic dependency over the rationals not taken into account).

s1=Together[TrigToExp[s]];
s2=s1/.{(-1)^Rational[p_,q_]->x1^(22*p/q),Sqrt[3]->x2, Complex[re_,im_]->re+im*x1^11};
polys={x1^22+1,x2^2-3,Numerator[s2]-result*Denominator[s2]}

Now we compute a Gröbner basis with respect to an efficient term ordering. We will use this to find a univariate polynomial for result via linear algebra.

vars={x1,x2,result};
Timing[gbdrlex=GroebnerBasis[polys,vars,MonomialOrder->DegreeReverseLexicographic];]

We now create a polynomial in result, of degree 20, with undetermined coefficients for all but the degree 0 term. As we know it must be nonzero for a minimal, hence irreducible polynomial, we may force it to be 1 and thus avoid an underdetermined system of linear equations. We reduce this polynomial with respect to the Gröbner basis just formed, then get the coefficients of the power products (these are in terms of our new variables).

deg=20;
newvars=Array[nv,deg];
newpoly=newvars.Table[result^j,{j,deg}]+1;
Timing[np1=Last[PolynomialReduce[newpoly1,gbdrlex,vars,MonomialOrder->DegreeReverseLexicographic]];]
coeffpolys = Map[Last,First[Internal`DistributedTermsList[np1,vars,MonomialOrder->DegreeReverseLexicographic]]];

Last we set these coefficients all to zero, solve for the indeterminates, use these to reform our univariate polynomial, and clear the denominators.

sol = First[Solve[coeffpolys==0, newvars]];
minpoly =adjustCoefficients[newpoly /. sol]

We check numerically that we obtained the correct polynomial.

Last[NRoots[minpoly==0, result]][[2]]==N[s]

Note that a direct computation either via the Mathematica built-in function RootReduce, or with GroebnerBasis using MonomialOrder->EliminationOrder, will not complete in reasonable time for this example. Morover one cannot take advantage of the special-case FGLM code provided above because this ideal is not in general position with respect to the variable of interest; there are 40 solutions but only 20 distinct ones for result. Alternatives include using a Gröbner walk (see for example [Amrheim, Gloor, and Küchlin 1997] for a description of how one obtains a fast walk) or a faster direct computation such as that discussed in [Faugére 1999]. An advantage to the method presented above, however, is that one need not compute a full lexicographic basis, but only obtain one polynomial therefrom after first finding a basis with respect to a (typically) cheaper term ordering.

Working with primitive elements in algebraic extensions of ℚ

Often when manipulating expressions that involve algebraic numbers it may be convenient to express them in terms of a primitive element [Becker, Weispfenning, and Kredel 1993, chapter 8; Dummit and Foote 1991 chapter 14; von zur Gathen and Gerhard 1999 chapter 22]. As is noted in that first reference, for example, one may express a primitive element as a generic linear combination (over the rationals) of the given algebraic numbers. For example, suppose (Subscript[x, 1]) and (Subscript[x, 2]) are distinct roots of the polynomial 9 - 8 x + 6(x^2)+ (x^4), and we wish to show that

expr=-8+6 (Subscript[x, 1])+((Subscript[x, 1])^3)+6 (Subscript[x, 2]) +((Subscript[x, 1])^2) (Subscript[x, 2])+(Subscript[x, 1]) ((Subscript[x, 2])^2)+((Subscript[x, 2])^3);

is zero. One tactic might be to evaluate the expression numerically but then one must find an a priori minimum precision required in order to determine that it vanishes. We give a different method below. First we create appropriate polynomials in (Subscript[x, 1]) and (Subscript[x, 2]).

poly1=(9-8*#1+6*#1^2+#1^4)&[(Subscript[x, 1])]
poly2=(9-8*#1+6*#1^2+#1^4)&[(Subscript[x, 2])]

Next we form an expression that says, in effect, that (Subscript[x, 1]) and (Subscript[x, 2]) are distinct (this is necessary because they have the same minimal polynomials). We do this by forcing the difference to be nonvanishing. For this we will require yet another new variable, to denote the reciprocal of that difference.

polyrecip=recip*((Subscript[x, 1])-(Subscript[x, 2]))-1;

Now we create a new polynomial in a new variable, (Subscript[x, 3]), which will be our primitive element. Assuming our choice of coefficients is generic (it is), this polynomial defines (Subscript[x, 3])as a primitive element for ℚ[(Subscript[x, 1]),(Subscript[x, 2])].

poly3 = (Subscript[x, 3]) - ((Subscript[x, 1])+5(Subscript[x, 2]))

We are now ready to form a Gröbner basis. We will eliminate the auxiliary variable recip, and use an ordering that makes the original variables higher than the primitive element.

Timing[gb=GroebnerBasis[{poly1,poly2,poly3,polyrecip},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])},recip];]

Now we may rewrite our expression in terms of (Subscript[x, 3]) alone. Clearly the result is explicitly zero if and only if the original expression was a representation of zero.

result = Last[PolynomialReduce[expr,gb,{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}]]

Implicitization and basis approximation.

One can use Gröbner bases to find the implicitization of a variety given parametrically. For example, suppose we have the surface, given parametrically by s and t, by the equations

x==((s^2)-1-(t^2))/((s^2)+(t^2)+1);
y==2 s/((s^2)+(t^2)+1);
z==2 s t/((s^2)+(t^2)+1);

We can clear denominators and add an auxiliary equation to prevent the denominator from vanishing. We then form a Gröbner basis with an ordering that eliminates the parameters and the auxiliary variable. This technique is developed in [Kalkbrenner 1990]. Detailed explanations may also be found in [Cox, Little, and O'Shea 1992, chapter 3, section 3; Becker, Weispfenning, and Kredel 1993, chapter 7, section 6; Adams and Loustaunau 1994, chapter 2, section 5 ].

den=(s^2)+(t^2)+1;
polys={1-(s^2)+(t^2)+den x,-2 s+den y,-2 s t+den z,-1+recip den};
GroebnerBasis[polys,{x,y,z},{s,t,recip},MonomialOrder->EliminationOrder]

The result is a basis for the set of implicit equations satisfied by the original parametric surface. There is only one in this example because we parametrized a hypersurface, that is, the codimension is one.

If the parametrization is proper, that is, points on the surface have one preimage in the parameter space, then one might wish to invert it to solve for parameters in terms of the surface points (this need can arise, for example, in the context of computer-aided geometric design). This inversion may be done as below.

GroebnerBasis[polys,{s,t,z},{recip},MonomialOrder->EliminationOrder,CoefficientDomain->RationalFunctions]

Note that we again get our implicit equation as well as polynomials that are linear in the parametrization variables. Hence we may solve for these, obtaining {s==(1+x)/y,t==z/y}.

Below is an example that required both implicitization and basis approximation methods in order to obtain a result. It was posted to the Usenet news group sci.math.symbolic in July 1997. We are given a parametricized plane curve of degree 11 in the parameter t.

polys={x-(110 (t^2)-495 (t^3)+1320 (t^4)-2772 (t^5)+5082 (t^6)-7590 (t^7)+8085 (t^8)-5555 (t^9)+2189 (t^10)-374 (t^11)),y-(22 t-110 (t^2)+330 (t^3)-1848 (t^5)+3696 (t^6)-3300 (t^7)+1650 (t^8)-550 (t^9)+88 (t^10)+22 (t^11))};

We want to find an implicit form. The development version of Mathematica will handle this provided one starts at 250 or so digits of precision (200 did not suffice). The poster was unable to obtain a result via direct computation, using Mathematica and at least one more specialized program [Unterkircher 1997]. So clearly the inexact method is of interest here. Results below are abbreviated for purposes of, well, brevity.

Short[Timing[poly250 = First[GroebnerBasis[
	polys, {x,y}, t, MonomialOrder->EliminationOrder,
	  CoefficientDomain->InexactNumbers[250]]]], 6]

The result contains decimals that seem short relative to their precision, which indicates that after rationalizing we may have small denominators. After playing a bit with Rationalize, it seemed that the common denominator would be 1024. And indeed we do get a polynomial with integer coefficients below.

Short[newpoly250 = First[fastRationalize[Expand[poly250*1024], {x,y}]],4]

This is in fact the correct exact implicit polynomial in {x,y} for the curve defined parametrically by polys. First evidence that this is so comes when one repeats the computation at 500 digits of precision, rationalizes, and gets the same polynomial. More importantly, one can prove that this is the implicit polynomial using the following facts.
(i) The polynomial is irreducible (Factor will show this)
(ii) Polynomial reduction with respect to a Gröbner basis will reduce  it to zero. 
As it happens, the input is already a Gröbner basis with respect to the lexicographic term order in the variables {x,y,t} (this is of course always the case for ideals defined by parametric equations). A simple reduction with respect to polys, in this term order, gives zero.

Last[PolynomialReduce[newpoly250,{x-(110 (t^2)-495 (t^3)+1320 (t^4)-2772 (t^5)+5082 (t^6)-7590 (t^7)+8085 (t^8)-5555 (t^9)+2189 (t^10)-374 (t^11)),y-(22 t-110 (t^2)+330 (t^3)-1848 (t^5)+3696 (t^6)-3300 (t^7)+1650 (t^8)-550 (t^9)+88 (t^10)+22 (t^11))},{x,y,t}]]

One can also do such an implicitization in two dimensions by computing a resultant e.g. via the Sylvester determinant method. Other techniques to handle this include the Gröbner walk [Amrheim, Gloor, and Küchlin 1997] and the F4 implementation of Gröbner bases [Faugére 1999] using an elimination term ordering; these results were privately communicated to the author by Oliver Gloor and Jean-Charles Faugére, respectively. The former uses a strategy, applicable to finding hypersurfaces, called "sudden death" in the reference. In essence it only computes as much of the basis conversion as it needs in order to obtain the polynomial that defines the hypersurface.

2.	Polynomial algebra with Gröbner bases

Several ways to invert a polynomial modulo powers of x

Given a polynomial p[x] in x with nonzero constant term, one can use Series to find the first n terms of the reciprocal. There are several other ways to do this. We could, for example, use PolynomialPowerMod (from the Mathematica standard package Algebra`PolynomialPowerMod`) with -1 as our power and x^n as the polynomial modulus (and we also can have an integer modulus). Below are some other methods. Note that relative timings for different methods may vary considerably depending on the types of coefficients used e.g. inexact numbers vs. exact numbers vs. symbolic entities, and on idiosyncracies in the coding. In the example below we have symbolic coefficients and work modulo 2.
It should be noted that PolynomialPowerMod, for negative powers, uses GroebnerBasis. We show a form of the inversion code it uses that also appears in [Lichtblau 1996]. This method of inverting one polynomial modulo another can be found in [Buchberger 1985]; extensions are in [Becker, Weispfenning, and Kredel 1993, chapter 7] and [Adams and Loustaunau 1993, chapter 2, section 6]. The example below came from a posting to the Usenet news group comp.soft-sys.math.mathematica.

PolynomialInverseMod[f_,g_,x_,p_Integer:0]:=Module[{finv,gb},gb=GroebnerBasis[{finv f-1,g},{finv,x},Modulus->p,CoefficientDomain->RationalFunctions];
PolynomialMod[finv/.First[Solve[gb[[2]]==0,finv]],p]]/;PrimeQ[p]||p===0

p[x]= 1+x a[1]+(x^2) a[2]+(x^3) a[3]+(x^16) a[16];

Timing[inv1= PolynomialInverseMod[p[x],x^16, x, 2];]

Normal[PolynomialMod[inv1*p[x], 2] + O[x]^15]

As an aside, note that polynomial reciprocals can be used to rationalize denominators. For example, if we have recip=1/(3+Sqrt[2]) and wish to write it as a polynomial in Sqrt[2] we might set x=Sqrt[2] and find the reciprocal of 3+x modulo the polynomial (x^2)-2.

GroebnerBasis[{recip*(3+x)-1, x^2-2}, {recip,x}]

We obtain the rationalized fraction by solving the second element for recip: (3-Sqrt[2])/7. By using a different term order (one that eliminates x) we could instead get the minimal polynomial for the algebraic number recip.

GroebnerBasis[{recip*(3+x)-1, x^2-2}, recip, x]

Here is another method to find a series reciprocal. We set up a system of equations by writing

recip = (Subscript[a, 0])+(Subscript[a, 1])*x + ... + (Subscript[a, n])*(x^n)

and then expanding recip*p[x], truncating powers greater than n-1, then setting the constant term to one and setting all terms of higher degree to zero, and solving for the coefficients. We do this with GroebnerBasis instead of Solve or LinearSolve in part because Gröbner bases are the topic of this note but also for reasons of efficiency that are explained below. The cost is a few more lines of code to handle the less convenient form of result in this case.

We can gain some speed using an iterated form of this method that is similar to Newton iterations. We start with one correct term (the constant term). We double that by obtaining the degree one term. We double again, obtaining the degree two and three terms, etc. When the above is rewritten in this fashion it runs much faster on our example.

seriesInvert2[poly_, degree_Integer, x_, mod_Integer:0] := Module[
	{iter=Ceiling[Log[2,degree]], j, result=1,
    corr, aa, unknowns, diff, coeffs, soln, reprule},
	For [j=0, j<iter, j++,
		diff = PolynomialQuotient[(poly*result-1), x^(2^j), x] +
		  O[x]^(2^j);
		unknowns = Array[aa, 2^j];
		corr = unknowns . Table[x^k, {k,0,2^j-1}];
		coeffs = CoefficientList[Normal[corr*poly + diff], x,
		  Modulus->mod];
		soln = Reverse[GroebnerBasis[coeffs, unknowns, Modulus->mod]];
		soln = -soln /. Thread[unknowns->0];
		reprule = Thread[unknowns->soln];
		corr = corr /. reprule;
		result += Expand[x^(2^j)*corr, Modulus->mod];
		];
	Collect[result, x]
	]

Timing[inv2 = seriesInvert2[p[x], 16, x, 2];]
Expand[inv2-inv1, Modulus->2]

Coefficient[inv2, x^15]

Why do we not explicitly set up a matrix for our linear equations? That is, why use GroebnerBasis instead of e.g. LinearSolve? It is a matter of efficiency. To use LinearSolve first we would need to convert explicit equations in variables to a matrix (this is not too bad, O((n^2)) complexity where n is the number of equations) and then we would use it. But solving this is going to be O((n^3)) (unless we use advanced options not yet generally available) whereas it need only be O((n^2)) because the system is triangular. To convince yourself that it is triangular, multiply out a few terms of recip*p[x] and you will see the constant term involves only one unknown, the degree one term involves that and another, the degree two term involves a third, etc. Now LinearSolve is inefficient precisely because it will not take advantage of all the zeroes in the upper triangle. GroebnerBasis, on the other hand, only sees terms that are represent in the polynomials (it has nothing equivalent to zeros in matrix elements for terms not present). Hence it is doing only the necessary work.

Another method might be a brute-force Taylor expansion.  Yet another way is to call Series on 1/p[x]. This is actually a bit slower than it might otherwise be because Series is not implemented over finite fields and hence must do slower operations over the rationals and later adjust for the modulus in a post-processing step.

Timing[inv3 = PolynomialMod[Normal[Series[1/p[x], {x,0,15}]], 2];]
Expand[inv3-inv1, Modulus->2]

Here is an explicit coding of the method used by Series. It also works by doubling the number of correct terms at each step.

seriesInvert4[poly_,deg_Integer,x_,mod_Integer:0]:=Module[{recip=1,wdeg,subpoly},For[wdeg=1,wdeg≤deg,wdeg*=2,subpoly=PolynomialMod[poly,{(x^wdeg)}];recip=2 recip-(recip^2) subpoly;recip=PolynomialMod[recip,{(x^wdeg),mod}];];Collect[recip,x]]
Timing[inv4=seriesInvert4[p[x],16,x,2];]
PolynomialMod[inv4 p[x],{(x^16),2}]

We explain briefly the above method. It "converges" quadratically in the sense that each step will compute twice as many terms as the previous step. Convergence is now measured in terms of degree rather than numerical proximity. It is actually just a form of Newton's method, working in an algebraic rather than analytic setting. We begin with the equation f(y) = 1/y - p where p=p(x) is a polynomial in some ring R[x] (or it could be a series in R[[x]]). The equation f(y)==0 implies that y = 1/p. We seek a truncated series solution (that is, a polynomial) for y. Suppose we have a solution, (Subscript[y, k]), that is correct modulo (x^k). We write 0=f[y]=f[(Subscript[y, k])] + f'[(Subscript[y, k])]*(y-(x^k))+((O[y-(x^k)])^2). Notice that f'[y] is just 1/(y^2). Solving for y above gives a value that is correct mod (x^(2*k)). In fact we simply get the familiar Newton's formula from calculus, but now the "error" is not measured by numerical length but rather by how small the degree of correction terms. (That is, lower degree menas larger error. One formalizes this notion using valuation rings.) After some algebra the equation above becomes (Subscript[y, 2k])=2(Subscript[y, k])-((Subscript[y, k])^2)p which is exactly what was used in the code above. A reference for this is chapter 16 of [Zippel 1993].  Note that for higher degree inversion when coefficients are numeric, one would typically want to use asymptotically fast methods for univariate polynomial multiplication.

Standard bases

For some purposes e.g. local algebra one wishes to have a term order wherein "smaller" degrees are considered "larger" terms. This is in essence the way one works with power series, in that the "main" terms, in contrast to polynomials, are those of smallest rather than largest degree. This generalization typically goes by the name of "standard bases". Moreover it encompasses what are called mixed orderings, wherein some variables are less than constants and others are greater, for purposes of term ordering. We treat the most basic case wherein all variables are ordered lower than constants. The simple algorithm demonstrated below by now belongs to computational folklore; it appears to have originated in [Lazard 1983]. The idea is that we first homogenize our polynomials with respect to a new variable and order the variables so that the new one is highest. We next compute a Gröbner basis with respect to this order, and last we dehomogenize. To handle mixed orders one would homogenize only with respect to the variables that are smaller than constants in the term ordering, and the new term order would put the homogenizing variable between these and the other variables.

The homogenization code below assumes that the input is expanded.

homogenize[poly_,vars_,new_]/;(Head[poly]=!=Plus && Head[poly]=!=List):=poly
homogenize[poly_,vars_,new_]:=Module[{degfunc,totdeg,j},degfunc=Apply[Plus,Map[Function[x,Exponent[#,x]&],vars]];
degfunc=Distribute[degfunc,Function,Plus];
totdeg=Max[Map[degfunc,Apply[List,poly]]];
Apply[Plus,Table[poly[[j]]*new^(totdeg-degfunc[poly[[j]]]),{j,Length[poly]}]]]
homogenize[polys_List,vars_, v_]:=Map[homogenize[#,vars, v]&,polys]

A useful application is the computation of generators of the tangent cone of a variety at a given point. This is discussed in [Cox, Little, and O'Shea 1992 chapter 9, section 7] albeit in the language of Gröbner (rather than standard) bases. We compute a standard basis and then take the initial forms of each member (these are the terms of lowest total degree in each polynomial). These generate the tangent cone. The example below assumes that the point of interest is the origin.

vars = {x,y,z};
polys = Expand[{x*y*z,x^2*y - 2*y^2*z, x*z+z(y^2-z^2)}];
hpolys = homogenize[polys,vars,w]

hgb = GroebnerBasis[hpolys, {w,x,y,z}, MonomialOrder->DegreeReverseLexicographic];
standardbasis = hgb/. w->1

Now we extract the initial forms. The code below assumes they are expanded polynomials (which of course they will be, if produced by GroebnerBasis).

initialForm[poly_Plus,vars_List] :=Module[
{t,newp},
newp = poly /. Thread[vars->t*vars];
Coefficient[newp,t,Exponent[newp,t,Min]] /. t->1
]
initialForm[poly_Plus,var_] := initialForm[poly,{var}]
initialForm[poly_,_] := poly

inits = Map[initialForm[#,vars]&,standardbasis]

This set of generators for the tangent cone can be simplified.

GroebnerBasis[inits,vars]

Note that it is not what we obtain by forming a basis from the initial forms of the original polynomials.

GroebnerBasis[Map[initialForm[#,vars]&,polys],vars]

Here is an example that is related to the topic of standard bases. Suppose we want to rewrite all elements in a polynomial, modulo some ideal, in such a way that we maximize total degrees. This might be useful e.g. as a preprocessing step prior to truncating high degrees. We take our polynomial and set of ideal generators as below.

poly=-(x^3)+9 x (y^2)+2 (y^3)+x (y^3)+6 (x^2) z+(x^3) z-3 y (z^2)+3 (z^3)+y (z^3);
polys=homogenize[{(x^3)+x (y^2)+(z^2)-2,x z+2 (y^2)-3 x y-1,(x^2)+2 y z-3},{x,y,z},w]

We will next obtain our Gröbner basis for the homogenized idea. While we might dehomogenize to recover the standard basis for the original set of polynomials, it is not necessary to do so in this example as we will continue with the homogenized basis.

hgb = GroebnerBasis[polys, {w,x,y,z}, MonomialOrder->DegreeReverseLexicographic]

Now we reduce poly. We must homogenize it, reduce using the homogenized basis, then dehomogenize the reductum.

Last[PolynomialReduce[homogenize[poly,{x,y,z},w], hgb, {w,x,y,z}, MonomialOrder->DegreeLexicographic]] /. w->1

Polynomial greatest common divisiors over extension fields

Polynomial greatest common divisor (gcd) computations are the most important manipulations in polynomial algebra after such basic functions as addition and expansion. One uses them to put rational functions over a common denominator, to do square-free factorization, and as a part of the implementation of full factorization. Moreover, when we find a Gröbner basis over a rational function field (that is, we have parameters in our coefficient field) then we must remove contents of coefficients, and this too involves a polynomial gcd extraction.

In this section we illustrate a Gröbner basis method for compting polynomial gcds. It is generally not competitive with other methods, hence usually falls in the category of "not quite practical." All the same it has considerable power, for example it is trivial to implement this over algebraic extensions of a ground field. It also can find good use in working over prime fields, as will be explained below.

The method we demonstrate is covered in [Adams and Loustaunau 1993, chapter 2, section 3] or [Cox, Little, and O'Shea 1992, chapter 4, section 3]. One first finds a polynomial least common multiple (lcm).  Say poly1 and poly2 are the polynomials. Letting w be a new indeterminate, one finds lcm as the intersection of the ideals formed by w*poly1, (1-w)*poly2, and the original variables (that is, we compute using a term order that eliminates w). If we want to work over a prime field ℤ /(p) (which we will denote in the sequel by (Subscript[ℤ, p])) then we simply compute this Gröbner basis in that field. If we want an algebraic extension of our ground field then we may augment by a new indeterminate and give the irreducible polynomial, over that field, for which that indeterminate is satisfied (giving an algebraic number as an extension to PolynomialGCD is essentially equivalent to giving a minimal polynomial). Finally, noting that the greatest common divisor, gcd,  is poly1*poly2/lcm, we have our algorithm to find gcd.

There is some subtlety to the variable ordering. It turns out that, in addition to having an elimination order that removes the auxiliary variable w, we must have the algebraic variable ordered lower than all others. Otherwise we may miss elements of minimal degree in the non-algebraic variable because terms in that variable appear to higher total degree.  For example, with one elimination variable, one polynomial variable, and one algebraic ring element, we are left with lexicographic order.

Below is the polynomial lcm code.

modPolynomialLCM[poly1_,poly2_,mod_Integer,ext_:0]/;mod==0||PrimeQ[mod]:=Module[{w,evars=Variables[ext],vars},vars=Complement[Variables[{poly1,poly2}],evars];Last[GroebnerBasis[{w poly1,(1-w) poly2,ext},Join[vars,evars],w,Modulus->mod]]/;Length[evars]≤1]

As noted, we have gcd =(poly1*poly2)/lcm. Now a new issue arises insofar as the quotient may not appear to divide exactly; it only need do so when we take into account the algebraic nature of the field extension generator. So
 we cannot use the Mathematica Together command. Instead we notice that

poly1*poly2 - lcm*q1 + ext*q2

for some {q1,q2} in our ring. Modulo the extension-defining polynomial ext and prime modulus we thus have  poly1*poly2/lcm = q1. We obtain q1 as the first element of the set of multipliers from PolynomialReduce[poly1*poly2, {lcm, ext}...]. Since generalized division is not unique it is possible to have ext divide some part of q1, that is, we must still reduce it. This all is put into the code below.

modPolynomialGCD[poly1_,poly2_,mod_Integer/;mod==0||PrimeQ[mod],ext_:0]:=Module[{lcm,rat,vars,evars=Variables[ext]},vars=Complement[Variables[{poly1,poly2}],evars];lcm=modPolynomialLCM[poly1,poly2,mod,ext];rat=PolynomialReduce[poly1 poly2,{lcm,ext},Join[vars,evars],Modulus->mod];Last[PolynomialReduce[rat[[1,1]],ext,Modulus->mod]]]

Below is a very simple example. We choose the algebraic extension specifically so that the gcd will be nontrivial; it will simply be the second polynomial. This is true both over the rationals and the integers modulo a prime.

g1=modPolynomialGCD[(x^2)-5,x-a,7,(a^2)-5]
g2=modPolynomialGCD[(x^2)-5,x-a,0,(a^2)-5]

Now we will make examples in {x,a} where a is an algebraic extension. The example we construct will have a nontrivial gcd only once we take into account the algebraic equation satisfied by this element.

randomPolynomial[deg_Integer,var_]:=Table[(var^j),{j,0,deg}].RandomInteger[{-10,10},deg+1]
randomBivariatePolynomial[a_,dega_,x_,degx_]:=Expand[RandomInteger[{1,10}] (x^degx)+(∑_j^degx)randomPolynomial[dega,a] (x^(j-1))]
SeedRandom[1111];
p1=randomBivariatePolynomial[a,2,x,5];
p2=randomBivariatePolynomial[a,2,x,5];
p3=randomBivariatePolynomial[a,2,x,5];
f=PolynomialReduce[Expand[p1 p3],(a^3)+37 a+51,{x,a},Modulus->97][[2]];
g=PolynomialReduce[Expand[p2 p3],(a^3)+37 a+51,{x,a},Modulus->97][[2]];
Timing[modPolynomialLCM[f,g,97,(a^3)+37 a+51]]
Timing[modPolynomialGCD[f,g,97,(a^3)+37 a+51]]

It should be mentioned that there is another simple approach to polynomial gcds via Gröbner bases. It is well known that the Gröbner basis of two univariate polynomials is their gcd. We can simply designate one variable as "main", and treat the rest as parameters by setting CoefficientDomain->RationalFunctions. The built-in GroebnerBasis code will not actually create any denominators as that would be quite inefficient. Furthermore it will keep numerators in simplest terms by removing common factors. So the idea is to remove the contents (which are in terms of the remaining variables), take their gcds, and then find the gcd of the primitive part. This can be thus be written as a recursive procedure (though the step of removing common factors, performed internally in GroebnerBasis, will of course rely on the built-in PolynomialGCD). As with the method above, this can also be used in modular fields and extensions thereof, simply by specifying a modulus and a minimal polynomial in a new variable to define the extension. This method of finding a polynomial gcd goes by the name of "polynomial primitive remainder sequence". A similar but possibly more efficient Gröbner basis method for extraction of greatest common divisors is discussed as algorithm 1 in [Sasaki and Suzuki 1992].

It turns out that this method for finding polynomial greatest common divisors can be used to advantage in some situations. When working in characteristic zero, say, one has recourse to a nice modular interpolation method for finding multivariate polynomial gcds, as described in [Zippel 1993]. It also may work in prime characteristic, but only when there are sufficiently many evaluation points (and no unlucky choices are made). A simple recourse is to use a subresultant polynomial remainder sequence [Collins 1967]. The problem is that this will often simply not run to completion in reasonable time. Another method is to ensure that enough points are available is to move to an extension field. A nice treatment of this technique, and of a more efficient variant, may be found in [Kaltofen and Monagan 1999]. This, however, is laborious to code. We describe a simple intermediate route. A reduced remainder sequence may be regarded as similar, computationally, to a Gröbner basis computation with lexicographic ordering imposed. But, as seen above, we can instead use a similar computation only with respect to a more desirable term ordering. Specifically we may use an elimination order that makes use of total degree in the remaining variables. We illustrate below using code that is quite similar to that already shown, but with the extension field functionality removed. It is not too difficult to merge these but we want to illustrate without added complications to the term ordering.

modPolynomialLCM[poly1_,poly2_,mod_Integer]/;mod==0||PrimeQ[mod]:=Module[{w,vars},vars=Variables[{poly1,poly2}]; Last[GroebnerBasis[{w*poly1,(1-w)*poly2},vars,w,MonomialOrder->EliminationOrder,Modulus->mod]]]

modPolynomialGCD[poly1_,poly2_,mod_Integer]/;mod==0||PrimeQ[mod]:=Module[{lcm,rat,vars},vars=Variables[{poly1,poly2}]; lcm=modPolynomialLCM[poly1,poly2,mod]; PolynomialReduce[poly1*poly2,{lcm},vars,Modulus->mod][[1,1]]]

The example below is taken from [Kaltofen and Monagan 1999] (where it is attributed to H. Dobbertin). While not blindingly fast, it does give a result in a few minutes. The standard modular algorithm [Zippel 1993] that is first attempted takes a bit under a second for the same example, for moduli sufficiently large that one may obtain enough evaluation points.

zz=(((b^2)+(a^2) (b^2)+(z^2)+(b^2) (z^2)+(a^2) (b^2) (z^2))/z);
aa=(((a^2)+(c^2)+(a^2) (c^2)+(b^2) (c^2)+(a^2) (b^2) (c^2))/a);
bb=(((b^2)+z+(b^2) z+(c^2) z+(b^2) (c^2) z)/b);
cc=((a+(c^2)+a (c^2)+a z+a (c^2) z)/c);
pp=(a^3) (b^2) (c^2) (z^2) ((bb^2) (cc^2) (zz^2)+(1+aa+(cc^2)) (((1+aa+zz))^2));
pp2=Expand[pp,Modulus->2];
dpp2da=PolynomialMod[D[pp2,a],2];

Timing[g2=modPolynomialGCD[pp2,dpp2da,2]]

Multivariate factorization via Gröbner bases

Generaly speaking there are nice algorithms to factor multivariate polynomials over various fields, that have complexity much better than Gröbner basis computations. All the same, there may be circumstances in which one might choose to use an approach via Gröbner bases. For example, if factoring over a prime field, it is known how to reduce from a general multivariate problem to a bivariate problem via homomorphic images [Kaltofen 1984]. Unfortunately one cannot easily reduce in this way to a one variable factorization problem because the univariate images will typically have more factors than the original polynomial. In other words, certain of the image factors must be grouped together in order to recover the multivariate factor. One algorithm, due to Kronecker, involves substituting sufficiently high powers of one variable for the others, factoring the univariate image, and trying all combinations of image factors until it recovers true factors [von zur Gathen and Gerhard 1999]. This has exponential complexity in the number of univariate factors. Various different methods have been devised in order to reconstruct more readily a bivariate factor from a univariate one; they typically use some flavor of Hensel lifting (see, for example, [Bernardin and Monagan 1997] and its several relevant references). The technique we employ below, using Gröbner bases, is a variant of an algorithm that appeared in the literature some time ago [Gianni and Trager 1985]. The description in the reference is a bit terse but it seems to require working in an extension field in order to guarantee that the hypotheses be satisfied (at least for sufficiently small fields). Moreover the author is unaware of any proof of this method in the literature. Indeed, the method and proof below grew out of an attempt to certify the correctness of the one in that paper.

We start with a bivariate polynomial f(x,y). We assume that it is square-free and has no factor involving only one variable (that is, it is primitive with respect to either). To enforce these conditions we may remove content and multiple factors via computations that involve polynomial greatest common divisors, which may be handled by the method presented above if need be. We then replace one variable, say x, by a low degree polynomial in q(y) (say the degree is n). We require that the total degree of the result be the same as that obtained by replacing x with a generic polynomial of degree n in y. This ensures that the image factors can be grouped in such a way that their product is the image of a factor of f(x,y). We also require that there be at least one square free image factor; the reason will be explained presently. These requirements can be enforced by the simple heuristic of trying several replacements until we find one that works.

We will illustrate with a simple example, and then explain the relevant theory. We will write a polynomial explicitly as a product, expand modulo 3, and then recover the factors (note the possibility that some of the original factors themselves factor nontrivially modulo 3).

poly = Expand[(2+x+(x^2) y+(y^2)) (1+2 y+(x^2) y+x (y^2)) (2+x+2 (x^2)+(x^2) y+(y^2)+2 x (y^2)+(y^3)), Modulus->3];

First we verify that it is primitive with respect to each variable.

FactorTerms[poly,y,Modulus->3]===poly
FactorTerms[poly,x,Modulus->3]===poly

We now form a univariate polynomial by substituting y+1 for x. Since the degree of our substitution in y is one, it suffices to check that the new degree in y agrees with that of substitution a random value times y for x.

xsub=y+1;
ypoly=Expand[poly/.x->xsub,Modulus->3];
Exponent[ypoly,y]===Exponent[poly/.x->RandomReal[] y,y]

We obtain the factorization of this univariate image.

yfax=FactorList[ypoly,Modulus->3]

Notice that the first factor is a constant, and all but one are square free. We will get rid of the powers, and drop the constant and the bad factor.

yfax = Map[First, Select[yfax, (!NumberQ[#[[1]]] && #[[2]]==1)&]]

We may now pick any of these to lift to a factor of poly. The idea, roughly, is to find a polynomial of smallest total degree in the ideal generated by poly and sufficiently high powers of the univariate image factor and the polynomial that defines the substitution equivalence. As we have square free image factors, we are guaranteed to obtain at least one nontrivial factor of poly in this way. We will show this for each of the good image factors.

First[GroebnerBasis[{poly,yfax[[1]]^9,(x-xsub)^9},{x,y},Modulus->3,MonomialOrder->DegreeReverseLexicographic]]

First[GroebnerBasis[{poly,yfax[[2]]^9,(x-xsub)^9},{x,y},Modulus->3,MonomialOrder->DegreeReverseLexicographic]]

First[GroebnerBasis[{poly,yfax[[3]]^3,(x-xsub)^3},{x,y},Modulus->3,MonomialOrder->DegreeReverseLexicographic]]

Notice that the first two actually give the same factor, up to a unit in (Subscript[ℤ, 3]). We can now divide out by all recovered factors and repeat the process (using a different substitution) with the factor that remains, to see if it has nontrivial factorization modulo 3.

We now explain why this works. Our set up is as follows. We begin with a bivariate polynomial f(x,y) in (Subscript[ℤ, p])[x,y] for some prime p. We assume we have found an irreducible polynomial in that ring, of the form q(x,y) = x-w(y), such that the image of f(x,y) modulo this irreducible satisfies the hypotheses presented above (in brief, it retains "full" degree and has a square free factor). We suppose that f=g*(Subscript[g, 2]) is a factorization of the bivariate, with g irreducible (we do not rule out the possibility that it might be equal to f(x,y), that is, there might be no nontrivial factorization). Modulo q(x,y) (that is, in the polynomial ring  (Subscript[ℤ, p])/q(x,y)) we suppose moreover that f(Subscript[≡, q])h*k, where h is irreducible, relatively prime to k, and divides g in this quotient ring, and neither h nor k is a unit (if the homomorphic image has only a trivial factorization, then the same is clearly true of f(x,y)). Note that h and k are univariate polynomials, since this quotient ring is isomorphic to (Subscript[ℤ, p])[y].

With these assumptions we may thus write

f(Subscript[≡, (q^(m+1))])(h+q*(Subscript[r, 1])+(q^2)*(Subscript[r, 2])+...+(q^m)*(Subscript[r, m])) * (k+q*(Subscript[s, 1])+(q^2)*(Subscript[s, 2])+...+(q^m)*(Subscript[s, m]))

This is just a statement of the Hensel Lemma [Zippel 1993; von zur Gathen and Gerhard 1999] for lifting factors in an 𝔪-adic setting (where 𝔪 is an ideal in our polynomial ring), as applied to the situation at hand; in this case, the ideal is generated by q.   We define factors of the lifted polynomials modulo powers of q:

(Subscript[h, 0])=h, (Subscript[h, 1])=h+q*(Subscript[r, 1]),(Subscript[h, 2])=h+q*(Subscript[r, 1])+(q^2)*(Subscript[r, 2]),...,  and likewise define (Subscript[k, 0])=k,(Subscript[k, 1])=k+q*(Subscript[s, 1]),....

We also define  ideals (Subscript[J, m])={f,(h^m),(q^m)}. With these definitions, we have the following results.

Proposition 1: (Subscript[h, m])∈(Subscript[J, m+1]) for all nonnegative integers m.

Proof:  We will proceed by induction.

(i) The base case is obvious since (Subscript[h, 0])=h and h∈(Subscript[J, 1]). To give some concrete feel for how the proof proceeds, we also do the next case explicitly.

(ii) Next case. By assumption we have

(h+q*(Subscript[r, 1]))*(k+q*(Subscript[s, 1])) = h*k + h*q*(Subscript[s, 1])+k*q*(Subscript[r, 1])+ O((q^2)) = f+O((q^2))

where we use O(v) to mean that the rest is divisible by v. Since f and (q^2) belong to (Subscript[J, 2]) we have

(Subscript[p, 1])= h*k + h*q*(Subscript[s, 1])+k*q*(Subscript[r, 1])∈(Subscript[J, 2])

We multiply (Subscript[p, 1]) by q and find that h*k*q∈(Subscript[J, 2]) (as the other terms clearly belong to that ideal). Also (h^2)*q ∈(Subscript[J, 2]) (because (h^2) is in the ideal). Hence for any polynomials {a,b} we have (a*h+b*k)*h*q∈(Subscript[J, 2]). Here is where we make use of our assumption that h is square free. Since {h,k}∈(Subscript[ℤ, p])[y], and each is a distinct factor of the univariate image of f modulo our homomorphism with h square free, it must be relatively prime to k. Hence by the Euclidean algorithm there exist polynomials {a,b}∈(Subscript[ℤ, p])[y] such that a*h+b*k=1. Thus we obtain h*q∈(Subscript[J, 2]). We use this to remove the middle term from (Subscript[p, 1]) and conclude that k*(h+q*(Subscript[r, 1])) = k*(Subscript[h, 1])∈(Subscript[J, 2]).

Since we also have (h^2) and h*q in (Subscript[J, 2]) we moreover obtain (h^2)+h*q=h*(Subscript[h, 1])∈(Subscript[J, 2]). Thus there are polynomials {a,b} such that (a*h+b*k)*(Subscript[h, 1])∈(Subscript[J, 2]). As above, we may conclude that (Subscript[h, 1])∈(Subscript[J, 2]). This gives the next case.

(iii) General induction step. We want to show that

(Subscript[h, m])∈(Subscript[J, m+1])={f,(h^(m+1)),(q^(m+1))} given that (Subscript[h, m-1])∈(Subscript[J, m])={f,(h^m),(q^m)}.

We break this into several sub-steps.

Claim (a): (h^m)*q, (h^(m-1))*(q^2),...h*(q^m)∈(Subscript[J, m+1]). We will prove this later.

Claim (b): q*(Subscript[h, m])∈(Subscript[J, m+1]). Proof: Since (Subscript[h, m])=(Subscript[h, m-1])+(q^m)*(Subscript[r, m]) and (Subscript[h, m-1])∈(Subscript[J, m])we have from induction that

(Subscript[h, m])=(a*f+b*(h^m)+c*(q^m))+(q^m)*(Subscript[r, m]). Now multiply by q and use claim (a) to demonstrate claim (b).

Claim (c): k*(Subscript[h, m])∈(Subscript[J, m+1]). Proof: By definition (Subscript[k, m])*(Subscript[h, m])∈(Subscript[J, m+1]) (because it is congruent to f modulo (q^(m+1))). Writing this product as

(Subscript[h, m])*(k+q*(Subscript[s, 1])+...+(q^m)*(Subscript[s, m])) and using claim (a) yields claim (c).

Claim (d): h*(Subscript[h, m])∈(Subscript[J, m+1]). Proof: Again we write (Subscript[h, m])=(a*f+b*(h^m)+c*(q^m))+(q^m)*(Subscript[r, m]) and this time multiply by h. We now use (from claim (a)) the fact that h*(q^m)∈(Subscript[J, m+1]) to conclude that h*(Subscript[h, m])∈(Subscript[J, m+1]).

Claim (e): (Subscript[h, m])∈(Subscript[J, m+1]). Proof: Using the previous two claims, we have (a*h+b*k)*(Subscript[h, m])∈(Subscript[J, m+1]) for all {a,b}. For reasons already noted, this suffices to conclude that (Subscript[h, m])∈(Subscript[J, m+1]).

It remains to prove claim (a). We have

(Subscript[k, m])*(Subscript[h, m])=(h+q*(Subscript[r, 1])+...+(q^m)*(Subscript[r, m]))*(k+q*(Subscript[s, 1])+...+(q^m)*(Subscript[s, m]))∈(Subscript[J, m+1])

We multiply by (q^m) to conclude that h*k*(q^m)∈(Subscript[J, m+1]) (because all other terms in this product are clearly in that ideal). Thus  (h^m)*k*(q^m)∈(Subscript[J, m+1]). But clearly (h^(m+1))*(q^m)∈(Subscript[J, m+1]), hence (a*h+b*k)*(h^m)*(q^m)∈(Subscript[J, m+1]) for all {a,b}, so we have (h^m)*(q^m)∈(Subscript[J, m+1]). Noting that (h^(m-1))*k*(q^m)∈(Subscript[J, m+1]) we obtain (a*h+b*k)*(h^(m-1))*(q^m)∈(Subscript[J, m+1])  for all {a,b} and so (h^(m-1))*(q^m)∈(Subscript[J, m+1]). We continue in this fashion and conclude that h*(q^m)∈(Subscript[J, m+1]).

By induction we already know that (Subscript[J, m])={f,(h^m),(h^(m-1))*q,...,h*(q^(m-1)),(q^m)}.

Moreover we just showed that (Subscript[J, m+1])={f,(h^(m+1)),(q^(m+1)),h*(q^m)}. But this is simply {f,(q^(m+1)),h*(Subscript[J, m])}, hence we have

(Subscript[J, m+1])={f,(q^(m+1)),(h^(m+1)),(h^(m-1))*q,...,h*(q^(m-1))}. This finishes claim (a) and thus finishes the proof of the proposition.

Proposition 2: The irreducible factor g∈(Subscript[J, m]) for all nonnegative integers m.

Proof: Note that proposition 1 made no claims of uniqueness of the lifted factors (Subscript[h, m]) in proposition 1 (one can attain uniqueness e.g. by imposing monicity requirements, as noted in the references, but we do not do so). Hence it applies to any lift of h as a factor of f modulo (q^(m+1)). As (Subscript[J, m]) is an ideal it thus contains all products that have such lifts as a factor. But clearly g is itself such a product, since its image is a product of h and possibly other image factors of f.

(Alternative proof: Recall our assumptions that g|f, g is irreducible, and h|g modulo q. Applying the Hensel Lemma to the factorization of g we obtain lifts (OverBar[(Subscript[h, m])]) for each nonnegative integer m, that is, (OverBar[(Subscript[h, m])])|g modulo (q^(m+1)). As g divides f we then see that (OverBar[(Subscript[h, m])])|f modulo (q^(m+1)). Now note that proposition 1 made no claims of uniqueness of the lifted factors (Subscript[h, m]) in proposition 1 (one can attain uniqueness e.g. by imposing monicity requirements, as noted in the references, but we do not do so). Hence it applies to any lift of h as a factor of f modulo (q^(m+1)). In particular it means that (OverBar[(Subscript[h, m])])∈(Subscript[J, m+1]). So (OverBar[(Subscript[h, m])]) is  a factor of g modulo (q^(m+1)) and it lies in (Subscript[J, m+1]). This implies that  g∈(Subscript[J, m+1]). The base case is trivial, hence g∈(Subscript[J, m]) for all nonnegative integers m.)

Note proposition 2 implies that (Subscript[J, m])={g,(h^m),(q^m)} since g∈(Subscript[J, m]) and g|f.

Proposition 3: (Subscript[J, m+1])={(Subscript[h, m]),(q^(m+1))} for all nonnegative integers m.

Proof: Recall we have (Subscript[h, m])=h+q*(Subscript[r, 1])+...+(q^m)*(Subscript[r, m]). Thus

((((Subscript[h, m])))^(m+1))=(((h+q*(Subscript[r, 1])+...+(q^m)*(Subscript[r, m])))^(m+1))=(h^(m+1))+s(h,q)

where terms in s(h,q) are of total degree at least m+1 in the polynomials h and q. Modulo (Subscript[h, m]) we have h=q*(Subscript[r, 1])+...+(q^m)*(Subscript[r, m]). Hence s(h,q)=O((q^(m+1))) modulo (Subscript[h, m]). This implies that ((((Subscript[h, m])))^(m+1))=(h^(m+1)) modulo {(Subscript[h, m]),(q^(m+1))}. Thus (Subscript[J, m+1])={g,(Subscript[h, m]),(q^(m+1))}. But by construction (Subscript[h, m]) divides g modulo (q^(m+1)), hence (Subscript[J, m+1])={(Subscript[h, m]),(q^(m+1))}.

Again, as no reference was made to which particular lift (Subscript[h, m]) of h we choose, proposition 3 holds for any such choice. We are now ready to finish the proof that the method is correct.

Proposition 4: Take m>??. Let (OverBar[h]) be an element of minimal total degree in (Subscript[J, m+1]). Then (OverBar[h])=g. This allows us to find g as above, using a degree-based ordering in our Gröbner basis computation.

Proof: To be written. As is the proposition itself.

We now make some remarks about the practical aspects of this method. First, while sufficient bounds for the power of h and q are given, it may be that a much lower power will suffice (moreover, experiments with the example presented above  demonstrate that the powers of h and q need not always be equal in order to get a correct factor in the lifting process). Given that the Gröbner basis computation is likely to be less strenuous at lower powers it may make sense to start with a low power, test the resulting polynomial to see whether it exactly divides f, and raise the lifting powers if not. Also note that "random" univariate polynomials of degree n over a prime field tend to have around log(n) factors [Knopfmacher and Knopfmacher 1993]. This gives reason to suspect that for most substitutions some image factors will have reasonably high degree relative to those of the corresponding true factors; such factors would typically require lower lifting powers.

We also remark that this method can be extended to the general multivariate case by substituting polynomials in the second variable, say, for the third, fourth, and so on. With high probability one obtains a bivariate polynomial whose factors are exactly he images of the factors of the original polynomial. This follows from [Kaltofen 1984; Zippel 1993].

The factorization algorithm as presented in [Gianni and Trager 1985] uses substititions of field elements for variables. This may necessitate moving to an extension field, in case the ground field does not have sufficiently many elements to satisfy the requirements imposed on the image polynomial. This can have two deleterious consequences. One is that polynomial arithmetic must now be done over the extension field. The other is that, unless the field extension is selected in a certain way (which typically involves using a large degree extension) one might get extra factors that need to be merged in order to recover actual factors over the ground field. How large an extension must be used to avoid this? If one uses a prime larger than the degree of the polynomial in every variable then one is guaranteed to have the same factorization pattern (see [Bernardin and Monagan 1997] and references therein). A quick proof of sufficiency is that merged factors would have to be conjugates in order for their product to lie in the ground ring. Using a prime degree extension guarantees that they must merge in groups at least as large as that prime. But it is chosen to exceed the degree in every variable, so there can be no merging; this forces all factors to lie in the base ring to begin with. In any case, the variation we presented above avoids extension fields but possibly uses substitutions that raise the total degree; this can make the univariate factorization and the lifting steps slower. Another  difference in the presentations is that they do not impose the same requirements for the image polynomials. The ones used here are simple to verify before proceeding to the Gröbner basis lifting step, which is why we chose them over those of the reference.

3.	Gröbner bases over modules and related computations

Computing a Gröbner basis over a vector space

We give a brief introduction to Gröbner bases over vector spaces and demonstrate a way to compute them. We will apply this method to finding conversion matrices and syzygies. The method for applying an ordinary Gröbner basis function to obtain a basis over a vector space is taken from [Lichtblau 1996]. More efficient variations, which require minor alterations to the Buchberger algorithm, may be found in [Adams and Loustaunau 1994; Becker, Weispfenning, and Kredel 1993]. As is shown in the references, this can be extended in a straightforward way to Gröbner bases of finitely generated modules over polynomial rings.

For purposes of exposition we will use a three dimensional vector space over a polynomial ring. So our vectors are ordered triples whose coordinates are polynomials. We will call the polynomial ring A. Then, as is often done in linear algebra, one could represent elements of this space as 3-tuples {a,b,c}, where each entry is in A. Instead we let {(Subscript[e, 1]), (Subscript[e, 2]), (Subscript[e, 3])} be three new "coordinate variables" and we represent elements as a*(Subscript[e, 1])+b*(Subscript[e, 2])+c*(Subscript[e, 3]), with {a, b, c} all elements of A.  It turns out that we can treat the new variables as we do all other variables, and just use GroebnerBasis to find a basis forover the vector space.

We want to form a Gröbner basis almost as usual, except that we never want to form S-polynomials between polynomials whose head terms live in different coordinates. We can enforce this simply by allowing them, but adding the polynomials (Subscript[e, j])*(Subscript[e, k]), for all pairs of indices {j,k}, so that all S-polynomials thus formed will immediately reduce to zero. This implies that all S-polynomials of interest must be linear in the coordinate variables, and this can only happen when the head terms of the two polynomials have the same coordinate.

There are various useful orders between the polynomial and coordinate variables. For our purposes we will want the "position over term" (POT) order (see [Adams and Loustaunau 1994, chapter 3, section 2]. That is, we order the coordinates and monomials in the first coordinate are always larger than monomials in later coordinates, etc., with ties broken by taking the largest monomial in the polynomial variables that resides in the given coordinate. This order in effect enforces something like Gaussian row reduction: if we regard the vectors as rows in a matrix, then we work first on the first column, then the second, etc.

There is some inefficiency in this method; having internal code to enforce the distinction between polynomial ring variables and coordinate variables would be better. All the same, we get tolerable results from this more primitive technique.

moduleGroebnerBasis[polys_, vars_, cvars_, opts___] := Module[
	{newpols, rels, len=Length[cvars], gb, j, k, ruls},
	rels = Flatten[Table[cvars[[j]]*cvars[[k]], {j,len}, {k,j,len}]];
	newpols = Join[polys, rels];
	gb = GroebnerBasis[newpols, Join[cvars, vars], opts];
	rul = Map[(# :> {})&, rels];
	gb = Flatten[gb /. rul];
	Collect[gb, cvars]
	]

We will see an application of this in the next section.

Finding a conversion matrix

Given a set of polynomials F={(Subscript[f, 1]), ..., (Subscript[f, n])} and a corresponding Gröbner basis G={(Subscript[g, 1]),... ,(Subscript[g, m])} with respect to some term order, we want to find the conversion matrix T such that G==T.F. Note that the reverse problem, finding S such that F==S.G, is easily handled by generalized division [Cox, Little, and O'Shea 1992, chapter 2, section 3], implemented in Mathematica as PolynomialReduce.

The method shown below is likely to be inefficient if there is more than a very small number of generators for F. This is because we will spend a lot of time reducing syzygies, which does not change the result but certainly influences speed. It does work nicely on a nontrivial example from [Cox, Little, and O'Shea 1992, chapter 2, section 9, exercise 13], shown below, most likely because in it F has only two generators. Note that much of [Adams and Loustaunau 1994, chapter 3] is devoted to applications in commutative algebra that require this matrix, so it is no useless exercise to find it. We begin by setting up the input. Regard the two generators as a column vector. Augment on the right with a 2x 2 identity matrix to obtain a 2x3 matrix. The identity matrix entries will keep track of multipliers of the original polynomials.

This is similar to methods that appear in [Caboaro and Traverso 1998]. The method below was shown at the 1998 IMACS-ACA conference, held about the same time as the ISSAC conference in which the referenced work was presented.

f1=(x^5)+(y^4)+(z^3)-1;
f2=(x^3)+(y^2)+(z^2)-1;
fF={f1,f2};
vars={x,y,z};
coords=Array[ee,3];
fmat={{f1,1,0},{f2,0,1}};
newfF=fmat.coords;
mgb=moduleGroebnerBasis[newfF,vars,coords];
mgb=Select[mgb,Coefficient[#1,coords[[1]]]=!=0&];
gG=(Coefficient[#1,coords[[1]]]&)/@mgb;
sS=First/@PolynomialReduce[fF,gG,vars];

Now we run a few sanity checks. We compare our result, gG, to ordinary GroebnerBasis. We check that we really get elements of gG as products of fF using the tT matrix under construction, and, vice versa, that we recover fF from gG by sS.

gb2 = GroebnerBasis[fF, vars];
gb2===gG
rul = {ee[1]->1,ee[2]->-f1,ee[3]->-f2};
Map[Expand[# /. rul]&, mgb] (* want zeroes *)
Expand[sS . gG - fF]

Let us now construct tT as a matrix. We use a simple tactic from the Mathematica standard package  LinearAlgebra`MatrixManipulation` function LinearEquationsToMatrices.

tT = Outer[D, mgb, Rest[coords]];
Expand[gG - tT.fF] (* want zeroes *)

This is almost identical to the extended gcd algorithm where one row reduces such a matrix over a Euclidean ring.In fact one can get the extended gcd using module Gröbner bases where the base ring is a Euclidean domain and we compute what is called a strong Gröbner basis over that ring.

Finding generating sets of syzygies using GroebnerBasis

Given a set of generators for an ideal one may require a set of generators for the syzygy module (this is useful for computing free resolutions of ideals: see [Adams and Loustaunau 1994, chapter 3, section 10]). We outline a way to find a syzygy module generating set in Mathematica.

We will assume that our generators are a Gröbner basis with respect to a given term order. If not then algorithms in [Adams and Loustaunau 1994, chapter 3, section 4] explain how to obtain the syzygies using the syzygies for a Gröbner basis and the transformation matrices to and from that basis that were discussed in the previous section.

One method is as follows (see e.g. [Adams and Loustaunau 1994, chapter 3, section 7] or [Becker, Weispfenning, and Kredel 1993, chapter 6, section 1]). We explicitly form S-polynomials of all pairs of elements in the basis. We record the multiplier terms used as vectors in an appropriate space (dimension=size of basis). We find the remaining multipliers using PolynomialReduce, because it outputs the generalized quotients as well as the reductum.

Another way to find these is similar to the way we found the conversion matrix in the last section. We first place our Gröbner basis polynomials, call them G={(Subscript[g, 1]),...,(Subscript[g, n])}, in a column vector. We form an n x n identity matrix, modify it in a way to be explained below, and place the result to the right of our column vector. We form a vector space Gröbner basis as before; this time we will be interested in rows whose first entries are zero. This technique is similar to ideas presented in

General theory  guarantees that the result of the module Gröbner basis will contain a generating set for the syzygy module. The difficulty is that in computing such a basis we may spend alot of time "reducing" the syzygy generating set. What would be useful is a term order such that one round of reduction of S-polynomials for the first column gives a Gröbner basis for the syzygies thus obtained. Such an ordering can be obtained as follows. Monomials in column one are larger than monomials in all other columns. Monomials in columns j and k, say, where {j,k}>1, are compared by adding the exponent vector of the lead monomial of (Subscript[g, j]) to the exponent vector of (Subscript[m, j]), likewise add the exponent vector of the lead monomial of (Subscript[g, k]) to the exponent vector of (Subscript[m, k]), and compare these using the term order for which G is a Gröbner basis. That this gives rise to a Gröbner basis for the syzygies is discussed in [Adams and Loustaunau 1994, chapter 3, section 7].

We do not have a simple way to describe this monomial order. Fortunately we can achieve it by the expedient of multiplying each column in the identity matrix by the lead monomial of the corresponding element in G. This will have the effect of adding the correct exponent vector for monomial comparisons in columns 2 through n+1 of the matrix thus formed, and when we are finished we will simply divide out the multiplier monomials we inserted into each column.

For simplicity it is assumed that G is a Gröbner basis with respect to lexicographic ordering in the code and example below.

syzygyBasis[gb_, vars_] := Module[
	{heads, nvars, len=Length[gb], hmonoms, nv, mat,
	  newpolys, fullvars, fvar, rvars, syzgb, rels, ruls}, 
	heads = Map[First[First[#]]&,
	  First[Internal`DistributedTermsList[gb, vars]]];
	nvars = Array[nv, len+1];
	fvar = First[nvars];
	rvars = Rest[nvars];
	hmonoms = Map[Apply[Times,vars^#]&, heads];
	mat = Transpose[Prepend[DiagonalMatrix[hmonoms], gb]];
	newpolys = mat . nvars;
	fullvars = Join[{fvar}, vars, rvars];
	rels = Flatten[Table[nvars[[j]]*nvars[[k]], {j,len+1}, {k,j,len+1}]];
	syzgb = GroebnerBasis[Join[newpolys,rels], fullvars];
	ruls = Map[(# :> {})&, rels];
	syzgb = Flatten[syzgb /. ruls];
	syzgb = Select[syzgb, FreeQ[#,fvar]&];
	syzgb = Outer[D, syzgb, rvars];
	Map[Together[#/hmonoms]&, syzgb]
	]

Timing[syz = syzygyBasis[gb, vars];]
Expand[syz . gb] (* want zeros *)

Computing the Hermite normal form of a polynomial matrix

Recall that a Gröbner basis for a set of linear polynomials over a given field 𝕂 is equivalent to a row-reduced matrix over that field. Suppose instead that we work over a Euclidean domain. The analogous concept to row-reduced form is the Hermite normal form (see, for example, [Hafner and McCurley 1991; Storjohann 1994; Storjohann and Labahn 1995]). Translated to polynomial ideal theory, this becomes what is called a "strong" Gröbner basis over the domain [Buchberger 1985; Kandri-Rody and Kapur 1988] of a set of linear polynomials. We will work with the particular domain of univariate polynomials over the rationals. In Mathematica one obtains such Gröbner bases by setting the option

CoefficientDomain->Polynomials[var]

where var is the variable in the base ring. To illustrate with examples we will use a generator of "random" univariate polynomials.

randomPolynomial[deg_Integer,var_]:=Table[(var^j),{j,0,deg}].RandomInteger[{-10,10},deg+1]

randomMatrix[degmax_,rows_,cols_,var_]:=Module[{deg},Table[deg=RandomInteger[{0,degmax}];randomPolynomial[deg,var],{rows},{cols}]]

We will now generate a random 5x8 matrix of polynomials of degree 2.

mat=randomMatrix[2,5,8,x];

Next we transform the rows of the matrix into polynomials using placeholder variables. Regard each column as defining a new variable. Each row then becomes a polynomial that is linear in these new variables. For example, we show the polynomial corresponding to the first row of mat.

len=Length[First[mat]];
newvars=Array[v,len];
First[mat].newvars

Now we form a module Gröbner basis of these polynomials over the ring of polynomials in x. Note that this is in some sense a trivial use of Gröbner basis because the list of non-placeholder variables empty.

Timing[mgb = moduleGroebnerBasis[mat.newvars, {}, newvars, CoefficientDomain->Polynomials[x]];]

Last, we convert back to matrix form. For brevity we show only the first row of the Hermite normal form.

hnf = Outer[D, Reverse[mgb], newvars];
First[hnf]

The reader might wonder why we needed a module Gröbner basis, given that we work with polynomials linear in the new variables. After all, head term power-products with respect to these variables are always either identical or distinct, and this is the situation in which one of Buchberger's criteria [Buchberger 1985] applies. The reason is that, when computing a strong basis over a Euclidean domain, only a weaker version of the criterion holds. This will be discussed in a future work on computation of strong Gröbner bases over such domains.

Note that we can also get a (non-unique) multiplier matrix as well as the Hermite normal form using tactics presented above for obtaining conversion matrices. An unfortunate drawback to this method for obtaining a polynomial Hermite normal form is that it is not terribly efficient. Some experimentation indicates moreover that small changes in degree or number of rows lead to disasterous speed hits. This is because no attempt has been made to control for intermediate polynomial swell. Experimentation indicates that the main problem is in growth of coefficients. Indeed, changing the GroebnerBasis call to use a prime modulus makes some of these harder examples again tractable. So one might modify this approach to use the Chinese remainder algorithm to reconstruct the polynomials from modular images. One might also use ideas in [Storjohann 1994; Storjohann and Labahn 1995] to make this computationally more efficient, though it might be preferable to code these from scratch rather than trying to retro-fit them within a Gröbner basis framework.

4.	Computing the "size" of some algebraic objects

There are problems for which one might be quite interested in some algebraic measure of size. Two important examples are finding the size of a finite solution set and finding the dimension of an ideal. The former can be recast as finding the dimension of a certain vector space, hence we count linearly independent generators. In the latter case it is algebraically independent generators that we count. We discuss ways to handle both.

Finding the size of a finite solution set to an algebraic system of equations

A common question is as follows: given some algebraic system with generic coefficients, how large is the solution set? Assuming there are finitely many solutions, this may be recast as a problem of finding the (finite) dimension of a certain vector space. Moreover it is easy to show that this must be constant for generic values of coefficients; a simple proof follows from the Shape Lemma, for example [Becker, Marinari, Mora and Traverso 1994]. In some cases one can simply compute a lexicographic Gröbner basis over a field of rational functions (so coefficients may vary). Quite commonly this is computationally too expensive. In such cases one may try any or all of the following tactics.

(i) Use a cheaper term order and count the number of terms in the "normal set" (this measures the vector space dimension). There is a nice discussion in [Corless 1996] of how to do this efficiently.
(ii) Use "random" (rather than generic parameter) coefficients, find the size, and be satisfied with a result that is only probabilistically correct.
(iii) Use method (ii) above several times with diffent sets of random coefficients.
(iv) If need be, do the computation in (ii) modulo some large prime. With high probability this will be a so-called "lucky" prime, and the vector space size will be the same as when working over the rationals [Traverso 1988; Winkler 1988].

Common situations for this sort of question to arise include enumerative geometry, computational geometry, and the physical sciences. For example, these tactics are heavily used in [Lichtblau 2000b] to handle a problem in enumerative geometry. We give a similar but simpler example below, taken from computational geometry. Given four lines in (ℝ^3), one wishes to determine all lines that intersect these four. One nice computational approach to this problem is provided in [Teller and Hohmeyer 2000]. It is not hard to show that generically the number must be finite. We will use tactics mentioned above to show that the number is two. It should be noted that we are solving for various parameters in the algebraic closure field, that is, ℂ. So for specific configurations one must still check whether the result lives in real space (also one must ascertain that the four given lines are generic, otherwise the solution set might be infinite).

One way to set up the problem is to write each line as an intersection of two planes, one parallel to the x axis and one parallel to the y axis. In other words, we are parametrizing by the z coordinate; generically this may be done. Write the fifth line also as such an intersection, but with undetermined coefficients. Then there are four points that each lie on one of the four lines and also on the fifth line. In this way we obtain 16 equations in 16 variables (12 coordinates for the 4 points and the 4 undetermined coefficients for the planes that give the fifth line). We also need an equation and variable to force that the planes giving the fifth line be independent (specifically, we force an appropriate determinant to be nonzero). All this is done below.

linepolys=Flatten[Table[{r[j]*x[j]+s[j]*z[j]+1,t[j]*y[j]+u[j]*z[j]+1},{j,4}]];
fifthlinepolys=Flatten[Table[{a*x[j]+b*z[j]+1,c*y[j]+d*z[j]+1},{j,4}]];
extra=Det[{{a,b},{c,d}}]*recip-1;
polys=Join[linepolys,fifthlinepolys,{extra}];
vars=Join[{a,b,c,d},Flatten[Table[{x[j],y[j],z[j]},{j,4}]]];

Some experiment indicates that a direct computation of an appropriate Gröbner basis is strenuous, at least for this formulation of the problem. Hence we will employ the tactic of using random values for the planes that give the four lines. When we find the basis we will eliminate the auxiliary variable that was used in the nondegeneracy condition.

params=Flatten[Table[{r[j],s[j],t[j],u[j]},{j,4}]];
reps=Thread[params->RandomInteger[{10,100},Length[params]]];
newpolys=polys/.reps;
Timing[gb=GroebnerBasis[newpolys,vars,recip]]

We see at once that there are two solutions.

Here is another way to set up the algebra. It may make more clear why we obtain two solutions. This time we will regard the four lines as given parametrically in terms of their z coordinates. We will begin again with generic parameters, and specialize for an example at the end.

line[j_]:={x-(a[j]*z+b[j]),y-(c[j]*z+d[j])};

We take two of these lines, and a point on the first. Then the set of lines through that point and intersecting the second line obviously gives a plane. We can find the implicit form of that plane by getting its normal, and this may be done by taking the cross product of any two vectors parallel to that plane. We obtain such a pair by taking two points on the second line, and subtracting from each our point on the first line. We use z[1] to denote the z coordinate of the first point, compute the x and y  coordinates in terms of that. For the two points on the second line we are free to choose any pair, so we take the ones with z coordinate of 0 and 1 respectively. So we have a set of planes, parametrized by z[1].

point[1]={a[1]*z[1]+b[1],c[1]*z[1]+d[1],z[1]};
point[2,1]={b[2],d[2],0};
point[2,2]={a[2]+b[2],c[2]+d[2],1};
vec[2,j_]:=point[2,j]-point[1];
normal[1,2]=Cross[vec[2,1],vec[2,2]];
plane[1,2]=Expand[normal[1,2].({x,y,z}-point[1])]

Next we take our third line and intersect it with such a plane. This determines a unique point on the third line.

polys=Append[line[3],plane[1,2]];
gb=GroebnerBasis[polys,{x,y,z}];
point[3]={x[3],y[3],z[3]}={x,y,z}/.First[Solve[gb==0,{x,y,z}]];

So for each point on the first line (that is, for each value of the parameter z[1]∈ℝ) we obtain a unique line intersecting both the second and third given lines, and moreover we have the coordinates of the intersection point on that third line. As the set of these intersecting lines gives a surface in (ℝ^3), we may parametrize it by

point[1] + t*(point[3]-point[1])

where t and z[1] are the parameters.

direction[1,3]=point[3]-point[1];
line[1,3]=point[1]+t*direction[1,3];
vars = {x,y,z};
surface = vars - line[1,3];

Clear[point]

As in a previous section, we will implicitize this surface. First we get a common denominator for the parametrization. Our polynomials will then consist of the numerators as well as a condition that the denominator not vanish. As before, this will give us yet another parameter to remove in the implicitization step.

ratfuns=Together[(vars-line[1,3])];
denom=Apply[PolynomialLCM,Denominator[ratfuns]];
surfpolys=Append[Numerator[ratfuns],recip*denom-1];

At this stage we resort to random rather than generic coefficients, to reduce the computational complexity.

params=Flatten[Table[{a[j],b[j],c[j],d[j]},{j,4}]];
reps=Thread[params->RandomInteger[{10,20},Length[params]]];
randomsurfpolys=surfpolys/.reps

Finally we are ready to get the implicit equation for this surface in (ℝ^3).

surface=First[GroebnerBasis[randomsurfpolys,vars,{t,z[1],recip},MonomialOrder->EliminationOrder]]

We have a quadric surface. Thus our fourth line will intersect it in two points, possibly with complex coordinates. This may help to explain why there are two lines through any four generic lines in three dimensional space.

Finding the dimension of an ideal

The dimension of an ideal in a polynomial ring over a field is, informally speaking, the number of indeterminates one needs to specify in order to obtain a finite point set. Technical definitions and relevant theory is contained in [Becker, Weispfenning, and Kredel 1993, chapters 6 section 3, chapter 7 section 5, and chapter 9 section 3]. In brief, given a polynomial ring in the indeterminates vars over a field 𝕂, and given an ideal I in that ring, we want to find what is called a maximally independent set of the indeterminates modulo I (this is similar to transcendence degree of a finite field extension). One can define this is in terms of the maximal length of chains of prime ideals , ordered by inclusion. Alternatively, as in the reference, we define this as the size of the largest subset S ⊂vars such that the elimination ideal 𝕂[S]⋂I is empty; it is shown that the definitions are equivalent.

For example, a "random" ideal generated by two polynomials in three variables over the rationals will contain no univariate polynomials, but there will be polynomials that involve all pairs of variables. Hence the dimension will be 1.

In order to apply the technology of Gröbner bases we really require some definition that may be stated with respect to a term order. To this end there is a notion of strong dimension that is more readily computed. It is a function of a specified term ordering for the ideal and may be found using the head terms of the Groebner bases for that ordering. For each head term we extract the set of variables it uses; any set of variables that contains none of these head term sets is strongly independent of the ideal with respect to that term ordering. It is moreover shown in the reference that this coincides with ordinary dimension. Thus what we need to find is a largest maximal subset of variables, S, that does not contain any of the head term variable subsets. Note that "largest maximal" is not redundant: largest refers to length, while maximal is in the sense of partial ordering by inclusion.

Here is a brief description of the code shown below. The function isIndependentSet will take a set of subsets that comprises the variable sets used in each head term. It will check whether a given subset of variables contains none of those subsets. This gives the strong independence criterion noted above.

The main function, getMaxIndependentSets, is based on the slightly more general DIMREC algorithm in the reference. As we are only interested in dimension and not in finding all possible maximal independent sets, we will take shortcuts. The idea is that we successively augment a given independent set with a new variable (we begin with the empty set which is vacuously independent). If this augmented set is also independent then call recursively to see if we can augment by more variables. If not, we add this one to a list of maximal (by inclusion ordering) independent sets. For efficiency we keep tabs on the largest thus far because we can take an early exit from the main loop whenever there is no hope of beating that maximum. While we do not necessarily find all maximal-by-inclusion independent subsets we will get at least one such with longest length, which suffices for our purposes.

The driver routine, krullDimension, finds the subsets of variables that appear in each head term as given by a particular term order.

firstContainsSecond[l1_,l2_]:=(Union[l1,l2]===l1)

isIndependentSet[set_,sets_]:=Map[!firstContainsSecond[set,#]&,sets]

getMaxIndependentSets[vars_,inset_,heds_,maxlen_,indx_,sets_]:=Module[{currentset,vlen=Length[vars],ilen=Length[inset],enlarged=False,newmax=maxlen,maxsets=sets},Do[If[ilen+vlen-i+1≤maxlen,Break[]];
currentset=Append[inset,vars[[i]]];
If[isIndependentSet[currentset,heds],{maxsets,enlarged,newmax}=getMaxIndependentSets[vars,currentset,heds,newmax,i+1,maxsets];
If[!enlarged,maxsets={maxsets,currentset};
newmax=Max[newmax,Length[currentset]];
enlarged=True;];],{i,indx,vlen}];
{maxsets,enlarged,newmax}]

idealDimension[ideal_,vars_]:=Module[{ord=MonomialOrder->DegreeReverseLexicographic,gb,pheads,heds,maxsets,el,mlen},
gb=GroebnerBasis[ideal,vars,ord];
pheads=Map[First[First[#]]&,First[Internal`DistributedTermsList[gb,vars,ord]]];
heds=Apply[And,(Map[(vars*#)&,(pheads/._?Positive->1)]/.0->Sequence[])];
{maxsets,el,mlen}=getMaxIndependentSets[vars,{},heds,0,1,{}];
mlen]

We illustrate with a pair of random examples. The first is an ideal generated by two polynomials in three variables so we expect the dimension to be one. The second is generated by three polynomials in five variables so the expected dimension is two.

polys1={-4+(x^2) y+3 x z,2-x y+(y^2)+z};
vars1={x,y,z};
Timing[krullDimension[polys1,vars1]]

polys2={-4+(x^2) y+3 w x z,-3+2 x-(w^2) x y+t (y^2)+t z,7+w (x^2) y+2 (t^2) x (z^2)-5 w y (z^2)};
vars2={t,w,x,y,z};
Timing[krullDimension[polys2,vars2]]

Summary and further directions

We covered computation and use of approximated Gröbner bases, simple FGLM-type conversion to lexicographic bases in zero dimensional ideals when the shape lemma holds, the tactic of obtaining an exact basis from an approximation that is sufficiently close to the exact result, and its use in an example involving approximate implicitization.
These approximatation techniques are quite valuable for industrial use, as, for example, in [Lichtblau 2000a]. The rationalization methods are perhaps just novelty items. Basis conversion via linear algebra is sometimes useful, though appears to be largely superceded by other techniques for fast Gröbner basis computation. Some of the newer competing technologies are noted below.

(i) P-adic lifting (find a basis mod p, use some flavor of Hensel-lifting to get a basis over the rationals. See [Winkler 1988] for more about this method.
(ii) Gröbner trace: use a modular computation as a "oracle" to avoid costly but ultimately wasteful S-polynomial reductions over the rationals [Traverso 1988].
(iii) Gröbner walk: a sophisticated basis-conversion algorithm from [Collart, Kalkbrenner, and Mall 1997].
(iv) Possibly some hybrid of (ii) and (iii) might be useful in cases where coefficient growth is a problem during various basis computations that occur in the walk algorithm.
(v) More efficient direct computation as described in [Faugére 1999]. This is very useful as a starting point for conversion methods such as the Gröbner walk.

The Gröbner walk has been implemented in Mathematica (using the built-in GroebnerBasis and related functions) and preliminary results look promising. We also plan to experiment with the trace method because for problems where intermediate coefficient growth is the bottleneck it holds great promise. One need only compute a basis for the first example in this paper, modulo a machine number prime, and compare to the time of computing a basis over the rationals, to see this. Moreover it may be the case that the trace method, used probabilistically, would give considerable speed boost when working over a rational function field. This might make tractable, for example, the computation of the generic case in the line-intersecting-four-lines problem.

Gröbner bases can be used to compute a myriad of things that are often approached in the computational mathematics literature by other means. For example, here we saw how to find series reciprocals and polynomial greatest common divisors. Small modifications to the input even allow us to use the garden-variety algorithm to extend GroebnerBasis itself, e.g. to find bases for modules. Applications include finding conversion matrices, finding syzygy module bases, and computing Hermite normal forms of matrices of univariate polynomials. Given the importance of conversion matrices to these various related applications in computational algebra, it would be useful to find them more efficiently using internal code; then it is merely a matter of using information that is available during the Gröbner basis computation but not thereafter. Such functionality may be implemented in a future version of Mathematica.

Acknowledgements

The author thanks Dan Grayson for inviting him to give two guest lectures to his beginning computational algebraic geometry class during the Spring 1998 semester. Much of the algorithmic material presented herein was developed for those lectures. Some material from section 3 was further refined and presented at the 1998 IMACS-ACA conference (Prague). He thanks Michael Trott for comments on a draft of this report.

References

W. Adams and P. Loustaunau (1994). An Introduction to Gröbner Bases. Graduate Studies in Mathematics 3. American Mathematical Society.

B. Amrheim, O. Gloor, and W. Küchlin (1997). On the walk. Theoretical Computer Science 187(1-2):179-202.

E. Becker, M. G. Marinari, T. Mora, and C. Traverso (1994). The shape of the Shape Lemma. Proceedings of the International Symposium on Symbolic and Algebraic Computation. 129-133. ACM Press.

T. Becker, V. Weispfenning, and H. Kredel (1993). Gröbner bases: A Computational Approach to Computer Algebra. Graduate Texts in Mathematics 141, Springer-Verlag.

L. Bernardin and M. Monagan (1997). Efficient multivariate factorization over finite fields. In: Applied Algebra, Algebraic Algorithms, and Error-Correcting Codes (AAECC-12). T. Mora and H. Mattson, eds. Lecture Notes in Computer Science 1255:15-28.

B. Buchberger (1985). Gröbner bases: An algorithmic method in polynomial ideal theory. In Multidimensional Systems Theory, chap 6. N. K. Bose, ed. D. Reidel Publishing Company.

B. Buchberger (1988). Applications of Gröbner bases in non-linear computational geometry. In: Proceedings of the International Symposium on Trends in Computer Algebra '87, R. Janssen, ed. Lecture Notes in Computer Science 296:52-80. Springer-Verlag.

M. Caboaro and C. Traverso (1998). Efficient algorithms for ideal operations (extended abstract). In Proceedings of the 1998 International Symposium on Symbolic and Algebraic Computation, O. Gloor, ed. 147-152. ACM Press.

S. Collart, M. Kalkbrenner, and D. Mall (1997). Converting bases with the Gröbner walk.. J. Symbolic Computation 24(3-4):465-469.

G. Collins (1967). Subresultants and reduced polynomial remiander sequences. Journal of the ACM 14(1):128-142.

R. Corless (1996). Editor's corner: Gröbner bases and matrix eigenproblems. Sigsam Bulletin: Communications in Computer Algebra 30(4):26-32.

D. Cox (1998). Introduction to Gröbner bases. In Proceedings of Symposia in Applied Mathematics 53, D. Cox and B. Sturmfels, eds. 1-24. ACM Press.

D. Cox, J. Little, and D. O'Shea (1992). Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Computer Algebra. Undergraduate Texts in Mathematics. Springer-Verlag.

D. Dummit and R. Foote (1991). Abstract Algebra. Prentice Hall.

J.-C. Faugére (1999). A new efficient algorithm for computing Gröbner bases ((Subscript[F, 4])). Journal of Pure and Applied Algebra 139:61-88. Elsevier Science.

J.-C. Faugére, P. Gianni, D. Lazard, and T. Mora (1993). Efficient change of ordering for Gröbner bases of zero-dimensional ideals. J. Symbolic Computation 16:329-344.

P. Gianni and B. Trager (1985). GCD's and factoring multivariate polynomials using Gröbner bases. In:Proceedings  EUROCAL '85, European Conference on Computer Algebra Vol. II, B. F. Caviness, ed. Lecture Notes in Computer Science 204:409-410. Springer-Verlag.

P. Gianni, B. Trager, and G. Zacharias (1988). Gröbner bases and primary decomposition of polynomial ideals. J. Symbolic Computation 6(2-3):149-167.

J. Hafner and K. McCurley (1991). Asymptotically fast triangularization of matrices over rings. SIAM Journal of Computing 20(6):1068-1083.

M. Kalkbrenner (1990). Implicitization of rational parametric curves and surfaces. In: Applied Algebra, Algebraic Algorithms, and Error-Correcting Codes (AAECC-8). S. Sakata, ed. Lecture Notes in Computer Science 508:249-259. Springer-Verlag.

E. Kaltofen (1984). Effective Hilbert irreducibility. In: Proceedings EUROSAM '84, International Symposium on Symbolic and Algebraic Computation, J. Fitch, ed. Lecture Notes in Computer Science 174:277-284.

E. Kaltofen and M. Monagan (1999). On the genericity of the modular polynomial GCD algorithm. In Proceedings of the 1999 International Symposium on Symbolic and Algebraic Computation, Sam Dooley, ed. ACM Press.

A. Kandri-Rody and D. Kapur (1988). Computing a Gröbner basis of a polynomial ideal over a Euclidean domain. J. Symbolic Computation 6(1):37-57.

A. Knopfmacher and J. Knopfmacher (1993). Counting irreducible factors of polynomials over a finite field. Discrete Mathematics 112:103-118. Elsevier Science Publishers.

D. Lazard (1983). Gröbner bases, Gaussian elimination, and resolution of systems of algebraic equations. In:  Proceedings EUROCAL '83, European Conference on Computer Algebra, J. A. van Hulzen, ed. Lecture Notes in Computer Science 162:146-156. Springer-Verlag.

D. Lichtblau (1996). Gröbner bases in Mathematica 3.0. The Mathematica Journal 6(4): 81-88.

D. Lichtblau (2000). Solving finite algebraic systems using numeric Gröbner bases and eigenvalues. In: Proceedings of the World Conference on Systemics, Cybernetics, and Informatics (SCI2000), Volume 10. (Concepts and Applications of Systemics, Cybernetics, and Informatics). M. Torres, J. Molero, Y. Kurihara, and A. David, eds.  555-560. International Institute of Informatics and Systemics.

D. Lichtblau (2000b). Finding cylinders through 5 points in (ℝ^3). Manuscript.

T. Sasaki and M. Suzuki (1992). Three new algorithms for multivariate polynomial GCD. J. Symbolic Computation 13:394-411.

M. Sofroniou and G. Spaletta. Precise Numerical Computation. To appear.

A. Storjohann (1994). Computation of Hermite and Smith Normal Forms of Matrices. Master's thesis, Department of Computer Science, University of Waterloo.

A. Storjohann and G. Labahn (1996). Preconditioning of rectangular polynomial matrices for efficient Hermite normal form computation. In Proceedings of the 1995 International Symposium on Symbolic and Algebraic Computation, A. H. M. Levelt, ed. 119-125. ACM Press.

S. Teller and M. Hohmeyer (2000). Determining the lines through four lines. To appear, Journal of Graphics Tools.

C. Traverso (1988). Gröbner trace algorithms. In: Proceedings of the International Symposium ISSAC '88, P. Gianni, ed. Lecture Notes in Computer Science 358:125-138. Springer-Verlag.

M. Trott (1997). Mathematica Solutions to the ISSAC '97 Challenge. SIGSAM Bulletin 31(4):2-35. Communications in Computer Algebra. ACM Press.

A. Unterkircher (1997). Private communication.

J. von zur Gathen and J. Gerhard (1999). Modern Computer Algebra. Cambridge University Press.

E. Weisstein (2000). Electronic mathWorld web site URL:
http://mathworld.wolfram.com/Antiprism.html

F. Winkler (1988). A p-adic approach to the computation of Gröbner bases. J. Symbolic Computation 6:287-304.

S. Wolfram (1999) The Mathematica Book (4th edition). Wolfram Media/Cambridge University Press.

R. Zippel (1993) Effective Polynomial Computation. Kluwer Academic Publishers.