Compute elements of n-torsion group of elliptic curve over finite field

edit

The following worked for me:

sage: E = EllipticCurve(GF(43), [0, 6])
....: [P.xy() for P in E if P.order() == 13]
....: 
[(13, 15),
 (13, 28),
 (26, 9),
 (26, 34),
 (27, 9),
 (27, 34),
 (33, 9),
 (33, 34),
 (35, 15),
 (35, 28),
 (38, 15),
 (38, 28)]

Also a slightly bigger case:

sage: E = EllipticCurve(GF(2027), [0, 1])
....: [P.xy() for P in E if P.order() == 13]
[(155, 48),
 (155, 1979),
 (579, 647),
 (579, 1380),
 (1304, 910),
 (1304, 1117),
 (1584, 379),
 (1584, 1648),
 (1732, 836),
 (1732, 1191),
 (1819, 992),
 (1819, 1035)]

For a "much bigger curve" (i.e. with more rational points) some improvements are necessary. For instance:

p = (3*10^7).next_prime()
E = EllipticCurve(GF(p), [0, 4])
ord = E.order()
print(f'E is the following curve:\n{E}')
print(f'E has order {ord} = {ord.factor()}')
print(f'E has the generator(s): {E.gens()}')

We get the following information on E:

E is the following curve:
Elliptic Curve defined by y^2 = x^3 + 4 over Finite Field of size 30000001
E has order 30010071 = 3 * 7 * 13 * 37 * 2971
E has the generator(s): ((14044277 : 14356696 : 1),)

Now

[P.xy() for P in E if P.order() == 13]

takes a looong time. But we can immediately get the elements of order $13$ as follows:

G = E.gens()[0]    # we already know there is one generator of order <ord>
n = ZZ(ord/13)
Q = n*G
[(k*Q).xy() for k in [1..12]]

This gives:

[(28289013, 19261067),
 (11842435, 11155846),
 (26389003, 19261067),
 (5321986, 10738934),
 (15676831, 11155846),
 (2480735, 11155846),
 (2480735, 18844155),
 (15676831, 18844155),
 (5321986, 19261067),
 (26389003, 10738934),
 (11842435, 18844155),
 (28289013, 10738934)]

Check:

sage: E( (2480735, 18844155) ).order()
13

Nowadays, you can use .division_field() to find the minimal extension where the full $\ell$-torsion is defined, and .torsion_basis() to find two points generating the $\ell$-torsion subgroup:

sage: F = E.division_field(13)
sage: EE = E.change_ring(F); EE
Elliptic Curve defined by y^2 = x^3 + 6 over Finite Field in t of size 43^6
sage: P,Q = EE.torsion_basis(13); P,Q
((41*t^5 + 23*t^4 + 27*t^3 + 11*t^2 + 9*t + 2 : 36*t^5 + 23*t^4 + 17*t^3 + 16*t^2 + 14*t + 11 : 1),
 (6*t^5 + 29*t^4 + 16*t^3 + 22*t^2 + 21*t + 30 : 32*t^5 + 28*t^4 + 3*t^3 + 24*t^2 + 11*t + 16 : 1))

The rest of the $\ell$-torsion can then be enumerated by taking linear combinations of the generating points:

sage: pts = {i*P+j*Q for i in range(13) for j in range(13)}
sage: len(pts)
169
sage: {13 * T for T in pts}
{(0 : 1 : 0)}