nf  = nfinit(a^2-5);
nf.disc
phi = (1+a)/2;
E   = ellinit([1,phi+1,phi,phi,0],nf);
N   = ellglobalred(E)[1]
tor = elltors(E) \\ Z/8Z
ellrootno(E)

om = E.omega
per = om[1][1]*om[2][1];
tam = elltamagawa(E)
bsd = (per*tam) / (tor[1]^2*sqrt(abs(nf.disc)))
ellbsd(E)
L1 = lfun(E,1)

E = ellinit([0,0,1,-7,6]);
ellidentify(E)
lfunorderzero(E)
ellratpoints(E,10)
G = ellgenerators(E)

R = matdet(ellheightmatrix(E,G))
ellbsd(E)*R
lfun(E,1,3)/3!

L = lfunsympow(E,2); L[2..5]
z = lfun(L,2)
om=E.omega[1]*imag(E.omega[2])
z*2*Pi*5077/(om*(2*Pi)^2)
ellmoddegree(E)

L=lfungenus2([x^2+x,x^3+1]);
L[2..5]
lfun(L,1)
lfunan(L,10)