\\ Global variable D: a prime which is 3 mod 4, so that -D
\\ is the discriminant of an imaginary quadratic field K.
D = 43;
w = quadgen(-D);
r =3
\\ The elliptic curve E is the CM curve of conductor
\\ D^2, with CM by K. Here D=11.
E = ellinit([0,0,1,-860,9707]);
\\ These are the periods of the elliptic curve
Omegap = E[15];
Omegam = E[16];
\\ The function teta(D,M) returns the theta series of weight 3
\\ attached to the discriminant -D. This is a form on
\\ Gamma_1(D), with odd character.
\\ It is stored as a vector v of length M^2,
\\ where v[n] is the nth
\\ fourier coefficient of the theta series.
\\
teta(D,M,\
res, u, q, a, b, nrm, cf ) = \
res = vector(M^2);\
for(a = -M,M,\
for(b = -M,M,\
u = a+ b*w;\
nrm = norm(u);\
if((nrm< M^2)& (nrm>0), res[nrm] = res[nrm]+ u^(r+1),);\
));\
return(res);
\\ The function evl(f,tau) evaluates the form f at the argument
\\ tau in the upper half-plane.
\\ It was mainly used as a debugging tool to check that the theta
\\ series computed by teta satisfies the required invariance
\\ properties.
evl(f,tau,\
q, res, M, j) = \
q = exp(2*Pi*I*tau);\
res=0;\
M = length(f);\
for(j=1,M,res=res+f[j]*q^j);\
return(res);
\\ prim(f,tau): Returns the vector whose j-th entry (0<= j <= r)
\\ is the integrals from (i.infinity) to tau
\\ of (z - tau-bar)^j f(z) dz, for f a modular form of weight r+2.
\\
prim(f,tau,\
res,M,n,q,cf) = \
res= vector(r+1);\
M = length(f);\
for(n=1,M, \
q = exp(2*Pi*I*n*tau);\
cf = q / (2*Pi*I*n);\
for(j=0,r,\
res[j+1] = res[j+1] + f[n]*cf;\
cf = (tau-conj(tau))^(j+1)* q/ (2*Pi*I*n) - (j+1)/(2*Pi*I*n)*cf ) );\
return(res);
\\ The function shp returns the Stark-Heegner point
\\ asssociated to the lattice generated by
\\ omega1 and omega2. This lattice is assumed to be
\\ positively oriented.
\\ We expect the resulting invariant to be a global
\\ point if [omega1,omega2] is the Neron lattice
\\ of an elliptic curve with CM by O_K.
shp(omega1,omega2,\
tau, vol,prims,res) = \
tau = omega1/omega2;\
prims = prim(f,tau);\
vol = (omega1*conj(omega2)-conj(omega1)*omega2);\
res = (2*Pi*I)^(r+1)* vol^(-r)* conj(omega2)^r * prims[r+1] ;\
return(res);
\\ Calculate the theta series.
f = teta(D,150)/2;
\\ [omega1, omega2] is an admissible basis for the Neron lattice of E.
omega1 = Omegam;\
omega2 = Omegap;\
\\ Compute the Stark-Heegner points
r =3
ell =2
n = 1
s = 22
epsl=-1
J = vector(ell+1)
P = vector(ell+1)
for(a=0, ell-1,\
J[a+1] = shp( (omega1 + a*omega2)/ell, omega2 ); \
JJ=J[a+1];\
for(s=0,42,\
J[a+1] = (JJ- s*omega1)/sqrt(-D);\
J[a+1] = n*J[a+1];\
P[a+1] = ellztopoint(E,J[a+1]);\
print("j = ", a);\
print("s= ",s);\
print(algdep(P[a+1][1],3));\
print("---------------------------");\
) )
J[ell+1] = epsl*shp( omega1, omega2/ell );
J[ell+1] = (J[ell+1]- epsl * ell^r *s*omega1)/sqrt(-D);
J[ell+1] = n*J[ell+1];
P[ell+1] = ellztopoint(E,J[ell+1]);
print("j = ", infty);
print(P[ell+1]);
print(algdep(P[ell+1][1],3));\
print("---------------------------");