\\ Global variable D: a prime which is 3 mod 4, so that -D \\ is the discriminant of an imaginary quadratic field K. D = 11 \\ The elliptic curve E is the CM curve of conductor \\ D^2, with CM by K. Here D=11. E = ellinit([0,-1,1,-7,10]) \\ 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(11,100); \\ 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(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(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 = vol^(-1)*(conj(omega2)*prim1(f,tau) - conj(omega1)*prim0(f,tau));\ return(4*Pi^2*res); \\ glob1 is the generator of E(Q) glob1 = ellpointtoz(E,[4,5]); \\ glob2 is the other generator of E(K) (or perhaps not...) glob2 = (1+sqrt(-D))/2*glob1; \\ [omega1, omega2] generate the Neron lattice of E. omega1 = Omegap*(D+sqrt(-D))/(2*D);\ omega2 = Omegap;\ \\ Compute the Stark-Heegner point associated to... print("Point 0"); J0 = 3*shp(omega1, 3*omega2 ); J0 = 3*J0; J0 = (J0+2*Omegam)/sqrt(-11); lindep([J0,glob1,glob2,Omegap,Omegam]) print("Point 1"); J1 = 3*shp(omega1+omega2, 3*omega2 ); J1 = 3*J1; J1 = (J1+2*Omegam)/sqrt(-11); lindep([J1,glob1,glob2,Omegap,Omegam]) print("Point 2"); J2 = 3*shp(omega1+2*omega2, 3*omega2 ); J2 = 3*J2; J2 = (J2+2*Omegam)/sqrt(-11); lindep([J2,glob1,glob2,Omegap,Omegam]) print("Point 3"); J3 = 3*shp(3*omega1, omega2 ); J3 = 3*J3; J3 = (J3+6*Omegam)/sqrt(-11); lindep([J3,glob1,glob2,Omegap,Omegam])