##################################################################### ##################################################################### # # IVA GROBNER PACKAGE FOR MAPLE V, RELEASE 5 # # 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 and V.5 # by C. D. Wensley, University of Wales, Bangor # (02/27/98 and 03/03/99) # # Further modifications by William Gryc, # Amherst College (08/13/99) # # 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 # gbr5.mpl. Then, if you start Maple from the directory # containing this file, you can load the file into Maple using # the command # # > read(`gbr5.mpl`): # # If the file resides in another directory, then the appropriate # path needs to be included in the read statement. # ##################################################################### ##################################################################### # # THE MAJOR COMMANDS # # This package has several major commands: # # 1) ring, which sets the term order and variables used by the # package. # # 2) div_alg (for "division algorithm"), which computes the # remainders AND quotients for the division algorithm. # # 3) slowbasis_gb, altbasis_gb and quickbasis_gb, which are # three versions of the basic Grobner basis algorithm. # # 4) min_gb and red_gb, which take a Grobner basis produced by # the three *_gb commands from 3) and make them minimal and # reduced respectively. # # 5) quot_mx (for "matrix of quotients"), which computes a matrix # telling you how to transform the Grobner basis into the # original basis. # # 6) mxgb (for "matrix grobner"), which computes a reduced # 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) quickbasis_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. # # Commands 2) - 6) use a monomial order which is set by the # user through the command ring(). The names of the predefined Maple # term orders (plex and tdeg) should not be used. The orders lex, # grevlex, grlex, [k,n] (elimination order), and [v1,...,vn] (matrix # order) are valid monomial orders (see ring() described below). # # These commands are described in more detail below. # ##################################################################### ##################################################################### # THE ring COMMAND # # The INPUT is: # # > ring(torder, varlist); # # where # torder is the order on the polynomials. The valid values # for torder are lex, grlex, grevlex, [k,n], and # [v1,...,vn]. [k,n] defines the kth elimination order # on n variables. Note that this means the number of # variables must equal n. [v1,...,vn] is a list of # lists, where vi is a n-length list used as a row # vector, and [v1,...,vn] defines a matrix order on n # variables. Note that the vi's must be linearly # independant and there must be n variables in varlist. # # varlist is the list of variables in the ring. # # The OUTPUT is: # # The term-order defined by the Grobner package. The output # is not terribly important in this command. The command's # main purpose is to set the monomial order for all other # commands which need a termorder. # ##################################################################### ##################################################################### # # THE grlex AND elimination COMMANDS # # Since they can be represented by matrix orders and are not # built into Maple, ring() sets grlex and elimination orders # using matrices. ring() uses the commands grlex() and # elimination() to get these matrices for grlex and elimination # orders. # # These commands can be used as part of the input of the # ring() command directly, although this is hardly necessary. # # 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: # # > elimination(k,n) # # where # n is the number of variables, # # k is the number of variables to be eliminated. # # The OUTPUT is: # # the vector list that specifies the kth elimination order on # n variables. # ##################################################################### # # EXAMPLES OF THE grlex AND elimination COMMANDS # # 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: # # > elimination(2,4); # # then the output is the vector list that specifies elimination # order that eliminates the first 2 of 4 variables: # # [[1,1,0,0],[0,0,1,1],[0,0,0,-1],[0,-1,0,0]] # ##################################################################### ##################################################################### # # THE div_alg COMMAND # # The INPUT is: # # > div_alg(f,[f1,...,fs]); # # where # f is the polynomial to be divided # # [f1,...,fs] is the list of input polynomials. # # The OUTPUT is: # # the remainder r and the quotients a1,...,as of the # division of f by f1,...,fs, with respect to the term # order set by ring(). In other words, # # f = a1 f1 + ... + as fs + r. # ##################################################################### # # EXAMPLE OF THE div_alg COMMAND # # If the input is: # # > ring(lex, [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); # # 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 slowbasis_gb, altbasis_gb, AND quickbasis_gb COMMANDS # # All these commands have the same format and do the same thing: # They all calculate a Grobner basis that is generally not # minimal nor reduced. They also print out the steps of the # Grobner basis calculation. If this is not desired, there is # an optional second argument, nosteps, that will supress this # printing. The difference between the three is how they use # Buchberger's algorithm. slowbasis_gb uses a naive version of # the algorithm, altbasis_gb uses a slightly more insightful version, # and quickbasis_gb uses a more advanced version. Since the input # and output are the same for all three commands, we will just use # quickbasis_gb as an example. # # The INPUT is: # # > quickbasis_gb([f1,...,fs]); # OR # > quickbasis_gb([f1,...,fs], nosteps); # # where # [f1,...,fs] is the list of input polynomials. # # nosteps is just the string "nosteps" (minus the quotes). # # The OUTPUT is: # # [g1,...,gt], a list of polynomial that form a Grobner # basis under the termorder set by ring() for . # If the second argument "nosteps" is not entered, the # steps of the calculation with be printed. # ##################################################################### ##################################################################### # # THE min_gb AND red_gb COMMANDS # # Three commands slowbasis_gb, altbasis_gb and quickbasis_gb # compute Grobner bases which in general are neither minimal nor # reduced. Hence there are companion commands min_gb and red_gb # to minimize and reduce them. These commands are described # below. # ##################################################################### ##################################################################### # # THE min_gb COMMAND # # The INPUT is: # # > min_mxgb([f1,...,fs]); # # where # [f1,...,fs] is the list of input polynomials which are # a Grobner basis in the termorder set by ring. # # The OUTPUT is: # # [[g1,...,gt]] # # where # [g1,...,gt] is a minimal Grobner basis of , # which in general is not reduced. # ##################################################################### ##################################################################### # # THE red_gb COMMAND # # The INPUT is: # # > red_gb([f1,...,fs]); # # where # [f1,...,fs] is the list of input polynomials which are # a minimal Grobner basis with respect to the termorder # set by ring(). # # The OUTPUT is: # # [g1,...,gt] # # where # [g1,...,gt] is a reduced Grobner basis of . # ##################################################################### ##################################################################### # # THE quot_mx COMMAND # # The INPUT is: # # > quot_mx([f1,...,fs],[g1,...,gt]); # # where # [f1,...,fs] is a list of polynomials # # [g1,...,gt] is a Grobner basis for the ideal # with respect to the term order set by ring(). # # The OUTPUT is: # # A matrix Q such that Q * [g1,...,gt]^T =[f1,...,fs]^T. # In other Q, the ith row of Q gives how fi is expressed # in terms of [g1,...,gt]. # ##################################################################### # # EXAMPLE OF THE quot_mx COMMAND # # If the input is: # # > ring(grevlex, [x,y]); # > quot_mx([x^2*y - 1,x*y^2 - x],[-y + x^2, y^2-1]); # # then the output is a 2 x 2 matrix which tells how to express # the original basis in terms of the Grobner basis: # # [ y 1 ] # [ ] # [ 0 x ] # # Thus # x^2*y - 1 = y*(-y+x^2) + 1*(y^2-1) # and # x*y^2 - x = 0*(-y+x^2) + x*(y^2-1) # ##################################################################### ##################################################################### # # THE mxgb COMMAND # # The INPUT is: # # > mxgb([f1,...,fs]); # # where # [f1,...,fs] is the list of input polynomials, and steps # of the Grobner basis calculation will be printed. # # OR # # > mxgb([f1,...,fs], nosteps); # # where # [f1,...,fs] is the list of input polynomials, and steps of # the Grobner basis calculation will no be printed. # # # The OUTPUT is: # # a list # # [[g1,...,gt], M] # # where [g1,...,gt] is the reduced Grobner basis of # with respect to the termorder set by # ring() and M is a t x s matrix giving the Grobner # basis as a combination of the fi's; i.e., # # [g1 ] [f1 ] # [g2 ] [f2 ] # [...] = M [...] # [...] [...] # [gt ] [fs ] # ##################################################################### # # EXAMPLE OF THE mxgb COMMAND # # If the input is: # # > ring(lex, [x,y,z]); # > Polys := [x^2+y^2+z^2-4, x^2+2*y^2-5, x*z-1]; # > gb := mxgb(Polys); # # 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], # [ # [ # # [ -1 1 0 ]] # [ ]] # [ 2z - z - x ]] # [ ]] # [ z^2 - 1/2 z^2 - 1/2 - 1/2 x *z ]] # # Also, the intermediate steps of the computation (the "unminimized" # basis and matrix step, and the "unreduced" basis and matrix step) # are also printed during computation, unless the optional second # argument "nosteps" is included. The command line for this option # would look like: # # > mxgb(Polys, nosteps); # ##################################################################### ##################################################################### # # THE quickbasis_mxgb, min_mxgb AND red_mxgb COMMANDS # # As explained above, the command mxgb operates by calling the # three commands: # 1) quickbasis_mxgb, which computes a Grobner basis (and associated # matrix) which is usually neither minimal nor reduced. # 2) min_mxgb, which takes a Grobner basis (and matrix) and # produces a minimal Grobner basis (and matrix). # 3) red_mxgb, which takes a minimal Grobner basis (and matrix) # and produces a reduced Grobner basis (and matrix). # Below, we give a brief description of these commands (though # here we don't include examples). # ##################################################################### ##################################################################### # # THE quickbasis_mxgb COMMAND # # The INPUT is: # # > quickbasis_mxgb([f1,...,fs]); # # where # [f1,...,fs] is the list of input polynomials. # # The OUTPUT is: # # a 3-element list: # # ["isgbasis",[g1,...,gt],coeff_mx] # # where # [g1,...,gt] is a Grobner basis of with # respect to the termorder set by ring(); the Grobner # basis is in general is neither minimal nor reduced, # # coeff_mx is the corresponding transformation matrix. # ##################################################################### ##################################################################### # # THE min_mxgb COMMAND # # The INPUT is: # # > min_mxgb("isgbasis", [f1,...,fs],coeff_mx); # # where # "isgbasis" is just a string # # [f1,...,fs] is the list of input polynomials which are # a grobner basis. # # coeff_mx is the corresponding transformation matrix for the # Grobner basis # # The OUTPUT is: # # a 3-element list: # # ["isminimal",[g1,...,gt],coeff_mx] # # where # [g1,...,gt] is a minimal Grobner basis of , # which in general is not reduced, # # coeff_mx is the corresponding transformation matrix. # ##################################################################### ##################################################################### # # THE red_mxgb COMMAND # # The INPUT is: # # > red_mxgb("isminimal",[f1,...,fs],coeff_mx); # # where # [f1,...,fs] is the list of input polynomials which are # a minimal Grobner basis with respect to the termorder # set by ring(). # # The OUTPUT is: # # a 2-element list # # [[g1,...,gt],coeff_mx] # # where # [g1,...,gt] is a reduced Grobner basis of , # # coeff_mx is the corresponding transformation matrix. # ##################################################################### ##################################################################### # # END OF DOCUMENTATION # ##################################################################### ##################################################################### # # HERE ARE STANDARD MAPLE PACKAGES USED BY gbr5. # ##################################################################### with(linalg); with(Groebner); with(Ore_algebra); ##################################################################### # # HERE IS THE CODE FOR mxgb. # ##################################################################### mxgb := proc(F::list) local basislist, minlist, redlist, temp1, temp2; basislist := quickbasis_mxgb(F, nosteps); if nargs = 1 then temp1 := basislist[2]; print(`Unminimized basis: `, temp1); temp2 := op(basislist[3]); print(`Unminimized matrix: `, temp2); fi; minlist := min_mxgb( basislist ); if nargs = 1 then temp1 := minlist[2]; print(`Minimized basis: `, minlist[2]); temp2 := op(minlist[3]); print(`Minimized matrix: `, temp2); fi; redlist := red_mxgb( minlist ); redlist[3] := op(redlist[3]); redlist[2..3]; end; ##################################################################### # # HERE IS THE CODE FOR _leadingterm. THIS FINDS LEADING TERMS # AND IS USED BY ALL OF THE GROBNER BASIS ROUTINES. # ##################################################################### _leadingterm := proc(f) local f1, mono; global morder; if type(op(morder), name) then ERROR(`Monomial Order not set. Use ring() to set it.`) fi; f1 := expand(f); mono := leadmon(f1, morder); mono; end; ##################################################################### # # HERE IS THE CODE FOR div_alg. # ##################################################################### div_alg := proc( g, Set ) 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 ); od; v := 0; h := g; while (h<>0) do lmg := _leadingterm( h ); 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 ); 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; ##################################################################### # # HERE IS THE CODE FOR quot_mx. # ##################################################################### quot_mx := proc( pols, bas ) 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 ); if ( d[1]<>0 ) then ERROR( `non-zero remainder` ); fi; for index2 to nb do Q[index1,index2] := d[2][index2]; od; od; op(Q); end; ##################################################################### # # HERE IS THE CODE FOR min_mxgb. THIS IS USED BY mxgb. # ##################################################################### min_mxgb := proc( gbrec ) local y, i, j, n, H, G, coeff_mx, numcol, minrec,lt; if not ( ( nargs = 1 ) and type( gbrec, list ) and ( nops( gbrec ) = 3 ) and type( gbrec[1], string ) and ( gbrec[1] = "isgbasis" ) ) then ERROR(`input not a matrixgrobner record` ); fi; G := gbrec[2]; coeff_mx := gbrec[3]; numcol := coldim( coeff_mx ); n:=nops(G); H:=G; i:=1; while i<=n do lt := _leadingterm( H[i] )[2]; j:=1; while j<=n do if i<>j then if divide(lt,_leadingterm( H[j] )[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; j:=j+1; od; i := i+1; od; minrec := [ "isminimal", H, coeff_mx ]; end; ##################################################################### # # HERE IS THE CODE FOR min_gb. # ##################################################################### min_gb := proc( gb ) local i, j, n, H, G, minrec,lt; if not isgbasis(gb) then ERROR(`input not a Grobner basis`); fi; G := gb; n:=nops(G); H:=G; i:=1; while i<=n do lt := _leadingterm( H[i] )[2]; j:=1; while j<=n do if i<>j then if divide(lt,_leadingterm( H[j] )[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; n:=n-1; j:=n; i := i-1; fi; fi; j:=j+1; od; i := i+1; od; minrec := H ; end; ##################################################################### # # HERE IS THE CODE FOR is_min_gbasis. THIS TESTS FOR MINIMALITY # AND IS USED BY red_gb. # ##################################################################### is_min_gbasis := proc(G::list) local i, j, n, lt; if not isgbasis(G) then RETURN(false); fi; n := nops(G); for i to n do for j to n do lt := _leadingterm(G[i])[2]; if i <> j then if divide(lt, _leadingterm(G[j])[2]) then RETURN(false); fi; fi; od; od; true; end; ##################################################################### # # HERE IS THE CODE FOR isgbasis. THIS TESTS FOR BEING A GROBNER # BASIS AND IS USED BY min_gb AND is_min_gbasis. # ##################################################################### isgbasis := proc(G::list) local i,j,n,r; global morder; n := nops(G); for i to n do for j to n do if i <> j then r := normalf(spoly(G[i], G[j], morder), G, morder); if r <> 0 then RETURN(false); fi; fi; od; od; true; end; ##################################################################### # # HERE IS THE CODE FOR quickbasis_mxgb. THIS IS USED BY mxgb. # ##################################################################### quickbasis_mxgb := proc(H::list) 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, divs; divs := 0; F:=[]; coeff_tab := table(sparse); for i to nops(H) do coeff_tab[i,i]:=1; od; F:=H; if nargs = 1 then print(`Current Basis: `.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]] ); lmg := _leadingterm( F[subscript[2]] ); 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])[2]) then ernie:=1; break; fi; od; if ernie = 0 then p := lcm(lmf,lmg); T := p*F[subscript[1]]/(lmf[1]*lmf[2]) - p*F[subscript[2]]/(lmg[1]*lmg[2]); Q := div_alg(simplify(T),F); 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 = 1 ) then print(`Added: `.r) fi; divs := divs + 1; 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, coeff_mx ]; end; ##################################################################### # # HERE IS THE CODE FOR quickbasis_gb. # ##################################################################### quickbasis_gb := proc(H::list) local x, y, setofpairs, z, p, l, F, r, i, n, lmf, lmg, ernie, numcol, u, T, Q, index1, index2, index3, subscript, gbrec, divs,localdivs; global morder; F:=[]; F:=H; divs := 0; localdivs := 0; if nargs = 1 then print(`Current Basis: `.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]] ); lmg := _leadingterm( F[subscript[2]] ); 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])[2]) then ernie:=1; break; fi; od; if ernie = 0 then p := lcm(lmf,lmg); T := p*F[subscript[1]]/(lmf[1]*lmf[2]) - p*F[subscript[2]]/(lmg[1]*lmg[2]); Q := normalf(simplify(T),F,morder); divs := divs + 1; localdivs := localdivs + 1; r := Q; if ( r <> 0 ) then for index3 to nops(F) do setofpairs:=[op(setofpairs),[index3,nops(F)+1]]; od; F := [op(F),r]; if ( nargs = 1 ) then print(`Added: `.r); print(`Local divisions: `.localdivs); fi; localdivs := 0; fi; fi; fi; od; if (nargs = 1) then print(`Total divisions performed: `. divs); fi; gbrec := F; end; ##################################################################### # # HERE IS THE CODE FOR slowbasis_gb. # ##################################################################### slowbasis_gb := proc(F::list) local G, Gprime, S, noneadded, p, q, printable, steps, divs,i,j,length, oldlength,localdivs; global morder; divs := 0; if nargs = 1 then print(`Current basis: `. F) fi; G := F; steps = 0; noneadded := false; oldlength := 1; while not noneadded do localdivs := 0; steps := steps + 1; printable := []; noneadded := true; Gprime := G; length := nops(Gprime); for i from 1 to length do p := Gprime[i]; for j from max(i+1, oldlength + 1) to length do q := Gprime[j]; S := normalf(spoly(p,q,morder), Gprime, morder); divs := divs + 1; localdivs := localdivs + 1; if S <> 0 then G := [op(G), S]; printable := [op(printable), S]; noneadded := false fi; od; od; oldlength := length; if nargs = 1 then print(`Added: `. printable); print(`Local divisions: `. localdivs); fi; od; if (nargs = 1) then print(`Total divisions performed: `. divs); fi; G; end; ##################################################################### # # HERE IS THE CODE FOR altbasis_gb. # ##################################################################### altbasis_gb := proc(F::list) local G, Gprime, S, noneadded, p, q, printable, steps, divs,i,j,length, oldlength,localdivs; global morder; divs := 0; if nargs = 1 then print(`Current basis: `.F) fi; G := F; steps = 0; noneadded := false; oldlength := 1; while not noneadded do localdivs := 0; steps := steps + 1; printable := []; noneadded := true; Gprime := G; length := nops(Gprime); for i from 1 to length do p := Gprime[i]; for j from max(i+1, oldlength + 1) to length do q := Gprime[j]; if is(p = q) then next fi; S := normalf(spoly(p,q,morder), G, morder); divs := divs + 1; localdivs := localdivs + 1; if S <> 0 then G := [op(G), S]; printable := [op(printable), S]; noneadded := false fi; od; od; oldlength := length; if nargs = 1 then print(`Added: `. printable); print(`Local divisions: `. localdivs) fi; od; if (nargs = 1) then print(`Total divisions performed: `. divs); fi; G; end; ##################################################################### # # HERE IS THE CODE FOR red_mxgb. THIS IS USED BY mxgb. # ##################################################################### red_mxgb := proc( minrec ) local x, Ca, t, i, j, k, l, m, n, a, o, Da, Db, r, J, J0, J1, 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]; grobmx := minrec[3]; 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) ] ); 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] )[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; ##################################################################### # # HERE IS THE CODE FOR red_gb. # ##################################################################### red_gb := proc(minrec) local x, t, i, j, k, l, m, n, a, o, Da, Db, r, J, J0, J1, index1, index; global morder; if not is_min_gbasis(minrec) then ERROR(`input not a minimal Grobner basis`); fi; J := minrec; n := nops(J); i:=1; if n<>1 then while (i <= n) do 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 := normalf( J[i], [ op(Da), op(Db) ], morder ); J := [ op(Da), r, op(Db) ]; i := i + 1; od; fi; J0 := [ ]; J1 := [ ]; for a to n do if (J[a]<>0) then J0 := [ op(J0), a ]; J1 := [ op(J1), J[a] ]; fi; od; n := nops( J0 ); J := J1; for a to n do x := _leadingterm( J[a] )[1]; o := J[a]/x; 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; J; end; ##################################################################### # # HERE IS THE CODE FOR ring. # ##################################################################### ring := proc(torder,varlist) local badform,i; global morder; badform := false; if is(torder, list) then if is(torder[1], list) then for i from 1 to (nops(torder)-1) do badform := is(nops(torder[i]) <> nops(torder[i+1])); if badform then break fi; od; badform := badform or is(nops(torder) <> nops(varlist)); if badform or is(nops(torder) <> nops(torder[1])) then ERROR(`Not a valid matrix`); else morder := termorder(poly_algebra(op(varlist)),'matrix' (torder, varlist)); fi; elif is(nops(torder) <> 2) or is(torder[1] >= torder[2]) or is(torder [1] <= 0) or is(torder[2]<>nops(varlist)) then ERROR(`Not a valid elimination form`); else morder := termorder(poly_algebra(op(varlist)),'matrix' (elimination(torder[1], torder[2]), varlist)); fi; elif is(torder = grlex) then morder := termorder(poly_algebra(op(varlist)),'matrix' (grlex(nops(varlist)), varlist)); elif is(torder = grevlex) then morder := termorder(poly_algebra(op(varlist)), tdeg(op(varlist))); elif is(torder = lex) then morder := termorder(poly_algebra(op(varlist)), plex(op(varlist))); elif (torder = current) then else ERROR(`First argument is not a valid monomial order`); fi; morder; end; ##################################################################### # # HERE IS THE CODE FOR grlex. THIS IS USED BY ring TO GENERATE # A MATRIX WHICH GIVES THE GRLEX ORDER. # ##################################################################### 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; ##################################################################### # # HERE IS THE CODE FOR elimination. THIS IS USED BY ring TO # GENERATE A MATRIX WHICH GIVES THE ELIMINATION ORDER. # ##################################################################### elimination := proc(k,n) local u,i,j; if k >= n or k = 0 then ERROR(`Not a valid elimination form.`); elif n=1 then [[1]]; else u := [0 $ i=1..n-j-1, -1, 0 $i=n-j+1..n]; [[1 $ i=1..k, 0 $i=k+1..n],[0 $ i=1..k, 1 $i=k+1..n],u $ j=0..n-k-2, u $ j =n-k..n-2] fi; end;