/* 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)
*/