/*** EXAMPLE: L-function of an elliptic curve over a quadratic field ***/

\\ This computes the L-function of the elliptic curve over
\\ K = Q(sqrt(37)) with everywhere good reduction given by
\\ E : y^2 + y = x^3 + 2x^2 - (19+8omeg)x + (28+11omeg)
\\ where omeg=(1+sqrt(37))/2.

read("computel");           \\ read the ComputeL package
                            \\ and set the default values
default(realprecision,28);  \\ set working precision; used throughout

\\ initialize L-function parameters

conductor=37^2;    \\ conductor for the exponential factor (not for E)
gammaV=[0,0,1,1];  \\ list of gamma-factors
weight=2;          \\ Functional equation is
sgn=1;             \\ Lstar(E,s)=sgn*Lstar(E,weight-s)

\\ Number of coefficients needed

print("Precision ",default(realprecision));
print("Need ",B0=cflength()," coefficients");

\\ Computation of the coefficients a_n of L(E/K,s) = sum a_n/n^s

\\ Square test in finite field F_p^2, p odd

{is_square_Fq(a,p)=
   \\ a=Mod(A,mu) where
   \\ mu irreducible degree 2 polynomial over F_p
   \\ A polynomial of degree <=1 over F_p
local(mu,b,B,bc,r);
if(a==0,return(0));
mu=component(a,1);
b=a^((p-1)/2);
B=component(b,2);
bc=Mod(subst(B,x,1-x),mu);
r=b*bc; \\ computes a^((p^2-1)/2)
if(r==1,1,-1);
}

\\ Computes a_p for prime p inert in K
\\ Remark : rustic algorithm, not te be used when p is big!

{ellap_inert(p)=
local(Nb,mu,t,X,r);
Nb=0;
mu=Mod(1,p)*x^2-Mod(1,p)*x-Mod(9,p); \\ minimal polynomial of omeg over F_p
t=Mod(Mod(1,p)*x,mu);
for(a=0,p-1,
for(b=0,p-1,
X=a+b*t;
r=X^3+2*X^2-(19+8*t)*X+Mod(28,p)+11*t;
if(p==2,Nb+=2*(r==0|r==1),
Nb+=1+is_square_Fq(1+4*r,p))));
p^2-Nb}

\\ Computes a_p for prime p split in K

{ellap_split(p)=
local(t,e1,e2);
t=polrootsmod(x^2-x-9,p);
e1=ellinit([Mod(0,p),Mod(2,p),Mod(1,p),-19-8*t[1],28+11*t[1]]);
e2=ellinit([Mod(0,p),Mod(2,p),Mod(1,p),-19-8*t[2],28+11*t[2]]);
[ellap(e1,p),ellap(e2,p)]
}

\\ Computes a_n up to some bound B

{coeffs(B)=
local(a,m,Lp,aP,f);
a=vector(B,i,0);
a[1]=1;
forprime(p=2,B,
m=floor(log(B)/log(p));
if(p==37,a[p]=-2;                 \\ p ramified in K
Lp=1/(1-a[p]*x+p*x^2)+O(x^(m+1)); \\ Euler factor at p
for(i=2,m,a[p^i]=polcoeff(Lp,i)));
if(kronecker(37,p)==1,aP=ellap_split(p);               \\ p split in K
Lp=1/((1-aP[1]*x+p*x^2)*(1-aP[2]*x+p*x^2))+O(x^(m+1)); \\ Euler factor at p
for(i=1,m,a[p^i]=polcoeff(Lp,i)));
if(kronecker(37,p)==-1,
if(m>=2,a[p^2]=ellap_inert(p);          \\ p inert in K
Lp=1/(1-a[p^2]*x^2+p^2*x^4)+O(x^(m+1)); \\ Euler factor at p
for(i=2,floor(m/2),a[p^(2*i)]=polcoeff(Lp,2*i)))));
for(n=2,B,
f=factor(n);
if(length(f~)>1,
a[n]=1;
for(i=1,length(f~),a[n]*=a[f[i,1]^f[i,2]])));
return(a)
}

\\ Computation of L(E/K,s)

coeff_vec=coeffs(B0);
initLdata("coeff_vec[k]"); \\ initialize L-function
print("Check functional equation : ",checkfeq());
print("L(E/Q(sqrt(37)),1) = ",L(1));