Maple 11 has been out for a while now so hopefully people have it. I thought I would write a short post detailing some of what was done in the area of Groebner bases. If you run the examples in Maple 10 and Maple 11, I would appreciate it if you could post the times and the specifications of your computer.
FGb
The most heavily marketed feature is that J.C. Faugere's FGb program has been integrated into Maple 11. This is a C library that is among the very fastest programs for computing Groebner bases. Magma can be faster in some situations and FGb in others due to different choices of strategy (homogenization, etc), but really, this is about the fastest software you can get right now. It supports computations with rational number coefficients (using a modular algorithm) and the integers mod p where p is less than half the wordsize of your computer. That is p < 65536 for a 32-bit machine or p < 2^32 for a 64-bit machine. FGb is called to compute "total degree" Groebner bases (tdeg in Maple) by the Groebner[Basis] command automatically if possible. FGb also supports elimination orders (lexdeg in Maple) but the default of the Groebner[Basis] command is to use a Groebner basis conversion algorithm for those. Here is an example showing a computation with FGb that you can try in both Maple 10 and Maple 11. Many problems that could not be computed in Maple 10 can now be computed easily with FGb.
cyclic6 := [x*y*z+y*z*t+z*t*u+t*u*v+u*v*x+v*x*y, x*y*z*t+y*z*t*u+z*t*u*v+t*u*v*x+u*v*x*y+v*x*y*z, x*y*z*t*u*v-1, x+y+z+t+u+v, x*y+y*z+z*t+t*u+u*v+x*v, x*y*z*t*u+y*z*t*u*v+z*t*u*v*x+t*u*v*x*y+u*v*x*y*z+v*x*y*z*t]:
TIMER := time(): G := Groebner[Basis](cyclic6, tdeg(x,y,z,t,u,v)): time()-TIMER;
Faster Groebner basis conversions (Groebner Walk and FGLM)
The Groebner basis conversion algorithms of Maple 11 were improved. These algorithms take Groebner bases for one term order and use them to compute a basis for another order. In particular, the Groebner[Basis] command uses them to compute lexicographic Groebner bases because these are often too hard to compute directly. There are some modest speedups in the FGLM algorithm, which uses linear algebra to compute the new basis using linear dependencies on the monomials with respect to the old basis. We still do not use a modular algorithm for this (Maple needs sparse modular linear algebra) however you can see that there has been some improvement from the example below. We can take the total degree basis computed for cyclic6 computed above and compute a lexicographic basis from it. The Groebner[Basis] command will do this automatically (in both Maple 10 and 11) when you ask for the lex basis.
TIMER := time(): L := Groebner[FGLM](G, tdeg(x,y,z,t,u,v), plex(x,y,z,t,u,v)): time()-TIMER;
The Groebner Walk algorithm received more substantial improvements because we missed some critical optimizations when we implemented the algorithm for Maple 10. This algorithm converts Groebner bases by computing and applying small transformations at every step of a "walk" across the Groebner cone.
circles := [(1-z)^2+(w+m-1)^2-m^2, x^2+y^2-m^2, (1-x-m)^2+w^2-m^2, (1-m-1/2*z)^2+(w-1/2-1/2*y-1/2*m)^2-m^2, 1/4*z^2+(1/2-1/2*y-1/2*m)^2-m^2]: G := Groebner[Basis](circles, tdeg(x,y,z,w,m)):
TIMER := time(): L := Groebner[Walk](G, tdeg(x,y,z,w,m), plex(x,y,z,w,m)): time()-TIMER;
Groebner[Solve]
The Groebner[Solve] command is a version of the Buchberger algorithm that factors intermediate polynomials and splits the polynomial system. It first appeared as the gsolve command in the old grobner package by Stephen Czapor, who had a lot of tricks to make it run fast. Later in Maple V R5 the Groebner package replaced grobner, but it was never able to duplicate the speed of this old code. We have tried, and if anyone out there has Maple V R4 I would appreciate it if you could post the times to solve the following systems using grobner[gsolve] and Maple 11.
caprasse := [y^2*z+2*x*y*t-2*x-z, 2*y*z*t+x*t^2-x-2*z, -x^3*z+4*x*y^2*z+4*x^2*y*t+2*y^3*t+4*x^2-10*y^2+4*x*z-10*y*t+2, -x*z^3+4*y*z^2*t+4*x*z*t^2+2*y*t^3+4*x*z+4*z^2-10*y*t-10*t^2+2]:
TIMER := time(): S := Groebner[Solve](caprasse, {x,y,z,t}): time()-TIMER;
discriminant4 := {y^2*z^2*t^2-4*x*z^3*t^2-4*y^3*t^3 +18*x*y*z*t^3-27*x^2*t^4-4*y^2*z^3+16*x*z^4+18*y^3*z*t -80*x*y*z^2*t-6*x*y^2*t^2+144*x^2*z*t^2-27*y^4 +144*x*y^2*z-128*x^2*z^2-192*x^2*y*t+256*x^3, -4*z^3*t^2+18*y*z*t^3-54*x*t^4+16*z^4-80*y*z^2*t-6*y^2*t^2 +288*x*z*t^2+144*y^2*z-256*x*z^2-384*x*y*t+768*x^2, 2*y*z^2*t^2-12*y^2*t^3+18*x*z*t^3-8*y*z^3+54*y^2*z*t -80*x*z^2*t-12*x*y*t^2-108*y^3+288*x*y*z-192*x^2*t, 2*y^2*z*t^2-12*x*z^2*t^2+18*x*y*t^3-12*y^2*z^2+64*x*z^3 +18*y^3*t-160*x*y*z*t+144*x^2*t^2+144*x*y^2-256*x^2*z, 2*y^2*z^2*t-8*x*z^3*t-12*y^3*t^2+54*x*y*z*t^2-108*x^2*t^3 +18*y^3*z-80*x*y*z^2-12*x*y^2*t+288*x^2*z*t-192*x^2*y}:
TIMER := time(): S := Groebner[Solve](discriminant4, {x,y,z,t}): time()-TIMER;
gerdt := {2*v*x-u*x+2*t*x, -2*z*x^2+10*y*x^2-20*x^3+7*t*z -35*t*y+70*t*x, 2*v*x^2-2*u*x^2+6*t*x^2-7*v*t+7*u*t-21*t^2, 2*z*y*x-6*y^2*x-z*x^2+13*y*x^2-5*x^3+14*v*x-28*t*x, z^2*x-3*z*y*x+5*z*x^2+14*v*x-28*t*x, -2*u*z*x-2*t*z*x +4*v*y*x+6*u*y*x-2*t*y*x-16*v*x^2-10*u*x^2+22*t*x^2+42*w*x, 28*v*z*x+8*u*z*x-20*t*z*x-88*v*y*x-24*u*y*x+68*t*y*x +156*v*x^2+40*u*x^2-132*t*x^2-252*w*x, -4*v*u*x+10*v*t*x +8*u*t*x-20*t^2*x+12*w*z*x-30*w*y*x+15*w*x^2, -2*v^2*x +v*u*x+2*v*t*x-2*u*t*x+4*t^2*x-6*w*z*x+12*w*y*x-6*w*x^2, 8*w*v*x-4*w*u*x+8*w*t*x, -2*y^3+4*z*y*x+5*y^2*x-6*z*x^2 -7*y*x^2+15*x^3+42*v*y-14*u*y-63*v*x+21*u*x-42*t*x+147*w, -9*z*x^3+45*y*x^3-135*x^4+14*u*y^2-14*t*y^2-28*u*z*x+70*t*z*x -14*u*y*x-196*t*y*x+28*u*x^2+602*t*x^2-294*v*u+98*u^2+294*v*t -98*u*t-147*w*z+735*w*y-2205*w*x, 6*v*x^3-9*u*x^3+36*t*x^3 -14*w*y^2-28*v*t*x+42*u*t*x-168*t^2*x+28*w*z*x+14*w*y*x -28*w*x^2+392*w*v-245*w*u+588*w*t}:
TIMER := time(): S := Groebner[Solve](gerdt, {x,y,z,t,u,v,w}): time()-TIMER;
gonnet := [a5*b0+a5*b1+a4*b3+a3*b4+2*a5*b4+a0*b5+2*a4*b5+b5+c5, a3*b0+a5*b0+a3*b1+a5*b1+a0*b3+a4*b3+a3*b4+a5*b4+a0*b5+a4*b5+b3+b5+c3+c5-1, a0*b0+a4*b0+a0*b1+a4*b1+a3*b2+a5*b2+a2*b3+a0*b4+a4*b4+a2*b5+b0+b1+b4+c0+c1+c4, a0*b1+a4*b1+a3*b2+a2*b3+b0+2*b1+b4+c1, a4*b0+a4*b1+a5*b2+a0*b4+2*a4*b4+a2*b5+b4+c4, a3*b0+2*a3*b1+a5*b1+a0*b3+a4*b3+a3*b4+2*b3+b5+c3, a2*b0+a2*b1+a0*b2+a4*b2+a2*b4+b2+c2, a5*b5, a5*b4+a4*b5, a4*b4, a5*b3+a3*b5, a5*b3+a3*b5+2*a5*b5, a3*b3+a5*b3+a3*b5+a5*b5, 2*a3*b3+a5*b3+a3*b5, a4*b2+a2*b4, a2*b2, a5*b1+a4*b3+a3*b4+b5, a4*b1+b4, a2*b1+b2]:
TIMER := time(): S := Groebner[Solve](gonnet, indets(gonnet)): time()-TIMER;
There are some other improvements in Groebner[Solve] as well (do not try in Maple 10). For example, instead of working with lists of and term orders, it is useful to have an object representing the system which keeps track of everything, such what Groebner bases have been computed. The PolynomialIdeals package (added in Maple 10) does this, and the Groebner[Solve] command can now use some of this extra information. In particular, it will check to see if a lexicographic basis has been computed and if so it will start from it. Here is an example of a fairly large system that is factored easily once a lexicographic basis is known.
# create an ideal
with(PolynomialIdeals):
filter9 := < m2*m4*m6-1/100, a*m4*b-7/500, a^2+m1^2-2/25, b^2+m7^2-37/50, m3^2+m5^2+m4^2+m2^2+m6^2-9401/10000, m4^2*m6^2+m2^2*m4^2+m3^2*m6^2+m2^2*m5^2+m3^2*m5^2+m2^2*m6^2-38589/1000000, m1*m3*m5*m7-m6*m1*m3*b+m2*m6*a*b-m2*a*m5*m7+81/10000, -m1*m2*m3*m4*b-a*m4*m5*m6*m7+a*m4*b*m6^2+a*m2^2*m4*b+39/25000, -m4^2*m7^2-m3^2*m7^2+2*m5*m6*b*m7-m2^2*m7^2-m5^2*m7^2 -m3^2*b^2-b^2*m6^2-m2^2*b^2+27173/40000 >:
L := Groebner[Basis](filter9, 'tord', order=plex):
IdealInfo:-KnownGroebnerBases(filter9);
infolevel[Groebner[Solve]] := 4:
TIMER := time(): S := Groebner[Solve](filter9): time()-TIMER;
Note that we used angled brackets to construct an ideal, and the output of the Groebner[Solve] command is a set of ideals. We have also used a new syntax that is available in the Groebner[Basis] command where it will choose the term order for you or order the variables heuristically. It returns a known Groebner basis if applicable, avoiding a computation.
Maple F4
Some people are aware of another implementation of F4 in the Maple application center. I wrote this implementation and it is not very good. It uses obscene amounts of memory and is slow. I wanted a modular algorithm but Maple doesn't have sparse modular linear algebra so I used dense linear algebra which can't scale. Maple 11 contains yet another implementation of F4, but this one is sparse. Why bother when there is FGb ? Well, it is for all of the cases that FGb doesn't handle. Large primes, weird term orders, polynomial coefficients, and non-commutative algebras are all supported by this code, and the sparse linear algebra is done in Maple using hash tables. Unfortunately we lost the modular algorithm but there is only so much one can do. The F4 algorithm is intrinsically faster than the Buchberger algorithm by an order of magnitude, and in Maple 11, it is always used by default. Here is an example of a non-commutative Weyl-algebra from the Maple test suite with polynomial coefficients in {a,b,c[1],...,c[n]}.
n := 3: S := add(x[i]*D[i],i=1..n):
A := Ore_algebra:-skew_algebra(comm={a,b,seq(c[i],i=1..n)}, seq(diff=[D[i],x[i]],i=1..n),polynom={seq(x[i],i=1..n)}):
P := Ore_algebra:-skew_product(S+a,S+b,A):
for i to n do l[i] := Ore_algebra:-skew_product(D[i],x[i]*D[i]+c[i]-1,A)-P end do:
M := Groebner:-MonomialOrder(A,lexdeg([seq(D[i],i=1..n)],[seq(x[i],i=1..n)])):
TIMER := time(): G := Groebner[Basis]([seq(l[i],i=1..n)],M)): time()-TIMER:
If we plug in numbers for the parameters and use an easier term order we can do the next problem as well. We use infolevel[F4] to print the sizes of the matrices, which are not that big. Most of the time is spent in the symbolic preprocessing, computing non-commutative products.
n := 4: S := add(x[i]*D[i],i=1..n):
A := Ore_algebra:-skew_algebra(comm={a,b,seq(c[i],i=1..n)}, seq(diff=[D[i],x[i]],i=1..n),polynom={seq(x[i],i=1..n)}):
P := Ore_algebra:-skew_product(S+a,S+b,A):
for i to n do l[i] := Ore_algebra:-skew_product(D[i],x[i]*D[i]+c[i]-1,A)-P end do:
M := Groebner:-MonomialOrder(A,tdeg(seq(D[i],i=1..n),seq(x[i],i=1..n))):
L := subs({a=1,b=2,seq(c[i]=-i, i=1..n)}, [seq(l[i],i=1..n)]):
infolevel[F4] := 5:
TIMER := time(): G := Groebner[Basis](L,M): time()-TIMER;
Fast Division of Multiple Polynomials
Maple 11 contains a new algorithm for multivariate division whose complexity is linear in the total number of terms encountered, regardless of whether one divides one polynomial or several polynomials by the Groebner basis. The standard division algorithm divides f by g by repeatedly subtracting multiples of g. If the quotient has n terms and g has m terms this uses O(n^2*m) operations to merge n*m terms when the polynomials are sparse. Geobuckets and other divide and conquer strategies reduce this to O(n*m*log(n)), but the algorithm in Maple 11 is actually O(n*m), however it is written in Maple code so it is "slow". To use the algorithm, make the first argument to Groebner[NormalForm] a list:
f := randpoly([x,y,z,t],degree=50,terms=1000):
G := [x^3+y^2+z+1, y^2+2*z+t, z^4-t^3+t+1, t^2-3]:
r := Groebner[NormalForm]([f], G, tdeg(x,y,z,t)):
Notice that the output is also a list. Of course you can input a list of several polynomials and get a list of several remainders. This code is now used by simplify/siderels as well, but for technical reasons it is called separately on each element if you give simplify a list. This code is also used by the Groebner[MultiplicationMatrix] command, which is now quite a bit faster. We compute the remainders of all the monomials (efficiently) at once.
G := Groebner[Basis](cyclic6, tdeg(x,y,z,t,u,v)):
ns, rv := Groebner[NormalSet](G, tdeg(x,y,z,t,u,v)):
TIMER := time(): M := Groebner[MultiplicationMatrix](v, ns, rv, G, tdeg(x,y,z,t,u,v)): time()-TIMER;
The first argument to MultiplicationMatrix can also now be a polynomial f as well. The eigenvalues of this matrix are the values of f on the variety (at the solutions) of the system. This post getting too long so I am going to split it up. In a future post I will try to focus a bit more on how to use the package including the new commands.
