The algorithm described in Section 3.3 of “The Hilbert series of symplectic quotients by the 2-torus” by Hans-Christian Herbig, Daniel Herden, and Christopher Seaton to compute the Hilbert series of the real regular functions on the symplectic quotient associated to a representation of the 2-torus.

herbighc@gmail.com
daniel herden@baylor.edu
seatonc@rhodes.edu

Definition of the Hilbert series function

Takes a single 2xn matrix (the weight matrix) as an argument, tests that it is generic and in standard form, and computes the Hilbert series as well as some expansions

hilbertSeries[mat_]:=Module[{a=mat,n,nn,d},
n=Length[a[[1]] ];
nn=Table[i,{i,1,n}];
(* define the determinant function *)
d[i_,j_]:=a[[1,i]]a[[2,j]]-a[[1,j]]a[[2,i]];

(* check to see if the matrix is in standard form and abort if not *)
If[Abs[ a[[1]] ]≠a[[1]]||Product[  a[[1,i]],{i,nn}]==0||DeleteDuplicates[   a[[1]]  ] ≠a[[1]],
Print["The matrix ",MatrixForm[a]," is not in standard form. The entries in the first row need to be distinct and positive. You can fix this by row operations or multiplying columns by -1."],

(* check to see if the matrix is effective and abort if not *)
(* define the list of all determinants as a sequence *)
dSeq=Flatten[Table[d[i,j],{i,nn},{j,nn}]]/.{stuff___}:>(stuff);
(* define the gcd of the determinants *)
gcd=GCD[dSeq];
(* define the rank of the weight matrix *)
rank=MatrixRank[a];
If[rank<2||gcd>1,
Print["The matrix ",MatrixForm[a]," is not effective; the rank is ",rank," and the gcd of the 2x2 minors is ",gcd],

(* check to see if the matrix is generic and abort if not *)
degenCheck=0;
i=1;
While[i≤n,
j=1;
While[j≤n,
k=1;
While[k≤n,
If[k≠i&&k≠j&&d[i,j]>0&&(  (d[i,j]+d[i,k]+d[j,k])==0   ||   (d[i,j]-d[i,k]-d[j,k])==0    ),
degenCheck=1;
iDeg=i;
jDeg=j;
kDeg=k;
];
k++];
j++];
i++];
If[degenCheck==1,
Print["The matrix ",MatrixForm[a]," is degenerate. When i = ",iDeg,", j = ",jDeg,", and k = ",kDeg,", we have d[i,j]+d[i,k]+d[j,k])= ",(d[iDeg,jDeg]+d[iDeg,kDeg]+d[jDeg,kDeg])," and d[i,j]-d[i,k]-d[j,k]= ",(d[iDeg,jDeg]-d[iDeg,kDeg]-d[jDeg,kDeg]),"."],
Clear[i,j,k];

(* here the matrix passes the checks, so we proceed to the computation *)
(* hilb is the final answer; once we compute each twice-shifted term, we add it to hilb *)
hilb=0;
(* we do the process for each i,j such that d[i,j] > 0 *)
i=1;
While[i≤n,
j=1;
While[j≤n,
If[d[i,j]>0,
(* define the term corresponding to i,j as a list *)
(* list is of the form:
{numerator constant, numerator monomial as list {sDeg,tDeg}, List of denom factors in the same form}
 *)
term={1,{0,0}, Join[ 
Table[{d[i,k], d[i,j]+d[j,k]},  {k,Complement[nn,{i,j}] }] ,
Table[{-d[i,k], d[i,j]-d[j,k]},  {k,Complement[nn,{i,j}] }] 
]};


(* if the t-power is negative, multiply the whole factor and numerator by a monomial so that the t-power is positive *)
counter=1;
While[counter<Length[term[[3]]   ],
If[term[[3,counter,2]]<0,
term[[3,counter,1]]=-term[[3,counter,1]];
term[[3,counter,2]]=-term[[3,counter,2]];
term[[1]]=-term[[1]];
term[[2]]=term[[2]]+term[[3,counter]];
];
counter++];
Clear[counter];

(* compute the t-degrees of the (old) numerator and denominator before shifting *)
tDegNum=term[[2,2]];
tDegDenom=Sum[term[[3,counter,2]],{counter,1,Length[term[[3]]  ]}];
(* compute the difference of t-degrees of the (new) numerator and (new) denominator after shifting *)
degBoundDiff=Floor[   (tDegNum-tDegDenom)/d[i,j]   ];

(* compute the manipulated (new) denominator as a list *)
newDenomList={};
counter1=1;
While[counter1≤Length[term[[3]] ],
counter2=1;
While[counter2<=GCD[  d[i,j],term[[3,counter1,2]]   ],
AppendTo[newDenomList,{term[[3,counter1,1]]d[i,j]/GCD[  d[i,j],term[[3,counter1,2]]  ], LCM[  d[i,j], term[[3,counter1,2]]   ]/d[i,j]}];
counter2++];
counter1++];
Clear[counter1,counter2];
(* compute the manipulated (new) denominator as a polynomial *)
newDenomPoly=Product[(1 - s^(newDenomList[[counter,1]]) t^(newDenomList[[counter,2]])),{counter,1,Length[newDenomList ]} ];

(* compute the t-degree of the manipulated (new) denominator *)
tDegDenomShifted=Sum[newDenomList[[counter,2]],{counter,1,Length[newDenomList  ]}];

(* compute the t-degree of the (new) numerator *)
tDegBound=tDegDenomShifted+degBoundDiff ;


(* apply the t-shift *)
(* first simplify the term to be shifted (to speed up the Series expansion *)
toExpand=FullSimplify[
(*numerator*)(term[[1]]*s^(term[[2,1]])*t^(term[[2,2]]))/(*denominator*)Product[(1 - s^(term[[3,counter,1]]) t^(term[[3,counter,2]])),{counter,1,Length[term[[3]] ]} 
]];
(* expand the term up to the degree bound and store the coefficients (functions in s) as a list termSeries *)
termSeries=CoefficientList[Series[toExpand,{t,0,tDegBound*d[i,j]}],t];
(* shift termSeries (i.e. apply (Subscript[U, dij])) *)
termSeriesShifted=Table[termSeries[[counter]],{counter,1,Length[termSeries],d[i,j]}];

(* express the shifted series as a polynomial, multiply by the shifted (new) numerator, and trucate at t-degree tDegBound *)
newNum=Normal[Series[Sum[termSeriesShifted[[counter]]*t^(counter-1),{counter,1,Length[termSeriesShifted]}]*newDenomPoly,{t,0,tDegBound}]   ];

(* t-shift is complete *)
(* now to perform the s-shift *)

(* express the result as a list of the form:
{numerator as a polynomial in s and t, List of denom factors in the same form {sDeg,tDeg} }
 *)
term={Expand[newNum],newDenomList};


(* if the s-power is negative, multiply the whole factor and numerator by a monomial so that the s-power is positive *)
counter=1;
While[counter<=Length[term[[2]]   ],
If[term[[2,counter,1]]<0,
term[[2,counter,1]]=-term[[2,counter,1]];
term[[2,counter,2]]=-term[[2,counter,2]];
term[[1]]=-term[[1]]*s^(term[[2,counter,1]])*t^(term[[2,counter,2]]);
];
counter++];
Clear[counter];


(* compute the s-degrees of the (old) numerator and denominator before shifting *)
sDegNum=Exponent[   term[[1]], s];
sDegDenom=Sum[term[[2,counter,1]],{counter,1,Length[term[[2]]  ]}];
(* compute the difference of s-degrees of the (new) numerator and (new) denominator after shifting *)
degBoundDiff=Floor[   (sDegNum-sDegDenom)/d[i,j]   ];


(* compute the manipulated (new) denominator as a list *)
newDenomList={};
counter1=1;
While[counter1≤Length[term[[2]] ],
counter2=1;
While[counter2<=GCD[  d[i,j],term[[2,counter1,1]]   ],
AppendTo[newDenomList,{LCM[  d[i,j], term[[2,counter1,1]]   ]/d[i,j],term[[2,counter1,2]]d[i,j]/GCD[  d[i,j],term[[2,counter1,1]]  ]}];
counter2++];
counter1++];
Clear[counter1,counter2];
(* compute the manipulated (new) denominator as a polynomial *)
newDenomPoly=Product[(1 - s^(newDenomList[[counter,1]]) t^(newDenomList[[counter,2]])),{counter,1,Length[newDenomList ]} ];

(* compute the s-degree of the manipulated (new) denominator *)
sDegDenomShifted=Sum[newDenomList[[counter,1]],{counter,1,Length[newDenomList  ]}];

(* compute the s-degree of the (new) numerator *)
sDegBound=sDegDenomShifted+degBoundDiff ;


(* apply the s-shift *)
(* first simplify the term to be shifted (to speed up the Series expansion *)
toExpand=FullSimplify[
(*numerator*)term[[1]]/(*denominator*)Product[(1 - s^(term[[2,counter,1]]) t^(term[[2,counter,2]])),{counter,1,Length[term[[2]] ]} ]];
(* first expand the term up to the degree bound and express the result as a Laurent polynomial termSeriesPoly. *)
termSeriesPoly=Normal[Series[toExpand,{s,0,sDegBound*d[i,j]}],Assumptions->{s>0,t>0}   ];


(* the result might be a Laurent polynomial, so figure out how much we have to shift to make it a polynomial *)
sMinDeg=Floor[  Min[Exponent[term,s,Min],0]/d[i,j]   ];
termSeries=CoefficientList[termSeriesPoly*s^(-sMinDeg*d[i,j]),  s];


(* shift termSeries (i.e. apply (Subscript[U, dij,s])) *)
termSeriesShifted=Table[termSeries[[counter]],{counter,1,Length[termSeries],d[i,j]}];


(* express the shifted series as a polynomial, multiply by the shifted (new) numerator, and trucate at s-degree tDegBound *)
newNum=Normal[Series[
Sum[termSeriesShifted[[counter]]*s^(counter-sMinDeg-1),{counter,1,Length[termSeriesShifted]}]*newDenomPoly,{s,0,sDegBound}]   ];


(* add the completed computation to hilb *)
hilb=hilb+((newNum/newDenomPoly)//.s->t)
(* the next bracket ends the If for d[i,j]>0 *)
];
(* the next bracket ends the computation for the chosen i,j. *)
j++];
i++];
Clear[i,j,k];

(* combine the terms to a single rational function and print *)
hilb=Together[hilb];

Print[hilb];
Print[Series[hilb,{t,0,40}]];
Print[Series[hilb,{t,1,8-2n}]]
(* ends the three If statements checking degenerate, effective, and standard form. *)
]]]
];

Examples

wMat={{1,2,3},{0,-1,-2}};
hilbertSeries[wMat]//Timing

wMat={{1,2,3,4},{0,-1,-2,-3}};
hilbertSeries[wMat]//Timing