\p 50 \\ Global variable D: a prime which is 3 mod 4, so that -D \\ is the discriminant of an imaginary quadratic field K. D = 7; \\ The elliptic curve E is the CM curve of conductor \\ D^2, with CM by K. Here D=7. E = ellinit([1,-1,0,-107,552]); \\ 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, q, a, b, nrm, cf ) = \ res = vector(M^2);\ a = 1;\ while(a<=2*M,\ if((a%2==1),b=1,b=2);\ while((a^2+D*b^2)<=4*M^2,\ nrm = (a^2+D*b^2)/4;\ cf = (a^2-D*b^2);\ res[nrm] = res[nrm] + cf;\ b = b+2);\ a = a+1);\ for(a=1,M, res[a^2] = res[a^2] + 2*a^2);\ for(b=1,floor(M/sqrt(D)), res[D*b^2] = res[D*b^2] - 2*D*b^2);\ return(res); f = teta(7,100)/2; \\ 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); \\ prim0(f,tau): Returns the integral from (i.infinity) to tau \\ of f(z) dz, for f a modular form of weight 3. \\ prim0(f,tau,\ res,M,n,q) = \ res=0;\ q = exp(2*Pi*I*tau);\ M = length(f);\ for(n=1,M, res=res+f[n]/n*q^n);\ res=res/(2*Pi*I);\ return(res); \\ prim1(f,tau): Returns the integral from (i.infinity) to tau \\ of z.f(z)dz, for f a modular form of weight 3. \\ prim1(f,tau,\ res,M,n,q) = \ res=0;\ q = exp(2*Pi*I*tau);\ M = length(f);\ for(n=1,M, res=res-f[n]/n^2*q^n);\ res=res/(2*Pi*I)^2;\ res=res+tau*prim0(f,tau);\ return(res); \\ per0(f,g): Returns the period integral of f(z)dz from \\ (i.infty) to g(i.infty), \\ Here g is taken to be a matrix in Gamma_0(D) with a(g) a square \\ modulo D. \\ per0(f,g,\ a,b,c,d,res) = \ a = g[1,1]; b=g[1,2];\ c = g[2,1]; d=g[2,2];\ res = prim0(f,a/c+ I/c);\ res = res - c*prim1(f,-d/c+I/c);\ res = res - d*prim0(f,-d/c+I/c);\ return(4*Pi^2*res); \\ per1(f,g): Returns the period integral of z.f(z)dz from \\ (i.infty) to g(i.infty), \\ Here g is taken to be a matrix in Gamma_0(D) with a(g) a square \\ modulo D. \\ per1(f,g,\ a,b,c,d,res) = \ a = g[1,1]; b=g[1,2];\ c = g[2,1]; d=g[2,2];\ res = prim1(f,a/c+ I/c);\ res = res - a*prim1(f,-d/c+I/c);\ res = res - b*prim0(f,-d/c+I/c);\ return(4*Pi^2*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,res) = \ tau = omega1/omega2;\ vol = (omega1*conj(omega2)-conj(omega1)*omega2);\ res =-(4*Pi^2)*vol^(-1)*(conj(omega2)*prim1(f,tau) - conj(omega1)*prim0(f,tau));\ return(res); \\ [omega1, omega2] generates the Neron lattice of E and is an admissible basis omega1 = Omegap;\ omega2 = (D - sqrt(-D))/2*Omegap;\ \\ Compute the Stark-Heegner point J = shp( omega1, omega2 ); \\ J = J/sqrt(-D); print("Stark-Heegner point of conductor 1"); lindep([J,omega1,omega2])