with(grobner); ##################################################################### ##################################################################### # # GROBNER PACKAGE FOR MAPLE V, RELEASE 3 # Modified by Allen Bell for use in Math 841 at UW-Milwaukee # ##################################################################### ##################################################################### # # HOW TO START # # 1) Copy this file under your favorite name, such as # grobner.sub or 841.sub, into the directory from which you # start Maple. # 2) Every time you start a Maple session and you want # to use this program you will have to compile the # program. # 3) To compile the program, first open a Maple session, # and at the Maple prompt, type # read `grobner.sub`: # or read `841.sub`: # (if stored on floppy, use read `a:\841.sub`:) # # NOTE: This package does not work on Release 2 of Maple V # ##################################################################### ##################################################################### # # THE THREE MAJOR COMMANDS # # This package has three major commands: # # 1) matrixgrobner, which computes a reduced grobner basis # together with a matrix telling you how to transform the # original basis into the grobner basis. # 2) grobnerbasis, which computes a grobner basis which in # general is neither minimal nor reduced. # 3) remainder, 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. # # These six commands are described in more detail below. # ##################################################################### # OTHER COMMANDS # # a) sp(f,g,varlist,termorder), which computes the # S-polynomial of f & g # b) _leadingterm(f,varlist,termorder), which computes # the leading term of g # c) _minimalgb(G,varlist,termorder), which computes a # minimal Groebner basis from a Groebner # basis G (given as a list) # d) _reducegrobner(G,varlist,termorder), which computes a # reduced Groebner basis from a minimal Groebner # basis G (given as a list) # ##################################################################### # # THE MATRIXGROBNER COMMAND # # The INPUT is: # matrixgrobner([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 # the reduced grobner basis [g1,...,gt] of and # an s x t matrix M 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]; # matrixgrobner(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: # # [ -1 1 0 ] # [ ] # [ 2 z - z - x ] # [ ] # [ z ^2 - 1/2 z^2 - 1/2 - 1/2 x *z ] # # # # [y ^2 - z^2 - 1, x + 2 z^3 - 3 z, 1/2 + z ^4 - 3/2 z ^2 ] # # The command "matrixgrobner(Polys,Var,lex(3))" would produce # the same output since there are three variables. # ##################################################################### ##################################################################### # # THE GROBNERBASIS COMMAND # # # The INPUT is: # grobnerbasis([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; # # The OUTPUT is # A grobner basis [g1,...,gt] of , which in # general is neither minimal nor reduced. # ##################################################################### # # 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]; # grobnerbasis(Polys,Var,lexic); # # then the output is a grobner basis G of Polys with respect to # lex order with x > y > z: # # [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] # # The command "grobnerbasis(Polys,Var,lex(3))" would produce the # same result since there are three variables. # ##################################################################### ##################################################################### # # THE REMAINDER COMMAND # # The INPUT is: # remainder(f,[f1,...,fs],[varlist],[vectorlist]); # # 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 REMAINDER # # 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]; # remainder(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 "remainder(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 remainder 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 remainder commands described # above. Thus, for example, the command # # remainder(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 remainder with respect grlex # order with z>y>x. # ##################################################################### ##################################################################### # # END OF DOCUMENTATION # ##################################################################### ##################################################################### # # THE PROGRAM FOR MATRIXGROBNER # ##################################################################### matrix.grobner:=proc(F,V,termorder) local grob, mingrob, redgrob; grob:=grobnerbasis(F,V,termorder); mingrob:=_minimalgb(grob,V,termorder); redgrob:=_reducegrobner(mingrob,V,termorder); redgrob; end; ##################################################################### ##################################################################### ##################################################################### # # THE PROGRAMS THAT FOLLOW ARE THE SUBROUTINES NEEDED TO RUN # # THE matrixgrobner PROGRAM # ##################################################################### ##################################################################### ##################################################################### # # 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); mono; fi; end; ####### ###### ## Procedure Allen Bell wrote for the S-polynomial ###### ###### with(grobner,leadmon); sp:=proc(f,g,V,U) local f1,g1,a,b,c,d; f1:=expand(f); g1:=expand(g); a:=_leadingterm(f1,V,U); b:=_leadingterm(g1,V,U); c:=lcm(a[2],b[2]); d:=c*f1/(a[1]*a[2])-c*g1/(b[1]*b[2]); simplify(d); end; ##################################################################### ##################################################################### # # FINDING REMAINDERS # # The INPUT is # remainder(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); remainder:=proc(g,Set,V,termorder) 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,termorder); od; v:=0; h:=g; while h<>0 do lmg := _leadingterm(h,V,termorder); for i to nops(Set) do f:=Set[i]; 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,termorder); fi; od; v:=simplify(v+lmg[1]*lmg[2]); h:=simplify(h-lmg[1]*lmg[2]); od; a:=convert(a,list); c:=[v,a]; c; end; ##################################################################### ##################################################################### # # FINDING A MINIMAL GROBNER BASIS # # The INPUT is # _minimalgb([grobner basis],[varlist],[vectorlist]); # # The input [grobner basis] 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]. # ##################################################################### with(linalg,submatrix); with(linalg,swaprow); _minimalgb:=proc(G,V,termorder) local y,i,j,n,H,Ga,Gb; global _grobM, _Ralph; 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,termorder)[2], _leadingterm(H[j],V,termorder)[2]) then Ga:=[H[1..i-1]]; Gb:=[H[i+1..n]]; for y from i to n-1 do _grobM:=swaprow(_grobM,y,y+1); od; _grobM:=submatrix(_grobM,1..n-1, 1.._Ralph); n:=n-1; if j with respect # to the term order defined by [vectorlist]. # #################################################################### with(grobner,leadmon); grobnerbasis:=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,tempmatrix; global _grobM,_Ralph; F:=[]; _grobM:=table(sparse); for i to nops(H) do _grobM[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 _Ralph:=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:=remainder(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 _grobM[n,l]:=-op(l,Q[2]); od; _grobM[n,subscript[1]]:=simplify(_grobM[n, subscript[1]]+p/(lmf[1]*lmf[2])); _grobM[n,subscript[2]]:=simplify(_grobM[n, subscript[2]]-p/(lmg[1]*lmg[2])); fi; fi; fi; od; if nops(F) > nops(H) then for x from _Ralph+2 to n do for y to _Ralph do for z from _Ralph+1 to x-1 do _grobM[x,y]:=simplify(_grobM[x,y]+ _grobM[x,z]*_grobM[z,y]); od; od; od; fi; tempmatrix := array(1..nops(F),1.._Ralph); for index4 to nops(F) do for index5 to _Ralph do tempmatrix[index4,index5] := _grobM[index4,index5]; od; od; _grobM := tempmatrix; F; end; #################################################################### #################################################################### # # FINDING A REDUCED GROBNER BASIS # # The INPUT is # _reducegrobner([minimal GB],[varlist],[vectorlist]); # # The OUTPUT is the reduced Grobner basis with respect to the term order # given by [vectorlist], and the matrix of transformation. # #################################################################### _reducegrobner:=proc(G,V,termorder) local x,Ca,t,i,j,k,l,n,a,o,Da,Db,r,J,_superM,index1,index2,tempmatrix; J:=[]; J:=G; n:=nops(J); 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)]; r:= remainder(J[i],[op(Da),op(Db)],V,termorder); 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; for a to n do x:=_leadingterm(J[a],V,termorder)[1]; o:=J[a]/x; for t to n do Ca[a,t]:=Ca[a,t]/x; od; J:=[op(1..a-1,J),o,op(a+1..n,J)]; od; tempmatrix := array(1..n,1..n); for index1 to n do for index2 to n do tempmatrix[index1,index2] := Ca[index1,index2]; od; od; _superM:=evalm(&*(tempmatrix,_grobM)); print(_superM); J; 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;