\\ Global variable D: a prime which is 3 mod 4, so that -D
\\ is the discriminant of an imaginary quadratic field K.
D = 67;
w = quadgen(-D);
\\ The elliptic curve E is the CM curve of conductor
\\ D^2, with CM by K. Here D=11.
E = ellinit([0,0,1,-7370,243528]);
\\ 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 point
print("Stark-Heegner point of conductor 1");
J = shp(omega1, omega2);
z1 = ellpointtoz(E,[201/4,-71/8]);
z2 = (1+sqrt(-D))/2* z1;
cc = lindep([J, z1,z2, omega1,omega2])
print(cc);
print("Modulo torsion, the Chow-Heegner point is this mutiple of the generator (201/4,-71/8)")
-(-cc[2]-cc[3]*w )/cc[1]/(1-2*w)