Enhanced number field selmer group capabilities by adding an iterator

and including this functionality for QQ.  Also added a function to
compute the preimage of a number field element to a subfield modulo
powers.

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
@@ -3821,10 +3822,7 @@ def _S_class_group_quotient_matrix(self, S):
 
     def selmer_group(self, S, m, proof=True):
         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:
 
@@ -3836,13 +3834,29 @@ def selmer_group(self, S, m, proof=True):
 
         OUTPUT:
 
-        A list of generators of `K(S,m)`.
+        A list of generators of `K(S,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 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]
+
+        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
             sage: K.selmer_group([K.ideal(2, -a+1)], 2)
             [2, -1]
             sage: K.selmer_group([K.ideal(2, -a+1), K.ideal(3, a+1)], 2)
@@ -3855,9 +3869,14 @@ def selmer_group(self, S, m, proof=True):
             [2, a + 1]
             sage: K.selmer_group([K.ideal(2, -a+1), K.ideal(3, a+1), 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]
 
         TESTS::
 
@@ -3874,7 +3893,6 @@ def selmer_group(self, S, m, proof=True):
         """
         units, clgp_gens = self._S_class_group_and_units(tuple(S), proof=proof)
         gens = []
-        from sage.rings.infinity import Infinity
         for unit in units:
             order = unit.multiplicative_order()
             if order == Infinity or order.gcd(m) != 1:
@@ -3905,6 +3923,54 @@ def selmer_group(self, S, m, proof=True):
                 gens.append(self(J.gens_reduced()[0]))
         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, 2, 4, 8, -1, -2, -4, -8]
+            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 = self.selmer_group(S=S, m=m, proof=proof)
+        f = lambda o: m if o is Infinity else o.gcd(m)
+        orders = [f(a.multiplicative_order()) for a in KSgens]
+        one = self.one_element()
+        from sage.misc.all import cartesian_product_iterator
+        for ev in cartesian_product_iterator([range(o) for o in orders]):
+            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):
         """
         List of all possible composite number fields formed from self and
@@ -4158,7 +4224,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])
@@ -7396,7 +7462,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
 
@@ -7631,7 +7697,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)
@@ -7668,7 +7738,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:
@@ -8332,8 +8405,8 @@ 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
+        The cyclotomic field `\QQ(\zeta_n)` coerces into the cyclotomic field
+        `\QQ(\zeta_m)` iff `n'|m`, where `n'` is the odd part of `n` if `4 \not
         | n` and `n'=n` otherwise.
 
         The morphism is consistant with the chosen embedding into `\CC`.
@@ -9823,12 +9896,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)

src/sage/rings/number_field/number_field_element.pyx

@@ -1633,6 +1633,26 @@ cdef class NumberFieldElement(FieldElement):
         else:
             raise ValueError, "%s not a %s-th root in %s"%(self, n, self._parent)
 
+    def is_nth_power(self, n):
+        r"""
+        Return True if self is an nth power in the given number field.
+
+        EXAMPLES::
+
+            sage: K.<a> = NumberField(x^4-7)
+            sage: K(7).is_nth_power(2)
+            True
+            sage: K(7).is_nth_power(4)
+            True
+            sage: K(7).is_nth_power(8)
+            False
+            sage: K((a-3)^5).is_nth_power(5)
+            True
+
+        ALGORITHM: Use PARI to factor `x^n` - ``self`` in `K`.
+        """
+        return len(self.nth_root(n, all=True))>0
+
     def __pow__(base, exp, dummy):
         """
         EXAMPLES::
@@ -3479,7 +3499,90 @@ cdef class NumberFieldElement(FieldElement):
         """
         return P.residue_symbol(self,m,check)
 
+    def descend_mod_power(self, K=QQ, d=2):
+        r"""
+        Return a list of elements of the subfield equal to self mod d'th powers.
+
+        INPUT:
 
+        - ``K`` (number field, default QQ) - a subfield of the parent number field.
+        - ``d`` (positive integer, default 2) - an integer at least 2
+
+        OUTPUT:
+
+        A list, possibly empty, of elements of ``K`` equal to self
+        modulo ``d``'th powers, i.e. the preimages of self under the
+        map $K^*/(K^*)^2 \rightarrow L^*/(L^*)^2$ where $L$ is the
+        parent of self.  A ValueError is raised if K does not embed
+        into $L$.
+
+        ALGORITHM:
+
+        All preimages must lie in the Selmer group $K(S,d)$ for a
+        suitable finite set of primes $S$, which reduces the question
+        to a finite set of possibilities.  We may take $S$ to be the
+        set of primes which ramify in $L$ together with those for
+        which the valuation of self is not divisible by $d$.
+
+        EXAMPLES:
+
+        A relative example::
+
+            sage: Qi.<i> = QuadraticField(-1)
+            sage: K.<zeta> = CyclotomicField(8)
+            sage: f = Qi.embeddings(K)[0]
+            sage: a = f(2+3*i) * (2-zeta)^2
+            sage: a.descend_mod_power(Qi,2)
+            [-3*i - 2, 2*i - 3]
+
+        An absolute example::
+
+            sage: K.<zeta> = CyclotomicField(8)
+            sage: K(1).descend_mod_power(QQ,2)
+            [1, -1, 2, -2]
+            sage: a = 17*K.random_element()^2
+            sage: a.descend_mod_power(QQ,2)
+            [17, -17, 34, -34]
+        """
+        if not self:
+            raise ValueError("must be nonzero")
+        L = self.parent()
+        if K is L:
+            return [self]
+
+        from sage.sets.set import Set
+
+        if K is QQ: # simpler special case avoids relativizing
+            # First set of primes: those which ramify in L/K:
+            S1 = L.absolute_discriminant().prime_factors()
+            # Second set of primes: those where self has nonzero valuation mod d:
+            S2 = Set([p.norm().support()[0]
+                      for p in self.support()
+                      if self.valuation(p)%d !=0])
+            S = S1 + [p for p in S2 if not p in S1]
+            return [a for a in K.selmer_group_iterator(S,d)
+                    if (self/a).is_nth_power(d)]
+
+        embs = K.embeddings(L)
+        if len(embs)==0:
+            raise ValueError("K = %s does not embed into %s" % (K,L))
+        f = embs[0]
+        LK = L.relativize(f, names='b')
+        # Unfortunately the base field of LK is not K but an
+        # isomorphic field, and we must make sure to use the correct
+        # isomorphism!
+        KK = LK.base_field()
+        h = [h for h in KK.embeddings(K) if f(h(KK.gen()))==L(LK(KK.gen()))][0]
+
+        # First set of primes: those which ramify in L/K:
+        S1 = LK.relative_discriminant().prime_factors()
+        # Second set of primes: those where self has nonzero valuation mod d:
+        S2 = Set([p.relative_norm().prime_factors()[0]
+                  for p in LK(self).support()
+                  if LK(self).valuation(p)%d !=0])
+        S = S1 + [p for p in S2 if not p in S1]
+        candidates = [h(a) for a in K.selmer_group_iterator(S,d)]
+        return [a for a in candidates if (self/f(a)).is_nth_power(d)]
 
 cdef class NumberFieldElement_absolute(NumberFieldElement):