(* A Mathematica program to evaluate the Hecke polynomial for a given 
prime, weight, level, and character.  The basic command is
      hecke[p_, wt_, lev_, ch_]
where p is a prime number, wt is a positive integer, lev is a
positive integer, and ch is a list whose length divides lev
and whose entries are just the values of a character modulo its
length. For example,
         hecke[3, 12, 5, {1,-1,-1,1,0}]
will return
    \!\(18385718256 + 341568\ x\^2 + x\^4\);
the roots are the values of a_3 for weight 12 level 5 modular form 
with the Legendre symbol mod 5 as the character.      
   Also, dim[wt,lev,ch] will return the dimension of the space.
Also, heckenew[p_,wt_,lev_,ch_] returns the characteristic polynomial
 of the appropriate newform space. Similarly, for dimnew[wt,lev,ch]
Note: This program does not work for p which divide the level.
 *)

crem[ls1_, ls2_] := 
  Block[{lcm = LCM[ls1[[2]], ls2[[2]]], t = 0, n = 0}, 
   While[t == 0, If[Mod[n, ls1[[2]]] == Mod[ls1[[1]],ls1[[2]]] && 
        Mod[n, ls2[[2]]] == Mod[ls2[[1]],ls2[[2]]], t = 1, 
		If[n > lcm, t = 1]]; n++]; n - 1]
 
psi[n_] := n*Block[{a = Divisors[n]}, 
    Product[If[PrimeQ[a[[m]]] == True, 1 + 1/a[[m]], 1], 
		{m, 1, Length[a]}]]
 
chival[ch_, n_, lev_] := If[Length[ch]==1,1,If[GCD[n, lev] == 1,
		 ch[[Mod[n, Length[ch]]]], 0]]

chival[ch_,n_]:=chival[ch,n,Length[ch]]

A3[n_,wt_,lev_,ch_]:=-
    Block[{aa=Divisors[n],laa,nchi=Length[ch],
	bb=Divisors[lev],
        lbb},laa=Length[aa];lbb=Length[bb];
      Sum[aa[[m1]]^(wt-1)*
          Sum[If[Mod[GCD[lev/nchi,n/aa[[m1]]-aa[[m1]]],
                GCD[bb[[m2]],lev/bb[[m2]]]]==0,
                EulerPhi[
                  GCD[bb[[m2]],lev/bb[[m2]]]]*
                      chival[ch,
                        crem[{aa[[m1]],bb[[m2]]},
				{n/aa[[m1]],lev/bb[[m2]]}],lev],
              0],{m2,1,lbb}]*
			If[m1==(laa+1)/2,1/2,1],{m1,1,(laa+1)/2}]]

A4[n_,wt_,lev_,ch_]:=If[wt==2&&Length[ch]==1,aaa=Divisors[n];
	Sum[If[GCD[lev,n/aaa[[m]]]==1,aaa[[m]],0],{m,1,Length[aaa]}],0]

A1[n_,wt_,lev_,ch_]:=If[IntegerQ[Sqrt[n]]==True,
	n^(wt/2-1)*chival[ch,Sqrt[n],lev]*(wt-1)/12*psi[lev],0]

mu[t_,ff_,n_,lev_,ch_]:=Block[{nf=GCD[lev,ff]},
	psi[lev]/psi[lev/nf]*Sum[If[Mod[x^2-t*x+n,lev*nf]==0,
		chival[ch,x,lev],0],{x,1,lev}]]

A2[n_,wt_,lev_,ch_]:=Block[{lim=If[IntegerQ[Sqrt[n]]==True,
	Sqrt[4*n]-1,Floor[Sqrt[4*n]]]},
	Simplify[-1/2*Sum[rho=(t+Sqrt[t^2-4*n])/2;
	rhob=(t-Sqrt[t^2-4*n])/2;aaaa=Divisors[t^2-4*n];
	(rho^(wt-1)-rhob^(wt-1))/(rho-rhob)*
	Sum[ff=Sqrt[aaaa[[mm]]];
	If[IntegerQ[ff]==True,hc[(t^2-4*n)/ff^2]*
	mu[t,ff,n,lev,ch],0],
		{mm,1,Length[aaaa]}],{t,-lim,lim}]]]

trace[n_,wt_,lev_,ch_]:=
	If[chival[ch,-1,lev]==(-1)^wt&&Mod[lev,Length[ch]]==0,
	A1[n,wt,lev,ch]+A2[n,wt,lev,ch]+
	A3[n,wt,lev,ch]+A4[n,wt,lev,ch],0]

Trunc[h_]:=h-Mod[h,1]

H[d_]:=Block[{a,b,c,det,hn,ww, ev=0, count=0, t={}},
                ev = Trunc[Sqrt[Abs[d]/3]];
                count=0;
                Do [If [Mod[b^2-d,4]==0,
                        det=(b^2-d)/4;
                        If[b>0, If[Mod[det,b]==0, 
                                t = Union[t,{{b, b, det/b}}]]];
                                Do [If [Mod[det, a]==0,
                                        t=Union[t,{{a, b, det/a}}];
                                        If[b>0,If[a!=det/a,
                                           t=Union[t,{{a, -b, det/a}}]]]],
                                           {a, b+1, Sqrt[det]} ]]
                        ,{b, 0, ev}];
ww=Sum[hn=t[[r]];If[hn[[1]]==hn[[3]],
	If[hn[[1]]==hn[[2]],1/3,If[hn[[2]]==0,1/2,1]],
		1],{r,1,Length[t]}];
                {t, ww}]

HH[d_]:=HH[d]=H[d][[2]]

QF[d_] := Block[{a, b, c, det, ev = 0, count = 0, t = {}}, 
   ev = Trunc[Sqrt[Abs[d]/3]]; count = 0; 
    Do[If[Mod[b^2 - d, 4] == 0, 
      det = (b^2 - d)/4; If[b > 0, 
        If[Mod[det, b] == 0, t = Union[t, {{b, b, det/b}}]]]; 
       Do[If[Mod[det, a] == 0, 
         t = Union[t, {{a, b, det/a}}]; 
          If[b > 0, If[a != det/a, t = Union[t, {{a, -b, det/a}}]]]], 
        {a, b + 1, Sqrt[det]}]], {b, 0, ev}]; {t, Length[t]}]

sqdiv[m_]:=Block[{div=Sqrt[Divisors[m]]},
	Union[Table[If[IntegerQ[div[[var]]]==True,div[[var]],1],
		{var,1,Length[div]}]]]

hc[d_]:=hh[d]=Block[{df=sqdiv[d]},
     Sum[MoebiusMu[df[[var]]]*HH[d/df[[var]]^2],{var,1,Length[df]}]]


dim[wt_,lev_,ch_]:=Block[{aa=Divisors[lev],nchi=Length[ch]},
	If[chival[ch,-1,lev]==(-1)^wt&&Mod[lev,Length[ch]]==0,
		(wt-1)/12*psi[lev]-((wt-1)/3-Floor[wt/3])*
	Sum[If[Mod[x^2+x+1,lev]==0,chival[ch,x,lev],0],{x,1,lev}]-
	((wt-1)/4-Floor[wt/4])*
	Sum[If[Mod[x^2+1,lev]==0,chival[ch,x,lev],0],{x,1,lev}]-
	1/2*Sum[If[Mod[lev/nchi,GCD[aa[[m]],lev/aa[[m]]]]==0,
	    EulerPhi[GCD[aa[[m]],lev/aa[[m]]]],0],{m,1,Length[aa]}]+
	If[wt==2&&nchi==1,1,0],0]]

tracenew[n_,wt_,lev_,ch_]:=Block[{dd=Divisors[lev]},
    Sum[mu2[lev/dd[[j]]]*trace[n,wt,dd[[j]],ch],
		{j,1,Length[dd]}]]

dimnew[wt_,lev_,ch_]:=Block[{dd=Divisors[lev]},
    Sum[mu2[lev/dd[[j]]]*dim[wt,dd[[j]],ch],
		{j,1,Length[dd]}]]

mu2[r_]:=Block[{dd=Divisors[r]},
    Sum[MoebiusMu[dd[[j]]]*MoebiusMu[r/dd[[j]]],{j,1,Length[dd]}]]

leg[p_]:=Table[JacobiSymbol[n,p],{n,1,p}]

pro[d_] := Product[x - a[n], {n, 1, d}]

pow[j_, d_] := Sum[a[n]^j, {n, 1, d}]

rep[n_] := rep[n] = 
   Table[2*Integrate[Sin[x]*Sin[(j + 1)*x]*2^n*Cos[x]^n,
 {x, 0, Pi}]/Pi, 
   {j, 0, n}]

U[m_]:=ChebyshevU[m,x]

Partitions[n_Integer] := Partitions[n,n]

Partitions[n_Integer,_] := {} /; (n<0)

Partitions[0,_] := { {} }

Partitions[n_Integer,1] := { Table[1,{n}] }

Partitions[_,0] := {}

Partitions[n_Integer,maxpart_Integer] :=
	Join[
		Map[(Prepend[#,maxpart])&,
		 Partitions[n-maxpart,maxpart]],
		Partitions[n,maxpart-1]
	]

weight[part_] := 
  (-1)^Length[part]*Product[part[[nn]], {nn, 1, Length[part]}]*
   Product[Count[part, nn]!, {nn, 1, Max[part]}]

t3[n_] := Block[{ppp},Table[ppp = Partitions[m]; 
    s[m] -> Sum[(-1)^m*Product[q[ppp[[j,l]]],
	 {l, 1, Length[ppp[[j]]]}]/
       weight[ppp[[j]]], {j, 1, Length[ppp]}], {m, 1, n}]]

hecke[p_,wt_,lev_,ch_] :=If[GCD[p,lev]==1,
	hhecke[p,wt,lev,ch] /. sub[p,wt,lev,ch],0]

heckef[p_,wt_,lev_,ch_]:=Factor[hecke[p,wt,lev,ch]]

heckenew[p_,wt_,lev_,ch_] :=If[GCD[p,lev]==1,
	hhecken[p,wt,lev,ch] /.
		 subn[p,wt,lev,ch],0]

heckenewf[p_,wt_,lev_,ch_]:=Factor[heckenew[p,wt,lev,ch]]

hhecke[p_,wt_,lev_,ch_] := Block[{deg = dim[wt,lev,ch],
		he,hec,heck},
    he = Sum[(-1)^j*s[j]*x^(deg - j), {j, 0, deg}];
	 hec = he /. t3[deg]; 
	   heck = hec /. sub1[p,wt,lev,ch];
	     Expand[heck] /. s[0] -> 1]

hhecken[p_,wt_,lev_,ch_] := Block[{deg = dimnew[wt,lev,ch],
		he,hec,heck},
    he = Sum[(-1)^j*s[j]*x^(deg - j), {j, 0, deg}];
	 hec = he /. t3[deg]; 
    heck = hec /. sub1n[p,wt,lev,ch];
      Expand[heck] /. s[0] -> 1]

sub1[p_,wt_,lev_,ch_]:=Block[{deg=dim[wt,lev,ch]},
	Table[q[j] -> Sum[rep[j][[i + 1]]*u[i]*
	(chival[ch,p,lev]*p^(wt-1))^((j-i)/2),
	 {i, 0, j}], {j, 1, deg}]]

sub1n[p_,wt_,lev_,ch_]:=Block[{deg=dimnew[wt,lev,ch]},
	Table[q[j] -> Sum[rep[j][[i + 1]]*u[i]*
	(chival[ch,p,lev]*p^(wt-1))^((j-i)/2),
	 {i, 0, j}], {j, 1, deg}]] 

sub[p_,wt_,lev_,ch_] := Block[{deg=dim[wt,lev,ch]},
		 Union[{u[0]->deg},
	Table[u[j] -> trace[p^j,wt,lev,ch], {j, 1, deg}]]]

subn[p_,wt_,lev_,ch_] := Block[{deg=dimnew[wt,lev,ch]},
		 Union[{u[0]->deg},
	Table[u[j] -> tracenew[p^j,wt,lev,ch], {j, 1, deg}]]]


sigma0[n_,ch_]:=Block[{aa=Divisors[n]},Sum[chival[ch,aa[[m]]],
	{m,1,Length[aa]}]]


lim=30

eis[ch_List]:=If[chival[ch,-1]==-1,1-2*Length[ch]/Sum[a*chival[ch,a],
			{a,1,Length[ch]}]*
		Sum[sigma0[n,ch]*q^n,{n,1,lim}]+O[q]^(lim+1),0]


eis[ch_List,lev_]:=If[chival[ch,-1]==-1,
	1-2*Length[ch]/Sum[a*chival[ch,a],{a,1,Length[ch]}]*
		Sum[sigma0[n,ch]*q^(n*lev),{n,1,lim/lev}]+
		O[q]^(lim+1),0]

eis[k_Integer]:=1-2*k/BernoulliB[k]*Sum[DivisorSigma[k-1,n]*q^n,
	{n,1,lim}]+O[q]^(lim+1)

eis[k_Integer,lev_]:=1-2*k/BernoulliB[k]*
	Sum[DivisorSigma[k-1,n]*q^(n*lev),{n,1,lim/lev}]+
		O[q]^(lim+1)



C3x10xchi2:=Block[{chi2={1,I,-I,-1,0},chi3={1,-I,I,-1,0},a,b,c,d},
	a=eis[chi2];b=eis[chi3];c=eis[chi2,2];d=eis[chi3,2];
	(2-I)/20*a*a*b+(-1+3*I)/8*a*a*d+(3-21*I)/40*b*b*d+
	(3+I)/20*a*b*c+(-4+3*I)/20*b*b*b]

C5x10xchi2x1:=Block[{chi2={1,I,-I,-1,0},chi3={1,-I,I,-1,0},a,b,c,d},
	a=eis[chi2];b=eis[chi3];c=eis[chi2,2];d=eis[chi3,2];
	(-1-7*I)/20*a^3*b^2+(1+3*I)/10*a*b^4+(3-6*I)/16*b^4*c+
	(-21+28*I)/80*a^2*b^2*c+(4+3*I)/16*a^3*b*d+
	(-18-9*I)/80*a*b^3*d]

C5x10xchi2x2:=Block[{chi2={1,I,-I,-1,0},chi3={1,-I,I,-1,0},a,b,c,d},
	a=eis[chi2];b=eis[chi3];c=eis[chi2,2];d=eis[chi3,2];
	(-13+9*I)/40*a^5+
	(1-I)/2*a^3*b^2+(-9+13*I)/40*a*b^4+
	(33+6*I)/40*a^4*c+
	(11-2*I)/16*b^4*c+
	(-23/16)*a^2*b^2*c+I/16*a^3*b*d+
	(-2-11*I)/80*a*b^3*d]

ell14:=Block[{ro=Exp[Pi*I/3],chi3,chi5,a,b,A,B},
	chi3={1,ro^2,ro,-ro,-ro^2,-1,0};
	chi5={1,-ro,-ro^2,ro^2,ro,-1,0};
		a=eis[chi3];A=eis[chi3,2];
		b=eis[chi5];B=eis[chi5,2];
	Simplify[1/7*a*b+2/7*A*B+
		(-3-3*I*Sqrt[3])/14*A*b+(-3+3*I*Sqrt[3])/14*a*B]]