/* Some simple utility functions for elliptic curves in GP Author: John Cremona (john.cremona@nottingham.ac.uk) Last edited: 2007-02-21 */ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ Negation of a point: \\ e is an ellinit, p a point on e; returns -p \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ ellneg(e,p)=ellsub(e,[0],p) /* Example: e=ellinit([0,0,1,-7,6]); p=[-2,3]; ellneg(e,p) returns [-2, -4] */ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ Member function to extract [a1,a2,a3,a4,a6] from a longer vector: \\ e is an ellinit; returns [a1,a2,a3,a4,a6] \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ e.ai=vecextract(e,[1,2,3,4,5]) \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ Inverse of a standard transformation [u,r,s,t]: \\ urst=[u,r,s,t] encoding a standard transformation of Weierstrass \\ equations (so u!=0); returns the inverse trabsformation \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ { ellinvert(urst)= local(u=urst[1],r=urst[2],s=urst[3],t=urst[4]); [1/u,-r/u^2,-s/u,(r*s-t)/u^3] } /* Example: e1=ellinit([1,2,3,4,6]); v=[u,r,s,t] \\generic e2=ellchangecurve(e1,v); e1==ellchangecurve(e2,ellinvert(v)) \\ =1 (true) P1=[-1,1]; ellisoncurve(e1,P1) \\ =1 (true) P2=ellchangepoint(P1,v) P1==ellchangepoint(P2,ellinvert(v)) \\ =1 (true) */ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ Composition of two standard transformations: \\ urst1, urts2 are two standard transformations; returns their composite \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ { ellcompose(urst1,urst2)= local(u1=urst1[1],r1=urst1[2],s1=urst1[3],t1=urst1[4], u2=urst2[1],r2=urst2[2],s2=urst2[3],t2=urst2[4]); [u1*u2,u1^2*r2+r1,u1*s2+s1,u1^3*t2+s1*u1^2*r2+t1] } /* Example: ellcompose([u,r,s,t],ellinvert([u,r,s,t])) == [1,0,0,0] e1=ellinit([1,2,3,4,6]); v1=[u1,r1,s1,t1] \\generic 1 v2=[u2,r2,s2,t2] \\generic 2 e2=ellchangecurve(e1,v1); e3=ellchangecurve(e2,v2); e3==ellchangecurve(e1,ellcompose(v1,v2)) */ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ \\ Roots of a polynomial! (not elliptic curve specific) \\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ { polratroots(pol)= fx=factor(pol); ans=[]; for(j=1,#(fx~),if(poldegree(fxj=fx[j,1])==1, ans=concat(ans,[-polcoeff(fxj,0)/polcoeff(fxj,1)]))); ans; } \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ Test of isomorphism between to curves e1,e2: either returns urst \\ such that ellchangecurve(e1,urst)==e2, or 0 \\ Requires field of definition such that issquare(), sqrt() and \\ (for the case j=0 only) factorization of x^3-a to get a cube root \\ \\ Do not expect this to work in characteristics 2,3!! \\ \\ We don't use the form issquare(u2,&u) since it does not work (for \\ example) over Z/pZ... \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ { elliso(e1,e2)= local(u6,u3,uu,u,u4,u2,r,s,t); \\ Check we have at least 13-vectors: if(#e1<13,e1=ellinit(e1,1)); if(#e2<13,e2=ellinit(e2,1)); \\ Check j-invariants equal: if(!(e1.j==e2.j),return(0)); \\ Special case where j==c4==0: if(e1.j==0,u6=e1.c6/e2.c6; if(!issquare(u6),return(0)); u3=polratroots(x^2-u6)[1]; uu=polratroots(x^3-u3); if(#uu==0,return(0),u=uu[1]), \\ Special case where j-1728==c6==0: if(e1.j==1728,u4=e1.c4/e2.c4; if(!issquare(u4),return(0)); u2=polratroots(x^2-u4)[1]; if(!issquare(u2),return(0)); u=polratroots(x^2-u2)[1], \\ Generic case, we have u^2: u2=(e1.c6*e2.c4)/(e2.c6*e1.c4); if(!issquare(u2),return(0)); u=polratroots(x^2-u2)[1])); \\Given u, solve for r,s,t: s=(u*e2.a1-e1.a1)/2; r=(u^2*e2.a2-e1.a2+s*e1.a1+s^2)/3; t=(u^3*e2.a3-e1.a3-r*e1.a1)/2; [u,r,s,t]; } compsq(e) = if((e.a1==0)&&(e.a3==0),e,ellchangecurve(e,[1/2,0,-e.a1/2,-e.a3/2])); \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ Linear combinations of points \\ e is an ellinit, pts a vector of points, v a vector of integers \\ with #v=#pts; returns sum(v[j]*pts[j]) on E \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ (May become obsolete when factorback() is adapted) { ellcomb(E,pts,v)= local(P); P=[0]; if(#pts==#v, for(i=1,#pts, if(v[i], P=elladd(E,P,ellpow(E,pts[i],v[i])))), print("Error in ellcomb: arguments 2 and 3 must have the same length!")); P } /* Example: e=ellinit([0,0,1,-7,6]); pts=[[-2,3],[-1,3],[0,2]]; ellcomb(e,pts,[-1,2,3]) returns [13961/100, -1649791/1000] */ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ The primes of bad reduction for an integral elliptic curve over Q \\ e is an ellinit; returns a vector of primes dividing the \\ discriminant (so these are bad reduction for this model but may not be \\ bad for the underlying elliptic curve if this is not a minimal model) \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ ellbadprimes(e)=factor(abs(e[12]))[,1]~ /* Example: ellbadprimes(ellinit([1,0,1,1,2])) returns [2,3,5] */ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ \\ Construct a vector of all points on e from a given vector of \\ (possible) x-coordinates); at most one from each x unless flag!=0 \\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ { ellpointsfromx(e,xi,flag=0)= local(yij,ans); yij=ellordinate(e,xi);maxny=if(flag,2,1); ans=vector(#xi,i,vector(min(#yij[i],maxny),j,[xi[i],yij[i][j]])); if(#ans,concat(ans),[]) } /* Example: e=ellinit([0,0,1,-7,6]); ellpointsfromx(e,vector(1004,j,j-4)) returns 18 integral points on e with x-coordinates between -3 and +1000 (these are all the integral points on e, up to sign, but that's another story! */ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\ \\ 2- and 3-torsion points (see ell_ff.gp for general case) \\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ { ell2torsion(e)= ellpointsfromx(e,polratroots(4*(x^3+e.a2*x^2+e.a4*x+e.a6)+(e.a1*x+e.a3)^2)); } { ell3torsion(e)= ellpointsfromx(e,polratroots(3*x^4+(e.a1^2+4*e.a2)*x^3+(3*e.a1*e.a3+6*e.a4)*x^2+(3*e.a3^2+12*e.a6)*x+e.a1^2*e.a6-e.a1*e.a3*e.a4+e.a2*e.a3^2+4*e.a2*e.a6-e.a4^2)); } /* Examples: e=ellinit([0,1,0,4,4]); ell3torsion(e) ell2torsion(e) returns [[0, 2], [0, -2]] [[-1, 0]] BUT e=ellinit(1.0*[0,0,0,-2,0]); ell2torsion(e) returns 4 points, since ellordinate(e,sqrt(2)) returns two values, both approximately zero! [Note added 2007-02-21: with Pari 2.4.1 it works correctly] e=ellinit(Mod(1,47)*[0,0,0,-2,0]); lift(ell2torsion(ellinit(Mod(1,47)*[0,0,0,-2,0]))) lift(ell3torsion(ellinit(Mod(1,47)*[0,0,0,-2,0]))) returns [[0, 0], [40, 0], [7, 0]] [[21, 30], [21, 17]] (we "lift"ed the results for easier readability) */