@@ -991,10 +991,18 @@ def abelian_group(self):
991991 self .gens .set_cache (gens )
992992 return AdditiveAbelianGroupWrapper (self .point_homset (), gens , orders )
993993
994- def torsion_basis (self , n ):
994+ def torsion_basis (self , n , * , algorithm = None ):
995995 r"""
996- Return a basis of the `n`-torsion subgroup of this elliptic curve,
997- assuming it is fully rational.
996+ Return a basis of the `n`-torsion subgroup of this elliptic curve.
997+
998+ INPUT:
999+
1000+ - ``n`` -- integer
1001+
1002+ - ``algorithm`` -- string (default: ``None``).
1003+ Currently available choices are ``"random"``, ``"structure"``,
1004+ and ``"divpoly"``. If ``algorithm`` is ``None``, the method
1005+ attempts to select the most suitable algorithm automatically.
9981006
9991007 EXAMPLES::
10001008
@@ -1017,7 +1025,7 @@ def torsion_basis(self, n):
10171025 sage: E.torsion_basis(23)
10181026 Traceback (most recent call last):
10191027 ...
1020- ValueError: curve does not have full rational 23-torsion
1028+ ValueError: curve does not have full rational 23-torsion, or very unlikely event
10211029 sage: F = E.division_field(23); F
10221030 Finite Field in t of size 101^11
10231031 sage: EE = E.change_ring(F)
@@ -1035,25 +1043,140 @@ def torsion_basis(self, n):
10351043 + 55*z11^5 + 23*z11^4 + 17*z11^3 + 90*z11^2 + 91*z11 + 68
10361044 : 1)
10371045
1046+ TESTS:
1047+
1048+ Check on random curves that all three algorithms return
1049+ equivalent results::
1050+
1051+ sage: while True:
1052+ ....: p = random_prime(100)
1053+ ....: e = randrange(1,4)
1054+ ....: E = choice(EllipticCurve(j=GF((p,e)).random_element()).twists())
1055+ ....: ds = E.abelian_group().invariants()
1056+ ....: if len(ds) < 2:
1057+ ....: continue
1058+ ....: ns = set(gcd(ds).divisors()) - {1}
1059+ ....: if ns: break
1060+ sage: n = choice(list(ns))
1061+ sage: assert n >= 2
1062+ sage: P1, Q1 = E.torsion_basis(n, algorithm='random')
1063+ sage: A1 = AdditiveAbelianGroupWrapper.from_generators([P1, Q1])
1064+ sage: P2, Q2 = E.torsion_basis(n, algorithm='structure')
1065+ sage: A2 = AdditiveAbelianGroupWrapper.from_generators([P2, Q2])
1066+ sage: P3, Q3 = E.torsion_basis(n, algorithm='divpoly')
1067+ sage: A3 = AdditiveAbelianGroupWrapper.from_generators([P3, Q3])
1068+ sage: assert A1 == A2
1069+ sage: assert A1 == A3
1070+ sage: assert A2 == A3
1071+
10381072 .. SEEALSO::
10391073
10401074 Use :meth:`~sage.schemes.elliptic_curves.ell_field.EllipticCurve_field.division_field`
1041- to determine a field extension containing the full `\ell `-torsion subgroup.
1075+ to determine a field extension containing the full `n `-torsion subgroup.
10421076
10431077 ALGORITHM:
10441078
1045- This method currently uses :meth:`abelian_group` and
1079+ If ``algorithm`` is ``"random"``, this method repeatedly
1080+ samples random points on the curve and distills a basis of the
1081+ `n`-torsion from them. This requires point counting.
1082+
1083+ If ``algorithm`` is ``divpoly``, this method uses division
1084+ polynomials to construct a basis of the `n`-torsion. The
1085+ complexity of this approach scales with the size of the prime
1086+ factors of `n`. This algorithm is usually much slower than
1087+ the others, but there might be situations in which it is
1088+ useful.
1089+
1090+ If ``algorithm`` is ``"structure"``, this method calls
1091+ :meth:`abelian_group` and
10461092 :meth:`AdditiveAbelianGroupWrapper.torsion_subgroup`.
1093+ Theoretically, this involves performing a superset of the work
1094+ of the ``"random"`` method, so it should never be the best
1095+ choice, but due to implementation details (PARI vs. Sage) this
1096+ approach can be faster than the ``"random"`` algorithm in
1097+ practice.
10471098 """
1048- # TODO: In many cases this is not the fastest algorithm.
1049- # Alternatives include factoring division polynomials and
1050- # random sampling (like PARI's ellgroup, but with a milder
1051- # termination condition). We should implement these too
1052- # and figure out when to use which.
1053- T = self .abelian_group ().torsion_subgroup (n )
1054- if T .invariants () != (n , n ):
1099+ n = ZZ (n )
1100+ if n < 2 :
1101+ raise ValueError ('computing an n-torsion basis only makes sense for n >= 2' )
1102+
1103+ if self .base_field ().characteristic ().divides (n ):
10551104 raise ValueError (f'curve does not have full rational { n } )
1056- return tuple (P .element () for P in T .gens ())
1105+
1106+ if algorithm is None :
1107+ if self .abelian_group .is_in_cache ():
1108+ algorithm = 'structure'
1109+ else :
1110+ algorithm = 'random'
1111+
1112+ if algorithm == 'random' :
1113+ # similar to AdditiveAbelianGroupWrapper.from_generators()
1114+ from sage .arith .misc import xlcm
1115+
1116+ N = self .cardinality ()
1117+ P = Q = self .zero ()
1118+ P ._order = ZZ (1 )
1119+
1120+ for step in range (999 ):
1121+ # check P,Q is a basis of the subgroup <P,Q> with ord(Q) | ord(P)
1122+ assert Q ._order .divides (P ._order )
1123+ # assert generic.has_order(P.weil_pairing(Q, P._order), Q._order, operation='*')
1124+
1125+ if Q ._order == n :
1126+ P *= P ._order // n
1127+ assert hasattr (P , '_order' )
1128+ break
1129+
1130+ cof = N .prime_to_m_part (n // Q ._order )
1131+ T = cof * self .random_point ()
1132+ T .set_order (multiple = N // cof , check = False )
1133+
1134+ # extend P using T as much as possible
1135+ m , k1 , k2 = xlcm (P ._order , T ._order )
1136+ m1 = P ._order // k1
1137+ m2 = T ._order // k2
1138+ P = m1 * P + m2 * T
1139+ P ._order = m
1140+
1141+ if not step :
1142+ continue
1143+
1144+ # remove the P component from T
1145+ l = generic .order_from_multiple (P .weil_pairing (T , P ._order ), P ._order , operation = '*' )
1146+ x = (l * T ).log (l * P )
1147+ T -= x * P
1148+ T .set_order (multiple = P ._order , check = False )
1149+
1150+ # for Q we only need the n-torsion part of T
1151+ T *= T ._order // n .gcd (T ._order )
1152+
1153+ # extend Q as much as possible
1154+ Q , m = generic .merge_points ((Q , Q ._order ), (T , T ._order ))
1155+ Q ._order = m
1156+
1157+ # remove the P component from Q
1158+ l = generic .order_from_multiple (P .weil_pairing (Q , P ._order ), P ._order , operation = '*' )
1159+ y = (l * Q ).log (l * P )
1160+ Q -= y * P
1161+ Q .set_order (multiple = P ._order , check = False )
1162+
1163+ else :
1164+ raise ValueError (f'curve does not have full rational { n } )
1165+
1166+ # assert generic.has_order(P.weil_pairing(Q, n), n, operation='*')
1167+
1168+ return P , Q
1169+
1170+ elif algorithm == 'divpoly' :
1171+ return super ().torsion_basis (n , algorithm = algorithm )
1172+
1173+ elif algorithm == 'structure' :
1174+ T = self .abelian_group ().torsion_subgroup (n )
1175+ if T .invariants () != (n , n ):
1176+ raise ValueError (f'curve does not have full rational { n } )
1177+ return tuple (P .element () for P in T .gens ())
1178+
1179+ raise ValueError (f'unknown algorithm { algorithm !r} )
10571180
10581181 def is_isogenous (self , other , field = None , proof = True ):
10591182 """
0 commit comments