Computation of the Hilbert series of the algebra of regular functions on a symplectic quotient by SU2

Hans-Christian Herbig (herbig@imf.au.dk), Daniel Herden (daniel_herden@baylor.edu), Christopher Seaton (seatonc@rhodes.edu)
A modification of the algorithm of Pedro Cayres (pedrocayres@poli.ufrj.br)
Aug 8, 2018

This computes the on-shell Hilbert series (the Hilbert Series of the ring of real regular invariant functions on the associated symplectic quotient). See arXiv:1809.07760 [math.SG] .

Define all functions and then run

compOn[d, m, e, n]

where
d= list of ODD dk, distinct
m= list of multiplicities of the elements of d
e= list of EVEN dk, distinct
n= list of multiplicities

The above weights are for the COTANGENT LIFTED representation, so all multiplicities must be even.

Note that this does not use the telescoping (it computes three sums), so it could be speeded up.

The function newFactor takes a polynomial p = poly in the variable z as input and factors out all factors of the form (1 - z^k) for positive integer k.  The output is a list with two entries; the first entry is a list indicating how many factors occur at each degree (e.g. {0, 1, 0, 2} corresponds to (1 - z)(1 - z^3)^2) and the second entry is the remaining factor of p which is not divisible by any (1 - z^k).

newFactor[p_]:=Module[{poly=p,k, polyFact},
(* Set up a blank table to record the factors. *)
polyFact=Table[0,{i,1, Exponent[poly,z]}];
(* set k = degree of the input polynomial.  We will search for factors starting with the highest possible degree k. *)
k=Exponent[poly,z];
While[k≥1,
(* If 1-z^k divides the polynomial, divide out this factor from poly, record this in the table polyFact, and check again. If not (or if the statement depends on the value of z, i.e. Mathematica interprets it as neither true nor false), increment k down. *)
If[PolynomialRemainder[poly, 1-z^k, z]==0,
polyFact[[k]]++;
poly=PolynomialQuotient[poly,1-z^k,z],
k=k-1,
k=k-1];];
(* Remove 0's at the end of the table polyFact. *)
k=Length[polyFact];
While[polyFact[[k]]==0,polyFact=Delete[polyFact,k];k=k-1;];
{polyFact,poly}]

The function manipPoly takes a list p = poly representing a product of factors of the form (1 - z^k) for positive integer k (identical to the first element of the output of newFactor, e.g. {0, 1, 0, 2} corresponds to (1 - z)(1 - z^3)^2) as well as a positive integer weight w = weight.  The output is a list with two entries; the first entry is a list indicating how many factors occur at each degree (e.g. {0, 1, 0, 2} corresponds to (1 - z)(1 - z^3)^2) and the second entry is the remaining factor of p which is not divisible by any (1 - z^k).  It changes each factor (1 - z^k) to (-1)^(gcd(k,w)+1)(1-z^(lcm (k,w)/w)^gcd(k,w).  The output is the resulting polynomial.

manipPoly[p_,weight_]:=Module[{poly=p,w=weight,sign,k,manPolyList,manPoly},
(* Set up a blank table manPolyList to record the factors of the new polynomial. *)
manPolyList=Table[0,{i,1,Length[poly]}];
(* As k changes from 1 to the highest degree of factors in poly, set the kth entry in manPolyList to be that of poly if gcd(k,w)=1 (no changes to the factors in this case).  If gcd(k,w) ≠ 1, then record the sign change (if any) and the new factor in manPolyList. *)
k=1;
While[k≤Length[poly],
If[GCD[k,w]==1,manPolyList[[k]]=manPolyList[[k]]+poly[[k]],
manPolyList[[LCM[w,k]/w]]=manPolyList[[LCM[w,k]/w]]+poly[[k]]*GCD[k,w];
];
k++];
(* Multiply the factors in manPolyList together to return a polynomial expression manPoly. *)
k=1;
manPoly=Product[(1-z^i)^manPolyList[[i]],{i,1,Length[manPolyList]}];
manPoly]

The function shiftSeries takes a rational function r = rat in z and a positive integer weight w = weight.  Letting D denote the degree of the denominator of rat, it computes the Maclaurin series of rat up to degree (D + 2)w.  It then shifts the series with weight a, i.e. replaces the series ∑(Subscript[c, k])(x^k) with the series ∑(Subscript[c, wk])(x^k).  The output is the shifted series (now up to degree (D + 2)) as a polynomial.

shiftSeries[rat_,weight_]:=Module[{r=rat,w=weight,k,poly,series},
series=CoefficientList[Series[r, {z, 0,w*(2+Exponent[Denominator[r],z])}],z];
(* A blank list to keep track of shifted series data. *)
newSeries={};
k=1;
(* Record in newSeries the coefficients of the MacLaurin series of rat in degrees divisible by w. *)
While[k≤Length[series],
newSeries=Append[newSeries, series[[k]]];
k+=w];
k=1;
(* Add the terms of the shifted series together in poly. *)
poly=0;
While[k≤(Exponent[Denominator[r],z]+2)&&k≤Length[newSeries],
poly=poly+newSeries[[k]]*z^(k-1);
k++];
poly]

Returns the multiplicities of all of the weights.

calcMult[degrees_, multiplicities_] := Module[{d = degrees, m = multiplicities, cm, i, j, cmI},
cm ={};
For[i = 0, i ≤ Max[d], i++,
cmI = 0;
For[j = 1, j ≤ Length[d],j++, If[Abs[Max[d] - 2 * i] ≤ d[[j]], cmI += m[[j]]];];
cm = Append[cm, cmI]];
cm]

Returns the parcel of the Hilbert Series associated to the term 1/(1 - tz^(d1 - 2 i1))^mi_i1.

calcPart[i1_,d1_,mi_,e1_,ni_,l_,lPole_] := Module[{a,indicesI,indicesJ,Phi,PhiList,PhiNum,PhiDen,PhiNew,PhiTilde,PhiTildeNum,PhiTildeDen ,PhiShiftedNum,PhiShiftedDen,PhiShiftedSum,i,j,m,aux,alShift},
indicesI=DeleteCases[Range[0,d1],i1];
indicesJ=Range[0,e1];
PhiTilde=Simplify[z^l(1-z^2)/((Product[(1-t*z^(d1-2*i))^mi[[i+1]],{i,indicesI}])*(Product[(1-t*z^(e1-2*j))^ni[[j+1]],{j,indicesJ}]))];
PhiTildeNum=Numerator[PhiTilde];
PhiTildeDen=Denominator[PhiTilde];
a=d1-2*i1;
(* computes the adjustment required to prevent the numerator of Phi from being a Laurent polynomial *)
alShift=If[lPole>0,Ceiling[lPole/Abs[a]+.5],0];
PhiList={};
For[i=0,i<mi[[i1+1]],i++,
m=mi[[i1+1]]-i;
(*Calculates Phi_i, a rational function that is multiplied by 1/(z(1 - tz^a)^(mi_i1-i)) in the integral*)
Phi=(-z^(-a))^i*(PhiTildeNum/.{t->z^(-a)})/(PhiTildeDen/.{t->z^(-a)})/Factorial[i];
PhiNum=Numerator[Phi];
PhiDen=Denominator[Phi];
(*Shifts the numerator and denominator of Phi_i*)
PhiShiftedDen=manipPoly[newFactor[PhiDen][[1]],-a];(**newFactor[PhiDen][[2]];*)
PhiShiftedNum=z^(-alShift)Normal[Series[shiftSeries[Phi*z^(alShift*Abs[a]),-a]*PhiShiftedDen,{z,0,Exponent[PhiShiftedDen,z] }]];
(*Saves the result in a list*)
PhiList=Append[PhiList,D[z^(m-1)*PhiShiftedNum/PhiShiftedDen,{z,m-1}]/Factorial[m-1]];
(* removed Simplify[ before D[ *)
PhiTildeNum=D[PhiTildeNum,{t,1}]*PhiTildeDen-D[PhiTildeDen,{t,1}]*PhiTildeNum;
(* removed Simplify[ before PhiTi *)
PhiTildeDen=PhiTildeDen^2];
(*Sum over all shifted Phi_i*)
PhiShiftedSum=Sum[PhiList[[i]],{i,1, mi[[i1+1]]}];
PhiShiftedSum]

This function computes the on-shell Hilbert series.

compOn[d_List,m_List,e_List,n_List]:=Module[{a,l,start,fin,d1,mi,e1,ni,lPole,PhiList,i1,hilbSerEven,hilbSerOdd,hilbSerN,hilbSer0,hilbSerP,h0,h2,h,i2,PsiList,PhiNew},
l=-2;
start=SessionTime[];
d1=Max[d];
mi=calcMult[d,m];
e1=Max[e];
ni=calcMult[e,n];
lPole=2;
PhiList={};
i1=Quotient[d1,2]+1;
While[i1≤d1,
PhiNew[z_] =calcPart[i1,d1,mi,e1,ni,l,lPole];
PhiList=Append[PhiList,PhiNew[z]];
i1++];
(* removed Simplify[ before Sum[ *)
hilbSerEven=Sum[PhiList[[i]],{i,1, d1-Quotient[d1,2]}];

PsiList={};
i2=Quotient[e1,2]+1;
While[i2≤e1,
PsiList=Append[PsiList,PsiNew=calcPart[i2,e1,ni,d1,mi,l,lPole]];
i2++];
(* removed Simplify[ before D[ *)
hilbSerOdd=Sum[PsiList[[i]],{i,1, e1-Quotient[e1,2]}];

hilbSerN = hilbSerEven + hilbSerOdd;

l=0;
start=SessionTime[];
d1=Max[d];
mi=calcMult[d,m];
e1=Max[e];
ni=calcMult[e,n];
lPole=0;
PhiList={};
i1=Quotient[d1,2]+1;
While[i1≤d1,
PhiNew[z_] =calcPart[i1,d1,mi,e1,ni,l,lPole];
PhiList=Append[PhiList,PhiNew[z]];
i1++];
(* removed Simplify[ before Sum[ *)
hilbSerEven=Sum[PhiList[[i]],{i,1, d1-Quotient[d1,2]}];

PsiList={};
i2=Quotient[e1,2]+1;
While[i2≤e1,
PsiList=Append[PsiList,PsiNew=calcPart[i2,e1,ni,d1,mi,l,lPole]];
i2++];
(* removed Simplify[ before Sum[ *)
hilbSerOdd=Sum[PsiList[[i]],{i,1, e1-Quotient[e1,2]}];

hilbSer0 = hilbSerEven + hilbSerOdd;


l=2;
start=SessionTime[];
d1=Max[d];
mi=calcMult[d,m];
e1=Max[e];
ni=calcMult[e,n];
lPole=0;
PhiList={};
i1=Quotient[d1,2]+1;
While[i1≤d1,
PhiNew[z_] =calcPart[i1,d1,mi,e1,ni,l,lPole];
PhiList=Append[PhiList,PhiNew[z]];
i1++];
(* removed Simplify[ before Sum[ *)
hilbSerEven=Sum[PhiList[[i]],{i,1, d1-Quotient[d1,2]}];

PsiList={};
i2=Quotient[e1,2]+1;
While[i2≤e1,
PsiList=Append[PsiList,PsiNew=calcPart[i2,e1,ni,d1,mi,l,lPole]];
i2++];
(* removed Simplify[ before Sum[ *)
hilbSerOdd=Sum[PsiList[[i]],{i,1, e1-Quotient[e1,2]}];

hilbSerP = hilbSerEven + hilbSerOdd;

h0[z_]:=hilbSer0;
h2[z_]:=hilbSerN+hilbSer0+hilbSerP;
h[z_]:=h0[z]-z^2h2[z]+z^4h2[z]-z^6 h0[z]//Together;
Print["The on-shell Hilbert Series is: ",Together[h[z]]];
deg=Sum[(d[[i]]+1)*m[[i]],{i,1,Length[d]}]+Sum[(e[[i]]+1)*n[[i]],{i,1,Length[e]}];
Print["The Laurent Series is: ",Series[h[z],{z,1,12-deg}]];
fin=SessionTime[];
Print["The duration of this evaluation was ",fin-start," seconds."]
];

One example application (V2+V3)

compOn[{3},{2},{2},{2}]