--p-adic Grobner bases newPackage("GrobnerValuations", Version => "0.9", Date => "6 September 2012", Authors => {{ Name => "Diane Maclagan", Email => "D.Maclagan@warwick.ac.uk", HomePage => "http://homepages.warwick.ac.uk/staff/D.Maclagan/"}, {Name => "Andrew Chan", Email => "a.j.chan@warwick.ac.uk", HomePage => "http://homepages.warwick.ac.uk/staff/A.J.Chan/"}}, Headline => "Grobner bases over fields with a valuation", DebuggingMode => false ) export { "pval", "Prime", "initForm", "Valuation", "normalForm", "ModPn", "grobnerVal" } --------------------------------------------------------------------------- -- CODE --------------------------------------------------------------------------- --compute the p-adic valuation of a rational number a pval = method( Options => {Prime => 2} ) pval(ZZ) := o -> a -> ( if a==0 then error("The rational number a must be nonzero"); if o.Prime<=0 then error("The prime number must be strictly positive"); val:=0; while (a % o.Prime)==0 do ( val=val+1; a=a//(o.Prime); ); return(val); ) pval(RingElement) := o -> a -> ( if liftable(a,ZZ)==false then error("pval does not support the ring element you have entered"); if o.Prime<=0 then error("The prime number must be strictly positive"); return(pval(lift(a,ZZ),Prime=>o.Prime)); ) pval(QQ) := o -> a -> ( if a==0 then error("The rational number a must be nonzero"); if o.Prime<=0 then error("The prime number must be strictly positive"); return(pval(numerator(a),Prime=>o.Prime)-pval(denominator(a),Prime=>o.Prime)); ) --compute smallest term of a polynomial f with respect to a weight vector --w and valuation pval --Note that my switching to the ring with weights seems rather cumbersome --using lots of maps and stuff, but couldn't find a better way --Outputs the lead mono smallTerm = method( Options => {Prime => 2, Valuation => pval} ) smallTerm(RingElement, List) := o ->(f,w)->( R := ring f; if #w != numgens(R) then error("The weight vector w has wrong size"); S := coefficientRing R; w' := -w; while any(w',alpha->alpha<=0) do w'= apply(w',alpha->alpha+1); T := newRing(R,Weights=>w); g := map(T, R); g' := map(R,T); f = g f; M := terms(f); I := toList(0..#M-1); E := exponents(f); C := flatten entries (coefficients f)#1; K := apply(I,alpha->flatten entries(o.Valuation(lift(C#alpha,S),Prime=>o.Prime)+matrix{w}*transpose(matrix{E#alpha}))); Kmin := min(K); L := sum apply(select(I,alpha->K#alpha==Kmin),beta->M#beta); --We have all these strange maps and stuff to that leadTerm takes into account the weights --Finding the smallest term wrt monomial order only, ignoring weights return(g' leadTerm(L)); ) --An intermediate step in the initial Form computation interForm = method( Options => {Prime => 2, Valuation => pval} ) interForm(RingElement, List) := o ->(f,w)->( R := ring f; k := coefficientRing R; M := terms(f); I := toList(0..#M-1); E := exponents(f); C := flatten entries (coefficients f)#1; K := apply(I,alpha->flatten entries(o.Valuation(lift(C#alpha,k),Prime=>o.Prime)+matrix{w}*transpose(matrix{E#alpha}))); Kmin := min(K); J := apply(select(I,alpha->K#alpha==Kmin),beta->M#beta); L := {}; for j from 0 to #J-1 do( A := coefficients(J#j); if k === QQ then ( c := lift((flatten entries(A#1))#0,k); d := lift(promote((o.Prime)^(-pval(c,Prime=>o.Prime))*c,ZZ/(o.Prime)),k) ) else ( c = lift((flatten entries(A#1))#0,ZZ); d = promote(lift(promote((o.Prime)^(-pval(c,Prime=>o.Prime))*c,ZZ/(o.Prime)),ZZ),k); ); L = L|{d*(flatten entries(A#0))#0}; ); return(sum L); ) --Compute the Leading Term of a polynomial. initForm = method( Options => {Prime => 2, Valuation => pval, ModPn =>false} ) --Compute the Leading Term of a polynomial (result in residue field) initForm(RingElement, List) := o ->(f,w)->( R := ring f; X := gens R; a := interForm(f,w,Prime=>o.Prime,Valuation=>o.Valuation); T := (ZZ/(o.Prime))(monoid R); h := map(T, R); b := h a; return(b); ) --Compute the Leading Term of an ideal (result in residue field) initForm(Ideal, List) := o ->(I,w)->( R := ring I; X := gens R; G := flatten entries grobnerVal(I,w,Prime=>o.Prime,Valuation=>o.Valuation,ModPn=>o.ModPn); g := #G; J := apply(G,alpha-> interForm(alpha,w,Prime=>o.Prime,Valuation=>o.Valuation)); T := (ZZ/(o.Prime))(monoid R); h := map(T, R); K := h matrix{J}; return(K); ) --compute the ecart of two polynomials f and g Ecart=(f,g)->( F := set(exponents(f)); G := exponents(g); n := #G; e := 0; i := 0; E := (while i < n do ( if (not(F #? (G#i))) then e=e+1; i=i+1); ); return(e); ) --compute polynomials in T whose starting terms divide starting term of f startDiv = method( Options => {Prime => 2, Valuation => pval} ) startDiv(List,RingElement,List) := o -> (w,f,T)->( R := ring f; n := numgens R; e := flatten(exponents(smallTerm(f,w,Prime=>o.Prime,Valuation=>o.Valuation))); M := {}; s := #T; i := 0; I := ( for a from 0 to (n-1) list a); S := ( while io.Prime,Valuation=>o.Valuation)) do i=i+1); for j from 0 to (s-1) do ( if all(I,x -> (S#j#x <= e#x)) then M=(M|{j}); ); return(for x in M list T#x); ) --Compute the S-polynomial of a pair of polynomials sPoly = method() sPoly(RingElement, RingElement, RingElement,RingElement) := (g,smallg,h,smallh)->( a:=(flatten entries monomials smallg)#0; b:=(flatten entries monomials smallh)#0; m:=lcm(a,b); coeffa:=coefficient(a,smallg); coeffb:=coefficient(b,smallh); S:=m*g*coeffb//(a*gcd(coeffa,coeffb))-m*h*coeffa//(b*gcd(coeffa,coeffb)); return(S); ) --Checks if one monomial divides the next one doesDivide = method() doesDivide(RingElement, RingElement) := (smallg, smallh)->( M := flatten exponents smallh-flatten exponents smallg; if all(M,alpha->alpha>=0) then return(true) else return(false); ) --compute the normal form of a homogeneous polynomial f with respect --to a list T of homogeneous polynomials. normalForm = method( Options => {Prime => 2, Valuation => pval} ) normalForm(RingElement,List,List) := o -> (f,T,w)->( R := ring f; S := coefficientRing R; n := numgens(R); P := T; smallP := apply(P,alpha->smallTerm(alpha,w,Prime=>o.Prime,Valuation=>o.Valuation)); Q := {}; Qpos := {}; s := #P; j := 0; q := f; if not(isHomogeneous(f)) then error("Polynomial(s) not homogeneous"); if not(all(T,a->isHomogeneous(a))) then error("Polynomial(s) not homogeneous"); H := new MutableHashTable; scan(0..#T-1, i -> H#(i,j) = 0); u := new MutableHashTable; u#0 = 1; r := new MutableHashTable; r#0 = 0; while q =!= 0_R do ( I := toList(0..(#P-1)); smallq := smallTerm(q,w,Prime=>o.Prime,Valuation=>o.Valuation); M := select(I,beta->all(flatten exponents smallq-flatten exponents smallP#(beta),alpha->alpha>=0)); --M is list indices of polys which divide q if #M==0 then ( r#(j+1)=r#j+smallq; q=q-smallq; u#(j+1)=u#j; for i from 0 to #T-1 do H#(i,j+1)=H#(i,j); ) else ( E:= apply(M,gamma->Ecart(q,P#gamma)); if min(E)>0 then ( Q=Q|{q}; P=P|{q}; smallP=smallP|{smallq}; Qpos=Qpos|{j}; s=#P); I' := toList(0..(#M-1)); minE := select(I',alpha->E#alpha==min(E)); delta:=M#(minE#0); g:=P#delta; --First bit is for working over the polynomial ring over QQ --and the second bit for working over ZZ/p^n if S===QQ then ( div:= smallq // smallP#delta; c:=lift(((flatten entries((coefficients(div))#1))#0),S); d:=(flatten entries((coefficients(div))#0))#0; ) else ( div= smallP#delta; coeffsmallq := lift((flatten entries((coefficients(smallq))#1))#0,ZZ); monosmallq := (flatten entries((coefficients(smallq))#0))#0; coeffdiv := lift((flatten entries((coefficients(div))#1))#0,ZZ); monodiv := (flatten entries((coefficients(div))#0))#0; c= coeffsmallq / coeffdiv; d= monosmallq // monodiv; ); a:=numerator c; b:=denominator c; rho:=b*q-a*d*g; alpha := lift((flatten entries((coefficients(smallq))#1))#0,S); if position(T,i->i==g) =!= null then ( u#(j+1)=b*u#(j); k:=position(T,omega->omega==g); (for i from 0 to #T-1 do ( if i =!= k then H#(i,j+1)=b*H#(i,j) else H#(i,j+1)=b*H#(i,j)+a*d); ); r#(j+1)=b*r#j; q=rho; ) else ( m:=position(Q,rho->rho==g); u#(j+1)=b*u#(j)-a*u#(Qpos#m); for i from 0 to #T-1 do H#(i,j+1)=b*H#(i,j)-a*H#(i,Qpos#m); r#(j+1)=b*r#j-a*r#(Qpos#m); q=rho; ); ); j=j+1; ); return(u#j,for i from 0 to #T-1 list H#(i,j),r#j); ) --Removed common factors from coefficients of polynomials --Note that we only use this when we add polynomials to the basis in the grobnerVal algorithm --We could also add this to the normalForm alg with some very small modifications to the algorithm removeCoeffFactors=(f)->( A := flatten entries (coefficients(f))#1; R := ring f; S := coefficientRing R; B := apply(A,alpha->lift(alpha,S)); d := gcd B; return(f // d); ) --Computes the reduced basis from a given basis reduce = method( Options => {Prime => 2, Valuation => pval} ) reduce(List,List) := o ->(G,w)->( Q := G; while G =!= {} do( h := G#0; G = delete(h,G); Q = delete(h,Q); k := last(normalForm(h,Q,w,Prime=>o.Prime,Valuation=>o.Valuation)); if k=!=0 then Q=Q|{k}; ); M:=apply(Q,j->removeCoeffFactors(j)); return(M); ) --BEGINNING OF GROBNER BASIS STUFF: --The following is code which computes Grobner Bases over QQ with --a valuation, in practice the p-adic one. --It is Buchberger based on work by Gebauer and Moller and taken from the --book Grobner Bases by Becker and Weispfenning --update updates the lists of polys in the Grobner basis and the S-pairs --to consider taking into consideration Buchberger's Criteria update = method() update(List,List,List) := (List1, List2, List3)->( Gold := List1#0; smallGold := List1#1; monomialsmallGold := apply(smallGold,alpha->(flatten entries monomials(alpha))#0); Bold := List2; h := List3#0; smallh := List3#1; monomialsmallh := (flatten entries monomials(smallh))#0; posh:= List3#2; C := for i from 0 to #Gold-1 list {posh,i}; D := {}; --This is while loop 1 while C =!= {} do( a := C#0; C = delete(a,C); g := Gold#(a#1); smallg := smallGold#(a#1); monomialsmallg := monomialsmallGold#(a#1); L := lcm(monomialsmallh,monomialsmallg); if L =!= monomialsmallh*monomialsmallg and ( all(C,alpha->doesDivide(lcm(monomialsmallh,monomialsmallGold#(alpha#1)),L)==false) and all(D,alpha->doesDivide(lcm(monomialsmallh,monomialsmallGold#(alpha#1)),L)==false) ) then D = D|{a}; ); --This is while loop 2 E := {}; while D =!= {} do ( b := D#0; D = delete(b,D); g = Gold#(b#1); smallg = smallGold#(b#1); monomialsmallg = monomialsmallGold#(b#1); L = lcm(monomialsmallh,monomialsmallg); if L =!= monomialsmallh*monomialsmallg then E = E|{b}; ); --This is while loop 3 Bnew := {}; while Bold =!= {} do ( c := Bold#0; Bold = delete(c,Bold); g1 := Gold#(c#0); smallg1 := smallGold#(c#0); monomialsmallg1 := monomialsmallGold#(c#0); g2 := Gold#(c#1); smallg2 := smallGold#(c#1); monomialsmallg2 := monomialsmallGold#(c#1); L = lcm(monomialsmallg1,monomialsmallg2); if doesDivide(monomialsmallh,L)==false or lcm(monomialsmallg1,monomialsmallh)==L or lcm(monomialsmallg2,monomialsmallh)==L then Bnew = Bnew|{c}; ); --This is while loop 3 Bnew = Bnew|E; Gnew := {}; smallGnew := {}; while Gold =!= {} do ( k := Gold#0; smallk := smallGold#0; monomialsmallk := monomialsmallGold#0; Gold = delete(k,Gold); smallGold = delete(smallk,smallGold); monomialsmallGold = delete(monomialsmallk,monomialsmallGold); if doesDivide(monomialsmallh,monomialsmallk)==false then ( Gnew = Gnew|{k}, smallGnew = smallGnew|{smallk} ); ); --Finalising things: Gnew = Gnew|{h}; smallGnew = smallGnew|{smallh}; return({Gnew,smallGnew},Bnew); ) --This is the new Grobner Val as following Becker/Weispfenning's GROBNERNEW2 gbVal = method( Options => {Prime => 2, Valuation => pval} ) gbVal(Ideal,List) := o -> (I,w)->( R:=ring I; F:=apply(I_*,j->removeCoeffFactors(j)); F =reduce(F,w,Prime=>o.Prime,Valuation=>o.Valuation); --reduce removes any unnecessary generators at this stage smallF := apply(F,alpha->smallTerm(alpha,w,Prime=>o.Prime,Valuation=>o.Valuation)); G := {}; smallG := {}; B := {}; while F =!= {} do( f := F#0; smallf := smallF#0; F = delete(f,F); smallF = delete(smallf,smallF); U := update({G,smallG},B,{f,smallf,#G}); G = U#0#0; smallG = U#0#1; B = U#1; ); numpairs := 0; while B =!= {} do ( monomialsmallG := apply(smallG,alpha->(flatten entries monomials(alpha))#0); A := apply(B,alpha->lcm(monomialsmallG#(alpha#0),monomialsmallG#(alpha#1))); --Choosing by min deg first degA := apply(A,alpha->degree alpha); minA := min degA; AA := select(A,alpha->degree alpha == minA); C := sum(AA); beta := (flatten entries monomials(smallTerm(C,w,Prime=>o.Prime,Valuation=>o.Valuation)))#0; gamma := position(A,i->i==beta); numpairs = numpairs +1; u := B#gamma; B = delete(u,B); g1 := G#(u#0); smallg1 := smallG#(u#0); g2 := G#(u#1); smallg2 := smallG#(u#1); S := sPoly(g1,smallg1,g2,smallg2); h := last(normalForm(S,G,w,Prime=>o.Prime,Valuation=>o.Valuation)); if h != 0_R or h!=0 then ( smallh := smallTerm(h,w,Prime=>o.Prime,Valuation=>o.Valuation); V := update({G,smallG},B,{h,smallh,#G}); G = V#0#0; smallG = V#0#1; B = V#1; ); ); --Compute a reduced minimal Grobner basis Q:=reduce(G,w,Prime=>o.Prime,Valuation=>o.Valuation); M:=(apply(G,j->j//(flatten entries((coefficients(smallTerm(j,w,Prime=>o.Prime,Valuation=>o.Valuation)))#1))#0)); return(matrix{M}); ) --The following code computed Grobner Bases over ZZ/p^n for some --suitably large n bigger than all valuations which will be seen. --We start with (I think it is) Buchberger's Second Criterion criterion = method( Options => {Prime => 2, Valuation => pval} ) criterion(List,List,List,List) := o -> (u,P,smallG,w)->( --u is {i,j} indexing the S-pair that we are considering --P is the list of S-pairs still to consider --smallG is the list of small terms of polynomials in list G --w is the weight vector i := u#0; j := u#1; R':= ring smallG#0; R:= ZZ(monoid R'); h:= map(R,R'); h':= map(R',R); L := h' lcm(h smallG#i,h smallG#j); altsmallG := delete(smallG#i,(delete(smallG#j,smallG))); A := delete(i,delete(j,toList(0..#smallG-1))); B := select(A,beta->all(flatten exponents L-flatten exponents smallG#beta,alpha->alpha>=0)); C := select(B, beta->not( member({i,beta},P)==true or member({j,beta},P)==true or member({beta,i},P)==true or member({beta,j},P)==true)); if C=={} then return(false) else return(true); ) buchbergerModPn = method( Options => {Prime => 2, Valuation => pval} ) buchbergerModPn(Ideal,List) := o -> (I,w)->( R:=ring I; R':= ZZ(monoid R); h:= map(R,R'); h':= map(R',R); G:=apply(I_*,j->removeCoeffFactors(j)); G =reduce(G,w,Prime=>o.Prime,Valuation=>o.Valuation); --reduce removes any unnecessary generators at this stage smallG := apply(G,alpha->smallTerm(alpha,w,Prime=>o.Prime,Valuation=>o.Valuation)); s:=#G; --P is a list of S-pairs remaining to be considered P:=flatten(i:=0;while i(flatten entries monomials(alpha))#0); A := apply(P,alpha->h lcm(h' monomialsmallG#(alpha#0),h' monomialsmallG#(alpha#1))); --A is a list of the lcm(lmfi),lm(fj)) for the S-pairs --We then consider the ones of the smallest degree first degA := apply(A,alpha->degree alpha); minA := min degA; AA := select(A,alpha->degree alpha == minA); C := sum(AA); beta := (flatten entries monomials(smallTerm(C,w,Prime=>o.Prime,Valuation=>o.Valuation)))#0; gamma := position(A,i->i==beta); u := P#gamma; a:=smallG#(u#0); b:=smallG#(u#1); c:=(flatten entries monomials a)#0; d:=(flatten entries monomials b)#0; m:=h lcm(h' c,h' d); P=delete(u,P); --Buchberger's Second Criterion: To remove just comment out the next line and the one further below witha bracket on it indicated if m =!=c*d and criterion(u,P,smallG,w,Prime=>o.Prime,Valuation=>o.Valuation)==false then( numpairs = numpairs +1; coeffa:=coefficient(c,a); coeffb:=coefficient(d,b); S:=m*(G#(u#0))*coeffb//(c*gcd(coeffa,coeffb))-m*(G#(u#1))*coeffa//(d*gcd(coeffa,coeffb)); r:=last(normalForm(S,G,w,Prime=>o.Prime,Valuation=>o.Valuation)); if r != 0_R or r!=0 then ( r=removeCoeffFactors(r); G=G|{r}; P=P|(i=0;while io.Prime,Valuation=>o.Valuation)}); --The bracket to delete if we do not want Buchberger's Second Criterion ); ); --We now compute a reduced, minimal basis with the following Q:=reduce(G,w,Prime=>o.Prime,Valuation=>o.Valuation); return(Q); ) --Compute the Leading Term of an ideal (result in residue field) --working Mod Pn initFormModPn = method( Options => {Prime => 2, Valuation => pval} ) initFormModPn(Ideal, List) := o ->(I,w)->( R := ring I; X := gens R; G := buchbergerModPn(I,w,Prime=>o.Prime,Valuation=>o.Valuation); g := #G; J := apply(G,alpha-> interFormModPn(alpha,w,Prime=>o.Prime,Valuation=>o.Valuation)); T := (ZZ/(o.Prime))(monoid R); h := map(T, R); K := h matrix{J}; return(K); ) --An intermediate step in the initial Form computation interFormModPn = method( Options => {Prime => 2, Valuation => pval} ) interFormModPn(RingElement, List) := o ->(f,w)->( R := ring f; S := coefficientRing R; g := smallTerm(f,w,Prime=>o.Prime,Valuation=>o.Valuation); m := (flatten entries((coefficients(g))#0))#0; return(m); ) --A subroutine in computing the Grobner basis mod P^n --Algorithm 4.2 from the paper preModPn = method( Options => {Prime => 2} ) preModPn(Ideal,List,Matrix) := o -> (I,w,J)->( L := matrix{I_*}; R := ring I; S := ring J; D := flatten unique apply(flatten entries J, alpha->degree alpha); h := map(R,S); J = h J; G := {}; for i from 0 to #D-1 do( d := D#i; M := matrix entries basis(d,I); preN := matrix entries basis(d,ideal J); monos := flatten entries matrix entries basis(d,R); N := flatten entries(J*preN); preA := L*M; --rows of preA are polynomials of deg d in I --we now need to trun this matrix into the coefficient matrix preA = transpose matrix entries (coefficients(preA,Monomials=>monos))#-1; --we now order cols so that the monos in initial ideal degree d come first positionsN := apply(N,alpha->position(monos,beta->beta==alpha)); positionsNotN := toList(0..#monos-1) - set positionsN; C := submatrix(preA,positionsN); A := C|submatrix(preA,positionsNotN); B := (inverse C)*A; F := matrix{monos}; F = matrix entries transpose (submatrix(F,positionsN)|submatrix(F,positionsNotN)); polys := flatten entries(B*F); H := select(flatten entries J,alpha->(degree(alpha))#0==d); F = flatten entries F; while H =!={} do ( h := H#0; j := position(F,alpha->alpha==h); H = delete(h,H); G = G|{polys#j}; ); ); return(G); ) --Computing the Grobner basis mod P^n --Algorithm 4.3 from the paper gbValModPn = method( Options => {Prime => 2, Valuation => pval} ) gbValModPn(Ideal,List) := o -> (I,w)->( R := ring I; S := (ZZ[]/(o.Prime)^100)(monoid R); X := flatten entries vars R; X = apply(toList(0..#X-1),alpha->X#alpha * (o.Prime)^(w#alpha)); f := map(R,R,X); h':= map(S,R); h := map(R,S); P := promote(ideal(o.Prime),R); I' := ideal apply(I_*,alpha->f alpha); I' = saturate(I',P); J := ideal apply((I')_*,alpha->h' alpha); onesVector := toList(numgens R :1); initialIdeal := initFormModPn(J,onesVector,Prime=>o.Prime,Valuation=>o.Valuation); G := preModPn(I,w,initialIdeal,Prime=>o.Prime); return(matrix{G}); ) --COMBINING THE TWO INTO ONE!!!! grobnerVal = method( Options => {Prime => 2, Valuation => pval, ModPn =>false} ) grobnerVal(Ideal,List) := o -> (I,w)->( if o.ModPn == false then ( return(gbVal(I,w,Prime=>o.Prime,Valuation=>o.Valuation))) else ( return(gbValModPn(I,w,Prime=>o.Prime, Valuation=>o.Valuation)); ); ) --------------------------------------------------------------------------- -- DOCUMENTATION --------------------------------------------------------------------------- beginDocumentation() needsPackage "SimpleDoc" debug SimpleDoc doc /// Key GrobnerValuations Headline Computes Grobner bases over a field with a valuation Description Text The goal of this package is to provide commands to compute Grobner bases of ideals over fields which are equipped with a valuation. It also contains tools to compute the initial ideal. The primary example is when the field is the p-adics. /// doc /// Key pval (pval,ZZ) (pval,QQ) (pval,RingElement) Headline Compute the p-adic valuation Usage pval a Inputs a:QQ Outputs :ZZ the p-adic valuation of a Description Text The $p$-adic valuation of rational number $a$ is the number $n$ for which $a=p^n u$ where $p$ does not divide $u$. Default valuation is the 2-adic. Example pval(2) pval(55) pval(2,Prime => 7) pval(49,Prime => 7) pval(60,Prime => 3) pval(1/2) pval(5/7) pval(3/7,Prime => 7) pval(128/3,Prime => 7) pval(128/3,Prime => 3) /// doc /// Key grobnerVal (grobnerVal,Ideal,List) Headline Compute a Grobner basis of an ideal over a rationals with the p-adic valuation Usage G = grobnerVal(I,w) Inputs I: Ideal w: List Outputs G: List Grobner basis of $I$ over a rationals with the p-adic valuation Description Text Computes the reduced Grobner Basis over a field with valuation and some weight vector $w$ of a given ideal $I$ by following Buchberger's Algorithm. Default valuation is the 2-adic Example QQ[x,y,z] grobnerVal(ideal{x-2*y,y-2*z,z-2*x},{1,1,1}) grobnerVal(ideal(7*x^4*z-18*y^5,2*x^3*z^4-3*y^7),{2,1,4},Prime=>199) grobnerVal(ideal(y^3+8*z^2*x,z^2+y*z),{1,7,3}) /// doc /// Key initForm (initForm,RingElement,List) (initForm,Ideal,List) Headline Compute the initial ideal of an ideal, or the initial form of a polynomial Usage I = initForm(v,w) Inputs v: RingElement or @ofClass Ideal@ w: List Outputs : RingElement the initial form of polynomial v with respect to weight vector w, or the initial ideal of an ideal. Output lives in the polynomial ring over the residue field Description Text The initial form of a polynomial is the monomial $ax^u$ which with respect to the weight vector $w$ and the $p$-adic valuation, has $pval(c)+w.u$ minimal amongst all monomials. It lives in the polynomial ring over the residue field. Example QQ[x,y,z] initForm(7*x*y+8*y^2+8*z^2+x*z,{1,2,1}) initForm(3*x*y+19*x^2+9*y^2,{1,1,1},Prime=>3) initForm(4*x^3*y+6*x*y+17,{3,2,1},Valuation=>pval) Text If the input is a homogenous ideal $I$, then the function first uses @TO grobnerVal@ to compute the Grobner basis of $I$ with respect to the selected valuation before computing the initial forms of this basis. The default valuation is the 2-adic. Example I=ideal(x+y+z,x^2+4*y^2+3*x*z,x*y*z+18*y^3+12*x*y^2) initForm(I,{1,3,7}) /// doc /// Key normalForm (normalForm,RingElement,List,List) Headline Compute the normal form of a polynomial with respect to a list of polynomials Usage r = normalForm(f,T,w) Inputs f: RingElement T: List w: List Outputs u: RingElement H: List r: RingElement the initial form of polynomial v with respect to weight vector w, or the initial ideal of an ideal. Output lives in the polynomial ring over the residue field Description Text Computes the normal form of a homogeneous polynomial over a field with valuations with respect to a given set of polynomials. If we a given a polynomial $f$ with list ${g_1,...,g_s}$ which is the list of polynomials with which we are trying to divide $f$ then the normal form is $u*f=h_1*g_1+...+h_s*g_s+r$. Default valuation is the 2-adic Example QQ[x,y,z] normalForm(13*x^3*y*z^2+12*y^6+8*x^6+y^2*z^4,{x*z+y^2*12,x^2+12*z^2+y*z},{7,2,1}) normalForm(x^10+y^10+z^10,{x+17*z,y-8*x,z-128*x},{1,3,2},Prime=>31) /// doc /// Key Prime [grobnerVal,Prime] [normalForm,Prime] [initForm,Prime] [pval,Prime] Headline Change the p-adic valuations prime Usage grobnerVal(...,Prime=>p) normalForm(...,Prime=>p) initForm(...,Prime=>p) pval(...,Prime=>p) Description Text This optional input allows us to change the p-adic valuation on the field. The default is the 2-adic. Below we show some examples of how to vary the prime in some examples. Example QQ[x,y,z] pval(2) pval(2,Prime => 7) initForm(3*x*y+19*x^2+9*y^2,{1,1,1}) initForm(3*x*y+19*x^2+9*y^2,{1,1,1},Prime=>3) normalForm(13*x^3*y*z^2+12*y^6+8*x^6+y^2*z^4,{x*z+y^2*12,x^2+12*z^2+y*z},{7,2,1}) normalForm(13*x^3*y*z^2+12*y^6+8*x^6+y^2*z^4,{x*z+y^2*12,x^2+12*z^2+y*z},{7,2,1},Prime=>7) /// doc /// Key Valuation [grobnerVal,Valuation] [normalForm,Valuation] [initForm,Valuation] Headline Change the field valuation Usage grobnerVal(...,Valuation=>val) normalForm(...,Valuation=>val) initForm(...,Valuation=>val) pval(...,Valuation=>val) Description Text This optional input allows us to change the valuation of the field. Currently this only the p-adic valuation is implemented, however the user may wish to define their own. The default is the 2-adic. We give examples of how it could in principle be used. Example QQ[x,y,z] initForm(3*x*y+19*x^2+9*y^2,{1,1,1},Valuation=>pval) normalForm(13*x^3*y*z^2+12*y^6+8*x^6+y^2*z^4,{x*z+y^2*12,x^2+12*z^2+y*z},{7,2,1},Valuation=>pval) /// doc /// Key ModPn [grobnerVal,ModPn] [initForm,ModPn] Headline Compute the Grobner Basis over Z/p^nZ for suitably large n Usage grobnerVal(...,ModPn=>false) initForm(...,ModPn=>false) Description Text This optional input allows us to do some of the Grobner basis computations in the polynomial ring over $\mathbb{Z}/p^N\mathbb{Z}$ for some suitably large $N$. Even though the algorithms theoretically terminate in finite time, the coefficients may get too large and thus require memory exceeding that of standard comuters. This can be remedied by doing some computations mod $p^N$ where $N$ is chosen to be larger than any valuation appearing in a reduced Grobner basis. The result is then lifted to the original polynomial ring. The default is not to implement ModPn as in general it will not be required. We implemented ModPn as we experienced extreme coefficient growth when computing Mustafin varieties. The example below is an example of one such computation of a Mustafin Variety: Example R=QQ[a..f] x1=matrix{{a},{b},{c}} x2=matrix{{d},{e},{f}} A1=matrix{{6,8,4},{16,12,8},{2,4,18}} A2=matrix{{16,4,5},{7,4,9},{6,7,3}} M=matrix{flatten entries(A1*x1),flatten entries(A2*x2)} I=minors(2,M) grobnerVal(I,{1,1,1,1,1,1},ModPn=>true) Text As a point of interest, it turned out that this particular example would terminate fine without having to work mod $p^N$. We do the computation here as a demonstration that we get the correct result working mod $p^N$. Example grobnerVal(I,{1,1,1,1,1,1},ModPn=>false) /// --doc /// -- Key -- "How to use this Package" -- Headline -- A brief explanation of how to use the package --/// --------------------------------------------------------------------------- -- TESTS ---------------------------------------------------------------------------