\\Pari package for calculating some coefs of Siegel cusp modular forms
\\R. E. Borcherds, version july 18 1996
\\Please send comments, bugs etc to reb@math.berkeley.edu
\\To use save this in a file called (say) siegel after editing
\\out any junk mail has added, and load it into pari using the command \rsiegel

\\This line seems to be needed to run this file under newer
\\versions of pari
default(compatible,3)

print1("Calculating stuff. Please wait..")

\\Root systems are represented as vectors with 48 components, 
\\where the first 24 give the number of An's, (1<=n<=24), 
\\the next 21 give the number of Dn's (4<=n<=24),
\\and the last 3 components give the number of E6's, E7's, and E8's. 

a(n)=if(n<=24, vector(48,x,x==n ), vector(48,x,0));
d(n)=\
	if(n<=24 && n>=4, vector(48,x,x==n+21),\
	if(n==2, a(1)*2,\
	if(n==3, a(3),\
	a(0))));
e(n)=\
	if(n>=6, vector(48,x,x==n+40),\
	if(n==1, a(1),\
	if(n==2, a(1),\
	if(n==3, a(2)+a(1),\
	if(n==4, a(4),\
	if(n==5, d(5),\
	a(0) ))))));

\\Order of the full automorphism group of connected root system j. 
autroot1=vector(48,x,\
	if(x==1,2,\
	if(x<=24,2*(x+1)!,\
	if(x==25,4!*2^3*6,\
	if(x<=45,(x-21)!*2^(x-21-1)*2,\
	if(x==46,6!*72*2,\
	if(x==47,7!*576,\
	if(x==48,8!*2*3*4*5*6*4*2*3,\
0))))))));

\\Order of full aut group of root system m
autroot(m)=prod(1,j=1,48,autroot1[j]^m[j]*m[j]!)

\\determinant of connected root system j.
detroot1=vector(48,x,\
	if(x<=24,x+1,\
	if(x<=45,4,\
	49-x)));
\\Order of determinant of root system m
detroot(m)=prod(1,j=1,48,detroot1[j]^m[j])

\\rank of root system
rankroot1=vector(48,x,if(x<=24,x,if(x<=45,x-21,x-40)));
rankroot(m)=sum(0,j=1,48,rankroot1[j]*m[j])

choose(n,m)=if(m>=0, prod(1,j=1,m,(n+1-j)/j), 0)

print1("a...")
niemeier(n)=[\
	a(0),\
	a(1)*24,\
	a(2)*12,\
	a(3)*8,\
	a(4)*6,\
	a(5)*4+d(4),\
	d(4)*6,\
	a(6)*4,\
	a(7)*2+d(5)*2,\
	a(8)*3,\
	a(9)*2+d(6),\
	d(6)*4,\
	a(11)+d(7)+e(6),\
	e(6)*4,\
	a(12)*2,\
	d(8)*3,\
	a(15)+d(9),\
	d(10)+e(7)*2,\
	a(17)+e(7),\
	d(12)*2,\
	a(24),\
	d(16)+e(8),\
	e(8)*3,\
	d(24)\
][n];

\\ embed(m,n) is the number of isometries of root system m
\\in root system n, for any root systems m,n. 
embed(m,n,  j,k,x)=\
j=0; for(k=1,48,if(m[k]!=0,j=k,0));\
if(j==0,1,\
sum(0,k=j,48,\
if(embed1[j,k] && n[k],\
n[k]*embed1[j,k]*embed(m-vector(48,x,x==j),n+complement1[j,k]-vector(48,x,x==k))\
+if(embed2[j,k] && n[k],\
n[k]*embed2[j,k]*embed(m-vector(48,x,x==j),n+complement2[j,k]-vector(48,x,x==k)),0)\
,0)))

\\embed1 and embed2 give the number of isometries from
\\root system j to root system k where 1<=j,k<=48 represent
\\connected dynkin diagrams. (embed2 is for when there are 2 orbits.)
print1("b...")
embed2=matrix(48,48,j,k,\
	autroot1[j]*\
	if(j==3 && k>=25 && k<=45, choose(k-21,3),\
	if(j==5 && k==47, 2^4*3^2*7,\
	if(j==7 && k==48, 2^6*3^3*5,\
0))));


\\complement 1 and complement 2 give the orthogonal complements of
\\root system j in root system k (1<=j,k<=48).
print1("c...")
complement2=matrix(48,48,j,k,\
	if(j==3 && k>=25 && k<=45, d(k-21-3),\
	if(j==5 && k==47, a(1),\
	if(j==7 && k==48, a(0),\
a(0)))));

print1("d...")
complement1=matrix(48,48,j,k,\
	if(j>=k,a(0),\
	if(k<=24,a(k-j-1),\
	if(k<=45 && j==1, d(k-21-2)+a(1),\
	if(k<=45 && j<=24,d(k-21-j-1),\
	if(k<=45, d(k-j),\
	if(k==46 && j==1, a(5),\
	if(k==46 && j==2, a(2)*2,\
	if(k==46 && j==3, a(1)*2,\
	if(k==46 && j==4, a(1),\
	if(k==46 && j==5, a(1),\
	if(k==46 && j==25, a(0),\
	if(k==46 && j==26, a(0),\
	if(k==46 && j==46, a(0),\
	if(k==46, a(0),\
	if(k==47 && j==1, d(6),\
	if(k==47 && j==2, a(5),\
	if(k==47 && j==3, a(3)+a(1),\
	if(k==47 && j==4, a(2),\
	if(k==47 && j==5, a(2),\
	if(k==47 && j==6, a(0),\
	if(k==47 && j==7, a(0),\
	if(k==47 && j==25, a(1)*3,\
	if(k==47 && j==26, a(1),\
	if(k==47 && j==27, a(1),\
	if(k==47 && j==46, a(0),\
	if(k==47 && j==47, a(0),\
	if(k==47, a(0),\
	if(k==48 && j==1, e(7),\
	if(k==48 && j==2, e(6),\
	if(k==48 && j==3, d(5),\
	if(k==48 && j==4, a(4),\
	if(k==48 && j==5, a(2)+a(1),\
	if(k==48 && j==6, a(1),\
	if(k==48 && j==7, a(1),\
	if(k==48 && j==8, a(0),\
	if(k==48 && j==25, d(4),\
	if(k==48 && j==26, a(3),\
	if(k==48 && j==27, a(1)*2,\
	if(k==48 && j==28, a(0),\
	if(k==48 && j==29, a(0),\
	if(k==48 && j==46, a(2),\
	if(k==48 && j==47, a(1),\
	if(k==48 && j==48, a(0),\
	if(k==48,a(0),\
a(0) )))))))))))))))))))))))))))))))))))))))))))));

print1("e...")
embed1=matrix(48,48,j,k,\
	autroot1[j]*\
	if(j>k,0,\
	if(k<=24,choose(k+1,j+1),\
	if(k<=45 && j<=24,choose(k-21,j+1)*2^j,\
	if(k<=45, choose(k-21,j-21),\
	if(k==46 && j==1, 2^2*3^2,\
	if(k==46 && j==2, 2^3*3*5,\
	if(k==46 && j==3, 2*3^3*5,\
	if(k==46 && j==4, 2^3*3^3,\
	if(k==46 && j==5, 2^2*3^2,\
	if(k==46 && j==25, 3^2*5,\
	if(k==46 && j==26, 3^3,\
	if(k==46 && j==46, 1,\
	if(k==46, 0,\
	if(k==47 && j==1, 3^2*7,\
	if(k==47 && j==2, 2^4*3*7,\
	if(k==47 && j==3, 2^2*3^2*5*7,\
	if(k==47 && j==4, 2^5*3^2*7,\
	if(k==47 && j==5, 2^4*3*7,\
	if(k==47 && j==6, 2^5*3^2,\
	if(k==47 && j==7, 2^2*3^2,\
	if(k==47 && j==25, 3^2*5*7,\
	if(k==47 && j==26, 2*3^3*7,\
	if(k==47 && j==27, 3^2*7,\
	if(k==47 && j==46, 2^2*7,\
	if(k==47 && j==47, 1,\
	if(k==47, a(0),\
	if(k==48 && j==1, 2^3*3*5,\
	if(k==48 && j==2, 2^5*5*7,\
	if(k==48 && j==3, 2^3*3^3*5*7,\
	if(k==48 && j==4, 2^7*3^3*7,\
	if(k==48 && j==5, 2^7*3^2*5*7,\
	if(k==48 && j==6, 2^8*3^3*5,\
	if(k==48 && j==7, 2^5*3^3*5,\
	if(k==48 && j==8, 2^6*3*5,\
	if(k==48 && j==25, 2*3^2*5^2*7,\
	if(k==48 && j==26, 2^3*3^3*5*7,\
	if(k==48 && j==27, 2^2*3^3*5*7,\
	if(k==48 && j==28, 2^3*3^3*5,\
	if(k==48 && j==29, 3^3*5,\
	if(k==48 && j==46, 2^5*5*7,\
	if(k==48 && j==47, 2^3*3*5,\
	if(k==48 && j==48, 1,\
	if(k==48,a(0),\
a(0) ))))))))))))))))))))))))))))))))))))))))))));

print("f...")
\\ From Conway and Sloane, page 413
denomoveraut =[\
	15570572852330496000,\
	31522712171959008000000,\
	312927932591898624000000,\
	437599241673834240000000,\
	180674574584719324741632,\
	52278522738634063872000,\
	1196560426451890500000,\
	8361079854908571648000,\
	2700612462901377024000,\
	225800767686574080000,\
	106690862731906252800,\
	19144966823230248000,\
	8082641116053504000,\
	373503391765504000,\
	834785957117952000,\
	156983146327507500,\
	33307587016704000,\
	4134535541136000,\
	3483146354688000,\
	67271626831500,\
	4173688995840,\
	271057837050,\
	63804560820,\
	24877125];
\\numerator and denominator of the mass for 24 dimensional lattices. 
num24=1027637932586061520960267;
den24=129477933340026851560636148613120000000;
mass24=num24/den24;

\\ Orders of automorphism groups of Niemeier lattices.
autnie(n)= den24/denomoveraut[n];


\\print("calculating matrix:")
\\innie(m,k)=embed(m,niemeier(k))
\\b=matrix(24,24,j,k,innie(if(j<=12,a(j-1),if(j<=15,e(j-7),d(j-12))),k));
\\print("inverting matrix")
\\cuspcoef=(1/b)*vvector(24,x,x==24);
cuspcoef= [1/152769576960, -1/3183476736, 1/591224832, -5/1146617856, 13/1990656000, -83/11943936000, 19/16307453952, 41/15676416000, -1/37158912000, -197/351151718400, 59/214990848000, -13/229323571200, -1/35831808000, 1/143327232000, 31/7685922816000, 37/11415217766400, -29/21069103104000, -1/7023034368000, 1/3511517184000, 53/4237899595776000, -1/1332620771328000, -1/3595793596416000, 1/10787380789248000, 1/2729207339679744000];


\\Coefficient of Siegel cusp form of genus 12, weight 12. 
cusp(m)=sum(0,j=1,24,cuspcoef[j]*embed(m,niemeier(j)))

eiscoef=denomoveraut/num24;
\\Coefficient of eisenstein series of weight 12 (const term=1)
e12(m)=sum(0,j=1,24,eiscoef[j]*embed(m,niemeier(j)))

aute82=autroot(e(8)*2);
halfautd16=autroot(d(16))/2;
mass16=1/aute82 + 1/halfautd16;
e8(m)=(embed(m,e(8)*2)/aute82 + embed(m,d(16))/halfautd16)/mass16

e4(m)=embed(m,e(8))
mass8=1/autroot(e(8));

\\Print a root system.
printroot(v)=for(j=1,48,printroot1(j,v[j]));
printroot1(j,k)=\
if(k>0,if(j<=24,print1("a",j),if(j<=45,print1("d",j-21),print1("e",j-40)));\
if(k>1,print1("^",k," "),print1(" "))  ,0);


print("a(n),d(n),e(n) root systems which can be added, e.g. a(3)*4+e(6)")
print("niemeier(j)    the root system of niemeier lattice j (1<=j<=24)")
print("autnie(j)      the order of the automorphism group of Niemeier lattice j")
print("autroot(m)     the order of the aut group of root system m")
print("detroot(m)     the determinant of root system m")
print("rankroot(m)    the rank of the root system m")
print("printroot(m)   will print the root system m")
print("embed(m,n)     the number of isometries from root system m to n")
print("cusp(m)        the coefficient corresponding to m of the cusp form of weight 12")
print("e12(m),e8(m),e4(m) coefficients of the Eisenstein series of weight 12,8,4")
print("mass24=num24/den24,mass16,mass8  M-S mass for lattices of dimension 24,16,8")