Eisenstein Fourier Coefficients Program
version 2.0
(written for Mathematica 10.0)

by David Yuen
yuen@lakeforest.edu
and Jerry Shurman
jerry@reed.edu

FlatGenusSymbol Program (execute this group)

Aux programs

The following routines are here because Mathematica 9 broke previous programs written in earlier versions.

ZeroMatrix[n_,m_]:=Table[0,{n},{m}];ZeroMatrix[n_]:=Table[0,{n},{n}];BlockMatrix[x___]:=ArrayFlatten[x];SubMatrix[a_,upperLeft_,dim_]:=Module[{},
If[a==={},Return[{}]];
If[(dim[[1]]===0)||(dim[[2]]===0),
{},
a[[Range[upperLeft[[1]],upperLeft[[1]]+dim[[1]]-1],Range[upperLeft[[2]],upperLeft[[2]]+dim[[2]]-1]]]
]
];
MyFactorInteger[n_]:=If[n===1,{},FactorInteger[n]];

Clear[Valuation];
Valuation[x_,p_]:=IntegerExponent[Numerator[x],p]-IntegerExponent[Denominator[x],p];

ElementaryMatrix1[n_,i_,j_]:=Module[{a},
a=IdentityMatrix[n];
a[[i,i]]=0;a[[j,j]]=0;
a[[i,j]]=1;a[[j,i]]=1;
a
]

ElementaryMatrix3[n_,i_,j_,alp_]:=Module[{a},
a=IdentityMatrix[n];
a[[i,j]]=alp;
a
]

MyReplaceBox[a_,b_,ii_,jj_]:=Module[{t=a,i,j},
Do[t[[ii-1+i,jj-1+j]]=b[[i,j]],
{i,1,Length[b]},{j,1,Length[b[[1]]]}];
t
];

MyDiagonalizeOddPrime[m_,p_]
MyDiagonalizeOddPrimeList[mx_,p_]

MyMinEntry[m_]:=Module[{ans,i,j,x},
ans={1,1};x=m[[1,1]];
Do[If[m[[i,i]]<x,x=m[[i,i]];ans={i,i}],{i,1,Length[m]}];
Do[If[m[[i,j]]<x,x=m[[i,j]];ans={i,j}],{i,1,Length[m]-1},{j,i+1,Length[m]}];
ans
];

MyDiagonalizeOddPrime[mx_,p_]:=Module[{m=mx,v,i,j,n=Length[mx],minv,minij,
Em, m2},
If[n==1,Return[mx]];
v=Table[Valuation[m[[i,j]],p],{i,1,n},{j,1,n}];
minij=MyMinEntry[v];
minv=v[[minij[[1]],minij[[2]]]];
(*Print[{v,minij,minv}];*)
If[minij[[1]]=!=1,
Em=ElementaryMatrix1[n,1,minij[[1]]];
m=Em.m.Em;
v=Table[Valuation[m[[i,j]],p],{i,1,n},{j,1,n}];
minij=MyMinEntry[v];
minv=v[[minij[[1]],minij[[2]]]];
(*Print[{v,minij,minv}];*)
];
If[minij[[2]]=!=1,
Em=ElementaryMatrix3[n,1,minij[[2]],1];
m=Em.m.Transpose[Em];
v=Table[Valuation[m[[i,j]],p],{i,1,n},{j,1,n}];
minij=MyMinEntry[v];
minv=v[[minij[[1]],minij[[2]]]];
(*Print[{v,minij,minv}];*)
];
(*Print[m//MatrixForm];*)
Do[
If[Valuation[m[[1,j]]/m[[1,1]],p]<0,
Print["ERROR 23525124626"];Return[Indeterminate]];
Em=ElementaryMatrix3[n,j,1,-m[[1,j]]/m[[1,1]]];
m=Em.m.Transpose[Em];
(*Print[m//MatrixForm];*)
,{j,2,n}
];
m2=SubMatrix[m,{2,2},{n-1,n-1}];
m2=MyDiagonalizeOddPrime[m2,p];
m=MyReplaceBox[m,m2,2,2];
m
];

MyDiagonalizeOddPrimeList[mx_,p_]:=Module[{m=mx,v,i,j,n=Length[mx],minv,minij,
Em, m2},
If[n==1,Return[Flatten[mx]]];
v=Table[Valuation[m[[i,j]],p],{i,1,n},{j,1,n}];
minij=MyMinEntry[v];
minv=v[[minij[[1]],minij[[2]]]];
(*Print[{v,minij,minv}];*)
If[minij[[1]]=!=1,
Em=ElementaryMatrix1[n,1,minij[[1]]];
m=Em.m.Em;
v=Table[Valuation[m[[i,j]],p],{i,1,n},{j,1,n}];
minij=MyMinEntry[v];
minv=v[[minij[[1]],minij[[2]]]];
(*Print[{v,minij,minv}];*)
];
If[minij[[2]]=!=1,
Em=ElementaryMatrix3[n,1,minij[[2]],1];
m=Em.m.Transpose[Em];
v=Table[Valuation[m[[i,j]],p],{i,1,n},{j,1,n}];
minij=MyMinEntry[v];
minv=v[[minij[[1]],minij[[2]]]];
(*Print[{v,minij,minv}];*)
];
(*Print[m//MatrixForm];*)
Do[
If[Valuation[m[[1,j]]/m[[1,1]],p]<0,
Print["ERROR 23525124626"];Return[Indeterminate]];
Em=ElementaryMatrix3[n,j,1,-m[[1,j]]/m[[1,1]]];
m=Em.m.Transpose[Em];
(*Print[m//MatrixForm];*)
,{j,2,n}
];
m2=SubMatrix[m,{2,2},{n-1,n-1}];
m2=MyDiagonalizeOddPrimeList[m2,p];
Join[{m[[1,1]]},m2]
];

GenusSymbolOdd[m_,p_]

GenusSymbolOdd[m_,p_]:=Module[
{mL,i,mV,mVU,ans,nfq,detfq,v,nd,e,q},
mL=MyDiagonalizeOddPrimeList[m,p];(*Print[mL];*)
mV=Table[Valuation[mL[[i]],p],{i,1,Length[mL]}];
mVU=Union[mV];ans={};
Do[nfq=0;detfq=1;v=mVU[[i]];q=p^v;
Do[
If[mV[[j]]==v,
nfq++;
detfq=detfq*mL[[j]]/q;
],
{j,1,Length[mV]}
];
nd=Numerator[detfq]*Denominator[detfq];(*Print[detfq];*)
If[Mod[nd,p]==0,Print["ERROR 4354514"];Return[Indeterminate]];
e=JacobiSymbol[nd,p];
AppendTo[ans,{q,e*nfq}];
,
{i,1,Length[mVU]}
];
ans
]

GenusSymbolTwo[m_]
fixed

Valuation2[x_]:=Valuation[x,2];

MyDiagonalizeTwoList[mx_]:=Module[{m=mx,v,i,j,n=Length[mx],minv,minij,
Em, m2,p=2,X,b,Y},
If[n==1,Return[{mx}]];
v=Table[Valuation[m[[i,j]],p],{i,1,n},{j,1,n}];
minij=MyMinEntry[v];
minv=v[[minij[[1]],minij[[2]]]];
(*Print[{v,minij,minv}];*)
If[minij[[1]]=!=1,
Em=ElementaryMatrix1[n,1,minij[[1]]];
m=Em.m.Em;
v=Table[Valuation[m[[i,j]],p],{i,1,n},{j,1,n}];
minij=MyMinEntry[v];
minv=v[[minij[[1]],minij[[2]]]];
(*Print[{v,minij,minv}];*)
];
If[minij[[2]]==1,
Do[
If[Valuation[m[[1,j]]/m[[1,1]],p]<0,
Print["ERROR 46726411134"];Return[Indeterminate]];
Em=ElementaryMatrix3[n,j,1,-m[[1,j]]/m[[1,1]]];
m=Em.m.Transpose[Em];
(*Print[m//MatrixForm];*)
,{j,2,n}
];
m2=SubMatrix[m,{2,2},{n-1,n-1}];
m2=MyDiagonalizeTwoList[m2];
Return[Join[{{{m[[1,1]]}}},m2]]
,
If[n==2,Return[{m}]];
If[minij[[2]]≠2,
Em=ElementaryMatrix1[n,2,minij[[2]]];
m=Em.m.Em
];
b=SubMatrix[m,{1,1},{2,2}];
X=SubMatrix[m,{3,1},{n-2,2}];
Y=-X.Inverse[b]; (*Check valences are nonnegative*)
If[Min[Valuation2/@Flatten[Y]]<0,
Print["ERROR 84581155"];Return[Indeterminate]];
Em=BlockMatrix[{{IdentityMatrix[2],Transpose[Y]},{Y,IdentityMatrix[n-2]}}];
m=Em.m.Transpose[Em];
m2=SubMatrix[m,{3,3},{n-2,n-2}];
m2=MyDiagonalizeTwoList[m2];
Return[Join[{ b },m2]]

];
Indeterminate

];

MyJacobiSymbol2[x_]:=Module[{a},
a=Mod[x,8];
Which[a==1,1,
a==3,-1,
a==5,-1,
a==7,1,
True,0]
]

GenusSymbolTwo[m_]:=Module[
{mL,i,mV,mVU,ans,nfq,detfq,v,nd,e,q,gr,newgr,oddEven,oddity,t,p=2},
mL=MyDiagonalizeTwoList[m];(*Print[mL];*)
mV=Table[Valuation[mL[[i,1,-1]],p],{i,1,Length[mL]}];
(*Use upper right hand corner for valence! *)
mVU=Union[mV];ans={};
Do[
nfq=0;detfq=1;v=mVU[[i]];q=p^v;gr={};
Do[
If[mV[[j]]==v,
nfq+=Length[mL[[j]]];
detfq=detfq*Det[mL[[j]]/q];
AppendTo[gr,mL[[j]]];
],
{j,1,Length[mV]}
];
oddEven=Union[Length/@gr];
If[oddEven=={1,2},
gr=q*MyFurtherReduce[gr/q];
oddEven={1};
];
nd=Numerator[detfq]*Denominator[detfq];(*Print[detfq];*)
If[Mod[nd,2]==0,Print["ERROR 35261134"];Return[Indeterminate]];
e=MyJacobiSymbol2[nd];
(*Print[{gr,oddEven,detfq,nd,e}];*)
If[oddEven=={1},
(*Type1*)
t=Sum[gr[[j,1,1]],{j,1,Length[gr]}];
t=t/q;
oddity=Mod[Numerator[t]*Denominator[t],8];
AppendTo[ans,{q,e*nfq,1,oddity}];

];
If[oddEven=={2},
(*Type2*)
AppendTo[ans,{q,e*nfq,2}];

];
,
{i,1,Length[mVU]}
];
ans
]

MyFurtherReduce[g_]:=Module[
{ans={},oddI,z,i,res},
Do[If[Length[g[[i]]]==1,oddI=i;Break[]],{i,1,Length[g]}];
z=g[[oddI,1,1]];
Do[
If[i≠oddI,
If[Length[g[[i]]]==1,
AppendTo[ans,g[[i]]],
res=MyTrick[z,g[[i]]];
ans=Join[ans,{{{res[[2]]}},{{res[[3]]}}}];
z=res[[1]];
]
],{i,1,Length[g]}
];
AppendTo[ans,{{z}}];
ans
];

MyTrick[z_,mx_]:=Module[
{a=mx[[1,1]],b=mx[[1,2]],c=mx[[2,2]]},
{c+z,(-b^2+a*c+a*z)/(c+z),(-(b^2*z)+a*c*z)/(-b^2+a*c+a*z)}
];

GenusSymbol[m_]
GenusSymbolCompact[m_]
FlatGenusSymbol[m_]

GenusSymbol[m_]:=Module[{li,p,i,ans},
li=FactorInteger[Det[m]][[All,1]];
li=Complement[li,{2}];
ans={{2,GenusSymbolTwo[m]}};
Do[
p=li[[i]];
AppendTo[ans,{p,GenusSymbolOdd[m,p]}];
,
{i,1,Length[li]}];
ans
];

GenusSymbolCompact[m_]:=Module[{li,part,i,j,ans={}},
li=GenusSymbol[m];
Do[
part=li[[i,2]];
Do[
If[part[[j,1]]≠1,AppendTo[ans,part[[j]]]]
,{j,1,Length[part]}
]
,
{i,1,Length[li]}];
ans
];

FlattenGenusSymbol[gs_]:=Module[{ans,n,i,j,a,b,p,li},
a=gs[[1]];
n=Apply[Plus,Abs[a[[2,All,2]]]];
ans={n,1};
li={};
For[i=1,i≤Length[gs],i++,
{p,a}=gs[[i]];b={};
For[j=1,j≤Length[a],j++,
If[a[[j,1]]≠1,
ans[[2]]*=a[[j,1]]^(Abs[a[[j,2]]]);
AppendTo[b,a[[j]]]
]
];
If[b=!={},
AppendTo[li,{p,Sort[b]}]
]
];
li=Reverse[Sort[li]];
ans=Join[ans,Flatten[li[[All,2]]]];
ans
]

FlatGenusSymbol[b_]:=FlattenGenusSymbol[GenusSymbol[b]];

GenusSymbolForPerlPreparation[m_]

UnFlatGS[li_]:=Module[{gs,j},
gs={li[[1]]};j=3;
While[j≤Length[li],
If[Mod[li[[j]],2]==0,
If[li[[j+2]]==1,
AppendTo[gs,{li[[j]],li[[j+1]],li[[j+2]],li[[j+3]]}];j+=4;,
AppendTo[gs,{li[[j]],li[[j+1]],li[[j+2]]}];j+=3;
]
,
AppendTo[gs,{li[[j]],li[[j+1]]}];j+=2;
]
];
gs
]

GenusSymbolForPerlPreparation[m_]:=UnFlatGS[FlatGenusSymbol[m]]

MoniEisensteinCoefficient Programs (execute this group)

Aux programs

MyFactorInteger[n_]:=If[n===1,{},FactorInteger[n]];

MyJacobiSymbolPrime[d_,p_]:=Module[{},
If[p===2,
Switch[Mod[d,8],
0,0,
1,1,
4,0,
5,-1,
_,Print["ERROR in MyJacobiSymbol: input not a fundamental discriminant"];Indeterminate
]
,
JacobiSymbol[d,p]
]
]

MyChiD[d_,a_]:=Module[{ans,i,f,p,e},
If[a===0,
If[d===1,1,0],
f=MyFactorInteger[a];
ans=1;
For[i=1,i≤Length[f],i++,
{p,e}=f[[i]];
If[p==-1,ans*=Sign[d]^e,
ans*=MyJacobiSymbolPrime[d,p]^e]
];
ans
]
]

BernoulliPolynomial[n_,x_]:=Module[{i},
Sum[Binomial[n,i]*BernoulliB[i]*x^(n-i),{i,0,n}]
];

TwistedBkD[k_,d_]:=Module[{a},
Abs[d]^(k-1)*Sum[MyChiD[d,a]*BernoulliPolynomial[k,a/Abs[d]],{a,1,Abs[d]}]
];

MyValence[x_,p_]:=IntegerExponent[Numerator[x],p]-IntegerExponent[Denominator[x],p];
MyFundDisc[d_]:=Module[{fL,i,s},
fL=MyFactorInteger[d];
s=Product[fL[[i,1]]^(Mod[fL[[i,2]],2]),{i,1,Length[fL]}];
(*s=square free part*)
If[Mod[s,4]===1,s,4s]
]

MyLFuncOneMinusK[k_,d_]:=-TwistedBkD[k,d]/k;

MonicEisensteinCoefficient[k_, gs_, fp_,x_]
ScaledEisensteinCoefficient[n_, k_, gs_, fp_,x_]
gs must be in flattened genus symbol format

EnkFactorEvenGenusGivenMatrix[n_,k_,B_]:=Module[{det2B,dB,fB,j},
(*detB=gs[[2]]/2^n; *)(*Remember B is half integral, really*)
det2B=Det[B];
dB=MyFundDisc[(-1)^(n/2)det2B];
fB=Sqrt[(-1)^(n/2)det2B/dB];
(*Print[(-1)^(n/2)det2B==dB*fB^2];*)
2^(n/2)/(Zeta[1-k]*Product[Zeta[1-(2k-2j)],{j,1,n/2}])*MyLFuncOneMinusK[(k-n/2),dB]
];
GetfBGivenMatrix[n_,k_,B_]:=Module[{det2B,dB,fB,j},
(*detB=gs[[2]]/2^n; *)(*Remember B is half integral, really*)
det2B=Det[B];
dB=MyFundDisc[(-1)^(n/2)det2B];
fB=Sqrt[(-1)^(n/2)det2B/dB];
fB
];
EnkFactorEvenGenusGS[n_,k_,gs_]:=Module[{det2B,dB,fB,j},
If[Mod[n,2]=!=0,Print["n not even in EnkFactorEvenGenusGS"];
Return[Indeterminate]];
(*detB=gs[[2]]/2^n; *)(*Remember B is half integral, really*)
det2B=gs[[2]];
dB=MyFundDisc[(-1)^(n/2)det2B];
fB=Sqrt[(-1)^(n/2)det2B/dB];
(*Print[(-1)^(n/2)det2B==dB*fB^2];*)
2^(n/2)/(Zeta[1-k]*Product[Zeta[1-(2k-2j)],{j,1,n/2}])*MyLFuncOneMinusK[(k-n/2),dB]
];
GetfBGivenGSEvenGenus[n_,k_,gs_]:=Module[{det2B,dB,fB,j},
(*detB=gs[[2]]/2^n; *)(*Remember B is half integral, really*)
det2B=gs[[2]];
dB=MyFundDisc[(-1)^(n/2)det2B];
fB=Sqrt[(-1)^(n/2)det2B/dB];
fB
];
EnkFactorOddGenusGS[n_,k_,gs_]:=Module[{j},
If[Mod[n,2]=!=1,Print["n not odd in EnkFactorOddGenusGS"];
Return[Indeterminate]];2^((n+1)/2)/(Zeta[1-k]*Product[Zeta[1-(2k-2j)],{j,1,(n-1)/2}])
];

MonicEisensteinCoeffient[k_,gs_,fp_,x_]:=Module[{n},
n=gs[[1]];
If[ Mod[n,2]==1,
MonicEisensteinCoeffientOddGenus[k,gs,fp,x],
MonicEisensteinCoeffientEvenGenus[k,gs,fp,x]
]
];
MonicEisensteinCoeffientEvenGenus[k_,gs_,fp_,x_]:=Module[
{i,myfp,efactor,fb,fs,oneCoeff,n},
n=gs[[1]];
If[Mod[k,2]=!=0,Print["weight must be even"];Return[Indeterminate]];
If[Mod[n,2]=!=0,Print["n must be even"];Return[Indeterminate]];
If[n=!=gs[[1]],Print["n does not match that in gs"];Return[Indeterminate]];
If[(Mod[n/2,4]==1)&&(k<n/2+3),Print["(Mod[n/2,4]==1)&&(k<n/2+3)"];Return[Indeterminate]];
If[(Mod[n/2,4]=!=1)&&(k<n/2+1),Print["(Mod[n/2,4]=!=1)&&(k<n/2+1)"];Return[Indeterminate]];

Clear[myfp];
For[i=1,i≤Length[fp],i++,
myfp[fp[[i,1]]]=fp[[i,2]];
];
efactor=EnkFactorEvenGenusGS[n,k,gs];
fb=GetfBGivenGSEvenGenus[n,k,gs];
fs=MyFactorInteger[fb];
oneCoeff=efactor*Product[myfp[fs[[i,1]]]/.{x->fs[[i,1]]^(k-n-1)},
{i,1,Length[fs]}];
oneCoeff
];
MonicEisensteinCoeffientOddGenus[k_,gs_,fp_,x_]:=Module[
{i,myfp,efactor,fb,fs,oneCoeff,det2B,n},
n=gs[[1]];
If[Mod[k,2]=!=0,Print["weight must be even"];Return[Indeterminate]];
If[Mod[n,2]=!=1,Print["n must be odd"];Return[Indeterminate]];
If[k≤n+1,Print["k>n+1 not true"];Return[Indeterminate]];
If[n=!=gs[[1]],Print["n does not match that in gs"];Return[Indeterminate]];

det2B=gs[[2]];
Clear[myfp];
For[i=1,i≤Length[fp],i++,
myfp[fp[[i,1]]]=fp[[i,2]];
];
efactor=EnkFactorOddGenusGS[n,k,gs];

If[Mod[det2B,2]==0,
fs=MyFactorInteger[det2B/2],
fs={MyFactorInteger[det2B]}
];
oneCoeff=efactor*Product[myfp[fs[[i,1]]]/.{x->fs[[i,1]]^(k-n-1)},
{i,1,Length[fs]}];
oneCoeff

];

(* NOTE Dilation factor. *)
ScriptEnkFactorEven[n_,k_]:=Module[{j},(2^(-n/2)*Zeta[1-k]*Product[Zeta[1-(2k-2j)],{j,1,n/2}])];ScaledEisensteinCoefficient[n_, k_, gs_, fp_,x_]:=Module[{},
ScriptEnkFactorEven[n,k]*MonicEisensteinCoeffient[k,gs,fp,x]
];

Example
Note: the Fp polynomials must be computed externally, say using King’s LISP program.

Example 6x6 matrix, weight 16

a=(({{2, 1, 1, 0, 1, 2}, {1, 4, 2, 2, 0, 1}, {1, 2, 4, 2, 0, 0}, {0, 2, 2, 4, 2, 2}, {1, 0, 0, 2, 4, 2}, {2, 1, 0, 2, 2, 8}}));Det[a]

GenusSymbol[a]

gs=FlatGenusSymbol[a]

inputToKingsProgram=GenusSymbolForPerlPreparation[a]

Output of King’s program yields
FpData[{6,{3, -1}, {4, -2, 1, 4}}]
={{3,1},{2,1+24*x^1+256*x^2+3072*x^3+16384*x^4}};

fp={{3,1},{2,1+24*x^1+256*x^2+3072*x^3+16384*x^4}};

6x6 matrix, weight 16

weight=16

monicCoeff=MonicEisensteinCoeffient[weight, gs, fp,x]

FactorInteger[Denominator[monicCoeff]]

ScaledEisensteinCoefficient[6, weight,gs, fp,x]

Same 6x6 matrix, weight 18

weight=18

monicCoeff=MonicEisensteinCoeffient[weight, gs, fp,x]

FactorInteger[Denominator[monicCoeff]]

ScaledEisensteinCoefficient[6, weight,gs, fp,x]

Same 6x6 matrix, weight 16, but as a coefficient in degree 8 Eisenstein series

weight=16

monicCoeff=MonicEisensteinCoeffient[weight, gs, fp,x]

FactorInteger[Denominator[monicCoeff]]

ScaledEisensteinCoefficient[8, weight,gs, fp,x]