\\Definition of an elliptic curve with coefficients [a1,a2,a3,a4,a6] in Q
(13:54) gp > E=ellinit([2,5,8,1136,531])
%1 = [2, 5, 8, 1136, 531, 24, 2288, 2188, -1295608, -54336, 1490400, -94122075952, 10026353442816/5882629747, [-0.4792552759371703523735735672, -2.760372362031414823813213216 - 33.67097412292706002213729956*I, -2.760372362031414823813213216 + 33.67097412292706002213729956*I]~, 0.6288018119558452216715109778, 0.3144009059779226108357554889 + 0.3242753052728374333438026774*I, -4.799518857659883600545936426 + 1.022223116 E-28*I, -2.399759428829941800272968213 - 7.471285875763482942024958904*I, 0.2039048995280950282462870983]
\\Discriminant of E
(14:14) gp > E.disc
%2 = -94122075952
\\j-invariant of E
(14:14) gp > E.j
%3 = 10026353442816/5882629747
\\Definition of the same curve mod 2003
(14:14) gp > E=ellinit([2,5,8,1136,531]*Mod(1,2003))
%4 = [Mod(2, 2003), Mod(5, 2003), Mod(8, 2003), Mod(1136, 2003), Mod(531, 2003), Mod(24, 2003), Mod(285, 2003), Mod(185, 2003), Mod(333, 2003), Mod(1748, 2003), Mod(168, 2003), Mod(1707, 2003), Mod(171, 2003), 0, 0, 0, 0, 0, 0]
(14:15) gp > E.disc
%5 = Mod(1707, 2003)
(14:15) gp > E.j
%6 = Mod(171, 2003)
(14:15) gp > E=ellinit([2,5,8,1136,531])
%7 = [2, 5, 8, 1136, 531, 24, 2288, 2188, -1295608, -54336, 1490400, -94122075952, 10026353442816/5882629747, [-0.4792552759371703523735735672, -2.760372362031414823813213216 - 33.67097412292706002213729956*I, -2.760372362031414823813213216 + 33.67097412292706002213729956*I]~, 0.6288018119558452216715109778, 0.3144009059779226108357554889 + 0.3242753052728374333438026774*I, -4.799518857659883600545936426 + 1.022223116 E-28*I, -2.399759428829941800272968213 - 7.471285875763482942024958904*I, 0.2039048995280950282462870983]
(14:15) gp > E.disc%2003
%8 = 1707
(14:15) gp > E.j%2003
%9 = 171
\\Trying to find a point on the curve, no solution for x = 1
(14:15) gp > polrootsmod(y^2+10*y+330,2003)
%10 = []~
\\There are two solutions for x = 1118 
(14:17) gp > polrootsmod(y^2+241*y+1017,2003)
%11 = [Mod(269, 2003), Mod(1493, 2003)]~
\\Discriminant of an elliptic curve in short Weierstrass form
(14:17) gp > E=ellinit([0,0,0,a,b])
%12 = [0, 0, 0, a, b, 0, 2*a, 4*b, -a^2, -48*a, -864*b, -64*a^3 - 432*b^2, -6912*a^3/(-4*a^3 - 27*b^2), 0, 0, 0, 0, 0, 0]
(14:21) gp > E.disc
%13 = -64*a^3 - 432*b^2
\\j-invariant of an elliptic curve in short Weierstrass form
(14:21) gp > E.j
%14 = -6912*a^3/(-4*a^3 - 27*b^2)
\\Isomorphic curve to E
(14:21) gp > Et=[0,0,0,4*a,8*b]
%15 = [0, 0, 0, 4*a, 8*b]
(14:27) gp > Et=ellinit([0,0,0,4*a,8*b])
%16 = [0, 0, 0, 4*a, 8*b, 0, 8*a, 32*b, -16*a^2, -192*a, -6912*b, -4096*a^3 - 27648*b^2, -6912*a^3/(-4*a^3 - 27*b^2), 0, 0, 0, 0, 0, 0]
(14:28) gp > E.j
%17 = -6912*a^3/(-4*a^3 - 27*b^2)
\\This curve has the same j-invariant it is isomorphic to E over \bar{K}. In case 2 is not a square in K, Et is not isomorphic to E
\over K but only over a quadratic extension of K   
(14:28) gp > Et.j
%18 = -6912*a^3/(-4*a^3 - 27*b^2)
\\Scalar multiplication can be done efficiently: take a random p
(15:02) gp > p=nextprime(random(2^200))
%31 = 1163309487432758728036238636255050284363371301768883524891401
\\Take a simple curve over Fp
(15:03) gp > E=ellinit([0,0,0,1,1]*Mod(1,p))
%32 = [Mod(0, 1163309487432758728036238636255050284363371301768883524891401), Mod(0, 1163309487432758728036238636255050284363371301768883524891401), Mod(0, 1163309487432758728036238636255050284363371301768883524891401), Mod(1, 1163309487432758728036238636255050284363371301768883524891401), Mod(1, 1163309487432758728036238636255050284363371301768883524891401), Mod(0, 1163309487432758728036238636255050284363371301768883524891401), Mod(2, 1163309487432758728036238636255050284363371301768883524891401), Mod(4, 1163309487432758728036238636255050284363371301768883524891401), Mod(1163309487432758728036238636255050284363371301768883524891400, 1163309487432758728036238636255050284363371301768883524891401), Mod(1163309487432758728036238636255050284363371301768883524891353, 1163309487432758728036238636255050284363371301768883524891401), Mod(1163309487432758728036238636255050284363371301768883524890537, 1163309487432758728036238636255050284363371301768883524891401), Mod(1163309487432758728036238636255050284363371301768883524890905, 1163309487432758728036238636255050284363371301768883524891401), Mod(412787237476140193819310483832437197677325300627668347542333, 1163309487432758728036238636255050284363371301768883524891401), 0, 0, 0, 0, 0, 0]
\\Find a point on E
(15:03) gp > ellordinate(E,9)
%41 = [Mod(852019353894870034995234465487017253976606804911790422115998, 1163309487432758728036238636255050284363371301768883524891401), Mod(311290133537888693041004170768033030386764496857093102775403, 1163309487432758728036238636255050284363371301768883524891401)]
(15:03) gp > P=[9,lift(%41[1])]
%42 = [9, 852019353894870034995234465487017253976606804911790422115998]
\\take a random scalar n
(15:04) gp > n=random(p)
%43 = 615685866273440222446493529976616807081324908238167371900531
\\compute [n]P
(15:04) gp > lift(ellpow(E,P,n))
%45 = [719813051721822900598504827836936046952829448871253428653792, 316638185574984655475300509640285363402140335232930391569668]
\\procedure to compute ell-division polynomial can be found at http://www.warwick.ac.uk/~masgaj/ftp/progs/pari/index.html
(15:11) gp > global(x,y);
elldivpol(e,n)=local(ans=elldivpol0(e,n));if(n%2==1,ans,ans*(2*y+e[1]*x+e[3]))
elldivpol0(e,n)=local(m,a1,a2,a3,a4,a6,t1,t2,f1,f2,psi24);a1=e[1];a2=e[2];a3=e[3];a4=e[4];a6=e[5];f1=x^3+a2*x^2+a4*x+a6;f2=a1*x+a3;n=abs(n);if(n==0,return(0));if(n==1,return(1));if(n==2,return(1));if(n==3,return(3*x^4+(a1^2+4*a2)*x^3+(3*a1*a3+6*a4)*x^2+(3*a3^2+12*a6)*x+a1^2*a6-a1*a3*a4+a2*a3^2+4*a2*a6-a4^2));if(n==4,return(2*x^6+(a1^2+4*a2)*x^5+(5*a1*a3+10*a4)*x^4+(10*a3^2+40*a6)*x^3+(10*a1^2*a6-10*a1*a3*a4+10*a2*a3^2+40*a2*a6-10*a4^2)*x^2+(a1^4*a6-a1^3*a3*a4+a1^2*a2*a3^2+8*a1^2*a2*a6-a1^2*a4^2-4*a1*a2*a3*a4-a1*a3^3-4*a1*a3*a6+4*a2^2*a3^2+16*a2^2*a6-4*a2*a4^2-2*a3^2*a4-8*a4*a6)*x+a1^3*a3*a6-a1^2*a3^2*a4+2*a1^2*a4*a6+a1*a2*a3^3+4*a1*a2*a3*a6-3*a1*a3*a4^2+2*a2*a3^2*a4+8*a2*a4*a6-a3^4-8*a3^2*a6-2*a4^3-16*a6^2));if(n%2==1,m=(n-1)/2;t1=elldivpol0(e,m+2)*elldivpol0(e,m)^3;t2=elldivpol0(e,m-1)*elldivpol0(e,m+1)^3;psi24=(4*f1+f2^2)^2;if(m%2==1,return(t1-psi24*t2),return(psi24*t1-t2)));m=n/2;t1=elldivpol0(e,m+2)*elldivpol0(e,m-1)^2;t2=elldivpol0(e,m-2)*elldivpol0(e,m+1)^2;elldivpol0(e,m)*(t1-t2);
(15:20) gp > e=ellinit([0,-1,1,0,0]);
(15:20) gp > for(j=1,5,print(j,": ",elldivpol(e,j)))
1: 1
2: (2*y + 1)
3: 3*x^4 - 4*x^3 + 3*x - 1
4: (4*y + 2)*x^6 + (-8*y - 4)*x^5 + (20*y + 10)*x^3 + (-20*y - 10)*x^2 + (8*y + 4)*x + (-2*y - 1)
5: 5*x^12 - 20*x^11 + 16*x^10 + 95*x^9 - 285*x^8 + 360*x^7 - 255*x^6 + 94*x^5 + 15*x^4 - 45*x^3 + 25*x^2 - 5*x
\\finding rational 5-torsion points
(15:20) gp > elldivpol(e,5)
%51 = 5*x^12 - 20*x^11 + 16*x^10 + 95*x^9 - 285*x^8 + 360*x^7 - 255*x^6 + 94*x^5 + 15*x^4 - 45*x^3 + 25*x^2 - 5*x
(15:20) gp > polrootsmod(%,2003)
%52 = [Mod(0, 2003), Mod(1, 2003)]~
(15:21) gp > e=ellinit([0,-1,1,0,0]*Mod(1,2003));
(15:21) gp > ellordinate(e,1)
%54 = [Mod(0, 2003), Mod(2002, 2003)]
(15:21) gp > P=[1,2002]
%55 = [1, 2002]
(15:21) gp > ellisoncurve(e,P)
%56 = 1
\\P is on e and is a 5-torsion point
(15:22) gp > ellpow(e,P,5)
%57 = [0]
\\in fact on this curve whatever the field, there will be K-rational points of 5-torsion because, the polynomial of 5-division is not
\\irreducible over Z
(15:22) gp > factor(5*x^12 - 20*x^11 + 16*x^10 + 95*x^9 - 285*x^8 + 360*x^7 - 255*x^6 + 94*x^5 + 15*x^4 - 45*x^3 + 25*x^2 - 5*x)
%58 = 
[x - 1 1]

[x 1]

[5*x^10 - 15*x^9 + x^8 + 96*x^7 - 189*x^6 + 171*x^5 - 84*x^4 + 10*x^3 + 25*x^2 - 20*x + 5 1]

\\n-division polynomial is of degree (n^2-1)/2
(15:24) gp > lift(elldivpol(e,25))
%60 = 25*x^312 + 1406*x^311 + 1963*x^310 + 1182*x^309 + 1214*x^308 + 1109*x^307 + 1192*x^306 + 154*x^305 + 1389*x^304 + 418*x^303 + 1177*x^302 + 518*x^301 + 1999*x^300 + 267*x^299 + 509*x^298 + 1714*x^297 + 990*x^296 + 878*x^295 + 1682*x^294 + 1593*x^293 + 1969*x^292 + 1837*x^291 + 357*x^290 + 613*x^289 + 1538*x^288 + 1183*x^287 + 232*x^286 + 1775*x^285 + 1784*x^284 + 311*x^283 + 1908*x^282 + 1253*x^281 + 91*x^280 + 1340*x^279 + 1857*x^278 + 408*x^277 + 473*x^276 + 1284*x^275 + 722*x^274 + 1809*x^273 + 1081*x^272 + 63*x^271 + 506*x^270 + 703*x^269 + 109*x^268 + 36*x^267 + 970*x^266 + 1059*x^265 + 1053*x^264 + 804*x^263 + 1113*x^262 + 830*x^261 + 1356*x^260 + 287*x^259 + 144*x^258 + 698*x^257 + 88*x^256 + 760*x^255 + 835*x^254 + 1991*x^253 + 1836*x^252 + 1998*x^251 + 1134*x^250 + 1467*x^249 + 1503*x^248 + 759*x^247 + 1384*x^246 + 1555*x^245 + 1959*x^244 + 852*x^243 + 1110*x^242 + 1262*x^241 + 142*x^240 + 356*x^239 + 1127*x^238 + 891*x^237 + 843*x^236 + 980*x^235 + 1390*x^234 + 733*x^233 + 1328*x^232 + 1528*x^231 + 1101*x^230 + 739*x^229 + 1158*x^228 + 722*x^227 + 1011*x^226 + 1631*x^225 + 1180*x^224 + 959*x^223 + 508*x^222 + 894*x^221 + 234*x^220 + 489*x^219 + 588*x^218 + 18*x^217 + 426*x^216 + 165*x^215 + 746*x^214 + 1127*x^213 + 291*x^212 + 671*x^211 + 752*x^210 + 1430*x^209 + 984*x^208 + 288*x^207 + 1633*x^206 + 1912*x^205 + 437*x^204 + 892*x^203 + 530*x^202 + 407*x^201 + 771*x^200 + 1894*x^199 + 1119*x^198 + 957*x^197 + 882*x^196 + 1271*x^195 + 1141*x^194 + 221*x^193 + 1836*x^192 + 742*x^191 + 230*x^190 + 150*x^189 + 801*x^188 + 1654*x^187 + 1975*x^186 + 428*x^185 + 694*x^184 + 1742*x^183 + 35*x^182 + 1883*x^181 + 1528*x^180 + 1178*x^179 + 1349*x^178 + 1259*x^177 + 777*x^176 + 1326*x^175 + 418*x^174 + 412*x^173 + 1402*x^172 + 142*x^171 + 1028*x^170 + 1883*x^169 + 1105*x^168 + 638*x^167 + 779*x^166 + 1910*x^165 + 1173*x^164 + 1987*x^163 + 753*x^162 + 1844*x^161 + 1795*x^160 + 22*x^159 + 136*x^158 + 1305*x^157 + 757*x^156 + 1210*x^155 + 703*x^154 + 1826*x^153 + 1980*x^152 + 1753*x^151 + 1086*x^150 + 264*x^149 + 1745*x^148 + 1372*x^147 + 1717*x^146 + 1562*x^145 + 369*x^144 + 1497*x^143 + 1825*x^142 + 1997*x^141 + 1810*x^140 + 866*x^139 + 1866*x^138 + 1658*x^137 + 1762*x^136 + 988*x^135 + 686*x^134 + 1833*x^133 + 677*x^132 + 961*x^131 + 723*x^130 + 820*x^129 + 1986*x^128 + 1581*x^127 + 1064*x^126 + 1955*x^125 + 1454*x^124 + 126*x^123 + 1491*x^122 + 709*x^121 + 1596*x^120 + 1531*x^119 + 156*x^118 + 579*x^117 + 1756*x^116 + 1718*x^115 + 1333*x^114 + 1945*x^113 + 437*x^112 + 1460*x^111 + 1729*x^110 + 1366*x^109 + 1327*x^108 + 1134*x^107 + 390*x^106 + 1944*x^105 + 111*x^104 + 1368*x^103 + 1359*x^102 + 809*x^101 + 59*x^100 + 113*x^99 + 759*x^98 + 796*x^97 + 172*x^96 + 1641*x^95 + 899*x^94 + 160*x^93 + 1516*x^92 + 1350*x^91 + 1992*x^90 + 791*x^89 + 1137*x^88 + 1714*x^87 + 310*x^86 + 61*x^85 + 64*x^84 + 330*x^83 + 1613*x^82 + 1081*x^81 + 1878*x^80 + 1381*x^79 + 195*x^78 + 202*x^77 + 572*x^76 + 539*x^75 + 1115*x^74 + 1257*x^73 + 346*x^72 + 422*x^71 + 1546*x^70 + 783*x^69 + 1093*x^68 + 436*x^67 + 1099*x^66 + 1023*x^65 + 74*x^64 + 1173*x^63 + 523*x^62 + 1796*x^61 + 201*x^60 + 1955*x^59 + 1706*x^58 + 1718*x^57 + 721*x^56 + 1979*x^55 + 1419*x^54 + 215*x^53 + 1849*x^52 + 286*x^51 + 1974*x^50 + 508*x^49 + 650*x^48 + 1737*x^47 + 1506*x^46 + 1153*x^45 + 1736*x^44 + 1865*x^43 + 424*x^42 + 1912*x^41 + 737*x^40 + 105*x^39 + 78*x^38 + 1474*x^37 + 1830*x^36 + 1460*x^35 + 617*x^34 + 440*x^33 + 1816*x^32 + 545*x^31 + 279*x^30 + 1573*x^29 + 708*x^28 + 1915*x^27 + 53*x^26 + 1627*x^25 + 274*x^24 + 1825*x^23 + 248*x^22 + 956*x^21 + 52*x^20 + 1101*x^19 + 730*x^18 + 1292*x^17 + 939*x^16 + 950*x^15 + 1575*x^14 + 1504*x^13 + 693*x^12 + 323*x^11 + 370*x^10 + 1691*x^9 + 453*x^8 + 1286*x^7 + 317*x^6 + 113*x^5 + 372*x^4 + 1566*x^3 + 1122*x^2 + 1978*x
(15:25) gp > (25^2-1)/2
%62 = 312
\\l-ogeneous curves have related j-invariants via the modular equation Phi_l
(15:25) gp > E=ellinit([0,0,0,1132,278]*Mod(1,2003))
%63 = [Mod(0, 2003), Mod(0, 2003), Mod(0, 2003), Mod(1132, 2003), Mod(278, 2003), Mod(0, 2003), Mod(261, 2003), Mod(1112, 2003), Mod(496, 2003), Mod(1748, 2003), Mod(168, 2003), Mod(1707, 2003), Mod(171, 2003), 0, 0, 0, 0, 0, 0]
(15:32) gp > Eh=ellinit([0,0,0,500,1005]*Mod(1,2003))
%64 = [Mod(0, 2003), Mod(0, 2003), Mod(0, 2003), Mod(500, 2003), Mod(1005, 2003), Mod(0, 2003), Mod(1000, 2003), Mod(14, 2003), Mod(375, 2003), Mod(36, 2003), Mod(982, 2003), Mod(744, 2003), Mod(515, 2003), 0, 0, 0, 0, 0, 0]
(15:32) gp > E.j
%65 = Mod(171, 2003)
(15:32) gp > Eh.j
%66 = Mod(515, 2003)
\\Phi2(E.j,Eh.j) = 0
(15:33) gp > Phi2=X^3+Y^3-X^2*Y^2+1488*(X*Y^2+X^2*Y)-162000*(X^2+Y^2)+40773375*X*Y+8748000000*(X+Y)-157464000000000
%67 = X^3 + (-Y^2 + 1488*Y - 162000)*X^2 + (1488*Y^2 + 40773375*Y + 8748000000)*X + (Y^3 - 162000*Y^2 + 8748000000*Y - 157464000000000)
(15:34) gp > subst(Phi2,X,E.j)
%68 = Y^3 + Mod(1114, 2003)*Y^2 + Mod(1514, 2003)*Y + Mod(1041, 2003)
(15:34) gp > subst(%,Y,Eh.j)
%69 = Mod(0, 2003)
(15:35) gp > polrootsmod(lift(%68),2003)
%70 = [Mod(515, 2003)]~