bnf = bnfinit(a^2+23);
bnr = bnrinit(bnf, 1);
bnr.clgp
Hecke = lfuncreate([bnr,[1]]); Hecke[2..5]

lfun(Hecke,0)
z=lfun(Hecke,0,1)
P=algdep(exp(z),3)
nfdisc(P)

mf = mfinit([23,1,-23],0);
S=mfeigenbasis(mf); #S
round(mfcoefs(S[1],100)[^1]-lfunan(Hecke,100),&e)
e

E=ellinit("40a1");
P=elldivpol(E,3)
Q=polresultant(P,y^2-elldivpol(E,2));
polgalois(Q)
R=nfsplitting(Q)

N=nfinit(R); G=galoisinit(N);
[T,o]=galoischartable(G); T~
o

minpoly(Mod(y^3+y, polcyclo(o,y)))

L = lfunartin(N,G,T[,4],o);
L[2..5]
lfuncheckfeq(L)

dT = galoischardet(G,T[,4],o)
dL = lfunartin(N,G,dT,o); dL[2..5]

S = lfunan(L,1000); SE = lfunan(E,1000);
Smod3 = round(real(S))+round(imag(S)/sqrt(2));
[(Smod3[i]-SE[i])%3|i<-[1..#Smod3],gcd(i,40)==1]

mfE=mfinit([40,2,1],0);
eiE = mfeigenbasis(mfE); #eiE
mfcoefs(eiE[1],100)[^1]==ellan(E,100)

mfL=mfinit([1800,1,-3],0);
eiL = mfeigenbasis(mfL); #eiL
round(subst(lift(mfcoefs(eiL[2],100)),y,sqrt(-2))[^1]-lfunan(L,100),&e)
e

mfgaloistype(mfL,eiL[2])
PR=polredabs(mfgaloisprojrep(mfL,eiL[2]))
PR==polredabs(nfsplitting(elldivpol(E,3)))