\p 100 \\ Global variable D: a prime which is 3 mod 4, so that -D \\ is the discriminant of an imaginary quadratic field K. D = 11; w = quadgen(-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]; gen = 2*Omegam+Omegap; \\ 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,50); \\ 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); \\ The function Lindep(p,x,y) finds a linear dependence \\ relation between p,x,y over O_K Lindep(p,x,y,\ x1,x2,y1,y2,res,cx,cy ) = \ x1 = x; x2 = (1+sqrt(-11))/2*x;\ y1 = y; y2 = (1+sqrt(-11))/2*y;\ res = lindep([p,x1,x2,y1,y2]);\ print(res);\ w=quadgen(-11);\ cx = res[2]+w*res[3];\ cy = res[4]+w*res[5];\ return([res[1],cx,cy]) ; \\ The function ratrec performs a rational recognition on \\ the real number x ratrec(x,\ p, res) = \ p= algdep(x,2);\ if(poldegree(p)==2,print("Warning: insufficient precision in ratrec"),);\ res = -polcoeff(p,0)/polcoeff(p,1) ; \\ The function ratrecpol performs a rational recognition on \\ a polynomial p with real coefficients. ratrecpol(f,\ d, res, j ) = \ d = poldegree(f);\ res=0;\ for(j=0,d,\ cf = ratrec(polcoeff(f,j));\ res=res+cf*x^j );\ return(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("computing point 0"); J0 = 3*shp(omega1, 7*omega2 ); J0= J0*7; J0 = (J0+2*Omegam)/sqrt(-11); J0 = (J0+Omegap+Omegam)/3; P0 = ellztopoint(E,J0); print("computing point 1"); J1 = 3*shp(omega1 + omega2, 7*omega2 ); J1 = J1*7; J1 = (J1+2*Omegam)/sqrt(-11); J1 = (J1+Omegap+Omegam)/3; P1 = ellztopoint(E,J1); print("computing point 2"); J2 = 3*shp(omega1 + 2*omega2, 7*omega2 ); J2 = J2*7; J2 = (J2+2*Omegam)/sqrt(-11); J2 = (J2+Omegap+Omegam)/3; P2 = ellztopoint(E,J2); print("computing point 3"); J3 = 3*shp(omega1 + 3*omega2, 7*omega2 ); J3 = J3*7; J3 = (J3+2*Omegam)/sqrt(-11); J3 = (J3+Omegap+Omegam)/3; P3 = ellztopoint(E,J3); print("computing point 4"); J4 = 3*shp(omega1 + 4*omega2, 7*omega2 ); J4 = J4*7; J4 = (J4+2*Omegam)/sqrt(-11); J4 = (J4+Omegap+Omegam)/3; P4 = ellztopoint(E,J4); print("computing point 5"); J5 = 3*shp(omega1 + 5*omega2, 7*omega2 ); J5 = J5*7; J5 = (J5+2*Omegam)/sqrt(-11); J5 = (J5+Omegap+Omegam)/3; P5 = ellztopoint(E,J5); print("computing point 6"); J6 = 3*shp(omega1 + 6*omega2, 7*omega2 ); J6 = J6*7; J6 = (J6+2*Omegam)/sqrt(-11); J6 = (J6+Omegap+Omegam)/3; P6 = ellztopoint(E,J6); print("Computing point Infinity"); J7 = 3* shp(7*omega1, omega2 ); J7 = J7*7; J7 = (J7+14*Omegam)/sqrt(-11); J7 = (J7+Omegap+Omegam)/3; J7=-J7; P7 = ellztopoint(E,J7); print(" "); print("Polynomial satisfied by the x coordinates"); fx = (x-P0[1])*(x-P1[1])*(x-P2[1])*(x-P3[1])*(x-P4[1])*(x-P5[1])*(x-P6[1])*(x-P7[1]) print(" ") print("This polynomial is close to"); fx = ratrecpol(fx) fx = fx*7^6*13^4; print(" ") print(" Factorisation = ") g = factor(fx) g=g[1,1] print(" ") print("Galois group of this polynomial") polgalois(g)