Trac #16456: Bug in descend_to method for elliptic curves

The function descend_to which was implemented in #9384 is incorrect.
The twisting parameter d in L* is only defined modulo {{{(L*)^2}}} and
one has to determine whether there is a representative which lies in K*,
which the implementation does not do, so some "positive" results are
missed.  Here is an example of this.
{{{
sage: K.<a> = NumberField(x^3-2)
sage: E = EllipticCurve('11a1').quadratic_twist(2)
sage: EK = E.change_ring(K)
sage: EK2 = EK.change_weierstrass_model((a,a,a,a+1))
sage: EK2.descend_to(QQ)
<None>
}}}
The heart of the problem being solved here is to determine whether
(given that j(E) is in K of course) the twisting parameter in
{{{L*/(L*)^2}}} is in the image of {{{K*/(K*)^2}}} or not, and the
implementation ignores this.  (For j=0 or 1728 the principle is the same
with squares replaced by 6th or 4th powers respectively.)

I can fix this for number fields (one can restrict from {{{K*/(K*)^2}}}
to a finite K(S,2) for an easily determined set S of primes) or for
finite fields, but it may not be possible to have this implemented for
arbitrary fields, which will cause a problem with the patching.

A second and independent bug is reported by Warren Moore:
{{{
sage: k.<i> = QuadraticField(-1)
sage: E = EllipticCurve(k,[0,0,0,1,0])
sage: E.descend_to(QQ) == None
True
}}}
This caused by a "naked Except" clause plus a call to f.preimage() for a
map f which has no attribute/method "preimage".

URL: http://trac.sagemath.org/16456
Reported by: cremona
Ticket author(s): John Cremona
Reviewer(s): Peter Bruin

src/sage/misc/binary_tree.pxd

@@ -3,7 +3,8 @@ include 'sage/ext/python.pxi'
 
 cdef struct binary_tree_node:
     int key
-    binary_tree_node *left, *right
+    binary_tree_node *left
+    binary_tree_node *right
     void *value
 
 #cdef binary_tree_node *BinaryTreeNode(int, object)

src/sage/rings/number_field/number_field.py

@@ -104,6 +104,7 @@
 from sage.misc.functional import is_odd
 from sage.misc.misc_c import prod
 from sage.categories.homset import End
+from sage.rings.all import Infinity
 
 import sage.rings.ring
 from sage.misc.latex import latex_variable_name
@@ -3819,50 +3820,99 @@ def _S_class_group_quotient_matrix(self, S):
         assert A[:c] == 1 and A[c:] == 0
         return Q[:c, :a]
 
-    def selmer_group(self, S, m, proof=True):
+    def selmer_group(self, S, m, proof=True, orders=False):
         r"""
-        Compute the Selmer group `K(S,m)`.  This is the subgroup of
-        `K^\times/(K^\times)^m` consisting of elements `a` such that
-        `K(\sqrt[m]{a})/K` is unramified at all primes of `K` outside
-        of `S`.
+        Compute the group `K(S,m)`.
 
         INPUT:
 
-        - ``S`` - A set of primes of self.
+        - ``S`` -- a set of primes of ``self``
 
-        - ``m`` - A positive integer.
+        - ``m`` -- a positive integer
 
-        - ``proof`` - If False, assume the GRH in computing the class group.
+        - ``proof`` -- if False, assume the GRH in computing the class group
+
+        - ``orders`` (default False) -- if True, output two lists, the
+          generators and their orders
 
         OUTPUT:
 
-        A list of generators of `K(S,m)`.
+        A list of generators of `K(S,m)`, and (optionally) their
+        orders as elements of `K^\times/(K^\times)^m`.  This is the
+        subgroup of `K^\times/(K^\times)^m` consisting of elements `a`
+        such that the valuation of `a` is divisible by `m` at all
+        primes not in `S`.  It fits in an exact sequence between the
+        units modulo `m`-th powers and the `m`-torsion in the
+        `S`-class group:
+
+        .. math::
+
+            1                                    \longrightarrow
+            O_{K,S}^\times / (O_{K,S}^\times)^m  \longrightarrow
+            K(S,m)                               \longrightarrow
+            \operatorname{Cl}_{K,S}[m]           \longrightarrow
+            0.
+
+        The group `K(S,m)` contains the subgroup of those `a` such
+        that `K(\sqrt[m]{a})/K` is unramified at all primes of `K`
+        outside of `S`, but may contain it properly when not all
+        primes dividing `m` are in `S`.
 
         EXAMPLES::
 
             sage: K.<a> = QuadraticField(-5)
             sage: K.selmer_group((), 2)
             [-1, 2]
-            sage: K.selmer_group([K.ideal(2, -a+1)], 2)
+
+        The previous example shows that the group generated by the
+        output may be strictly larger than the 'true' Selmer group of
+        elements giving extensions unramified outside `S`, since that
+        has order just 2, generated by `-1`::
+
+            sage: K.class_number()
+            2
+            sage: K.hilbert_class_field('b')
+            Number Field in b with defining polynomial x^2 + 1 over its base field
+
+        When `m` is prime all the orders are equal to `m`, but in general they are only divisors of `m`::
+
+            sage: K.<a> = QuadraticField(-5)
+            sage: P2 = K.ideal(2, -a+1)
+            sage: P3 = K.ideal(3, a+1)
+            sage: K.selmer_group((), 2, orders=True)
+            ([-1, 2], [2, 2])
+            sage: K.selmer_group((), 4, orders=True)
+            ([-1, 4], [2, 2])
+            sage: K.selmer_group([P2], 2)
             [2, -1]
-            sage: K.selmer_group([K.ideal(2, -a+1), K.ideal(3, a+1)], 2)
-            [2, a + 1, -1]
-            sage: K.selmer_group((K.ideal(2, -a+1),K.ideal(3, a+1)), 4)
+            sage: K.selmer_group((P2,P3), 4)
             [2, a + 1, -1]
-            sage: K.selmer_group([K.ideal(2, -a+1)], 3)
+            sage: K.selmer_group((P2,P3), 4, orders=True)
+            ([2, a + 1, -1], [4, 4, 2])
+            sage: K.selmer_group([P2], 3)
             [2]
-            sage: K.selmer_group([K.ideal(2, -a+1), K.ideal(3, a+1)], 3)
+            sage: K.selmer_group([P2, P3], 3)
             [2, a + 1]
-            sage: K.selmer_group([K.ideal(2, -a+1), K.ideal(3, a+1), K.ideal(a)], 3)  # random signs
+            sage: K.selmer_group([P2, P3, K.ideal(a)], 3)  # random signs
             [2, a + 1, a]
+
+        Example over `\QQ` (as a number field)::
+
             sage: K.<a> = NumberField(polygen(QQ))
             sage: K.selmer_group([],5)
             []
+            sage: K.selmer_group([K.prime_above(p) for p in [2,3,5]],2)
+            [2, 3, 5, -1]
+            sage: K.selmer_group([K.prime_above(p) for p in [2,3,5]],6, orders=True)
+            ([2, 3, 5, -1], [6, 6, 6, 2])
 
         TESTS::
 
             sage: K.<a> = QuadraticField(-5)
-            sage: S = K.selmer_group([K.ideal(2, -a+1), K.ideal(3, a+1), K.ideal(a)], 3)
+            sage: P2 = K.ideal(2, -a+1)
+            sage: P3 = K.ideal(3, a+1)
+            sage: P5 = K.ideal(a)
+            sage: S = K.selmer_group([P2, P3, P5], 3)
             sage: S == [2, a + 1, a] or S == [2, a + 1, -a]
             True
 
@@ -3879,21 +3929,26 @@ def selmer_group(self, S, m, proof=True):
             sage: S = K.selmer_group((), 4)
             sage: all(4.divides(x.valuation(p)) for x in S)
             True
-
         """
         units, clgp_gens = self._S_class_group_and_units(tuple(S), proof=proof)
         gens = []
-        from sage.rings.infinity import Infinity
+        ords = []
         for unit in units:
             order = unit.multiplicative_order()
-            if order == Infinity or order.gcd(m) != 1:
+            if order == Infinity:
                 gens.append(unit)
+                ords.append(m)
+            else:
+                m1 = order.gcd(m)
+                if m1!= 1:
+                    gens.append(unit)
+                    ords.append(m1)
         card_S = len(S)
         if card_S != 0:
             from sage.matrix.constructor import Matrix
             H = self.class_group()
-            ords = [g.order() for g in H.gens()]
-            pari_ords = pari(ords).Col()
+            gen_ords = [g.order() for g in H.gens()]
+            pari_ords = pari(gen_ords).Col()
             Sords = [H(s).order() for s in S]
             MS = Matrix(ZZ, [H(s).exponents() for s in S]).transpose()
             pari_MS = pari(MS)
@@ -3913,7 +3968,57 @@ def selmer_group(self, S, m, proof=True):
                     Spart = prod([S[i] ** (exps[i] % Sords[i]) for i in range(card_S)])
                     J *= Spart
                 gens.append(self(J.gens_reduced()[0]))
-        return gens
+                ords.append(d)
+        if orders:
+            return gens, ords
+        else:
+            return gens
+
+    def selmer_group_iterator(self, S, m, proof=True):
+        r"""
+        Return an iterator through elements of the finite group `K(S,m)`.
+
+        INPUT:
+
+        - ``S`` -- a set of primes of ``self``
+
+        - ``m`` -- a positive integer
+
+        - ``proof`` -- if False, assume the GRH in computing the class group
+
+        OUTPUT:
+
+        An iterator yielding the distinct elements of `K(S,m)`.  See
+        the docstring for :meth:`NumberField_generic.selmer_group` for
+        more information.
+
+        EXAMPLES::
+
+            sage: K.<a> = QuadraticField(-5)
+            sage: list(K.selmer_group_iterator((), 2))
+            [1, 2, -1, -2]
+            sage: list(K.selmer_group_iterator((), 4))
+            [1, 4, -1, -4]
+            sage: list(K.selmer_group_iterator([K.ideal(2, -a+1)], 2))
+            [1, -1, 2, -2]
+            sage: list(K.selmer_group_iterator([K.ideal(2, -a+1), K.ideal(3, a+1)], 2))
+            [1, -1, a + 1, -a - 1, 2, -2, 2*a + 2, -2*a - 2]
+
+        Examples over `\QQ` (as a number field)::
+
+            sage: K.<a> = NumberField(polygen(QQ))
+            sage: list(K.selmer_group_iterator([], 5))
+            [1]
+            sage: list(K.selmer_group_iterator([], 4))
+            [1, -1]
+            sage: list(K.selmer_group_iterator([K.prime_above(p) for p in [11,13]],2))
+            [1, -1, 13, -13, 11, -11, 143, -143]
+        """
+        KSgens, ords = self.selmer_group(S=S, m=m, proof=proof, orders=True)
+        one = self.one_element()
+        from sage.misc.all import cartesian_product_iterator
+        for ev in cartesian_product_iterator([range(o) for o in ords]):
+            yield prod([p**e for p,e in zip(KSgens, ev)], one)
 
     def composite_fields(self, other, names=None, both_maps=False, preserve_embedding=True):
         """
@@ -4168,7 +4273,7 @@ def composite_fields(self, other, names=None, both_maps=False, preserve_embeddin
                             self_to_F = self.hom([a])
                             other_to_F = other.hom([(~self.hom([F(a_in_F)]))(F(b_in_F))])
                             F = self
-                            k = sage.rings.infinity.Infinity
+                            k = Infinity
                             i -= 1
                         elif r.degree() == n:
                             other_to_F = other.hom([b])
@@ -7406,7 +7511,7 @@ def places(self, all_complex=False, prec=None):
             R = sage.rings.real_double.RDF
             C = sage.rings.complex_double.CDF
 
-        elif prec == sage.rings.infinity.Infinity:
+        elif prec == Infinity:
             R = sage.rings.all.AA
             C = sage.rings.all.QQbar
 
@@ -7641,7 +7746,11 @@ def relativize(self, alpha, names, structure=None):
                 raise ValueError("Co-domain of morphism must be self")
             L = alpha.domain()
             alpha = alpha(L.gen()) # relativize over phi's domain
-            f = L.defining_polynomial() # = alpha.minpoly()
+            if L is QQ:
+                from sage.rings.all import polygen
+                f = polygen(QQ)
+            else:
+                f = L.defining_polynomial() # = alpha.minpoly()
         else:
             # alpha must be an element coercible to self
             alpha = self(alpha)
@@ -7678,7 +7787,10 @@ def relativize(self, alpha, names, structure=None):
         if structure is None:
             from sage.rings.number_field.structure import RelativeFromAbsolute
             structure = RelativeFromAbsolute(self, alpha)
-        return L.extension(f, names[0], structure=structure)
+        if L is QQ:
+            return L.extension(f, names[0])
+        else:
+            return L.extension(f, names[0], structure=structure)
 
     # Synonyms so that terminology appropriate to relative number fields
     # can be applied to an absolute number field:
@@ -8342,14 +8454,17 @@ def _latex_(self):
 
     def _coerce_map_from_(self, K):
         r"""
-        The cyclotomic field `\Q(\zeta_n)` coerces into the cyclotomic field
-        `\Q(\zeta_m)` iff `n'|m`, where `n'` is the odd part of `n` if `4 \not
-        | n` and `n'=n` otherwise.
+        Return a coercion map from `K` to ``self``, or None.
+
+        The cyclotomic field `\QQ(\zeta_n)` coerces into the
+        cyclotomic field `\QQ(\zeta_m)` if and only if `n' \mid m`,
+        where `n'` is the odd part of `n` if `4 \nmid n` and `n' = n`
+        otherwise.
 
-        The morphism is consistant with the chosen embedding into `\CC`.
+        The morphism is consistent with the chosen embedding into `\CC`.
 
-        If `K` is not a cyclotomic field, the normal coercion rules for number
-        fields are used.
+        If `K` is not a cyclotomic field, the normal coercion rules
+        for number fields are used.
 
         EXAMPLES::
 
@@ -9833,12 +9948,12 @@ def refine_embedding(e, prec=None):
 
     # We first compute all the embeddings at the new precision:
     if sage.rings.real_mpfr.is_RealField(RC) or RC is RDF:
-        if prec==sage.rings.infinity.Infinity:
+        if prec == Infinity:
             elist = K.embeddings(sage.rings.qqbar.AA)
         else:
             elist = K.real_embeddings(prec)
     else:
-        if prec==sage.rings.infinity.Infinity:
+        if prec == Infinity:
             elist = K.embeddings(sage.rings.qqbar.QQbar)
         else:
             elist = K.complex_embeddings(prec)