\p 70
r  = 1

\\ Global variable D: a prime which is 3 mod 4, so that -D
\\ is the discriminant of an imaginary quadratic field K.
D = 19;
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,-38,90]);

\\ 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 = 1
ell =11
n = 1 
s = 1 
epsl=1
J = vector(ell+1)
P = vector(ell+1)
for(a=0, ell-1,\
 J[a+1] = shp( (omega1 + a*omega2)/ell, omega2 ); \
 J[a+1] = (J[a+1]- s*omega1)/sqrt(-D);\
 J[a+1] = n*J[a+1];\
 P[a+1] = ellztopoint(E,J[a+1]);\
 print("j = ", a);\
 print(P[a+1]);\
 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("---------------------------");