##################################################################### ##################################################################### # # GROBNER PACKAGE FOR MAPLE V, RELEASE 4 # # by Albert Lin and Philippe Loustaunau, # George Mason University # # Minor modifications and corrections by David Cox, # Amherst College (04/21/99) # # Modifications to work under Maple V.4 # by C. D. Wensley, University of Wales, Bangor # (02/27/98) # # All comments, suggestions, and bug reports should be sent to: # # David A. Cox, Amherst College # # email: dac@cs.amherst.edu # ##################################################################### ##################################################################### # # HOW TO START # # Suppose that the file containing this Maple package is called # grobnerR4.sub. Then, if you start Maple from the directory # containing this file, you can load the file into Maple using # the command # # > read(`grobnerR4.sub`): # # If the file resides in another directory, then the appropriate # path needs to be included in the read statement. # ##################################################################### ##################################################################### # # THE THREE MAJOR COMMANDS # # This package has two major commands: # # 1) mxgb, which computes a reduced grobner basis # together with a matrix telling you how to transform the # original basis into the grobner basis. # mxgb is the union of three separate commands: # (i) basis_mxgb which computes a grobner basis which in # general is neither minimal nor reduced; # (ii) min_mxgb which converts a basis to a minimal basis; # (iii) red_mxgb which reduces a minimal basis. # # 2) div_alg, which computes the remainders AND quotients for # the division algorithm. # # These commands use a monomial order which is specified by the # user as a list of vectors. The names of the predefined Maple # term orders (plex and tdeg) should not be used. However, # three commands (lex, grlex, grevlex) are provided that make it # easier to use some of the more common term orders (see below). # # These commands are described in more detail below. # ##################################################################### ##################################################################### # # THE MATRIX GROBNER COMMAND # # The INPUT is: # mxgb( [f1,...,fs], [varlist], [vectorlist] ); # # [f1,...,fs] is the list of input polynomials. # # [varlist] is the list of variables that appear # in f1, ... , fs; e.g., [x,y,z]. # # [vectorlist] is the list of vectors that determines # the term order; e.g., if we have 3 # variables, say [x,y,z], then # [[1,0,0],[0,1,0],[0,0,1]] # gives the lexicographic order with # x < y < z. # # However, if [varlist] is [y,x,z], and # [vectorlist] is still # [[1,0,0],[0,1,0],[0,0,1]], # then the order is lexicographic but with # y < x < z. # # Any other term order can be obtained # using a list of vectors; e.g., if [x,y,z] # is the [varlist], then the total degree # lexicographic order with x>y>z is given by # [[1,1,1],[1,0,0],[0,1,0]]; # the total degree reverse lexicographic with # x>y>z is given by # [[1,1,1],[0,0,-1],[0,-1,0]]. # # If an ordering is used often, the user can # ``give'' it a name; for instance, # lexic:=[[1,0,0],[0,1,0],[0,0,1]]; # and then type "lexic" instead of the list of # vectors. However, the user should not use # the names of the predefined term orders in # Maple (plex for lexicographic, and tdeg for # total degree lexicographic). # # To make it easier to specify term orders, # we have provided three functions lex, # grlex and grevlex, whose input is the # number of variables and whose output is # the appropriate list of vectors. For # example, the three lists given above are # lex(3), grlex(3) and grevlex(3) respectively. # These commands are described below in more # detail. # # The OUTPUT is a list # [ [g1,...,gt], M ] # where [g1,...,gt] is the reduced grobner of and # M is an s x t matrix giving the grobner basis as a # combination of the fi's; i.e., # # [f1 ] # [f2 ] # [g1,...,gt] = M [...] # [...] # [fs ] # ##################################################################### # # EXAMPLE OF MATRIXGROBNER # # If the input is: # # > lexic := [ [1,0,0], [0,1,0], [0,0,1] ]; # > Var := [x,y,z]; # > Polys := [ x^2+y^2+z^2-4, x^2+2*y^2-5, x*z-1 ]; # > mxgb( Polys, Var, lexic ); # # then the output is the reduced grobner basis G of Polys # with respect to the lexicographic ordering with x > y > z, # and the matrix M of transformation such that G=M*Polys: # # [ [ y^2 - z^2 - 1, x + 2z^3 - 3z, 1/2 + z^4 - 3/2 z^2 ], coeff_mx ] # # and then the following input displays the transformation matrix: # # > M := gb[2]: print( M ); # # [ -1 1 0 ] # [ ] # [ 2z - z - x ] # [ ] # [ z^2 - 1/2 z^2 - 1/2 - 1/2 x *z ] # # The command > mxgb( Polys, Var, lex(3) ); would produce # the same output since there are three variables. # ##################################################################### ##################################################################### # # THE GROBNERBASIS COMMAND # # # The INPUT is: # > basis_gbmx( [f1,...,fs], [varlist], [vectorlist] ); # # where [f1,...,fs] is the list of input polynomials. # # [varlist] is the list of variables that appear # in f1, ... , fs; e.g., [x,y,z]. # # [vectorlist] is the list of vectors that determines # the term order; # # The OUTPUT is a 5-element list: # [ `isgbasis`, [g1,...,gt], [varlist], [vectorlist], coeff_mx ] # # where [g1,...,gt] is a grobner basis of , # which in general is neither minimal nor reduced, # # coeff_mx is the corresponding transformation matrix. # ##################################################################### # # EXAMPLE OF GROBNERBASIS # # If the input is: # # > lexic := [ [1,0,0], [0,1,0], [0,0,1] ]; # > Var := [x,y,z]; # > Polys := [ x^2+y^2+z^2-4, x^2+2*y^2-5, x*z-1 ]; # > bgb := basis_mxgb( Polys, Var, lexic ); # # then the output is a grobner basis G of Polys with respect to # lex order with x > y > z: # # [ isgbasis, # [x^2+y^2+z^2-4, x^2+2*y^2-5, x*z-1, # -y^2+z^2+1, x+2*z^3-3z, -1-2*z^4+3*z^2], # [x,y,z], [[1,0,0],[0,1,0],[0,0,1]], mxgb ; # # The trabsformation matrix may be displayed using: # # > bM := bgb[5]: print( bM ); # # The command > basis_mxgb( Polys, Var, lex(3) ); # would produce the same result since there are three variables. # ##################################################################### ##################################################################### # # THE DIVISION ALGORITHM COMMAND # # The INPUT is: # > div_alg( f, [f1,...,fs], [varlist], [vectorlist] ); # where # f is the polynomial to be divided # # [f1,...,fs] is the list of input polynomials. # # [varlist] is the list of variables that appear # in f1, ... , fs; e.g., [x,y,z]. # # [vectorlist] is the list of vectors that determines # the term order; # # The OUTPUT is # the remainder r and the quotients a1,...,as of the # division of f by f1,...,fs, # i.e. f = a1 f1 + ... + as fs + r. # ##################################################################### # # EXAMPLE OF DIVION ALGORITHM # # If the input is: # # > lexic := [ [1,0,0], [0,1,0], [0,0,1] ]; # > Var := [x,y,z]; # > Polys := [ x^2+y^2+z^2-4, x^2+2*y^2-5, x*z-1 ]; # > div_alg( x^3+2*y*z, Polys, Var, lexic ); # # then the output is the remainder r and quotients [a1,a2,a3] # of the division of x^3 + 2*y*z on division by Polys: # # [ -x*y^2 + 4*x + 2*y*z - z, [x, 0, -z] ] # # So the remainder is -x*y^2 + 4*x + 2*y*z - z # and the quotients are x, 0 and -z. # The command > div_alg( x^3+2*y*z, Polys, Var, lex(3) ); # would produce the same output since there are three variables. # ##################################################################### ##################################################################### # # THE LEX, GRLEX AND GREVLEX COMMANDS # # Since it can be awkward to manually type in the vector lists # that give the standard monomial orders, we have included three # functions to generate the appropriate vector list. # # These commands can be used as part of the input of the # matrixgrobner, grobnerbasis and div_alg commands. The # examples show how this is done. # # The INPUT is # lex(n), where n is the number of variables. # # The OUTPUT is # the vector list that specifies lexicographic order for # n variables. # # The INPUT is # grlex(n), where n is the number of variables. # # The OUTPUT is # the vector list that specifies graded lexicographic # order for n variables. # # The INPUT is # grevlex(n), where n is the number of variables. # # The OUTPUT is # the vector list that specifies graded reverse # lexicographic order for n variables. # ##################################################################### # # EXAMPLES OF LEX, GRLEX AND GREVLEX # # If the input is: # # lex(4); # # then the output is the vector list that specifies lex order # for 4 variables: # # [[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]] # # If the input is: # # grlex(4); # # then the output is the vector list that specifies lex order # for 4 variables: # # [[1,1,1,1],[1,0,0,0],[0,1,0,0],[0,0,1,0]] # # If the input is: # # grevlex(4); # # then the output is the vector list that specifies lex order # for 4 variables: # # [[1,1,1,1],[0,0,0,-1],[0,0,-1,0],[0,-1,0,0]] # # An important observation is that these commands, with the # appropriate number of variables, can be used as input for the # matrixgrobner, grobnerbasis and div_alg commands described # above. Thus, for example, the command # # div_alg( x^7+y^7+z^7, [x^2+y*z, y^3+y^2*z^2], [z,y,x], grlex(3) ); # # would compute the appropriate div_alg with respect grlex # order with z > y > x and find no reduction possible. # ##################################################################### ##################################################################### # # END OF DOCUMENTATION # ##################################################################### ##################################################################### # # THE PROGRAM FOR MATRIXGROBNER # ##################################################################### mxgb := proc( F, V, termorder ) local basislist, minlist, redlist; basislist := basis_mxgb( F, V, termorder ); minlist := min_mxgb( basislist ); redlist := red_mxgb( minlist ); redlist[2..3]; end; ##################################################################### ##################################################################### print( `Loading functions: lex, grlex, grevlex, div_alg, _leadingterm` ); print( ` mxgb (which comprises: basis_mxgb, min_mxgb, red_mxgb)` ); print( `and quot_mx (which expresses initial polys in terms of basis)` ); print( `WARNING: do not use the orders "plex" or "tdeg" with these functions!` ); ##################################################################### # # THE PROGRAMS THAT FOLLOW ARE THE SUBROUTINES NEEDED TO RUN mxgb # ##################################################################### ##################################################################### # # FINDING LEADING TERMS # # The INPUT is # _leadingterm(f,[varlist],[vectorlist]); # # The OUTPUT is the leading monomial of f with respect to # the term order given by the vectors in [vectorlist]. # ##################################################################### with(grobner,leadmon); _leadingterm := proc(f,V,U) local f1,h,mono,t,a,i,j,k,l,m,n,P,r,_Wa; f1:=expand(f); if type(f1,monomial) then leadmon(f1,V,plex); else n:=nops(U); _Wa:=convert(U,array); m:=nops(V); P:=array(1..m,1..nops(f1)); r:=array(1..n,1..nops(f1)); for i to nops(f1) do t[i]:=op(i,f1); od; for j to nops(f1) do for k to m do P[k,j]:=degree(t[j],V[k]); od; od; r:=evalm(&*(_Wa,P)); a:=1; for l from 2 to nops(f1) do for h to n do if r[h,l]r[h,a] then a:=l; break; fi; od; od; mono:=leadmon(simplify(op(a,f1)),V,plex); # print(r, a, mono); mono; fi; end; ##################################################################### # # FINDING REMAINDERS # # The INPUT is # div_alg( f, [f1,f2,...,fs], [varlist], [vectorlist]); # # The OUPUT is the remainder r and quotients a1,...,as of the # division of f by f1,...,fs; # i.e., f = a1 fa + ... + as fs + r. # ##################################################################### with(grobner,leadmon); div_alg := proc( g, Set, V, ord ) local h, i, a, v, f, lmf, lmg, b, c, d; a := array(1..nops(Set)); for b to nops(Set) do f := Set[b]; a[b] := 0; lmf[b] := _leadingterm( f, V, ord ); od; v := 0; h := g; while (h<>0) do lmg := _leadingterm( h, V, ord ); for i to nops(Set) do f := Set[i]; if (f <> 0) then d := degree( denom( simplify( lmg[2]/lmf[i][2] ) ) ); if d=0 then a[i] := simplify( a[i]+lmg[1]*lmg[2]/(lmf[i][1]*lmf[i][2]) ); h := simplify( h-f*lmg[1]*lmg[2]/(lmf[i][1]*lmf[i][2]) ); if h<>0 then i:=0 fi; lmg := _leadingterm( h, V, ord ); fi; else a[i]:=0; fi; od; v := simplify( v+lmg[1]*lmg[2] ); h := simplify( h-lmg[1]*lmg[2] ); od; a := convert( a, list ); c := [v,a]; end; ##################################################################### # # DIVISION ALGORITHM MATRIX # # The INPUT is # > quot_mx( [f1,f2,...,fs], [g1,...,gt], [varlist], [vectorlist]); # where the fj are in the ideal with gbasis [g1,...,gt], # so that the div_alg remainders are all zero. # # The OUPUT is the matrix of quotients # Q = [ [a11,...,a1t], ..., [as1,...,ast] ] # i.e. Q.[g1,...,gt] = [f1,...fs] # ##################################################################### quot_mx := proc( pols, bas, V, ord ) local np, nb, index1, index2, d, Q; np := nops( pols ); nb := nops( bas ); Q := matrix( np, nb, 0 ); for index1 to np do d := div_alg( pols[index1], bas, V, ord ); if ( d[1]<>0 ) then ERROR( `non-zero remainder` ); fi; for index2 to nb do Q[index1,index2] := d[2][index2]; od; od; Q; end; ##################################################################### # # FINDING A MINIMAL GROBNER BASIS # # The INPUT is # > min_mxgb( grobner basis list from basis_mxgb ) :- # # [ `isgbasis`, [g1,...,gt], [varlist], [vectorlist], coeff_mx ] # # The input [g1,...,gt] must be a grobner basis with respect # to the term order defined by [vectorlist]. # # The OUTPUT is a minimal grobner basis with respect to the # term order defined by [vectorlist]: # # [ `ismingbasis`, [g1,...,gt], [varlist], [vectorlist], coeff_mx ] # ##################################################################### with(linalg,submatrix); with(linalg,swaprow); with(linalg,coldim); min_mxgb := proc( gbrec ) local y, i, j, n, H, G, V, ord, coeff_mx, numcol, minrec; if not ( ( nargs = 1 ) and type( gbrec, list ) and ( nops( gbrec ) = 5 ) and type( gbrec[1], string ) and ( gbrec[1] = `isgbasis` ) ) then ERROR(`input not a matrixgrobner record` ); fi; G := gbrec[2]; V := gbrec[3]; ord := gbrec[4]; coeff_mx := gbrec[5]; numcol := coldim( coeff_mx ); n:=nops(G); H:=G; i:=1; while i<=n do; j:=1; while j<=n do; if i<>j then if divide(_leadingterm( H[i], V, ord )[2], _leadingterm( H[j], V, ord )[2]) then if (i=1) then H := H[2..n]; elif (i=n) then H := H[1..(n-1)]; else H := [ op( H[1..(i-1)] ), op( H[(i+1)..n] ) ]; fi; for y from i to n-1 do coeff_mx:=swaprow( coeff_mx, y, y+1 ); od; coeff_mx:=submatrix( coeff_mx, 1..n-1, 1..numcol ); n:=n-1; j:=n; i:=i-1; fi; fi; fi; j:=j+1; od; i:=i+1; od; minrec := [ `isminimal`, H, V, ord, coeff_mx ]; end; ##################################################################### # # FINDING A GROBNER BASIS # # The INPUT is # > basis_mxgb( [f1,..,fs], [varlist], [vectorlist] ); # # The OUTPUT is a 5-element list: # [ `isgbasis`, [g1,...,gt], [varlist], [vectorlist], coeff_mx ] # # where [g1,...,gt] is a grobner basis for # with respect to the term order defined by [vectorlist], # which in general is neither minimal nor reduced, # and coeff_mx is the corresponding transformation matrix. # #################################################################### with(grobner,leadmon); basis_mxgb := proc( H, Vl, ord) local x, y, setofpairs, z, p, l, F, r, i, n, lmf, lmg, ernie, u, T, Q, index1, index2, index3, index4, index5, subscript, coeff_tab, coeff_mx, numcol, gbrec; F:=[]; coeff_tab := table(sparse); for i to nops(H) do coeff_tab[i,i]:=1; od; F:=H; if nargs > 3 then print(F) fi; setofpairs:=[]; for index1 to nops(H)-1 do for index2 from index1+1 to nops(H) do setofpairs:=[ op(setofpairs), [index1, index2] ]; od; od; while nops(setofpairs) <>0 do numcol := nops(H); subscript := setofpairs[1]; setofpairs := [ op(2..nops(setofpairs), setofpairs) ]; lmf := _leadingterm( F[subscript[1]], Vl, ord ); lmg := _leadingterm( F[subscript[2]], Vl, ord ); if ( gcd( lmf[2], lmg[2] ) <> 1 ) then ernie := 0; for u to subscript[1]-1 do if divide(lcm(lmf[2],lmg[2]),_leadingterm(F[u],Vl,ord)[2]) then ernie:=1; break; fi; od; if ernie = 0 then p := lcm(lmf[2],lmg[2]); T := p*F[subscript[1]]/(lmf[1]*lmf[2]) - p*F[subscript[2]]/(lmg[1]*lmg[2]); Q := div_alg(simplify(T),F,Vl,ord); r := Q[1]; if ( r <> 0 ) then for index3 to nops(F) do setofpairs:=[op(setofpairs),[index3,nops(F)+1]]; od; F := [op(F),r]; if ( nargs > 3 ) then print(F) fi; n := nops(F); for l to n-1 do coeff_tab[n,l] := - op(l,Q[2]); od; coeff_tab[n,subscript[1]] := simplify( coeff_tab[n,subscript[1]] + p/(lmf[1]*lmf[2]) ); coeff_tab[n,subscript[2]] := simplify( coeff_tab[n,subscript[2]] - p/(lmg[1]*lmg[2]) ); fi; fi; fi; od; if ( nops(F) > nops(H) ) then for x from numcol+2 to n do for y to numcol do for z from numcol+1 to x-1 do coeff_tab[x,y] := simplify( coeff_tab[x,y] + coeff_tab[x,z]*coeff_tab[z,y] ); od; od; od; fi; coeff_mx := array( 1..nops(F), 1..numcol ); for index4 to nops(F) do for index5 to numcol do coeff_mx[ index4, index5 ] := coeff_tab[ index4, index5 ]; od; od; gbrec := [ `isgbasis`, F, Vl, ord, coeff_mx ]; end; #################################################################### # # FINDING A REDUCED GROBNER BASIS # # The INPUT is # > red_mxgb( minimal Grobner basis list from min_mxgb ); # # The OUTPUT is a 5-element list: # [ `isredgbasis`, [g1,...,gt], [varlist], [vectorlist], coeff_mx ] # # where [g1,...,gt] is the reduced Grobner basis for # with respect to the term order defined by [vectorlist], # and coeff_mx is the corresponding transformation matrix. # #################################################################### with(linalg,multiply); with(linalg,matrix); with(linalg,submatrix); red_mxgb := proc( minrec ) local x, Ca, t, i, j, k, l, m, n, a, o, Da, Db, r, J, J0, J1, V, ord, coeff_mx, index1, index2, redmx, grobmx, grob1; if not ( ( nargs = 1 ) and type( minrec[1], string ) and ( minrec[1] = `isminimal` ) ) then ERROR(`input not a minimal matrixgrobner record` ); fi; J := minrec[2]; V := minrec[3]; ord := minrec[4]; grobmx := minrec[5]; n := nops(J); m := coldim( grobmx ); grob1 := submatrix( grobmx, 1..n, 1..m ); i:=1; Ca:=table(sparse); if n<>1 then while (i <= n) do # Da:=[op(1..i-1,J)]; # Db:=[op(i+1..n,J)]; if (i=1) then Da := [ ]; Db := J[2..n]; elif (i=n) then Da := J[1..(n-1)]; Db := [ ]; else Da := J[1..(i-1)]; Db := J[(i+1)..n]; fi; r := div_alg( J[i], [ op(Da), op(Db) ], V, ord ); J := [ op(Da), r[1], op(Db) ]; for j to n do for k to i-1 do Ca[i,j]:=simplify(Ca[i,j]-r[2][k]*Ca[k,j]); od; od; for l from i+1 to n do Ca[i,l]:=simplify(Ca[i,l]-r[2][l-1]); od; Ca[i,i]:=simplify(Ca[i,i]+1); i:=i+1; od; else Ca[1,1]:=1; fi; redmx := matrix( n, n, 0 ); for index1 to n do for index2 to n do redmx[ index1, index2 ] := Ca[ index1, index2 ]; od; od; J0 := [ ]; J1 := [ ]; for a to n do if (J[a]<>0) then J0 := [ op(J0), a ]; J1 := [ op(J1), J[a] ]; fi; od; redmx := submatrix( redmx, J0, J0 ); grob1 := submatrix( grob1, J0, 1..m ); n := nops( J0 ); J := J1; for a to n do x := _leadingterm( J[a], V, ord )[1]; o := J[a]/x; for t to n do redmx[a,t] := redmx[a,t]/x; od; if (a=1) then J := [ o, op(2..n,J) ]; elif (a=n) then J := [ op(1..n-1,J), o ]; else J := [ op(1..a-1,J), o, op(a+1..n,J) ]; fi; od; coeff_mx := multiply( redmx, grob1 ); [ `isreduced`, J, coeff_mx ]; end; ##################################################################### # # GENERATING VECTOR LISTS FOR LEX, GRLEX AND GREVLEX # # The INPUT is lex(n), grlex(n) or grevlex(n), and the OUTPUT is # the vector list that generates the appropriate order when # there are n variables. # ##################################################################### lex := proc(n) local u,i,j; if n=1 then [[1]]; else u := [0 $ i=1..j-1,1,0 $ i=j+1..n]; [[1,0 $ i=1..n-1],u $ j=2..n-1,[0 $ i=1..n-1,1]] fi; end; grlex := proc(n) local u,i,j; if n=1 then [[1]]; else u := [0 $ i=1..j-1,1,0 $ i=j+1..n]; [[1 $ i=1..n],u $ j=1..n-1] fi; end; grevlex := proc(n) local u,i,j; if n=1 then [[1]]; else u := [0 $ i=1..n-j-1,-1,0 $ i=1..j]; [[1 $ i=1..n],[0 $ i=1..n-1,-1], u $ j=1..n-2] fi; end;