Trac #34174: insufficient precision in scaling elliptic curves over n…

…umber fields by units

<this was first reported on the sage-devel google group:
https://groups.google.com/g/sage-devel/c/s0B7OqpB0KU/m/18eHSiRWAAAJ ;
as requested I'm opening this ticket>

I am running into an issue with computing isogeny classes of elliptic
curves over number fields.

Here is what I entered:

{{{
K.<a> = QuadraticField(4569)
myJ = 46969655/32768
E = EllipticCurve(j=K(myJ))
C = E.isogeny_class()
}}}

This raises a `KeyError` in the first instance, which seems to be
handled via a direct call to `IsogenyClass_EC_NumberField`, but this
then raises a `ValueError: Cannot convert infinity or NaN to Sage
Integer`. See the above linked discussion for the full stacktrace.

This error has occurred on sage-9.4 on a system whose `uname -a` is
`Linux LEGENDRE 5.4.0-91-generic #102-Ubuntu SMP Fri Nov 5 16:31:28 UTC
2021 x86_64 x86_64 x86_64 GNU/Linux`, as well as on sage-9.7-beta5 on
`Linux Barinder 5.4.0-121-generic #137-Ubuntu SMP Wed Jun 15 13:33:07
UTC 2022 x86_64 x86_64 x86_64 GNU/Linux`.

Changing the field K to many other quadratic fields does not yield any
error. However I did also notice this error with many j-invariants in
`QuadraticField(6537)`, so it somehow seems to be base-field dependent.

Two issues are fixed in this ticket: (1) computing the isogeny class
could fail for a curve defined by a non-integral model; (2) precision
was not handled well in scaling equations by units.

URL: https://trac.sagemath.org/34174
Reported by: gh-BarinderBanwait
Ticket author(s): John Cremona
Reviewer(s): David Lowry-Duda

build/pkgs/configure/checksums.ini

@@ -1,4 +1,4 @@
 tarball=configure-VERSION.tar.gz
-sha1=7df04dc0c2e0d7cdfd767dfb86dd3a60d3af1767
-md5=f8e92133f5447c39196fc45a026ed49d
-cksum=1963708627
+sha1=32693d220d5b1c5d743b33e2a5ab8b1cbd07d077
+md5=ee385184c3a968a0dcf6534c538ff25d
+cksum=4278547201

build/pkgs/configure/package-version.txt

@@ -1 +1 @@
-e3ababd22689e7cc6e89a01b1e8bd6d770e7b1cf
+74ea2f1b4b1eb592f2f9b1aaba35cc3bc09d5af3

src/sage/schemes/elliptic_curves/ell_number_field.py

@@ -775,27 +775,44 @@ def _scale_by_units(self):
             sage: E1 = E.scale_curve(u^5)
             sage: E1._scale_by_units().ainvs() == E.ainvs()
             True
+
+        TESTS:
+
+        See :trac:`34174`.  This used to raise an error due to insufficient precision::
+
+            sage: K.<a> = QuadraticField(4569)
+            sage: j = 46969655/32768
+            sage: E = EllipticCurve(j=K(j))
+            sage: C = E.isogeny_class()
         """
         K = self.base_field()
         r1, r2 = K.signature()
         if r1 + r2 == 1:  # unit rank is 0
             return self
 
-        prec = 1000  # lower precision works badly!
-        embs = K.places(prec=prec)
         degs = [1]*r1 + [2]*r2
         fu = K.units()
-        from sage.matrix.all import Matrix
-        U = Matrix([[e(u).abs().log()*d for d,e in zip(degs,embs)] for u in fu])
-        A = U*U.transpose()
-        Ainv = A.inverse()
-
         c4, c6 = self.c_invariants()
-        c4s = [e(c4) for e in embs]
-        c6s = [e(c6) for e in embs]
+
+        from sage.matrix.all import Matrix
         from sage.modules.free_module_element import vector
-        v = vector([(x4.abs().nth_root(4)+x6.abs().nth_root(6)).log()*d for x4,x6,d in zip(c4s,c6s,degs)])
-        es = [e.round() for e in -Ainv*U*v]
+
+        prec = 1000 # initial value, will be increased if necessary
+        ok = False
+        while not ok:
+            embs = K.places(prec=prec)
+            c4s = [e(c4) for e in embs]
+            c6s = [e(c6) for e in embs]
+
+            U = Matrix([[e(u).abs().log()*d for d,e in zip(degs,embs)] for u in fu])
+            v = vector([(x4.abs().nth_root(4)+x6.abs().nth_root(6)).log()*d for x4,x6,d in zip(c4s,c6s,degs)])
+            w = -(U*U.transpose()).inverse()*U*v
+            try:
+                es = [e.round() for e in w]
+                ok = True
+            except ValueError:
+                prec *= 2
+
         u = prod([uj**ej for uj,ej in zip(fu,es)])
         return self.scale_curve(u)
 

src/sage/schemes/elliptic_curves/gal_reps_number_field.py

@@ -468,7 +468,8 @@ def _non_surjective(E, patience=100):
 
 
 def Frobenius_filter(E, L, patience=100):
-    r""" Determine which primes in L might have an image contained in a
+    r"""
+    Determine which primes in L might have an image contained in a
     Borel subgroup, by checking of traces of Frobenius.
 
     .. NOTE::
@@ -526,7 +527,7 @@ def Frobenius_filter(E, L, patience=100):
         sage: [len(E.isogenies_prime_degree(l)) for l in [2,3]]
         [1, 1]
     """
-    E = _over_numberfield(E)
+    E = _over_numberfield(E).global_integral_model()
     K = E.base_field()
 
     L = list(set(L)) # Remove duplicates from L and makes a copy for output
@@ -1138,7 +1139,8 @@ def Billerey_P_l(E, l):
 
     INPUT:
 
-    - ``E`` -- an elliptic curve over a number field `K`
+    - ``E`` -- an elliptic curve over a number field `K`, given by a
+      global integral model.
 
     - ``l`` -- a rational prime
 
@@ -1171,7 +1173,8 @@ def Billerey_B_l(E,l,B=0):
 
     INPUT:
 
-    - ``E`` -- an elliptic curve over a number field `K`
+    - ``E`` -- an elliptic curve over a number field `K`, given by a
+      global integral model.
 
     - ``l`` (int) -- a rational prime
 
@@ -1211,7 +1214,8 @@ def Billerey_R_q(E, q, B=0):
 
     INPUT:
 
-    - ``E`` -- an elliptic curve over a number field `K`
+    - ``E`` -- an elliptic curve over a number field `K`, given by a
+      global integral model.
 
     - ``q`` -- a prime ideal of `K`
 
@@ -1247,7 +1251,7 @@ def Billerey_R_q(E, q, B=0):
 
 
 def Billerey_B_bound(E, max_l=200, num_l=8, small_prime_bound=0, debug=False):
-    """
+    r"""
     Compute Billerey's bound `B`.
 
     We compute `B_l` for `l` up to ``max_l`` (at most) until ``num_l``
@@ -1258,7 +1262,8 @@ def Billerey_B_bound(E, max_l=200, num_l=8, small_prime_bound=0, debug=False):
 
     INPUT:
 
-    - ``E`` -- an elliptic curve over a number field `K`.
+    - ``E`` -- an elliptic curve over a number field `K`, given by a
+      global integral model.
 
     - ``max_l`` (int, default 200) -- maximum size of primes l to check.
 
@@ -1364,7 +1369,8 @@ def Billerey_R_bound(E, max_l=200, num_l=8, small_prime_bound=None, debug=False)
 
     INPUT:
 
-    - ``E`` -- an elliptic curve over a number field `K`.
+    - ``E`` -- an elliptic curve over a number field `K`, given by a
+      global integral model.
 
     - ``max_l`` (int, default 200) -- maximum size of rational primes
       l for which the primes q above l are checked.
@@ -1512,6 +1518,18 @@ def reducible_primes_Billerey(E, num_l=None, max_l=None, verbose=False):
         sage: E = EllipticCurve_from_j(K(2268945/128)).global_minimal_model() # c.f. [Sut2012]
         sage: reducible_primes_Billerey(E)
         [7]
+
+    TESTS:
+
+    Test that this function works with non-integral models (see :trac:`34174`)::
+
+        sage: K.<a> = QuadraticField(4569)
+        sage: j = 46969655/32768
+        sage: E = EllipticCurve(j=K(j))
+        sage: EK = E.change_ring(K)
+        sage: C = EK.isogeny_class(minimal_models=False)
+        sage: len(C)
+        4
     """
     #verbose=True
     if verbose:
@@ -1525,8 +1543,12 @@ def reducible_primes_Billerey(E, num_l=None, max_l=None, verbose=False):
 
     K = E.base_field()
     DK = K.discriminant()
-    ED = E.discriminant().norm()
-    B0 = ZZ(6*DK*ED).prime_divisors()  # TODO: only works if discriminant is integral
+
+    # We replace E by an integral model if necessary, since this
+    # function and the helper functions need this:
+    E1 = E.global_integral_model()
+    ED = E1.discriminant().norm()
+    B0 = ZZ(6*DK*ED).prime_divisors()
 
     # Billeray's algorithm will be faster if we tell it to ignore
     # small primes; these can be tested using the naive algorithm.
@@ -1535,16 +1557,16 @@ def reducible_primes_Billerey(E, num_l=None, max_l=None, verbose=False):
         print("First doing naive test of primes up to {}...".format(max_l))
 
     max_small_prime = 200
-    OK_small_primes = reducible_primes_naive(E, max_l=max_small_prime, num_P=200, verbose=verbose)
+    OK_small_primes = reducible_primes_naive(E1, max_l=max_small_prime, num_P=200, verbose=verbose)
     if verbose:
         print("Naive test of primes up to {} returns {}.".format(max_small_prime, OK_small_primes))
 
-    B1 = Billerey_B_bound(E, max_l, num_l, max_small_prime, verbose)
+    B1 = Billerey_B_bound(E1, max_l, num_l, max_small_prime, verbose)
     if B1 == [0]:
         if verbose:
             print("...  B_bound ineffective using max_l={}, moving on to R-bound".format(max_l))
 
-        B1 = Billerey_R_bound(E,max_l, num_l, max_small_prime, verbose)
+        B1 = Billerey_R_bound(E1,max_l, num_l, max_small_prime, verbose)
         if B1 == [0]:
             if verbose:
                 print("... R_bound ineffective using max_l={}",format(max_l))
@@ -1559,7 +1581,7 @@ def reducible_primes_Billerey(E, num_l=None, max_l=None, verbose=False):
         print("... combined bound = {}".format(B))
 
     num_p = 100
-    B = Frobenius_filter(E, B, num_p)
+    B = Frobenius_filter(E1, B, num_p)
     if verbose:
         print("... after Frobenius filter = {}".format(B))
     return B