(* This is the Groebner Package for version 3.0 of Mathematica *) (* Created by Susan Goldstine, Will Gryc and David Cox of Amherst College. Last updated August 3, 1999. *) (* INSTALLATION To use this package, you should copy this file under a name, such as groebner30.m, into the directory from which you start Mathematica. Every time you start a Mathematica session you will need to load this package. This can be done by giving the command << groebner30.m at the Mathematica prompt. The first part of this file is the documentation that explains how to use the package. *) (* GENERAL DESCRIPTION This package does a variety of commands related to Groebner bases and can be used with the text "Ideal, Varieties, and Algorithms." The user should be warned that the package was written for teaching purposes and hence very slow compared to the GroebnerBasis and the PolynomialReduce commands which are part of Mathematica 3.0. This package is designed to work over the rational numbers Q. Hence the commands in the package cannot be used for polynomials whose coefficients lie in finite fields or function fields. *) (* CONTENTS OF PACKAGE The major commands in the Groebner package are as follows: MonOrder: sets the monomial order. PRemainder: computes the remainder in the division algorithm. PQuotient: computes the quotients in the division algorithm. Buchberger: a naive version of the Buchberger algorithm. AltBuchberger: a bit more discriminating version of the Buchberger algorithm. QuickBuchberger: a more efficient Buchberger algorithm. The package also includes commands to do the following: SPoly: computes S-polynomials. ReduceGroebner: reduces Groebner basis to a reduced one. GroebnerQ: tests to see if you have a Groebner basis. IdealQ: solves the ideal membership problem. RadicalQ: solves the radical membership problem. FiniteQ: determines if there are finitely many solutions. VSDimension: determines the spanning cosets and dimension of a quotient ring. Radical: determines generators for the radical of a zero-dimensional ideal. RadicalIdealQ: determines if a zero-dimensional ideal is radical. ConvertGB: converts a Groebner basis of a zero-dimensional ideal from one monomial and variable order to a Lex or Grlex Groebner basis with another variable order. IdealofPoints: determines a Groebner basis for the ideal of all polynomials which vanish on a given finite set of points. *) (* DESCRIPTIONS OF THE COMMANDS The descriptions of MonOrder and PRemainder should be read carefully. For each of the commands described below, on-line documentation is available through the usual Mathematica help facility. Any parameters listed in parentheses are optional. MonOrder: This command is used to set the monomial order used in all other commands in the package. The allowable monomial orders are denoted as follows: Lex (lexicographic order) Grlex (graded lexicographic order) Grevlex (graded reverse lexicographic order) {k,n} (elimination order, eliminating the first k of n variables) {{a11,...,a1n},...,{an1,...,ann}} (matrix order for n variables) Note the use of capital letters in the names of the first three of these orders. Thus, to change the monomial order to Grlex, you would issue the command MonOrder[Grlex] Furthemore, the command MonOrder[] will return the current monomial order. This is useful if you've forgotten what the order is. The default monomial order is Lex. PRemainder: The format for this command is PRemainder[f,{f1,...,fs},varlist] where f is the polynomial to be divided and {f1,...,fs} is the list of polynomials to divide by. The output will be the remainder of f on division by f1,...fs. The "varlist" is a list of variables. You should be aware that the order of the variables in varlist is important since it (together with what you specified in MonOrder) determines the monomial order. PQuotient: The format is PQuotient[f,{f1,...,fs},varlist] where f is the polynomial to be divided and {f1,...,fs} is the list of polynomials to divide by. The output will be a list of quotients {a1,...,as} of f on division by f1,...fs. The "varlist" is a list of variables that determines the order of the variables in the current monomial order. Buchberger (or B for short): The format is Buchberger[{f1,...,fs},varlist, options...] or B[{f1,...,fs},varlist, options...] where {f1,...,fs} is a list of polynomials and, "varlist" determines the order of the variables used in the current monomial order. The output of this command is a Groebner basis for the ideal , and it also prints out the number of times that PRemainder is performed. The algorithm employed is a naive version of the Buchberger algorithm, with the one efficiency that no S-polynomial that has already been checked is checked again. The output need not be a reduced Groebner basis. The steps of the computations will be printed. If this is not desired, the option PrintSteps should be set to False. AltBuchberger (or AB for short): The format is AltBuchberger[{fi,...,fs}, varlist, options...] or AB[{f1,...,fs}, varlist, options...] This command returns a Groebner basis for the ideal , with respect to the monomial order set by MonOrder and the order of variables given by varlist. This commands differs only slightly to Buchberger, in that it uses the latest constructed basis for all S-polynomial divisions. The output need not be a reduced Groebner basis. The steps of the computations will be printed. If this is not desired, the option PrintSteps should be set to False. QuickBuchberger (or QB for short): The format is QuickBuchberger[{f1,...,fs},varlist, options...] or QB[{f1,...,fs},varlist, options...] This command uses a more efficient version of the Buckberger algorithm to compute a Groebner basis for the ideal . The output need not be a reduced Groebner basis. As usual, varlist determines the order of the variables used in the current monomial order. The steps of the computations will be printed. If this is not desired, the option PrintSteps should be set to False. MathematicaBuchberger (or MB for short): The format is MathematicaBuchberger[{f1,...,fs},varlist] or MB[{f1,...,fs},varlist] This command directly uses Mathematica's GroebnerBasis command to find a reduced Groebner basis for the ideal with respect to the current monomial order. As usual, varlist determines the order of the variables. SPoly: The format is SPoly[poly1,poly2,varlist] This command returns the S-polynomial of poly1 and poly2, and "varlist" determines the order of the variables used in the current monomial order. ReduceGroebner: The format is ReduceGroebner[{f1,...,fs},varlist,(torder)] This command takes a Groebner basis {f1,...,fs} and produces a reduced Groebner basis for the same ideal. If specified, torder is the monomial ordering. If left unspecified, the monomial ordering is determined by what is set by MonOrder. Be sure not to change the monomial order or "varlist" from when you produced the Groebner basis. GroebnerQ: The format is GroebnerQ[{f1,...,fs},varlist,(torder)] This command returns either True or False, depending on whether {f1,...,fs} forms a Groebner basis under the current monomial ordering or torder, if it is specified, with "varlist" determining the order of the variables used in the monomial order. IdealQ: The format is IdealQ[f,{f1,...,fs},varlist] This is an implementation of the Ideal Membership Algorithm. The command returns either True or False, depending on whether the polynomial f is in the ideal , with "varlist" determining the order of the variables used in the current monomial order. RadicalQ: The format is RadicalQ[f,{f1,...,fs},varlist] This is an implementation of the Radical Membership Algorithm. The command returns either True or False, depending on whether the polynomial f is in the radical of the ideal , with "varlist" determining the order of the variables used in the current monomial order. FiniteQ: The format is FiniteQ[{f1,...,fs},varlist] This is an implementation of the Finiteness Algorithm. The command returns True or False, depending on whether the set of equations specified by f1=...=fs=0 has a finite number of solutions over the complex numbers. VSDimension: The format is VSDimension[{f1,...,fs},varlist,(torder)] This command returns the basis remainders and the dimension (as a vector space) of the polynomial ring modulo the ideal . The output is {Infinity, Infinity} if the quotient is infinite dimensional. If specified, torder sets the monomial order, and if unspecified, the monomial order is the current order set by MonOrder. Radical: The format is Radical[{f1,...fs}, varlist] This command returns generators of the radical of the zero-dimensional ideal defined by . The returned polynomials do NOT necessarily form a Groebner basis. RadicalIdealQ: The format is RadicalIdealQ[{f1,...,fs}, varlist] This command returns True if the zero-dimensional ideal generated by {f1,...,fs} is radical, and False otherwise. ConvertGB: The format is ConvertGB[{f1,...,fs}, varlist, oldvarlist, options...] In this command, {f1,...,fs} is a Groebner basis of a zero-dimensional ideal for oldvarlist and the current monomial order set by MonOrder. This order specification can be changed by setting the option OldOrderType to the appropriate order. The other options are OrderType and PrintSteps. When OrderType is set to Grlex, ConvertGB returns a Groebner basis with respect to varlist and grlex that generates the same ideal as {f1,...,fs}. If set to Lex or left unset, ConvertGB returns a Groebner basis with respect to varlist and lex. In both cases, the remainder monomials with respect to the monomial ordering are also returned. If PrintSteps is set to True, every time the command adds a polynomial to the constructed Groebner basis, it will print that polynomial and print the remainder monomials added since the last polynomial was added, if any were added. IdealofPoints: The format is IdealofPoints[{P1,...,Ps}, varlist, (torder)] This command returns the Groebner basis of the ideal of all polynomials which vanish on points {P1,...,PS}. Each point is given as an n-dimensional vector. If left unspecified, torder is the monomial order set by MonOrder. Note that there is also on-line documentation for each of these commands. For example, to find out about the "Buchberger" command, type ?Buchberger at the Mathematica prompt. END OF DOCUMENTATION *) BeginPackage["Buchberger`"] MonOrder::usage = "MonOrder[] returns the current monomial order; MonOrder[order] changes the current monomial ordering to order. The valid arguments for MonOrder are Lex, Grlex (graded lex), Grevlex (graded reverse lex), {k, n} (the k-th elimination order of n variables), and {v1, v2,..., vn} (a matrix order for n variables, where v1,...,vn are rows of an n x n matrix)." Lex::usage = "Lex is an argument of MonOrder which changes the monomial ordering to lex order." LexLT::usage = "LexLT[poly, varlist] returns the leading term in lex order of the polynomial poly. The variable varlist determines the ordering of the variables." Grlex::usage = "Grlex is an argument of MonOrder which changes the monomial ordering to graded lex order." GrlexLT::usage = "GrlexLT[poly, varlist] returns the leading term in grlex order of the polynomial poly. The variable varlist determines the ordering of the variables." Grevlex::usage = "Grevlex is an argument of MonOrder which changes the monomial ordering to graded reverse lex order." GrevlexLT::usage = "GrevlexLT[poly, varlist] returns the leading term in grevlex order of the polynomial poly. The variable varlist determines the ordering of the variables." MatrixLT::usage = "MatrixLT[poly, varlist] returns the leading term, using the matrix order previously set by the MonOrder command. The variable varlist determines the ordering of the variables. " ElimMatrix::usage = "ElimMatrix[k, n] returns the n x n matrix which gives the kth-elimination order on n variables. Note that k must be strictly less than n." PQuotient::usage = "PQuotient[f,{f1,...,fs},varlist] returns the list of quotient polynomials of f divided by {f1,...,fs}, i.e. {a1,...,as} where a1*f1 +...+ as*fs + r = f. The variable varlist determines the order of variables used in the current monomial order." PRemainder::usage = "PRemainder[f,{f1,...,fs}, varlist] returns the remainder of f divided by {f1,...,fs}. The variable varlist determines the order of the variables used in the current monomial order." SPoly::usage = "SPoly[poly1,poly2, varlist] returns the S-polynomial of poly1 and poly2. The variable varlist determines the order of the variables used in the current monomial order." Buchberger::usage = "Buchberger[{f1,...,fs},varlist,options...] returns a Groebner basis for the ideal and prints out the number of times that PRemainder is performed. The algorithm is made as efficient as it can be without major structural changes by assuring that no S-polynomial that has already been checked is checked again. The variable varlist determines the order of the variables used in the current monomial order. The steps of the computations will be printed. If this is not desired, the option PrintSteps should be set to False." B::usage = "B[{f1,...,fs}, varlist, options...] is the same command as Buchberger, but with an abbreviated name." AltBuchberger::usage ="AltBuchberger[{f1,..,fs}, varlist, options...] returns a Groebner basis for the ideal , with respect to the monomial order set by MonOrder and the order of variables given by varlist. This commands differs only slightly to Buchberger, in that it uses the latest constructed basis for all S-polynomial divisions. The steps of the computations will be printed. If this is not desired, the option PrintSteps should be set to False" AB::usage = "AB[{f1,...,fs}, varlist, options...] is the same command as AltBuchberger, but with an abbreviated name." QuickBuchberger::usage = "QuickBuchberger[{f1,...,fs},varlist,options...] uses a highly modified, more efficient version of the Buchberger algorithm to compute a Groebner basis for the ideal . The variable varlist determined the order of the variables used in the current monomial ordering. The steps of the computations will be printed. If this is not desired, the option PrintSteps should be set to False." QB::usage = "QB[{f1,...,fs}, varlist,options...] is the same command as QuickBuchberger, but with an abbreviated name." MathematicaBuchberger::usage = "MathematicaBuchberger[{f1,...,fs}, varlist] uses the built-in GroebnerBasis command to find the Groebner basis of {f1,...,fs} relative to varlist and the monomial order determined by MonOrder[]." MB::usage = "MB[{f1,...fs}, varlist] is the same command as MathematicaBuchberger, but with an abbreviated name." GroebnerQ::usage = "GroebnerQ[{f1,...,fs},varlist,(torder)] determines whether {f1,...,fs} forms a Groebner basis under the order torder, or the current monomial ordering if torder is unspecified. The variable varlist determines the order of the variables used in torder." ReduceGroebner::usage = "ReduceGroebner[basis,varlist,(torder)] takes a Groebner basis and produces a reduced Groebner basis for the same ideal. The variable varlist determines the order of the variables used in the current monomial ordering. If left unspecified, torder is the order determined by MonOrder[]." IdealQ::usage = "IdealQ[f,{f1,...,fs}, varlist] determines whether the polynomial f is in the ideal . The variable varlist determines the order of the variables used in the current monomial ordering." RadicalQ::usage = "RadicalQ[f,{f1,...,fs}, varlist] determines whether the polynomial f is in the radical of the ideal . The argument varlist determines the order of the variables used in the current monomial ordering." FiniteQ::usage = "FiniteQ[{f1,...,fs}, varlist] determines whether the set of equations specified by f1=...=fs=0 has a finite number of solutions over the complex numbers." VSDimension::usage = "VSDimension[{f1,...,fs},varlist, (torder)] returns the monomial basis and the dimension (as a vector space) of the polynomial ring modulo the ideal . If the quotient ring is of infinite dimension, {Infinity, Infinity} is returned. If left unspecified, order is the order determined by MonOrder. Note: input torder as you would for MonOrder (i.e., input the 3rd elimination order on 5 variables {3, 5}, as you would for MonOrder[])." ConvertGB::usage = "ConvertGB[{f1,...,fs}, varlist, oldvarlist, options...] takes a Groebner basis of a zero-dimensional ideal computed using the variable order of oldvarlist and the monomial order set by MonOrder to make a Lex Groebner basis relative to the varlist. The option OldOrderType can be set to any monomial order, and should be set to the order of which {f1,...,fs} is a Groebner basis, if this order is different than the order set by MonOrder. You specify OldOrderType as you would for MonOrder (i.e., input the 3rd elimination order on 5 variables {3, 5}, as you would for MonOrder[]). The other options for ConvertGB are OrderType, which can be set to Grlex for a Groebner basis conversion to Grlex, and PrintSteps which, if set to True, will print out the current Groebner basis and remainder basis after the algorithm places a new element in the constructed Groebner basis." Radical::usage ="Radicial[{f1,...,fs}, varlist] returns generators, NOT necessarily a Groebner basis, of the radical of the zero-dimensional ideal generated by {f1,...,fs}." RadicalIdealQ::usage="RadicalBasisQ[{f1,...,fs}, varlist] returns True if the zero-dimensional ideal generated by {f1,...,fs} is radical, and false if it is not." IdealofPoints::usage = "IdealofPoints[{P1,...,Ps}, varlist, (torder)] returns a Groebner basis of the ideal of all polynomials vanishing on points P1,...,Ps. Each point given is an n-dimensional vector. If left unspecified, torder is the order currently given by MonOrder[]." OldOrderType::usage = "Option for ConvertGB that tells the command what in what order the basis given as a parameter is a Groebner basis. Default is the order set by MonOrder." OrderType::usage ="Option for ConvertGB that selects the type of Groebner basis returned by ConvertGB: Lex or Grlex. Default is Lex." PrintSteps::usage = "Option for Buchberger, AltBuchberger, QuickBuchberger, and ConvertGB that determines if the function prints out intermediate steps. Default is False for ConvertGB and True for the others." Options[ConvertGB] = {OldOrderType-> MonOrder[], OrderType -> Lex, PrintSteps -> False} Options[Buchberger] ={PrintSteps -> True} Options[AltBuchberger] = {PrintSteps -> True} Options[QuickBuchberger] = {PrintSteps -> True} Begin["`Private`"] OrderList = {Lex,LexLT,Lexicographic} MonOrder[order_:Automatic]:= Module[{validtype=True, badform = False, i}, If [order===Lex, OrderList[[1]] = Lex; OrderList[[2]] = LexLT; OrderList[[3]] = Lexicographic, If [order===Grlex, OrderList[[1]] = Grlex; OrderList[[2]] = GrlexLT; OrderList[[3]] = DegreeLexicographic, If [order===Grevlex, OrderList[[1]] = Grevlex; OrderList[[2]] = GrevlexLT; OrderList[[3]] = DegreeReverseLexicographic, If [Length[order] != 0, If [Length[order[[1]]] == 0, If[((Length[order] != 2) || (order[[1]] >= order[[2]])), Print["Error: Not a valid Elimination form"], OrderList[[1]] = order; OrderList[[2]] = MatrixLT; OrderList[[3]] = ElimMatrix[order[[1]], order[[2]]]], Do [ If [Length[order[[i-1]]] != Length[order[[i]]], badform = True], {i, 2, Length[order]}]; If [(Length[order] != Length[order[[1]]]) || badform, Print["Error: Matrix must be in n x n form"], If [Det[order]===0, Print["Error: Matrix rows are linearly dependent."], OrderList[[1]] = order; OrderList[[2]] = MatrixLT; OrderList[[3]] = order]]], validtype = False]]]]; If [(!validtype) && (order=!=Automatic), Print["Error: ", order," is not a valid monomial order."]]; OrderList[[1]] ] ElimMatrix[k_, n_] := Module[{elim ={}, i, j, temp}, Do [ temp = {}; If [i==1, Do[If [j <= k, AppendTo[temp, 1], AppendTo[temp, 0]], {j, 1, n}], If [i==2, Do[If [j <= k, AppendTo[temp, 0], AppendTo[temp, 1]], {j, 1, n}], If [((i > 2) && (i <= (n-k+1))), Do[If [j == (n + 3 - i), AppendTo[temp, -1],AppendTo[temp, 0]], {j, 1, n}], If [i > (n-k+1), Do[If [j == (n + 2 - i), AppendTo[temp, -1],AppendTo[temp, 0]], {j, 1, n}], Print ["Error!!!"]]]]]; AppendTo[elim, temp], {i, 1, n}]; elim ] MatrixLT[polys_List,varlist_] := Map [MatrixLT[#,varlist]&,polys] MatrixLT[poly_, varlist_] := Module[{IVAtempvar = varlist[[1]]}, Divide[MonomialList[(IVAtempvar*poly), varlist, MonomialOrder -> OrderList[[3]]][[1]], IVAtempvar] ] LexLT[polys_List,varlist_]:= Map [LexLT[#,varlist]&,polys] LexLT[poly_, varlist_] := Module[{IVAtempvar = varlist[[1]]}, Divide[MonomialList[(IVAtempvar*poly), varlist, MonomialOrder -> Lexicographic][[1]], IVAtempvar] ] GrlexLT[polys_List,varlist_]:= Map [GrlexLT[#,varlist]&,polys] GrlexLT[poly_, varlist_] := Module[{IVAtempvar = varlist[[1]]}, Divide[MonomialList[(IVAtempvar*poly), varlist, MonomialOrder -> DegreeLexicographic][[1]], IVAtempvar] ] GrevlexLT[polys_List,varlist_]:= Map [GrevlexLT[#,varlist]&,polys] GrevlexLT[poly_, varlist_] := Module[{IVAtempvar = varlist[[1]]}, Divide[MonomialList[(IVAtempvar*poly), varlist, MonomialOrder -> DegreeReverseLexicographic][[1]], IVAtempvar] ] PRemainder[poly_, polylist_List, varlist_List] := Module[{}, PolynomialReduce[poly, polylist, varlist, MonomialOrder -> OrderList[[3]]][[2]] ] PQuotient[poly_, polylist_List, varlist_List] := Module[{}, PolynomialReduce[poly, polylist, varlist, MonomialOrder -> OrderList[[3]]][[1]] ] MonoLCM [mono1_,mono2_,vars_]:= Module[{loop}, Product[ vars[[loop]]^Max[ Exponent[mono1,vars[[loop]]], Exponent[mono2,vars[[loop]]] ], {loop,1,Length[vars]} ] ] SPoly [poly1_,poly2_,varlist_]:= Module[{ltfunc,vars=varlist,lt1,lt2,lcm}, ltfunc = OrderList[[2]]; lt1 = ltfunc[poly1,vars]; lt2 = ltfunc[poly2,vars]; lcm = MonoLCM[lt1,lt2,vars]; Simplify[(lcm/lt1)*poly1 - (lcm/lt2)*poly2] ] B[polylist_, varlist_,opts___Rule] := Buchberger[polylist, varlist,opts] Buchberger [polylist_,varlist_,opts___Rule]:= Module[{basis=polylist,basisprime,vars=varlist,breakloop, length=Length[polylist],oldlength=1,spoly,spolyrem, loop1,loop2,step=0,printtable,divloop,divisions=0,steps}, steps = PrintSteps/.{opts}/.Options[Buchberger]; breakloop = False; While[!breakloop, step++; divloop = 0; basisprime = basis; Do [spoly=SPoly[basisprime[[loop1]], basisprime[[loop2]],vars]; spolyrem=PRemainder[spoly,basisprime,vars]; divloop++; If[!(spolyrem===0),AppendTo[basis,spolyrem]], {loop1,1,length}, {loop2,Max[loop1+1,oldlength+1],length}]; oldlength = length; length = Length[basis]; printtable = Table[basis[[loop1]], {loop1,oldlength+1,length}]; If[steps, If[printtable==={}, Print["Step ",step,": nothing added."], Print["Step ",step,": ",printtable," added."]]; If[divloop===1, Print[" 1 division performed."], Print[" ",divloop," divisions performed."]]]; divisions = divisions + divloop; breakloop = (oldlength===length) ]; If[steps, Print [divisions," divisions total."]]; basis ] AltBuchberger [polylist_,varlist_,opts___Rule]:= Module[{basis=polylist,basisprime,vars=varlist,breakloop, length=Length[polylist],oldlength=1,spoly,spolyrem, loop1,loop2,step=0,printtable,divloop,divisions=0,steps}, steps = PrintSteps/.{opts}/.Options[AltBuchberger]; breakloop = False; While[!breakloop, step++; divloop = 0; basisprime = basis; Do [spoly=SPoly[basisprime[[loop1]], basisprime[[loop2]],vars]; spolyrem=PRemainder[spoly,basis,vars]; divloop++; If[!(spolyrem===0),AppendTo[basis,spolyrem]], {loop1,1,length}, {loop2,Max[loop1+1,oldlength+1],length}]; oldlength = length; length = Length[basis]; printtable = Table[basis[[loop1]], {loop1,oldlength+1,length}]; If[steps, If[printtable==={}, Print["Step ",step,": nothing added."], Print["Step ",step,": ",printtable," added."]]; If[divloop===1, Print[" 1 division performed."], Print[" ",divloop," divisions performed."]]]; divisions = divisions + divloop; breakloop = (oldlength===length) ]; If[steps, Print [divisions," divisions total."]]; basis ] AB[polylist_,varlist_,opts___Rule] := AltBuchberger[polylist, varlist, opts] QB[polylist_, varlist_,opts___Rule] := QuickBuchberger[polylist, varlist, opts] QuickBuchberger[polylist_,varlist_,opts___Rule]:= Module[{ltfunc,vars=varlist,blist,basis=polylist, length=Length[polylist],indi,indj,indk, lti,ltj,ltk,lcm,criterion,divloop=0,step=1, spoly,spolyrem,loop,divisions=0, steps}, steps = PrintSteps/.{opts}/.Options[QuickBuchberger]; ltfunc = OrderList[[2]]; blist = Apply[Union,Table[{indi,indj}, {indi,1,length},{indj,indi+1,length}]]; While [!(blist==={}), indi = blist[[1]][[1]]; indj = blist[[1]][[2]]; lti = ltfunc[basis[[indi]],varlist]; ltj = ltfunc[basis[[indj]],varlist]; lcm = MonoLCM[lti,ltj,vars]; If [!NumberQ[lcm/(lti*ltj)], criterion = False; Do [If[!MemberQ[blist,{indi,indk}] &&!MemberQ[blist,{indk,indi}] &&!MemberQ[blist,{indj,indk}] &&!MemberQ[blist,{indk,indj}], ltk = ltfunc[basis[[indk]], varlist]; criterion = NumberQ[ Denominator[lcm/ltk]] ]; If[criterion, Break[]], {indk,1,length}]; If [!criterion, spoly=SPoly[basis[[indi]],basis[[indj]],vars]; spolyrem=PRemainder[spoly,basis,vars]; divisions++; divloop++; If [!(spolyrem===0), length++; AppendTo[basis,spolyrem]; If[steps, Print["Step ", step,": ", spolyrem, " added."]; Print[" ",divloop," divisions performed."]]; step++; divloop = 0; Do [AppendTo[blist,{loop,length}], {loop,length-1}] ] ] ]; blist = Drop[blist,1]; ]; If[steps, Print["Step ", step,": nothing added."]; Print[" ", divloop, " divisions performed."]; Print[divisions," divisions total."]]; basis ] MathematicaBuchberger[polylist_List, varlist_List] := GroebnerBasis[polylist, varlist, MonomialOrder -> OrderList[[3]]] MB[polylist_, varlist_] := GroebnerBasis[polylist, varlist, MonomialOrder -> OrderList[[3]]] GroebnerQ[polylist_,varlist_, torder_:Automatic]:= Module[{vars=varlist,length=Length[polylist],spoly,spolyrem, loop1,loop2,groeb=True, order, temp}, If [torder =!= Automatic, temp = MonOrder[]; MonOrder[torder]]; Do [spoly = SPoly[polylist[[loop1]],polylist[[loop2]],vars]; spolyrem = PRemainder[spoly,polylist,vars]; If[!(spolyrem===0),groeb=False; Break[]], {loop1,1,length},{loop2,loop1+1,length}]; If [torder =!= Automatic, MonOrder[temp]]; groeb ] ReduceGroebner[basis_,varlist_, torder_:Automatic]:= Module[{ltfunc,vars=varlist,length=Length[basis],basisprime=basis, divlist,loop=1,loop2,rem,lt,temp, oldorder}, oldorder = MonOrder[]; If [torder =!= Automatic, MonOrder[torder]]; While [(length>1)&&(loop<=length), divlist = Drop[basisprime,{loop,loop}]; rem = PRemainder[basisprime[[loop]],divlist,vars]; If [rem===0, basisprime = divlist; length--, basisprime[[loop]] = rem;loop++]; ]; ltfunc = OrderList[[2]]; Do [lt = ltfunc[basisprime[[loop2]],vars]; basisprime[[loop2]] = Simplify[basisprime[[loop2]]* (MonoLCM[lt,1,vars]/lt)], {loop2,1,length}]; MonOrder[oldorder]; basisprime ] IdealQ[poly_,polylist_,varlist_]:= Module[{rem,basis,oldorder}, oldorder = MonOrder[]; MonOrder[Grevlex]; rem = PRemainder[poly,polylist,varlist]; If [!(rem===0), basis = MB[polylist,varlist]; rem = PRemainder[poly,basis,varlist] ]; MonOrder[oldorder]; rem===0 ] RadicalQ[poly_, polylist_, varlist_] := Module[{listprime, varlistprime, extraly}, listprime = Append[polylist, 1 - extraly*poly]; varlistprime = Append[varlist, extraly]; GroebnerBasis[listprime, varlistprime, MonomialOrder -> DegreeReverseLexicographic] === {1} ] FiniteQ[polylist_,varlist_]:= Module [{basis, ltbasis, finite, i, varprime}, basis = GroebnerBasis[polylist, varlist, MonomialOrder -> DegreeReverseLexicographic]; ltbasis = GrevlexLT[basis, varlist]; finite = True; Do [finite = finite && Apply[Or, Map[SameQ[{varlist[[i]]}, Variables[#]]&, ltbasis]], {i, 1, Length[varlist]}]; finite ] VSDimension[polylist_,varlist_, order_:Automatic]:= Module [{oldorder, ltfunc, vars=varlist, basis, ltbasis, boundlist, freelist={}, poslist, dim, i, j, k, remlist, remlistprime, bound, elt, loop, loop2, endloop, var}, oldorder = MonOrder[]; If [order === Automatic, MonOrder[oldorder], MonOrder[order]]; ltfunc = OrderList[[2]]; endloop = Length[vars]; basis = GroebnerBasis[polylist,vars, MonomialOrder -> OrderList[[3]]]; ltbasis = ltfunc[basis, vars]; boundlist = Table[Infinity,{endloop}]; Do [freelist = Map[SameQ[{vars[[loop]]}, Variables[#]]&, ltbasis]; poslist = Flatten[Position[freelist,True]]; If [!(poslist==={}), boundlist[[loop]] = Exponent[ltbasis[[poslist[[1]]]],vars[[loop]]] ], {loop,1,endloop}]; MonOrder[oldorder]; If [Apply[Or,Map[(#===Infinity)&,boundlist]], {Infinity,Infinity}, If[Apply[Or,Map[(#==0)&,boundlist]], {0,0}, remlist = {1}; Do [var = vars[[i]]; bound = boundlist[[i]]; remlistprime = Table[remlist[[k]], {k, Length[remlist]}]; Do [elt = remlistprime[[j]]; l = 1; While[l <= (bound -1), elt = var * elt; If[!Apply[Or, Map[NumberQ[Denominator[elt/#]]&,ltbasis]], AppendTo[remlist, elt], Break[]]], {j, 1, Length[remlistprime]}], {i, 1, endloop}]; {remlist, Length[remlist]}]] ] RadicalIdealQ[polys_List, vars_] := Module[{i, grbasis, var, varsprime, grembasis, monorem, mono, blex, bvect, coor, addtoreduced, tempmatrix, monocor, dim, reduced, temppoly, j, notinideal, isradical = True}, If [!(FiniteQ[polys, vars]), Print["Error: Polynomials given do not form a zero-dimensional ideal"], grbasis = GroebnerBasis[polys, vars, MonomialOrder -> DegreeReverseLexicographic]; grembasis = VSDimension[grbasis, vars, Grevlex][[1]]; Do[var = vars[[i]]; varsprime=Drop[vars,{i,i}]; AppendTo[varsprime, var]; mono = 1; bvect = {}; dim = 0; blex = {}; While[True, monorem = PolynomialReduce[mono, grbasis, vars, MonomialOrder -> DegreeReverseLexicographic][[2]]; coor = GetCoor[monorem, grembasis]; tempmatrix = Transpose[Append[bvect, coor]]; tempmatrix = Transpose[RowReduce[tempmatrix]]; monocoor = tempmatrix[[Length[bvect] + 1]]; If [dim === Length[grembasis], addtoreduced = True, addtoreduced = (monocoor[[Length[bvect] + 1]] == 0)]; If [addtoreduced, temppoly = mono; Do[temppoly = temppoly - (monocoor[[j]]*blex[[j]]), {j, 1, Length[blex]}]; reduced = temppoly; Break[], AppendTo[blex, mono]; AppendTo[bvect, coor]; dim = dim + 1; mono = mono*var]]; deriv = D[reduced, var]; If[Exponent[PolynomialGCD[reduced,deriv],var] >0, isradical = False; Break[]], {i, 1, Length[vars]}]; isradical] ] Radical[polys_List, vars_] := Module[{i, grbasis, rbasis, var, varsprime, grembasis, monorem, mono, blex, bvect, coor, addtoreduced, tempmatrix, monocor, dim, reduced, temppoly, j}, If [!(FiniteQ[polys, vars]), Print["Error: Polynomials given do not form a zero-dimensional ideal"], grbasis = GroebnerBasis[polys, vars, MonomialOrder -> DegreeReverseLexicographic]; grembasis = VSDimension[grbasis, vars, Grevlex][[1]]; rbasis = polys; Do[var = vars[[i]]; varsprime=Drop[vars,{i,i}]; AppendTo[varsprime, var]; mono = 1; bvect = {}; dim = 0; blex = {}; While[True, monorem = PolynomialReduce[mono, grbasis, vars, MonomialOrder -> DegreeReverseLexicographic][[2]]; coor = GetCoor[monorem, grembasis]; tempmatrix = Transpose[Append[bvect, coor]]; tempmatrix = Transpose[RowReduce[tempmatrix]]; monocoor = tempmatrix[[Length[bvect] + 1]]; If [dim === Length[grembasis], addtoreduced = True, addtoreduced = (monocoor[[Length[bvect] + 1]] == 0)]; If [addtoreduced, temppoly = mono; Do[temppoly = temppoly - (monocoor[[j]]*blex[[j]]), {j, 1, Length[blex]}]; reduced = temppoly; Break[], AppendTo[blex, mono]; AppendTo[bvect, coor]; dim = dim + 1; mono = mono*var]]; deriv = D[reduced, var]; gcdpoly = PolynomialGCD[reduced, deriv]; If[Exponent[gcdpoly,var] > 0, AppendTo[rbasis, Expand[Simplify[Divide[reduced, gcdpoly]]]]], {i, 1, Length[vars]}]; rbasis] ] ConvertGB[polys_List, vars_List, origvars_List, opts___Rule] := Module[{daprint, datype, daold}, daold = OldOrderType/.{opts}/.Options[ConvertGB]; datype = OrderType/.{opts}/.Options[ConvertGB]; daprint = PrintSteps/.{opts}/.Options[ConvertGB]; If [daprint === True || daprint === False, If [datype === Grlex, ConvertGBGrlex[polys, vars, origvars, daprint, daold], If[datype === Lex, ConvertGBLex[polys, vars, origvars, daprint, daold], Print[datype, " is not a valid order to convert to. ", "Please use Lex or Grlex"]]], Print[daprint, " is not a boolean value. Please use True or False."]] ] ConvertGBLex[polys_List, vars_List, origvars_List, shouldprint_, origor_:Automatic]:= Module[{glex = {}, blex ={}, mono = 1, bvect = {}, addtoglex, order, grembasis, monorem, firstcoor, monocoor, tempmatrix, i, temppoly, k, ltglex, monoprime, j, dim = 0, oldorder, addremainder = {}}, oldorder = MonOrder[]; If [origor =!= Automatic, MonOrder[origor]; order = origor, order = oldorder]; If [!GroebnerQ[polys, origvars, order], MonOrder[oldorder]; Print["Error: The set ",polys," does not define a ",order, " Groebner basis with respect to ",origvars], grembasis = VSDimension[polys, origvars, order][[1]]; If [grembasis === Infinity, MonOrder[oldorder]; Print["Dimension of Quotient Ring is infinite; base cannot be changed."], While[True, monorem = PolynomialReduce[mono, polys, origvars, MonomialOrder -> OrderList[[3]]][[2]]; firstcoor = GetCoor[monorem, grembasis]; If[Length[bvect] == 0, addtoglex = False, tempmatrix = Append[bvect, firstcoor]; tempmatrix = Transpose[tempmatrix]; tempmatrix = Transpose[RowReduce[tempmatrix]]; monocoor = tempmatrix[[(Length[bvect] + 1)]]; If [dim === Length[grembasis], addtoglex = True, addtoglex = (monocoor[[(Length[bvect] + 1)]] == 0)]]; If[addtoglex, temppoly = mono; Do[temppoly = temppoly - (monocoor[[i]]*blex[[i]]), {i, 1, Length[blex]}]; AppendTo[glex, temppoly]; If[shouldprint, If[addremainder =!= {}, Print["Added to remainders: ", addremainder]]; Print["Added to Groebner basis: ", temppoly]; addremainder = {}]; If[Variables[LexLT[temppoly, vars]] == {vars[[1]]}, Break[]], (*termination step*) AppendTo[blex, mono]; AppendTo[addremainder, mono]; AppendTo[bvect, firstcoor]; dim = dim +1]; (*Next Monomial Step*) mono = mono * vars[[Length[vars]]]; k = Length[vars]; If[glex =!= {}, ltglex = LexLT[glex, vars]; While[Apply[Or, Map[NumberQ[Denominator[mono/#]]&,ltglex]], k = k - 1; monoprime = 1; Do [ monoprime = monoprime * vars[[j]]^(Exponent[mono, vars[[j]]]), {j, 1, k - 1}]; mono = monoprime * vars[[k]]^(Exponent[mono, vars[[k]]] + 1)] ] ] ]]; MonOrder[oldorder]; {glex, blex} ] ConvertGBGrlex[polys_List, vars_List, origvars_List, shouldprint_, origor_:Automatic] := Module[{glex = {}, blex ={}, mono = 1, bvect = {}, addtoglex, order, grembasis, monorem, firstcoor, monocoor, tempmatrix, vsresult, i, temppoly, k, ltglex, monoprime, j, dim = 0, oldorder, fixdim, boxset = False, box, boxprime, addremainder = {}}, oldorder = MonOrder[]; If [origor =!= Automatic, MonOrder[origor]; order = origor, order = oldorder]; If [!GroebnerQ[polys, origvars, order], MonOrder[oldorder]; Print["Error: The set ",polys," does not define a ",order, " Groebner basis with respect to ",origvars], vsresult = VSDimension[polys, origvars, order]; grembasis = vsresult[[1]]; fixdim = vsresult[[2]]; If [grembasis === Infinity, MonOrder[oldorder]; Print["Dimension of Quotient Ring is infinite; base cannot be changed."], While[True, monorem = PolynomialReduce[mono, polys, origvars, MonomialOrder -> OrderList[[3]]][[2]]; firstcoor = GetCoor[monorem, grembasis]; If[Length[bvect] == 0, addtoglex = False, tempmatrix = Append[bvect, firstcoor]; tempmatrix = Transpose[tempmatrix]; tempmatrix = Transpose[RowReduce[tempmatrix]]; monocoor = tempmatrix[[(Length[bvect] + 1)]]; If [dim === Length[grembasis], addtoglex = True, addtoglex = (monocoor[[(Length[bvect] + 1)]] == 0)]]; If[addtoglex, temppoly = mono; Do[temppoly = temppoly - (monocoor[[i]]*blex[[i]]), {i, 1, Length[blex]}]; AppendTo[glex, temppoly], (*termination step moved to next monomial*) AppendTo[addremainder, mono]; AppendTo[blex, mono]; AppendTo[bvect, firstcoor]; dim = dim +1]; k = Length[vars]; If[addtoglex && shouldprint, If[addremainder =!= {}, Print["Added to remainders: ", addremainder]]; Print["Added to Groebner basis: ", temppoly]; addremainder = {}]; (*Next Monomial Step*) If[glex == {}, mono = NextGrlexMon[mono, vars], ltglex = GrlexLT[glex, vars]; If [dim == fixdim && !boxset && FiniteQ[ltglex, vars], box = VSDimension[ltglex, vars, Grlex][[1]]; boxset = True]; If [boxset, If[Length[box] == fixdim, Break[]; (*termination step here*)], boxprime = {}; Do[Print["4"]; If[!NumberQ[Denominator[box[[j]]/mono]], AppendTo[boxprime, box[[j]]]], {j, 1, Length[box]}]; box = boxprime]; mono = NextGrlexMon[mono, vars]; While[Apply[Or, Map[NumberQ[Denominator[mono/#]]&,ltglex]], mono = NextGrlexMon[mono, vars]]; ] ] ]]; MonOrder[oldorder]; {glex, blex} ] NextGrlexMon[mono_, vars_List] := Module[{k, j, monoprime, degree, i}, k = Length[vars]; monoprime = 0; While[k > 1, If [(Exponent[(mono/vars[[k]]), vars[[k]]] >= 0), monoprime = 1; degree = 0; Do[degree = degree + Exponent[mono, vars[[i]]], {i, 1, Length[vars]}]; Do[If[(j < k - 1), degree = degree - Exponent[mono, vars[[j]]]; monoprime = monoprime*vars[[j]]^(Exponent[mono, vars[[j]]]), If[(j == k - 1), degree = degree - Exponent[mono, vars[[j]]] - 1; monoprime = monoprime*vars[[j]]^(Exponent[mono, vars[[j]]] + 1), If[j == Length[vars], monoprime = monoprime*vars[[j]]^degree]]], {j, 1, Length[vars]}]; Break[], k = k - 1]]; If [monoprime === 0, monoprime = vars[[Length[vars]]]^(Exponent[mono, vars[[1]]] + 1)]; monoprime ] (*This is to get the coordinates of a vector known to be in the span of the basis. This function is used in the ConvertGB functions*) GetCoor[expr_, basis_List] := Module[{coors, i, j, k, temp, extra, num}, coors = Table[0, {Length[basis]}]; Do[If[NumberQ[basis[[i]]], num = 0; exprprime = expr + extra; Do [If[NumberQ[exprprime[[k]]], num = num + exprprime[[k]]], {k, 1, Length[exprprime]}]; coors[[i]] = num/basis[[i]], temp = Coefficient [expr, basis[[i]]] + extra; Do[If[NumberQ[temp[[j]]], coors[[i]] = temp[[j]]], {j, 1, Length[temp]}] ], {i, 1, Length[basis]}]; coors ] IdealofPoints[points_List, vars_List, torder_:Automatic] := Module[{G = {}, Q = {}, S = {}, L = {1}, M = {}, min, eval, i, j, k, C, temp, vect, set, zerovect, Lprime, temptable, temporder, ltfunct, multable, ltG, iprime, jprime, varlist=vars}, If[!Apply[And, Map[SameQ[Length[vars], Length[#]]&, points]], Print["Points must be the same dimension as the length of the ", "variable list."]; Return[]]; temporder = MonOrder[]; If [torder =!= Automatic, MonOrder[torder]]; ltfunc = OrderList[[2]]; zerovect = Table[0, {Length[points]}]; While [L =!= {}, (*C2*) min = MonomialList[Sum[L[[i]], {i, Length[L]}], MonomialOrder -> OrderList[[3]]][[Length[L]]]; Do[iprime = i; If[L[[i]] === min, Break[]], {i, 1, Length[L]}]; L = Drop[L, {iprime}]; (*C3*) eval = {}; If[L==={1}, eval = Table[1, {Length[points]}], Do[temp = min; Do[temp = temp/.vars[[k]] -> points[[j]][[k]], {k, 1, Length[vars]}]; AppendTo[eval,temp], {j, 1, Length[points]}]]; vect = eval; C = {}; If[M =!= {}, Do[temp = M[[i]]; Do[jprime = j; If[temp[[j]] =!= 0, AppendTo[C,vect[[j]]/temp[[j]]]; Break[]], {j, 1, Length[points]}]; vect = vect - C[[i]]*temp, {i, 1, Length[M]}]]; (*C4 and C5*) temp = min; If [Length[S] > 0, Do[temp = temp - C[[i]]*S[[i]], {i, 1, Length[S]}]]; If[vect === zerovect, (*C4*) AppendTo[G, Expand[temp]]; Lprime = {}; Do[If[!NumberQ[Denominator[L[[i]]/min]], AppendTo[Lprime, L[[i]]]], {i, 1, Length[L]}]; L = Lprime, (*C5*) AppendTo[M, vect]; AppendTo[S, temp]; AppendTo[Q, min]; temptable = Table[min*vars[[j]], {j, 1, Length[vars]}]; ltG = ltfunct[G, vars]; Do[If[L === {} || Apply[And, Map[!NumberQ[Denominator[ temptable[[i]]/#]]&,L]], If[ltG === {} || Apply[And, Map[!NumberQ[Denominator[ temptable[[i]]/#]]&,ltG]], AppendTo[L, temptable[[i]]]]], {i, 1, Length[temptable]}] ] ]; MonOrder[temporder]; {G, Q} ] End[] Protect [OrderList,PrintSteps,OrderType,MonOrder,ElimMatrix,MonoLCM, Lex,LexLT,Grlex,GrlexLT,Grevlex,GrevlexLT,MatrixLT,PQuotient, PRemainder,SPoly,Buchberger,B,AltBuchberger,AB,QuickBuchberger, QB,MathematicaBuchberger,MB,GroebnerQ,ReduceGroebner,IdealQ, RadicalQ,FiniteQ,VSDimension,Radical,RadicalIdealQ,ConvertGB, ConvertGBLex,ConvertGBGrlex,IdealofPoints,OldOrderType] EndPackage[]