def duke_imamoglu_lift(self, f, f_k, precision, half_integral_weight = False) :
        """
        INPUT:
        
        - ``half_integral_weight``   -- If ``False`` we assume that `f` is the Fourier expansion of a
                                        Jacobi form. Otherwise we assume it is the Fourier expansion
                                        of a half integral weight elliptic modular form.
        """
        
        if half_integral_weight :
            coeff_index = lambda d : d
        else :
            coeff_index = lambda d : ((d + (-d % 4)) // 4, (-d) % 4)
        
        coeffs = dict()
        
        for t in precision.iter_positive_forms() :
            dt = (-1)**(precision.genus() // 2) * t.det()
            d = fundamental_discriminant(dt)
            eps = Integer(isqrt(dt / d))
    
            coeffs[t] = 0 
            for a in eps.divisors() :
                d_a = abs(d * (eps // a)**2)
                 
                coeffs[t] = coeffs[t] + a**(f_k - 1) * self._kohnen_phi(a, t) \
                                        * f[coeff_index(d_a)]

        return coeffs

 def __pow__(self, n):
     """
     Returns this species to the power n. This uses a binary
     exponentiation algorithm to perform the powering.
     
     EXAMPLES::
     
         sage: X = species.SingletonSpecies()
         sage: (X^2).generating_series().coefficients(4)
         [0, 0, 1, 0]
         sage: X^1 is X
         True
         sage: A = X^32
         sage: A.digraph()
         Multi-digraph on 6 vertices
     """
     from sage.rings.all import Integer
     import operator
     n = Integer(n)
     if n <= 0:
         raise ValueError, "only positive exponents are currently supported"
     digits = n.digits(2)
     squares = [self]
     for i in range(len(digits)-1):
         squares.append(squares[-1]*squares[-1])
     return reduce(operator.add, (s for i,s in zip(digits, squares) if i != 0))

    def _element_constructor_(self, x) :
        """
        TESTS::
            sage: from psage.modform.fourier_expansion_framework.monoidpowerseries.monoidpowerseries_basicmonoids import *
            sage: from psage.modform.fourier_expansion_framework.monoidpowerseries.monoidpowerseries_ring import MonoidPowerSeriesRing_generic
            sage: mps = MonoidPowerSeriesRing_generic(QQ, NNMonoid(False))
            sage: h = mps(1) # indirect doctest
            sage: h = mps(mps.monoid().zero_element())
            sage: h = mps.zero_element()
            sage: K.<rho> = CyclotomicField(6)
            sage: mps = MonoidPowerSeriesRing_generic(K, NNMonoid(False))
            sage: h = mps(rho)
            sage: h = mps(1)
        """
        if isinstance(x, int) :
            x = Integer(x)
            
        if isinstance(x, Element) :
            P = x.parent()

            if P is self.coefficient_domain() :
                return self._element_class( self, {self.monoid().zero_element(): x},
                                                self.monoid().filter_all() )
            elif self.coefficient_domain().has_coerce_map_from(P) :
                return self._element_class( self, {self.monoid().zero_element(): self.coefficient_domain()(x)},
                                                self.monoid().filter_all() )
            elif P is self.monoid() :
                return self._element_class( self, {x: self.base_ring().one_element},
                                                self.monoid().filter_all() )
                
        return MonoidPowerSeriesAmbient_abstract._element_constructor_(self, x)

    def __pow__(self, n):
        r"""
        Returns this species to the power `n`.

        This uses a binary exponentiation algorithm to perform the
        powering.

        EXAMPLES::

            sage: One = species.EmptySetSpecies()
            sage: X = species.SingletonSpecies()
            sage: X^2
            Product of (Singleton species) and (Singleton species)
            sage: X^5
            Product of (Singleton species) and (Product of (Product of
            (Singleton species) and (Singleton species)) and (Product
            of (Singleton species) and (Singleton species)))

            sage: (X^2).generating_series().coefficients(4)
            [0, 0, 1, 0]
            sage: (X^3).generating_series().coefficients(4)
            [0, 0, 0, 1]
            sage: ((One+X)^3).generating_series().coefficients(4)
            [1, 3, 3, 1]
            sage: ((One+X)^7).generating_series().coefficients(8)
            [1, 7, 21, 35, 35, 21, 7, 1]

            sage: x = QQ[['x']].gen()
            sage: coeffs = ((1+x+x+x**2)**25+O(x**10)).padded_list()
            sage: T = ((One+X+X+X^2)^25)
            sage: T.generating_series().coefficients(10) == coeffs
            True
            sage: X^1 is X
            True
            sage: A = X^32
            sage: A.digraph()
            Multi-digraph on 6 vertices

        TESTS::

            sage: X**(-1)
            Traceback (most recent call last):
            ...
            ValueError: only positive exponents are currently supported
        """
        from sage.rings.all import Integer
        import operator
        n = Integer(n)
        if n <= 0:
            raise ValueError("only positive exponents are currently supported")
        digits = n.digits(2)
        squares = [self]
        for i in range(len(digits) - 1):
            squares.append(squares[-1] * squares[-1])
        return reduce(operator.mul, (s for i, s in zip(digits, squares)
                                     if i != 0))

    def __init__(self, abvar, p):
        """
        Create a `p`-adic `L`-series.

        EXAMPLES::

            sage: J0(37)[0].padic_lseries(389)
            389-adic L-series attached to Simple abelian subvariety 37a(1,37) of dimension 1 of J0(37)
        """
        Lseries.__init__(self, abvar)
        p = Integer(p)
        if not p.is_prime():
            raise ValueError("p (=%s) must be prime"%p)
        self.__p = p

def validate_label(label):
    parts = label.split('.')
    if len(parts) != 3:
        raise ValueError("it must be of the form g.q.iso, with g a dimension and q a prime power")
    g, q, iso = parts
    try:
        g = int(g)
    except ValueError:
        raise ValueError("it must be of the form g.q.iso, where g is an integer")
    try:
        q = Integer(q)
        if not q.is_prime_power(): raise ValueError
    except ValueError:
        raise ValueError("it must be of the form g.q.iso, where g is a prime power")
    coeffs = iso.split("_")
    if len(coeffs) != g:
        raise ValueError("the final part must be of the form c1_c2_..._cg, with g=%s components"%(g))
    if not all(c.isalpha() and c==c.lower() for c in coeffs):
        raise ValueError("the final part must be of the form c1_c2_..._cg, with each ci consisting of lower case letters")

    def apply_TT(self, j):
        """
        Apply the matrix `TT=[-1,-1,0,1]` to the `j`-th Manin symbol.

        INPUT:

        - ``j`` - (int) a symbol index

        OUTPUT: see documentation for apply()

        EXAMPLE::

            sage: from sage.modular.modsym.manin_symbol_list import ManinSymbolList_gamma0
            sage: m = ManinSymbolList_gamma0(5,8)
            sage: m.apply_TT(4)
            [(38, 1)]
            sage: [m.apply_TT(i) for i in xrange(10)]
            [[(37, 1)],
            [(41, 1)],
            [(39, 1)],
            [(40, 1)],
            [(38, 1)],
            [(36, 1)],
            [(31, -1), (37, 1)],
            [(35, -1), (41, 1)],
            [(33, -1), (39, 1)],
            [(34, -1), (40, 1)]]
        """
        k = self._weight
        i, u, v = self._symbol_list[j]
        u, v = self.__syms.normalize(-u-v,u)
        if (k-2-i) % 2 == 0:
            s = 1
        else:
            s = -1
        z = []
        a = Integer(i)
        for j in range(i+1):
            m = self.index((k-2-i+j, u, v))
            z.append((m, s * a.binomial(j)))
            s *= -1
        return z

    def apply_T(self, j):
        """
        Apply the matrix `T=[0,1,-1,-1]` to the `j`-th Manin symbol.

        INPUT:

        - ``j`` - (int) a symbol index

        OUTPUT: see documentation for apply()

        EXAMPLE::

            sage: from sage.modular.modsym.manin_symbol_list import ManinSymbolList_gamma0
            sage: m = ManinSymbolList_gamma0(5,8)
            sage: m.apply_T(4)
            [(3, 1), (9, -6), (15, 15), (21, -20), (27, 15), (33, -6), (39, 1)]
            sage: [m.apply_T(i) for i in xrange(10)]
            [[(5, 1), (11, -6), (17, 15), (23, -20), (29, 15), (35, -6), (41, 1)],
            [(0, 1), (6, -6), (12, 15), (18, -20), (24, 15), (30, -6), (36, 1)],
            [(4, 1), (10, -6), (16, 15), (22, -20), (28, 15), (34, -6), (40, 1)],
            [(2, 1), (8, -6), (14, 15), (20, -20), (26, 15), (32, -6), (38, 1)],
            [(3, 1), (9, -6), (15, 15), (21, -20), (27, 15), (33, -6), (39, 1)],
            [(1, 1), (7, -6), (13, 15), (19, -20), (25, 15), (31, -6), (37, 1)],
            [(5, 1), (11, -5), (17, 10), (23, -10), (29, 5), (35, -1)],
            [(0, 1), (6, -5), (12, 10), (18, -10), (24, 5), (30, -1)],
            [(4, 1), (10, -5), (16, 10), (22, -10), (28, 5), (34, -1)],
            [(2, 1), (8, -5), (14, 10), (20, -10), (26, 5), (32, -1)]]
        """
        k = self._weight
        i, u, v = self._symbol_list[j]
        u, v = self.__syms.normalize(v,-u-v)
        if (k-2) % 2 == 0:
            s = 1
        else:
            s = -1
        z = []
        a = Integer(k-2-i)
        for j in range(k-2-i+1):
            m = self.index((j, u, v))
            z.append((m, s * a.binomial(j)))
            s *= -1
        return z

def xi_degrees(n,p=2, reverse=True):
    r"""
    Decreasing list of degrees of the xi_i's, starting in degree n.

    INPUT:

    - `n` - integer
    - `p` - prime number, optional (default 2)
    - ``reverse`` - bool, optional (default True)

    OUTPUT: ``list`` - list of integers

    When `p=2`: decreasing list of the degrees of the `\xi_i`'s with
    degree at most n.

    At odd primes: decreasing list of these degrees, each divided by
    `2(p-1)`.

    If ``reverse`` is False, then return an increasing list rather
    than a decreasing one.

    EXAMPLES::

        sage: sage.algebras.steenrod.steenrod_algebra_bases.xi_degrees(17)
        [15, 7, 3, 1]
        sage: sage.algebras.steenrod.steenrod_algebra_bases.xi_degrees(17, reverse=False)
        [1, 3, 7, 15]
        sage: sage.algebras.steenrod.steenrod_algebra_bases.xi_degrees(17,p=3)
        [13, 4, 1]
        sage: sage.algebras.steenrod.steenrod_algebra_bases.xi_degrees(400,p=17)
        [307, 18, 1]
    """
    from sage.rings.all import Integer
    if n <= 0: return []
    N = Integer(n*(p-1) + 1)
    l = [int((p**d-1)/(p-1)) for d in range(1,N.exact_log(p)+1)]
    if not reverse:
        return l
    l.reverse()
    return l

    def ap(self, p):
        """
        Return a list of the eigenvalues of the Hecke operator `T_p`
        on all the computed eigenforms.  The eigenvalues match up
        between one prime and the next.

        INPUT:

        - ``p`` - integer, a prime number

        OUTPUT:

        - ``list`` - a list of double precision complex numbers

        EXAMPLES::

            sage: n = numerical_eigenforms(11,4)
            sage: n.ap(2) # random order
            [9.0, 9.0, 2.73205080757, -0.732050807569]
            sage: n.ap(3) # random order
            [28.0, 28.0, -7.92820323028, 5.92820323028]
            sage: m = n.modular_symbols()
            sage: x = polygen(QQ, 'x')
            sage: m.T(2).charpoly('x').factor()
            (x - 9)^2 * (x^2 - 2*x - 2)
            sage: m.T(3).charpoly('x').factor()
            (x - 28)^2 * (x^2 + 2*x - 47)
        """
        p = Integer(p)
        if not p.is_prime():
            raise ValueError("p must be a prime")
        try:
            return self._ap[p]
        except AttributeError:
            self._ap = {}
        except KeyError:
            pass
        a = Sequence(self.eigenvalues([p])[0], immutable=True)
        self._ap[p] = a
        return a

def hilbert_class_polynomial(D, algorithm=None):
    r"""
    Returns the Hilbert class polynomial for discriminant `D`.

    INPUT:

    - ``D`` (int) -- a negative integer congruent to 0 or 1 modulo 4.

    - ``algorithm`` (string, default None).

    OUTPUT:

    (integer polynomial) The Hilbert class polynomial for the
    discriminant `D`.

    ALGORITHM:

    - If ``algorithm`` = "arb" (default): Use Arb's implementation which uses complex interval arithmetic.

    - If ``algorithm`` = "sage": Use complex approximations to the roots.

    - If ``algorithm`` = "magma": Call the appropriate Magma function (if available).

    AUTHORS:

    - Sage implementation originally by Eduardo Ocampo Alvarez and
      AndreyTimofeev

    - Sage implementation corrected by John Cremona (using corrected precision bounds from Andreas Enge)

    - Magma implementation by David Kohel

    EXAMPLES::

        sage: hilbert_class_polynomial(-4)
        x - 1728
        sage: hilbert_class_polynomial(-7)
        x + 3375
        sage: hilbert_class_polynomial(-23)
        x^3 + 3491750*x^2 - 5151296875*x + 12771880859375
        sage: hilbert_class_polynomial(-37*4)
        x^2 - 39660183801072000*x - 7898242515936467904000000
        sage: hilbert_class_polynomial(-37*4, algorithm="magma") # optional - magma
        x^2 - 39660183801072000*x - 7898242515936467904000000
        sage: hilbert_class_polynomial(-163)
        x + 262537412640768000
        sage: hilbert_class_polynomial(-163, algorithm="sage")
        x + 262537412640768000
        sage: hilbert_class_polynomial(-163, algorithm="magma") # optional - magma
        x + 262537412640768000

    TESTS::

        sage: all([hilbert_class_polynomial(d, algorithm="arb") == \
        ....:      hilbert_class_polynomial(d, algorithm="sage") \
        ....:        for d in range(-1,-100,-1) if d%4 in [0,1]])
        True

    """
    if algorithm is None:
        algorithm = "arb"

    D = Integer(D)
    if D >= 0:
        raise ValueError("D (=%s) must be negative"%D)
    if not (D%4 in [0,1]):
         raise ValueError("D (=%s) must be a discriminant"%D)

    if algorithm == "arb":
        import sage.libs.arb.arith
        return sage.libs.arb.arith.hilbert_class_polynomial(D)

    if algorithm == "magma":
        magma.eval("R<x> := PolynomialRing(IntegerRing())")
        f = str(magma.eval("HilbertClassPolynomial(%s)"%D))
        return IntegerRing()['x'](f)

    if algorithm != "sage":
        raise ValueError("%s is not a valid algorithm"%algorithm)

    from sage.quadratic_forms.binary_qf import BinaryQF_reduced_representatives
    from sage.rings.all import RR, ZZ, ComplexField
    from sage.functions.all import elliptic_j

    # get all primitive reduced quadratic forms, (necessary to exclude
    # imprimitive forms when D is not a fundamental discriminant):

    rqf = BinaryQF_reduced_representatives(D, primitive_only=True)

    # compute needed precision
    #
    # NB: [http://arxiv.org/abs/0802.0979v1], quoting Enge (2006), is
    # incorrect.  Enge writes (2009-04-20 email to John Cremona) "The
    # source is my paper on class polynomials
    # [http://hal.inria.fr/inria-00001040] It was pointed out to me by
    # the referee after ANTS that the constant given there was
    # wrong. The final version contains a corrected constant on p.7
    # which is consistent with your example. It says:

    # "The logarithm of the absolute value of the coefficient in front
#.........这里部分代码省略.........

def check_prime(K,P):
    r"""
    Function to check that `P` determines a prime of `K`, and return that ideal.

    INPUT:

    - ``K`` -- a number field (including `\QQ`).

    - ``P`` -- an element of ``K`` or a (fractional) ideal of ``K``.

    OUTPUT:

    - If ``K`` is `\QQ`: the prime integer equal to or which generates `P`.

    - If ``K`` is not `\QQ`: the prime ideal equal to or generated by `P`.

    .. note::

       If `P` is not a prime and does not generate a prime, a TypeError is raised.

    EXAMPLES::

        sage: from sage.schemes.elliptic_curves.ell_local_data import check_prime
        sage: check_prime(QQ,3)
        3
        sage: check_prime(QQ,QQ(3))
        3
        sage: check_prime(QQ,ZZ.ideal(31))
        31
        sage: K.<a>=NumberField(x^2-5)
        sage: check_prime(K,a)
        Fractional ideal (a)
        sage: check_prime(K,a+1)
        Fractional ideal (a + 1)
        sage: [check_prime(K,P) for P in K.primes_above(31)]
        [Fractional ideal (5/2*a + 1/2), Fractional ideal (5/2*a - 1/2)]
        sage: L.<b> = NumberField(x^2+3)
        sage: check_prime(K, L.ideal(5))
        Traceback (most recent call last):
        ..
        TypeError: The ideal Fractional ideal (5) is not a prime ideal of Number Field in a with defining polynomial x^2 - 5
        sage: check_prime(K, L.ideal(b))
        Traceback (most recent call last):
        TypeError: No compatible natural embeddings found for Number Field in a with defining polynomial x^2 - 5 and Number Field in b with defining polynomial x^2 + 3
    """
    if K is QQ:
        if P in ZZ or isinstance(P, integer_types + (Integer,)):
            P = Integer(P)
            if P.is_prime():
                return P
            else:
                raise TypeError("The element %s is not prime" % (P,) )
        elif P in QQ:
            raise TypeError("The element %s is not prime" % (P,) )
        elif is_Ideal(P) and P.base_ring() is ZZ:
            if P.is_prime():
                return P.gen()
            else:
                raise TypeError("The ideal %s is not a prime ideal of %s" % (P, ZZ))
        else:
            raise TypeError("%s is neither an element of QQ or an ideal of %s" % (P, ZZ))

    if not is_NumberField(K):
        raise TypeError("%s is not a number field" % (K,) )

    if is_NumberFieldFractionalIdeal(P) or P in K:
        # if P is an ideal, making sure it is an fractional ideal of K
        P = K.fractional_ideal(P)
        if P.is_prime():
            return P
        else:
            raise TypeError("The ideal %s is not a prime ideal of %s" % (P, K))

    raise TypeError("%s is not a valid prime of %s" % (P, K))

    def galois_action(self, t, N):
        r"""
        Suppose this cusp is `\alpha`, `G` a congruence subgroup of level `N`
        and `\sigma` is the automorphism in the Galois group of
        `\QQ(\zeta_N)/\QQ` that sends `\zeta_N` to `\zeta_N^t`. Then this
        function computes a cusp `\beta` such that `\sigma([\alpha]) = [\beta]`,
        where `[\alpha]` is the equivalence class of `\alpha` modulo `G`.

        This code only needs as input the level and not the group since the
        action of galois for a congruence group `G` of level `N` is compatible
        with the action of the full congruence group `\Gamma(N)`.


        INPUT:

           - `t` -- integer that is coprime to N

           - `N` -- positive integer (level)

        OUTPUT:

           - a cusp


        .. WARNING::

            In some cases `N` must fit in a long long, i.e., there
            are cases where this algorithm isn't fully implemented.

        .. NOTE::

            Modular curves can have multiple non-isomorphic models over `\QQ`.
            The action of galois depends on such a model. The model over `\QQ`
            of `X(G)` used here is the model where the function field
            `\QQ(X(G))` is given by the functions whose fourier expansion at
            `\infty` have their coefficients in `\QQ`. For `X(N):=X(\Gamma(N))`
            the corresponding moduli interpretation over `\ZZ[1/N]` is that
            `X(N)` parametrizes pairs `(E,a)` where `E` is a (generalized)
            elliptic curve and `a: \ZZ / N\ZZ \times \mu_N \to E` is a closed
            immersion such that the weil pairing of `a(1,1)` and `a(0,\zeta_N)`
            is `\zeta_N`. In this parameterisation the point `z \in H`
            corresponds to the pair `(E_z,a_z)` with `E_z=\CC/(z \ZZ+\ZZ)` and
            `a_z: \ZZ / N\ZZ \times \mu_N \to E` given by `a_z(1,1) = z/N` and
            `a_z(0,\zeta_N) = 1/N`.
            Similarly `X_1(N):=X(\Gamma_1(N))` parametrizes pairs `(E,a)` where
            `a: \mu_N \to E` is a closed immersion.

        EXAMPLES::

            sage: Cusp(1/10).galois_action(3, 50)
            1/170
            sage: Cusp(oo).galois_action(3, 50)
            Infinity
            sage: c=Cusp(0).galois_action(3, 50); c
            50/67
            sage: Gamma0(50).reduce_cusp(c)
            0

        Here we compute the permutations of the action for t=3 on cusps for
        Gamma0(50). ::

            sage: N = 50; t=3; G = Gamma0(N); C = G.cusps()
            sage: cl = lambda z: exists(C, lambda y:y.is_gamma0_equiv(z, N))[1]
            sage: for i in range(5): print i, t^i, [cl(alpha.galois_action(t^i,N)) for alpha in C]
            0 1 [0, 1/25, 1/10, 1/5, 3/10, 2/5, 1/2, 3/5, 7/10, 4/5, 9/10, Infinity]
            1 3 [0, 1/25, 7/10, 2/5, 1/10, 4/5, 1/2, 1/5, 9/10, 3/5, 3/10, Infinity]
            2 9 [0, 1/25, 9/10, 4/5, 7/10, 3/5, 1/2, 2/5, 3/10, 1/5, 1/10, Infinity]
            3 27 [0, 1/25, 3/10, 3/5, 9/10, 1/5, 1/2, 4/5, 1/10, 2/5, 7/10, Infinity]
            4 81 [0, 1/25, 1/10, 1/5, 3/10, 2/5, 1/2, 3/5, 7/10, 4/5, 9/10, Infinity]

        TESTS:

        Here we check that the galois action is indeed a permutation on the
        cusps of Gamma1(48) and check that :trac:`13253` is fixed. ::

            sage: G=Gamma1(48)
            sage: C=G.cusps()
            sage: for i in Integers(48).unit_gens():
            ...     C_permuted = [G.reduce_cusp(c.galois_action(i,48)) for c in C]
            ...     assert len(set(C_permuted))==len(C)

        We test that Gamma1(19) has 9 rational cusps and check that :trac:`8998`
        is fixed. ::

            sage: G = Gamma1(19)
            sage: [c for c in G.cusps() if c.galois_action(2,19).is_gamma1_equiv(c,19)[0]]
            [2/19, 3/19, 4/19, 5/19, 6/19, 7/19, 8/19, 9/19, Infinity]


        REFERENCES:

            - Section 1.3 of Glenn Stevens, "Arithmetic on Modular Curves"

            - There is a long comment about our algorithm in the source code for this function.

        AUTHORS:

            - William Stein, 2009-04-18

        """
#.........这里部分代码省略.........

def discriminants_with_bounded_class_number(hmax, B=None, proof=None):
    """
    Return dictionary with keys class numbers `h\le hmax` and values the
    list of all pairs `(D, f)`, with `D<0` a fundamental discriminant such
    that `Df^2` has class number `h`.  If the optional bound `B` is given,
    return only those pairs with fundamental `|D| \le B`, though `f` can
    still be arbitrarily large.

    INPUT:

    - ``hmax`` -- integer
    - `B` -- integer or None; if None returns all pairs
    - ``proof`` -- this code calls the PARI function ``qfbclassno``, so it
      could give wrong answers when ``proof``==``False``.  The default is
      whatever ``proof.number_field()`` is.  If ``proof==False`` and `B` is
      ``None``, at least the number of discriminants is correct, since it
      is double checked with Watkins's table.

    OUTPUT:

    - dictionary

    In case `B` is not given, we use Mark Watkins's: "Class numbers of
    imaginary quadratic fields" to compute a `B` that captures all `h`
    up to `hmax` (only available for `hmax\le100`).

    EXAMPLES::

        sage: v = sage.schemes.elliptic_curves.cm.discriminants_with_bounded_class_number(3)
        sage: list(v)
        [1, 2, 3]
        sage: v[1]
        [(-3, 3), (-3, 2), (-3, 1), (-4, 2), (-4, 1), (-7, 2), (-7, 1), (-8, 1), (-11, 1), (-19, 1), (-43, 1), (-67, 1), (-163, 1)]
        sage: v[2]
        [(-3, 7), (-3, 5), (-3, 4), (-4, 5), (-4, 4), (-4, 3), (-7, 4), (-8, 3), (-8, 2), (-11, 3), (-15, 2), (-15, 1), (-20, 1), (-24, 1), (-35, 1), (-40, 1), (-51, 1), (-52, 1), (-88, 1), (-91, 1), (-115, 1), (-123, 1), (-148, 1), (-187, 1), (-232, 1), (-235, 1), (-267, 1), (-403, 1), (-427, 1)]
        sage: v[3]
        [(-3, 9), (-3, 6), (-11, 2), (-19, 2), (-23, 2), (-23, 1), (-31, 2), (-31, 1), (-43, 2), (-59, 1), (-67, 2), (-83, 1), (-107, 1), (-139, 1), (-163, 2), (-211, 1), (-283, 1), (-307, 1), (-331, 1), (-379, 1), (-499, 1), (-547, 1), (-643, 1), (-883, 1), (-907, 1)]
        sage: v = sage.schemes.elliptic_curves.cm.discriminants_with_bounded_class_number(8, proof=False)
        sage: [len(v[h]) for h in v]
        [13, 29, 25, 84, 29, 101, 38, 208]

    Find all class numbers for discriminant up to 50::

        sage: sage.schemes.elliptic_curves.cm.discriminants_with_bounded_class_number(hmax=5, B=50)
        {1: [(-3, 3), (-3, 2), (-3, 1), (-4, 2), (-4, 1), (-7, 2), (-7, 1), (-8, 1), (-11, 1), (-19, 1), (-43, 1)], 2: [(-3, 7), (-3, 5), (-3, 4), (-4, 5), (-4, 4), (-4, 3), (-7, 4), (-8, 3), (-8, 2), (-11, 3), (-15, 2), (-15, 1), (-20, 1), (-24, 1), (-35, 1), (-40, 1)], 3: [(-3, 9), (-3, 6), (-11, 2), (-19, 2), (-23, 2), (-23, 1), (-31, 2), (-31, 1), (-43, 2)], 4: [(-3, 13), (-3, 11), (-3, 8), (-4, 10), (-4, 8), (-4, 7), (-4, 6), (-7, 8), (-7, 6), (-7, 3), (-8, 6), (-8, 4), (-11, 5), (-15, 4), (-19, 5), (-19, 3), (-20, 3), (-20, 2), (-24, 2), (-35, 3), (-39, 2), (-39, 1), (-40, 2), (-43, 3)], 5: [(-47, 2), (-47, 1)]}
    """
    # imports that are needed only for this function
    from sage.structure.proof.proof import get_flag
    import math
    from sage.misc.functional import round

    # deal with input defaults and type checking
    proof = get_flag(proof, 'number_field')
    hmax = Integer(hmax)

    # T stores the output
    T = {}

    # Easy case -- instead of giving error, give meaningful output
    if hmax < 1:
        return T

    if B is None:
        # Determine how far we have to go by applying Watkins's theorem.
        v = [largest_fundamental_disc_with_class_number(h) for h in range(1, hmax+1)]
        B = max([b for b,_ in v])
        fund_count = [0] + [cnt for _,cnt in v]
    else:
        # Nothing to do -- set to None so we can use this later to know not
        # to do a double check about how many we find.
        fund_count = None
        B = Integer(B)

    if B <= 2:
        # This is an easy special case, since there are no fundamental discriminants
        # this small.
        return T

    # This lower bound gets used in an inner loop below.
    from math import log
    def lb(f):
        """Lower bound on euler_phi."""
        # 1.79 > e^gamma = 1.7810724...
        if f <= 1: return 0  # don't do log(log(1)) = log(0)
        return f/(1.79*log(log(f)) + 3.0/log(log(f)))

    for D in range(-B, -2):
        D = Integer(D)
        if is_fundamental_discriminant(D):
            h_D = D.class_number(proof)
            # For each fundamental discriminant D, loop through the f's such
            # that h(D*f^2) could possibly be <= hmax.  As explained to me by Cremona,
            # we have h(D*f^2) >= (1/c)*h(D)*phi_D(f) >= (1/c)*h(D)*euler_phi(f), where
            # phi_D(f) is like euler_phi(f) but the factor (1-1/p) is replaced
            # by a factor of (1-kr(D,p)*p), where kr(D/p) is the Kronecker symbol.
            # The factor c is 1 unless D=-4 and f>1 (when c=2) or D=-3 and f>1 (when c=3).
            # Since (1-1/p) <= 1 and (1-1/p) <= (1+1/p), we see that
            #     euler_phi(f) <= phi_D(f).
            #
            # We have the following analytic lower bound on euler_phi:
#.........这里部分代码省略.........

    def an_numerical(self, prec = None,
                         use_database=True, proof=None):
        r"""
        Return the numerical analytic order of `Sha`, which is
        a floating point number in all cases.

        INPUT:

        - ``prec`` - integer (default: 53) bits precision -- used
          for the L-series computation, period,  regulator, etc.
        - ``use_database`` - whether the rank and generators should
          be looked up in the database if possible. Default is ``True``
        - ``proof`` - bool or ``None`` (default: ``None``, see proof.[tab] or
          sage.structure.proof) proof option passed
          onto regulator and rank computation.

        .. note::

            See also the :meth:`an` command, which will return a
            provably correct integer when the rank is 0 or 1.

        .. WARNING::

            If the curve's generators are not known, computing
            them may be very time-consuming.  Also, computation of the
            L-series derivative will be time-consuming for large rank and
            large conductor, and the computation time for this may
            increase substantially at greater precision.  However, use of
            very low precision less than about 10 can cause the underlying
            PARI library functions to fail.

        EXAMPLES::

            sage: EllipticCurve('11a').sha().an_numerical()
            1.00000000000000
            sage: EllipticCurve('37a').sha().an_numerical()
            1.00000000000000
            sage: EllipticCurve('389a').sha().an_numerical()
            1.00000000000000
            sage: EllipticCurve('66b3').sha().an_numerical()
            4.00000000000000
            sage: EllipticCurve('5077a').sha().an_numerical()
            1.00000000000000

        A rank 4 curve::

            sage: EllipticCurve([1, -1, 0, -79, 289]).sha().an_numerical()  # long time (3s on sage.math, 2011)
            1.00000000000000

        A rank 5 curve::

            sage: EllipticCurve([0, 0, 1, -79, 342]).sha().an_numerical(prec=10, proof=False)  # long time (22s on sage.math, 2011)
            1.0

        See :trac:`1115`::

            sage: sha=EllipticCurve('37a1').sha()
            sage: [sha.an_numerical(prec) for prec in xrange(40,100,10)]  # long time (3s on sage.math, 2013)
            [1.0000000000,
            1.0000000000000,
            1.0000000000000000,
            1.0000000000000000000,
            1.0000000000000000000000,
            1.0000000000000000000000000]
        """
        if prec is None:
            prec = RealField().precision()
        RR = RealField(prec)
        prec2 = prec+2
        RR2 = RealField(prec2)
        try:
            an = self.__an_numerical
            if an.parent().precision() >= prec:
                return RR(an)
            else: # cached precision too low
                pass
        except AttributeError:
            pass
        # it's critical to switch to the minimal model.
        E = self.Emin
        r = Integer(E.rank(use_database=use_database, proof=proof))
        L = E.lseries().dokchitser(prec=prec2)
        Lr= RR2(L.derivative(1,r))  # L.derivative() returns a Complex
        Om = RR2(E.period_lattice().omega(prec2))
        Reg = E.regulator(use_database=use_database, proof=proof, precision=prec2)
        T = E.torsion_order()
        cp = E.tamagawa_product()
        Sha = RR((Lr*T*T)/(r.factorial()*Om*cp*Reg))
        self.__an_numerical = Sha
        return Sha

#.........这里部分代码省略.........
    factors are small, this can be efficient even if the modulus itself
    is large::

        sage: sorted(solve_mod([x^2 == 41], 10^20))
        [(4538602480526452429,), (11445932736758703821,), (38554067263241296179,),
        (45461397519473547571,), (54538602480526452429,), (61445932736758703821,),
        (88554067263241296179,), (95461397519473547571,)]

    We solve a simple equation modulo 2::

        sage: x,y = var('x,y')
        sage: solve_mod([x == y], 2)
        [(0, 0), (1, 1)]

    .. warning::

       The current implementation splits the modulus into prime
       powers, then naively enumerates all possible solutions
       (starting modulo primes and then working up through prime
       powers), and finally combines the solution using the Chinese
       Remainder Theorem.  The interface is good, but the algorithm is
       very inefficient if the modulus has some larger prime factors! Sage
       *does* have the ability to do something much faster in certain
       cases at least by using Groebner basis, linear algebra
       techniques, etc. But for a lot of toy problems this function as
       is might be useful. At least it establishes an interface.


    TESTS:

    Make sure that we short-circuit in at least some cases::

        sage: solve_mod([2*x==1], 2*next_prime(10^50))
        []

    Try multi-equation cases::

        sage: x, y, z = var("x y z")
        sage: solve_mod([2*x^2 + x*y, -x*y+2*y^2+x-2*y, -2*x^2+2*x*y-y^2-x-y], 12)
        [(0, 0), (4, 4), (0, 3), (4, 7)]
        sage: eqs = [-y^2+z^2, -x^2+y^2-3*z^2-z-1, -y*z-z^2-x-y+2, -x^2-12*z^2-y+z]
        sage: solve_mod(eqs, 11)
        [(8, 5, 6)]

    Confirm that modulus 1 now behaves as it should::

        sage: x, y = var("x y")
        sage: solve_mod([x==1], 1)
        [(0,)]
        sage: solve_mod([2*x^2+x*y, -x*y+2*y^2+x-2*y, -2*x^2+2*x*y-y^2-x-y], 1)
        [(0, 0)]


    """
    from sage.rings.all import Integer, Integers, crt_basis
    from sage.symbolic.expression import is_Expression
    from sage.misc.all import cartesian_product_iterator
    from sage.modules.all import vector
    from sage.matrix.all import matrix

    if not isinstance(eqns, (list, tuple)):
        eqns = [eqns]
    eqns = [eq if is_Expression(eq) else (eq.lhs()-eq.rhs()) for eq in eqns]
    modulus = Integer(modulus)
    if modulus < 1:
        raise ValueError("the modulus must be a positive integer")
    vars = list(set(sum([list(e.variables()) for e in eqns], [])))
    vars.sort(key=repr)

    if modulus == 1: # degenerate case
        ans = [tuple(Integers(1)(0) for v in vars)]
        return ans

    factors = modulus.factor()
    crt_basis = vector(Integers(modulus), crt_basis([p**i for p,i in factors]))
    solutions = []

    has_solution = True
    for p,i in factors:
        solution =_solve_mod_prime_power(eqns, p, i, vars)
        if len(solution) > 0:
            solutions.append(solution)
        else:
            has_solution = False
            break


    ans = []
    if has_solution:
        for solution in cartesian_product_iterator(solutions):
            solution_mat = matrix(Integers(modulus), solution)
            ans.append(tuple(c.dot_product(crt_basis) for c in solution_mat.columns()))

    # if solution_dict == True:
    # Relaxed form suggested by Mike Hansen (#8553):
    if solution_dict:
        sol_dict = [dict(zip(vars, solution)) for solution in ans]
        return sol_dict
    else:
        return ans

def check_prime(K,P):
    r"""
    Function to check that `P` determines a prime of `K`, and return that ideal.

    INPUT:

    - ``K`` -- a number field (including `\QQ`).

    - ``P`` -- an element of ``K`` or a (fractional) ideal of ``K``.

    OUTPUT:

    - If ``K`` is `\QQ`: the prime integer equal to or which generates `P`.

    - If ``K`` is not `\QQ`: the prime ideal equal to or generated by `P`.

    .. note::

       If `P` is not a prime and does not generate a prime, a TypeError is raised.

    EXAMPLES::

        sage: from sage.schemes.elliptic_curves.ell_local_data import check_prime
        sage: check_prime(QQ,3)
        3
        sage: check_prime(QQ,ZZ.ideal(31))
        31
        sage: K.<a>=NumberField(x^2-5)
        sage: check_prime(K,a)
        Fractional ideal (a)
        sage: check_prime(K,a+1)
        Fractional ideal (a + 1)
        sage: [check_prime(K,P) for P in K.primes_above(31)]
        [Fractional ideal (5/2*a + 1/2), Fractional ideal (5/2*a - 1/2)]
    """
    if K is QQ:
        if isinstance(P, (int,long,Integer)):
            P = Integer(P)
            if P.is_prime():
                return P
            else:
                raise TypeError("%s is not prime"%P)
        else:
            if is_Ideal(P) and P.base_ring() is ZZ and P.is_prime():
                return P.gen()
        raise TypeError("%s is not a prime ideal of %s"%(P,ZZ))

    if not is_NumberField(K):
        raise TypeError("%s is not a number field"%K)

    if is_NumberFieldFractionalIdeal(P):
        if P.is_prime():
            return P
        else:
            raise TypeError("%s is not a prime ideal of %s"%(P,K))

    if is_NumberFieldElement(P):
        if P in K:
            P = K.ideal(P)
        else:
            raise TypeError("%s is not an element of %s"%(P,K))
        if P.is_prime():
            return P
        else:
            raise TypeError("%s is not a prime ideal of %s"%(P,K))

    raise TypeError("%s is not a valid prime of %s"%(P,K))

def victor_miller_basis(k, prec=10, cusp_only=False, var='q'):
    r"""
    Compute and return the Victor Miller basis for modular forms of
    weight `k` and level 1 to precision `O(q^{prec})`.  If
    ``cusp_only`` is True, return only a basis for the cuspidal
    subspace.

    INPUT:

    - ``k`` -- an integer

    - ``prec`` -- (default: 10) a positive integer

    - ``cusp_only`` -- bool (default: False)

    - ``var`` -- string (default: 'q')

    OUTPUT:

        A sequence whose entries are power series in ``ZZ[[var]]``.

    EXAMPLES::

        sage: victor_miller_basis(1, 6)
        []
        sage: victor_miller_basis(0, 6)
        [
        1 + O(q^6)
        ]
        sage: victor_miller_basis(2, 6)
        []
        sage: victor_miller_basis(4, 6)
        [
        1 + 240*q + 2160*q^2 + 6720*q^3 + 17520*q^4 + 30240*q^5 + O(q^6)
        ]

        sage: victor_miller_basis(6, 6, var='w')
        [
        1 - 504*w - 16632*w^2 - 122976*w^3 - 532728*w^4 - 1575504*w^5 + O(w^6)
        ]

        sage: victor_miller_basis(6, 6)
        [
        1 - 504*q - 16632*q^2 - 122976*q^3 - 532728*q^4 - 1575504*q^5 + O(q^6)
        ]
        sage: victor_miller_basis(12, 6)
        [
        1 + 196560*q^2 + 16773120*q^3 + 398034000*q^4 + 4629381120*q^5 + O(q^6),
        q - 24*q^2 + 252*q^3 - 1472*q^4 + 4830*q^5 + O(q^6)
        ]

        sage: victor_miller_basis(12, 6, cusp_only=True)
        [
        q - 24*q^2 + 252*q^3 - 1472*q^4 + 4830*q^5 + O(q^6)
        ]
        sage: victor_miller_basis(24, 6, cusp_only=True)
        [
        q + 195660*q^3 + 12080128*q^4 + 44656110*q^5 + O(q^6),
        q^2 - 48*q^3 + 1080*q^4 - 15040*q^5 + O(q^6)
        ]
        sage: victor_miller_basis(24, 6)
        [
        1 + 52416000*q^3 + 39007332000*q^4 +