The Groebner Package

by Frederic Chyzak, Algorithms Project, INRIA, France

Frederic.Chyzak@inria.fr

The Groebner package extends the functionality of the grobner package into several directions, namely:

calculations over more complicated ground fields are now possible with a single package;

new term orders are available, in view of new applications and improved efficiency;

calculations of common invariants of ideals and varieties are available;

facilities for change of orderings are available.

> with(Groebner);

Some syntactical differences exist between the grobner and the Groebner packages. These differences are motivated by the extension of functionality, particularly in the number of term orderings. Specifically, a call

grobner[gbasis](G,[x,y],plex);

becomes

Groebner[gbasis](G,plex(x,y));

while the call

grobner[gbasis](G,{x,y});

with being the standard has no counterpart in Groebner. For instance,

> grobner[gbasis]([x^2-2*x*z+5,x*y^2+y*z^3,3*y^2-8*z^3],[y,x,z],plex);

becomes

> Groebner[gbasis]([x^2-2*x*z+5,x*y^2+y*z^3,3*y^2-8*z^3],plex(y,x,z));

Moreover, for complicated calculations, the Groebner package requires data structures that are created by the Ore_algebra package.

> with(Ore_algebra);

Ground fields

The Groebner package supports transcendental and algebraic extensions of the field of rational numbers and finite fields. This includes algebraic number fields. Here is a sample computation with complex numbers in .

> G:=[x^2+y^2-I,sqrt(2)*y+sqrt(2)*(x-sqrt(2))-1];

> GB:=gbasis(G,tdeg(x,y));

For computations with modular integers, one has to declare the corresponding algebra of polynomials. This is dealt with by the Ore_algebra package. (Ore polynomials extend commutative polynomials, and the general mechanism to deal with general algebras is in this package.)

> A:=poly_algebra(x,y,z,characteristic=5);

> G:=[x^2-2*x*z+5,x*y^2+y*z^3,3*y^2-8*z^3]:

> A[normalizer](G);

> T:=termorder(A,plex(x,y,z)):

> gbasis(G,T);

Term orderings

The grobner package allows two types of term orderings:

, for "efficient" normal forms,

, for "slow" triangularizations (system solving).

There was a need for specialized orderings for the purpose of efficient elimination and for term orderings defined by a matrix for the purpose of generality. The Groebner package implements both orderings, the and types respectively, together with a more practical for weighted orderings. User-defined orderings are also available.

As an example, here is a calculation that is performed far faster using than using . Assume we want to find the optima of

> phi:=x^3+2*x*y*z-z^2;

on the sphere

> const:=x^2+y^2+z^2-1;

By Lagrange method, we compute

> G:={seq(diff(phi,v)-diff(const,v)*t,v=[x,y,z]),const};

and we eliminate by the use of an appropriate term order

> GB:=gbasis(G,lexdeg([t],[x,y,z]));

> remove(has,GB,t);

In the general case, we would have to compute numerical estimates and substitute into . Here, we can go further and compute all the possible values for , by elimination of and .

> GB2:=gbasis(%,lexdeg([x,y],[z]));

> op(remove(has,GB2,{x,y}));

> sol_z:=solve(%,z);

Substituting into the system yields the possible tuples where the optima take place.

> map(op,[seq(map(`union`,[solve(convert(subs(z=i,GB2),set),{x,y})],{z=i}),i=sol_z)]);

Finally, we get the minimum and maximum values.

> min(op(map(subs,%,phi))),max(op(map(subs,%,phi)));

Similar calculations using need more than an hour.

As another example, the Groebner package allows for the computation of standard bases. We define the leading term of a polynomial as its minimum term instead of its maximum term by using a term order rather than a term order:

> L:=1+x+y+z+x^2+y^2+z^2;

> leadterm(L,plex(x,y,z));

> leadterm(L,plex_min(x,y,z));

Then Groebner bases can be computed for homogeneous ideals. Here is an example in .

> G:=[y*w^2-z^3,x*w^3-z^4];

Compare the usual Groebner basis

> gbasis(G,plex(w,x,y,z));

which allows for divisions by decreasing powers:

> reduce(w^5*x^3*y^2*z,%,plex(w,x,y,z));

and the standard basis

> gbasis(G,plex_min(w,x,y,z));

which allows for divisions by increasing powers:

> reduce(w^5*x^3*y^2*z,%,plex_min(w,x,y,z));

Ideal invariants

The dimension of a variety (zero for isolated points, one for lines, two for planes) is the Hilbert dimension of the corresponding annihilating ideal. The Groebner package allows one to compute this Hilbert dimension, together with the Hilbert polynomial and Hilbert series.

As a first example, here is a zero-dimensional ideal, corresponding to finitely many points in the variety.

> G:=[3*y^2-8*z^3,x^2-2*z*x+5,3*y^3+8*x*y^2]:

> hilbertdim(G,tdeg(y,x,z));

> hilbertseries(G,tdeg(y,x,z),s);

> subs(s=1,%);

> hilbertpoly(G,tdeg(y,x,z),s);

As another example, here is an ideal of positive dimension.

> G:=[x^3*y^2+3*x^2*y^2+y^3+1]:

> hilbertdim(G,tdeg(x,y));

> hilbertseries(G,tdeg(x,y),s);

> series(%,s,10);

> hilbertpoly(G,tdeg(y,x),s);

FGLM algorithm

The FGLM algorithm was introduced to perform changes of orderings for zero-dimensional ideals.

> G:=[x+y+z,x*y+y*z+z*x,x*y*z-1]:

> GB:=gbasis(G,tdeg(x,y,z));

> hilbertdim(GB,tdeg(x,y,z));

The implementation is a fully user-customizable generalization:

1. choice of the normal form:

> NFproc:=proc(t,NF,T) local r,s; global GB;
r:=reduce(t,GB,tdeg(x,y,z),s);
NF[t]:=r/s
end:

2. computation of dependencies between those normal forms:

> FDproc:=proc(M,NF)
local ind_lst,t,eta,sol,zero,sys,rel;
ind_lst:=map(op,[indices(NF)]);
for t in M do
ind_lst:=remove(divide,ind_lst,t)
od;
zero:=expand(numer(normal(add(eta[t]*NF[t],
t=ind_lst))));
sys:={coeffs(zero,{x,y,z})};
sol:=[readlib(`solve/linear`)(sys,
{seq(eta[t],t=ind_lst)})];
rel:=subs(sol[1],add(eta[t]*t,t=ind_lst));
`if`(rel=0,FAIL,
collect(primpart(numer(subs(map(n->n=1,
map2(op,1,select(evalb,sol[1]))),
rel)),{x,y,z}),{x,y,z},
distributed,factor))
end:

3. choice of the termination condition:

> TMproc:=proc(border,monoideal,TOrd)
border<>[]
end:

> A:=poly_algebra(x,y,z):

> T:=termorder(A,plex(z,y,x)):

> FGLM:=[fglm(NFproc,FDproc,TMproc,T)]:

> GB2:=map(op,[entries(FGLM[2])]);

Groebner bases for skew algebras

The functionality of the Groebner package is also available for skew polynomial inputs. Whereas in the commutative case the ideals are two-sided, in the skew case we focus on left-sided ideal. This corresponds to skew polynomials acting on functions on the left.

As a first application, we derive a proof that Legendre polynomials satisfy

> Sum(orthopoly[P](n,z)*u^n,n=0..infinity)=1/sqrt(1-2*z*u+u^2);

Legendre polynomials are defined by the following set of operators in the Ore algebra built over a differential operator and a shift operator .

> G:={(1-z^2)*Dz^2-2*z*Dz+n*(n+1),
(n+2)*Sn^2-(2*n+3)*z*Sn+(n+1),
(1-z^2)*Dz*Sn+(n+1)*z*Sn-(n+1)};

It follows that the summand in the identity we want to prove satisfies

> G:=map(numer@normal,subs(Sn=Sn/u,G)) union {u*Du-n};

in the new algebra . (Since the index of summation is , we consider polynomials in this variable to allow for its elimination.)

> A:=skew_algebra(diff=[Dz,z],diff=[Du,u],shift=[Sn,n],polynom=n):

> T:=termorder(A,lexdeg([n],[Dz,Du,Sn])):

> GB:=gbasis(G,T);

The previous equations are satisfied by the summand. Setting in them yields equations satisfied by the sum. A final Groebner basis calculation returns these equations

> GB:=gbasis(subs(Sn=1,remove(has,GB,n)),T):

> collect(GB,{Dz,Du},distributed,factor);

Solving this simple system of PDE's yields the announced closed form.

The same method applies to perform definite integration. Here, we show that the integral

> Int((1+x*y+2*y^2)*exp(x^2*y^2/(1+2*y^2))/y^n/y/(1+2*y^2)^(2/3),
y=-infinity..infinity);

satisfies a recurrence of order 6. With a more refined method also based on Groebner basis calculations, one could obtain a second order recurrence so as to get a proof that the integral is an integral representation for Hermite polynomials. These polynomials are known to Maple as .

They satisfy the system

> G:={Sn*y-1,
3*(1+x*y+2*y^2)*(1+2*y^2)^2*y*Dy+24*n*y^6+8*y^6+12*n*x*y^5+16*x*y^5
-12*x^2*y^4+20*y^4+36*n*y^4+12*n*x*y^3+8*x*y^3-6*x^3*y^3+14*y^2
-6*x^2*y^2+18*n*y^2+3*n*x*y+3+3*n};

in the Ore algebra

> A:=skew_algebra(diff=[Dy,y],shift=[Sn,n],comm=x,polynom=y):

built on a differential operator and a shift operator . (Since the variable of integration is , we consider polynomials in this variable to allow for its elimination.)

The elimination is obtained by using a suitable term order.

> T:=termorder(A,lexdeg([y],[Dy,Sn])):

> GB:=gbasis(G,T);

The above equations are also satisfied by the integrand. Setting in them yields an equation satisfied by the integral.

> collect(primpart(subs(Dy=0,op(remove(has,GB,y))),Sn),Sn,factor);

The order is easily reduced by 1 to get the announced 6 order equation.

> collect(subs(n=n-1,%/Sn),Sn,factor);

Groebner bases for modules

Groebner bases for modules are available in the case of K[x[1],...,x[p]]-modules of the form K[x[1],...,x[p],s[1],...,s[q]], where the 's and the 's need not commute. Those module structures are obtained by disallowing the multiplications by the 's along all calculations.

As an example of a Groebner basis for the -module , we simulaneously compute a Groebner basis for an ideal generated by , and we compute each polynomial in as a combination of the initial polynomials.

The initial basis is:

> G:=[x+y+z,x*y+y*z+z*x,x*y*z-1];

Let us introduce the module generated by

> MG:=[seq(s^(i-1)+G[i]*s^nops(G),i=1..nops(G))];

in the polynomial algebra

> A:=poly_algebra(x,y,z,s):

To refrain from any multiplication by in the calculations, we use a third argument in the call to termorder.

> T:=termorder(A,lexdeg([s],[x,y,z]),[s]):

> MGB:=collect(gbasis(MG,T),s);

Up to 0's, the Groebner basis for is

> GB:=map(coeff,MGB,s,3);

as checked by

> gbasis(G,tdeg(x,y,z));

Furthermore, we get the following relations between the elements of :

> seq(coeff(MGB[i],s,3)=subs(dummy='G',add(coeff(MGB[i],s,j)*dummy[j+1],j=0..2)),i=1..3);

And we get the Groebner basis for in terms of the elements of :

> seq(coeff(MGB[i],s,3)=subs(dummy='G',add(coeff(MGB[i],s,j)*dummy[j+1],j=0..2)),i=4..6);

or, in matrix notation:

> matrix(3,1,[seq(coeff(MGB[i],s,3),i=4..6)])=
matrix(3,3,[seq(subsop(seq(j+1=coeff(MGB[i],s,j),j=0..2),[0,0,0]),
i=4..6)])&*matrix(3,1,G);