def submodule_generated_by_images(self, M):
        """
        Return the submodule of this ambient modular symbols space
        generated by the images under all degeneracy maps of M. The space M
        must have the same weight, sign, and group or character as this
        ambient space.

        EXAMPLES::

            sage: ModularSymbols(6, 12).submodule_generated_by_images(ModularSymbols(1,12))
            Modular Symbols subspace of dimension 12 of Modular Symbols space of dimension 22 for Gamma_0(6) of weight 12 with sign 0 over Rational Field
        """
        S = self.zero_submodule()
        if self.level() % M.level() == 0:
            D = arith.divisors(self.level() // M.level())
        elif M.level() % self.level() == 0:
            D = arith.divisors(M.level() // self.level())
        else:
            D = []
        for t in D:
            d = M.degeneracy_map(self, t)
            if d.codomain() != self:
                raise ArithmeticError, "incompatible spaces of modular symbols"
            S += d.image()

        if self.is_full_hecke_module(compute=False):
            S._is_full_hecke_module = True

        return S

def find_product_decomposition(g, k, lmbda=1):
    r"""
    Try to find a product decomposition construction for difference matrices.

    INPUT:

    - ``g,k,lmbda`` -- integers, parameters of the difference matrix

    OUTPUT:

    A pair of pairs ``(g1,lmbda),(g2,lmbda2)`` if Sage knows how to build
    `(g1,k,lmbda1)` and `(g2,k,lmbda2)` difference matrices and ``False``
    otherwise.

    EXAMPLES::

        sage: from sage.combinat.designs.difference_matrices import find_product_decomposition
        sage: find_product_decomposition(77,6)
        ((7, 1), (11, 1))
        sage: find_product_decomposition(616,7)
        ((7, 1), (88, 1))
        sage: find_product_decomposition(24,10)
        False
    """
    for lmbda1 in divisors(lmbda):
        lmbda2 = lmbda//lmbda1

        # To avoid infinite loop:
        # if lmbda1 == lmbda, then g1 should not be g
        # if lmbda2 == lmbda, then g2 should not be g
        if lmbda1 == lmbda:
            if lmbda2 == lmbda:
                div = divisors(g)[1:-1]
            else:
                div = divisors(g)[:-1]
        else:
            if lmbda2 == lmbda:
                div = divisors(g)[1:]
            else:
                div = divisors(g)

        for g1 in div:
            g2 = g//g1
            if g1 > g2:
                break
            if (difference_matrix(g1,k,lmbda1,existence=True) and
                difference_matrix(g2,k,lmbda2,existence=True)):
                return (g1,lmbda1),(g2,lmbda2)

    return False

def test_Hecke_relations(a,b,C):
    r"""Testing Hecke relations for the Fourier coefficients in C



    INPUT:
    -''C'' -- dictionary of complex (Fourier coefficients)
    -''a'' -- integer
    -''b'' -- integer

    OUTPUT:
    -''diff'' -- real : |C(a)C(b)-C(ab)| if (a,b)=1

    EXAMPLE::

    
    sage: S=MaassWaveForms(Gamma0(1))
    sage: R=mpmath.mpf(9.53369526135355755434423523592877032382125639510725198237579046413534)
    sage: Y=mpmath.mpf(0.85)
    sage: C=coefficients_for_Maass_waveforms(S,R,Y,10,20,12)
    sage: d=test_Hecke_relations(C,2,3); mppr(d)
    '9.29e-8'
    sage: C=coefficients_for_Maass_waveforms(S,R,Y,30,50,20)
    sage: d=test_Hecke_relations(C,2,3); mppr(d)
    '3.83e-43'
    
    
    """
    c=gcd(Integer(a),Integer(b))
    lhs=C[0][a]*C[0][b]
    rhs=mpmath.mpf(0)
    for d in divisors(c):
        rhs=rhs+C[0][Integer(a*b/d/d)]
    return abs(rhs-lhs)

    def _coset_reduction_data_second_coord(G):
        """
        Compute data used for determining the canonical coset
        representative of an element of SL_2(Z) modulo G. This
        function specifically returns data needed for the second part
        of the reduction step (the second coordinate).

        INPUT:
            self

        OUTPUT:
            a dictionary v with keys the divisors of N such that v[d]
            is the subgroup {h in H : h = 1 (mod N/d)}.

        EXAMPLES::

            sage: G = GammaH(240,[7,239])
            sage: G._coset_reduction_data_second_coord()
            {1: [1], 2: [1], 3: [1], 4: [1], 5: [1, 49], 6: [1], 48: [1, 191], 8: [1], 80: [1, 7, 49, 103], 10: [1, 49], 12: [1], 15: [1, 49], 240: [1, 7, 49, 103, 137, 191, 233, 239], 40: [1, 7, 49, 103], 20: [1, 49], 24: [1, 191], 120: [1, 7, 49, 103, 137, 191, 233, 239], 60: [1, 49, 137, 233], 30: [1, 49, 137, 233], 16: [1]}
            sage: G = GammaH(1200,[-1,7]); G
            Congruence Subgroup Gamma_H(1200) with H generated by [7, 1199]
            sage: K = G._coset_reduction_data_second_coord().keys() ; K.sort()
            sage: K == divisors(1200)
            True
        """
        H = G._list_of_elements_in_H()
        N = G.level()
        v = { 1: [1] , N: H }
        for d in [x for x in divisors(N) if x > 1 and x < N ]:
            N_over_d = N // d
            v[d] = [x for x in H if x % N_over_d == 1]
        return v

def number_of_Gamma0_NFCusps(N):
    """
    Returns the total number of orbits of cusps under the action of the
    congruence subgroup `\\Gamma_0(N)`.

    INPUT:

    - ``N`` -- a number field ideal.

    OUTPUT:

    ingeter -- the number of orbits of cusps under Gamma0(N)-action.

    EXAMPLES::

        sage: k.<a> = NumberField(x^3 + 11)
        sage: N = k.ideal(2, a+1)
        sage: from sage.modular.cusps_nf import number_of_Gamma0_NFCusps
        sage: number_of_Gamma0_NFCusps(N)
        4
        sage: L = Gamma0_NFCusps(N)
        sage: len(L) == number_of_Gamma0_NFCusps(N)
        True
    """
    k = N.number_field()
    # The number of Gamma0(N)-sub-orbits for each Gamma-orbit:
    from sage.rings.arith import divisors
    s = sum([len(list((d+N/d).invertible_residues_mod(k.unit_group().gens()))) \
                                                for d in divisors(N)])
    # There are h Gamma-orbits, with h class number of underlying number field.
    return s*k.class_number()

def AllCusps(N):
    r"""
    Return a list of CuspFamily objects corresponding to the cusps of
    `X_0(N)`.

    INPUT:

    -  ``N`` - (integer): the level


    EXAMPLES::

        sage: AllCusps(18)
        [(Inf), (c_{2}), (c_{3,1}), (c_{3,2}), (c_{6,1}), (c_{6,2}), (c_{9}), (0)]
    """
    try:
        N = ZZ(N)
        assert N>0
    except TypeError:
        raise TypeError, "N must be an integer"
    except AssertionError:
        raise AssertionError, "N must be positive"
    c = []
    for d in divisors(N):
        n = num_cusps_of_width(N, d)
        if n == 1:
            c.append(CuspFamily(N, d))
        elif n > 1:
            for i in xrange(n):
                c.append(CuspFamily(N, d, label=str(i+1)))
    return c

    def _b_power_k(self, k):
        r"""
        An expression involving moebius inversion in the powersum generators.

        For a positive value of ``k``, this expression is

        .. MATH::

            \frac{1}{k} \sum_{d|k} \mu(d/k) p_d.

        INPUT:

        - ``k`` -- a positive integer

        OUTPUT:

        - an expression in the powersum basis of the symmetric functions

        EXAMPLES::

            sage: st = SymmetricFunctions(QQ).st()
            sage: st._b_power_k(1)
            p[1]
            sage: st._b_power_k(2)
            -1/2*p[1] + 1/2*p[2]
            sage: st._b_power_k(6)
            1/6*p[1] - 1/6*p[2] - 1/6*p[3] + 1/6*p[6]

        """
        if k == 1:
            return self._p([1])
        if k > 0:
            return ~k * self._p.sum(moebius(k/d)*self._p([d])
                                    for d in divisors(k))

def TD_find_product_decomposition(k,n):
    r"""
    Attempts to find a factorization of `n` in order to build a `TD(k,n)`.

    If Sage can build a `TD(k,n_1)` and a `TD(k,n_2)` such that `n=n_1\times
    n_2` then a `TD(k,n)` can be built (from the function
    :func:`transversal_design`). This method returns such a pair of integers if
    it exists, and ``None`` otherwise.

    INPUT:

    - ``k,n`` (integers) -- see above.

    .. SEEALSO::

        :func:`TD_product` that actually build a product

    EXAMPLES::

        sage: from sage.combinat.designs.orthogonal_arrays import TD_find_product_decomposition
        sage: TD_find_product_decomposition(6, 84)
        (7, 12)

        sage: TD1 = designs.transversal_design(6, 7)
        sage: TD2 = designs.transversal_design(6, 12)
        sage: from sage.combinat.designs.orthogonal_arrays import TD_product
        sage: TD = TD_product(6, TD1, 7, TD2, 12)
    """
    from sage.rings.arith import divisors
    for n1 in divisors(n)[1:-1]: # we ignore 1 and n
        n2 = n//n1
        if transversal_design(k, n1, existence = True) and transversal_design(k, n2, existence = True):
            return n1,n2
    return None

def AllCusps(N):
    r"""
    Return a list of CuspFamily objects corresponding to the cusps of
    `X_0(N)`.

    INPUT:

    -  ``N`` - (integer): the level


    EXAMPLES::

        sage: AllCusps(18)
        [(Inf), (c_{2}), (c_{3,1}), (c_{3,2}), (c_{6,1}), (c_{6,2}), (c_{9}), (0)]
        sage: AllCusps(0)
        Traceback (most recent call last):
        ...
        ValueError: N must be positive
    """
    N = ZZ(N)
    if N <= 0:
        raise ValueError("N must be positive")

    c = []
    for d in divisors(N):
        n = num_cusps_of_width(N, d)
        if n == 1:
            c.append(CuspFamily(N, d))
        elif n > 1:
            for i in xrange(n):
                c.append(CuspFamily(N, d, label=str(i+1)))
    return c

    def _Weyl_law_consts(self):
        r"""
        Compute constants for the Weyl law on self._G

        OUTPUT:

        - tuple of real numbers

        EXAMPLES::


            sage: M=MaassWaveForms(MySubgroup(Gamma0(1)))
            sage: M._Weyl_law_consts()
            (0, 2/pi, (log(pi) - log(2) + 2)/pi, 0, -2)
        """
        import mpmath
        pi=mpmath.fp.pi
        ix=Integer(self._G.index())
        nc=Integer(len(self._G.cusps()))
        if(self._G.is_congruence()):
            lvl=Integer(self._G.level())
        else:
            lvl=0
        n2=Integer(self._G.nu2())
        n3=Integer(self._G.nu3())
        c1=ix/Integer(12)
        c2=Integer(2)*nc/pi
        c3=nc*(Integer(2)-ln(Integer(2))+ln(pi))/pi
        if(lvl<>0):
            A=1
            for q in divisors(lvl):
                num_prim_dc=0
                DG=DirichletGroup(q)
                for chi in DG.list():
                    if(chi.is_primitive()):
                        num_prim_dc=num_prim_dc+1
                for m in divisors(lvl):
                    if(lvl % (m*q) == 0   and m % q ==0 ): 
                        fak=(q*lvl)/gcd(m,lvl/m)
                        A=A*Integer(fak)**num_prim_dc        
            c4=-ln(A)/pi
        else:
            c4=Integer(0)
        # constant term
        c5=-ix/144+n2/8+n3*2/9-nc/4-1
        return (c1,c2,c3,c4,c5)

    def reduce_basis(self, long_etas):
        r"""
        Produce a more manageable basis via LLL-reduction.

        INPUT:


        - ``long_etas`` -  a list of EtaGroupElement objects (which
          should all be of the same level)


        OUTPUT:


        - a new list of EtaGroupElement objects having
          hopefully smaller norm


        ALGORITHM: We define the norm of an eta-product to be the
        `L^2` norm of its divisor (as an element of the free
        `\ZZ`-module with the cusps as basis and the
        standard inner product). Applying LLL-reduction to this gives a
        basis of hopefully more tractable elements. Of course we'd like to
        use the `L^1` norm as this is just twice the degree, which
        is a much more natural invariant, but `L^2` norm is easier
        to work with!

        EXAMPLES::

            sage: EtaGroup(4).reduce_basis([ EtaProduct(4, {1:8,2:24,4:-32}), EtaProduct(4, {1:8, 4:-8})])
            [Eta product of level 4 : (eta_1)^8 (eta_4)^-8,
            Eta product of level 4 : (eta_1)^-8 (eta_2)^24 (eta_4)^-16]
        """
        N = self.level()
        cusps = AllCusps(N)
        r = matrix(ZZ, [[et.order_at_cusp(c) for c in cusps] for et in long_etas])
        V = FreeModule(ZZ, r.ncols())
        A = V.submodule_with_basis([V(rw) for rw in r.rows()])
        rred = r.LLL()
        short_etas = []
        for shortvect in rred.rows():
            bv = A.coordinates(shortvect)
            dict = {}
            for d in divisors(N):
                dict[d] = sum( [bv[i]*long_etas[i].r(d) for i in xrange(r.nrows())])
            short_etas.append(self(dict))
        return short_etas

def p1NFlist(N):
    """
    Returns a list of the normalized elements of `\\mathbb{P}^1(R/N)`, where
    `N` is an integral ideal.

    INPUT:

    -  ``N`` - integral ideal (the level or modulus).

    EXAMPLES::

        sage: k.<a> = NumberField(x^2 + 23)
        sage: N = k.ideal(3)
        sage: from sage.modular.modsym.p1list_nf import p1NFlist, psi
        sage: len(p1NFlist(N))==psi(N)
        True
    """
    k = N.number_field()

    L = [MSymbol(N, k(0),k(1), check=False)]
    #N.residues() = iterator through the residues mod N
    L = L+[MSymbol(N, k(1), r, check=False) for r in N.residues()]

    from sage.rings.arith import divisors
    for D in divisors(N):
        if not D.is_trivial() and D!=N:
            #we find Dp ideal coprime to N, in inverse class to D
            if D.is_principal():
                Dp = k.ideal(1)
                c = D.gens_reduced()[0]
            else:
                it = k.primes_of_degree_one_iter()
                Dp = it.next()
                while not Dp.is_coprime(N) or not (Dp*D).is_principal():
                    Dp = it.next()
                c = (D*Dp).gens_reduced()[0]
            #now we find all the (c,d)'s which have associated divisor D
            I = D + N/D
            for d in (N/D).residues():
                if I.is_coprime(d):
                    M = D.prime_to_idealM_part(N/D)
                    u = (Dp*M).element_1_mod(N/D)
                    d1 = u*d + (1-u)
                    L.append(MSymbol(N, c, d1, check=False).normalize())
    return L

    def cardinality(self):
        r"""
        Return the number of integer necklaces with the evaluation ``content``.

        The formula for the number of necklaces of content `\alpha`
        a composition of `n` is:

        .. MATH::

            \sum_{d|gcd(\alpha)} \phi(d)
            \binom{n/d}{\alpha_1/d, \ldots, \alpha_\ell/d},

        where `\phi(d)` is the Euler `\phi` function.

        EXAMPLES::

            sage: Necklaces([]).cardinality()
            0
            sage: Necklaces([2,2]).cardinality()
            2
            sage: Necklaces([2,3,2]).cardinality()
            30
            sage: Necklaces([0,3,2]).cardinality()
            2

        Check to make sure that the count matches up with the number of
        necklace words generated.

        ::

            sage: comps = [[],[2,2],[3,2,7],[4,2],[0,4,2],[2,0,4]]+Compositions(4).list()
            sage: ns = [ Necklaces(comp) for comp in comps]
            sage: all( [ n.cardinality() == len(n.list()) for n in ns] )
            True
        """
        evaluation = self._content
        le = list(evaluation)
        if not le:
            return 0

        n = sum(le)

        return sum(euler_phi(j)*factorial(n/j) / prod(factorial(ni/j)
                    for ni in evaluation) for j in divisors(gcd(le))) / n

def summand(part, n):
    """
    Create the summand used in the Harrison count for a given partition.

    Args:
        part (tuple): A partition of `n` represented as a tuple.
        n (int): The integer for which `part` is a partition.

    Returns:
        int: The summand corresponding to the partition `part` of `n`.
    """

    t = 1
    count = list(cycle_count(part, n)) + (factorial(n)-n)*[0]
    for i in range(1,n+1):
        for j in range(1,n+1):
            s = sum([d*(count[d-1]) for d in divisors(lcm(i,j))])
            t = t*(s**(count[i-1]*count[j-1]*gcd(i,j)))
    t = t*factorial(n)/(prod(factorial(count[d-1])*(d**(count[d-1])) for d in range(1,n+1)))
    return t

    def _find_cusps(self):
        r"""
        Return an ordered list of inequivalent cusps for self, i.e. a
        set of representatives for the orbits of self on
        `\mathbb{P}^1(\QQ)`.  These are returned in a reduced
        form; see self.reduce_cusp for the definition of reduced.

        ALGORITHM:
            Uses explicit formulae specific to `\Gamma_0(N)`: a reduced cusp on
            `\Gamma_0(N)` is always of the form `a/d` where `d | N`, and `a_1/d
            \sim a_2/d` if and only if `a_1 \cong a_2 \bmod {\rm gcd}(d,
            N/d)`.

        EXAMPLES::

            sage: Gamma0(90)._find_cusps()
            [0, 1/45, 1/30, 1/18, 1/15, 1/10, 1/9, 2/15, 1/6, 1/5, 1/3, 11/30, 1/2, 2/3, 5/6, Infinity]
            sage: Gamma0(1).cusps()
            [Infinity]
            sage: Gamma0(180).cusps() == Gamma0(180).cusps(algorithm='modsym')
            True
        """
        N = self.level()
        s = []

        for d in arith.divisors(N):
            w = arith.gcd(d, N//d)
            if w == 1:
                if d == 1:
                    s.append(Cusp(1,0))
                elif d == N:
                    s.append(Cusp(0,1))
                else:
                    s.append(Cusp(1,d))
            else:
                for a in xrange(1, w):
                    if arith.gcd(a, w) == 1:
                        while arith.gcd(a, d//w) != 1:
                            a += w
                        s.append(Cusp(a,d))
        return sorted(s)

    def nregcusps(self):
        r"""
        Return the number of orbits of regular cusps for this subgroup. A cusp is regular
        if we may find a parabolic element generating the stabiliser of that
        cusp whose eigenvalues are both +1 rather than -1. If G contains -1,
        all cusps are regular.

        EXAMPLES::

            sage: GammaH(20, [17]).nregcusps()
            4
            sage: GammaH(20, [17]).nirregcusps()
            2
            sage: GammaH(3212, [2045, 2773]).nregcusps()
            1440
            sage: GammaH(3212, [2045, 2773]).nirregcusps()
            720

        AUTHOR:

        - Jordi Quer
        """
        if self.is_even():
            return self.ncusps()

        N = self.level()
        H = self._list_of_elements_in_H()

        c = ZZ(0)
        for d in [d for d in divisors(N) if d**2 <= N]:
            Nd = lcm(d,N//d)
            Hd = set([x%Nd for x in H])
            if Nd - 1 not in Hd:
                summand = euler_phi(d)*euler_phi(N//d)//(2*len(Hd))
                if d**2==N:
                    c = c + summand
                else:
                    c = c + 2*summand
        return c

def lift2smallest_field2(a):
    """
    INPUT: a is an element of a finite field GF(q) OUTPUT: the element
    b of the smallest subfield F of GF(q) for which F(b)=a.

    EXAMPLES::

        sage: from sage.coding.code_constructions import lift2smallest_field2
        sage: FF.<z> = GF(3^4,"z")
        sage: a = z^40
        sage: lift2smallest_field2(a)
        (2, Finite Field of size 3)
        sage: FF.<z> = GF(2^4,"z")
        sage: a = z^15
        sage: lift2smallest_field2(a)
        (1, Finite Field of size 2)

    .. warning::

       Since coercion (the FF(b) step) has a bug in it, this
       *only works* in the case when you *know* F is a prime field.

    AUTHORS:

    - David Joyner
    """
    FF = a.parent()
    q = FF.order()
    if q.is_prime():
        return a,FF
    p = q.factor()[0][0]
    k = q.factor()[0][1]
    for d in divisors(k):
        F = GF(p**d,"zz")
        for b in F:
            if FF(b) == a:
                return b, F

def find_product_decomposition(k,n):
    r"""
    Look for a factorization of `n` in order to build an `OA(k,n)`.

    If Sage can build a `OA(k,n_1)` and a `OA(k,n_2)` such that `n=n_1\times
    n_2` then a `OA(k,n)` can be built by a product construction (which
    correspond to Wilson's construction with no truncated column). This
    function look for a pair of integers `(n_1,n_2)` with `n1 \leq n_2`, `n_1
    \times n_2 = n` and such that both an `OA(k,n_1)` and an `OA(k,n_2)` are
    available.

    INPUT:

    - ``k,n`` (integers) -- see above.

    OUTPUT:

    A pair ``f,args`` such that ``f(*args)`` is an `OA(k,n)` or ``False`` if no
    product decomposition was found.

    EXAMPLES::

        sage: from sage.combinat.designs.orthogonal_arrays_recursive import find_product_decomposition
        sage: f,args = find_product_decomposition(6, 84)
        sage: args
        (6, 7, 12, ())
        sage: _ = f(*args)
    """
    from sage.rings.arith import divisors
    for n1 in divisors(n)[1:-1]: # we ignore 1 and n
        n2 = n//n1  # n2 is decreasing along the loop
        if n2 < n1:
            break
        if orthogonal_array(k, n1, existence=True) and orthogonal_array(k, n2, existence=True):
            return simple_wilson_construction, (k,n1,n2,())
    return False

    def additive_lift(self, forms, weight, with_character = False, is_integral = False) :
        """
        Borcherds additive lift to hermitian modular forms of
        degree `2`. This coinsides with Gritsenko's arithmetic lift after
        using the theta decomposition.
        
        INPUT:
            - ``forms``          -- A list of functions accepting an integer and
                                    returning a q-expansion.
            - ``weight``         -- A positive integer; The weight of the lift.
            - ``with_character`` -- A boolean (default: ``False``); Whether the
                                    lift has nontrivial character.
            - ``is_integral``    -- A boolean (default: ``False``); If ``True``
                                    use rings of integral q-expansions over `\Z`.
        
        ALGORITHME:
            We use the explicite formulas in [D].

        TESTS::
            sage: from hermitianmodularforms.hermitianmodularformd2_fegenerators import HermitianModularFormD2AdditiveLift
            sage: HermitianModularFormD2AdditiveLift(4, [1,0,0], -3, 4).coefficients()
            {(2, 3, 2, 2): 720, (1, 1, 1, 1): 27, (1, 0, 0, 2): 270, (3, 3, 3, 3): 2943, (2, 1, 1, 3): 2592, (0, 0, 0, 2): 9, (2, 2, 2, 2): 675, (2, 3, 2, 3): 2160, (1, 1, 1, 2): 216, (3, 0, 0, 3): 8496, (2, 0, 0, 3): 2214, (1, 0, 0, 3): 720, (2, 1, 1, 2): 1080, (0, 0, 0, 1): 1, (3, 3, 2, 3): 4590, (3, 1, 1, 3): 4590, (1, 1, 1, 3): 459, (2, 0, 0, 2): 1512, (1, 0, 0, 1): 72, (0, 0, 0, 0): 1/240, (3, 4, 3, 3): 2808, (0, 0, 0, 3): 28, (3, 2, 2, 3): 4752, (2, 2, 2, 3): 1350}
            sage: HermitianModularFormD2AdditiveLift(4, [0,1,0], -3, 6).coefficients()
            {(2, 3, 2, 2): -19680, (1, 1, 1, 1): -45, (1, 0, 0, 2): -3690, (3, 3, 3, 3): -306225, (2, 1, 1, 3): -250560, (0, 0, 0, 2): 33, (2, 2, 2, 2): -13005, (2, 3, 2, 3): -153504, (1, 1, 1, 2): -1872, (3, 0, 0, 3): -1652640, (2, 0, 0, 3): -295290, (1, 0, 0, 3): -19680, (2, 1, 1, 2): -43920, (0, 0, 0, 1): 1, (3, 3, 2, 3): -948330, (3, 1, 1, 3): -1285290, (1, 1, 1, 3): -11565, (2, 0, 0, 2): -65520, (1, 0, 0, 1): -240, (0, 0, 0, 0): -1/504, (3, 4, 3, 3): -451152, (0, 0, 0, 3): 244, (3, 2, 2, 3): -839520, (2, 2, 2, 3): -108090}
        """
        if with_character and self.__D % 4 != 0 :
            raise ValueError( "Characters are only possible for even discriminants." )

        ## This will be needed if characters are implemented
        if with_character :
            if (Integer(self.__D / 4) % 4) in [-2,2] :
                alpha = (-self.__D / 4, 1/2)
            else :
                alpha = (-self.__D / 8, 1/2)
        
        #minv = 1/2 if with_character else 1
        
        R = self.power_series_ring()
        q = R.gen(0)
            
        (vv_expfactor, vv_basis) = self._additive_lift_vector_valued_basis()
        
        vvform = dict((self._reduce_vector_valued_index(k), R(0)) for (k,_) in self._semireduced_vector_valued_indices_with_discriminant_offset(1))

        for (f,b) in zip(forms, vv_basis) :
            ## We have to apply the scaling of exponents to the form
            f = R( f(self._qexp_precision()) ).add_bigoh(self._qexp_precision()) \
                 .subs({q : q**vv_expfactor})
            
            if not f.is_zero() :
                for (k,e) in b.iteritems() :
                    vvform[k] = vvform[k] + e * f
        
        ## the T = matrix(2,[*, t / 2, \bar t / 2, *] th fourier coefficients of the lift
        ## only depends on (- 4 * D * det(T), eps = gcd(T), \theta \cong t / eps)
        ## if m != 1 we consider 2*T
        maass_coeffs = dict()

        ## TODO: use divisor dictionaries
        if not with_character :
            ## The factor for the exponent of the basis of vector valued forms
            ## and the factor D in the formula for the discriminant are combined
            ## here 
            vv_expfactor = vv_expfactor // (- self.__D)
            for eps in self._iterator_content() : 
                for (theta, offset) in self._semireduced_vector_valued_indices_with_discriminant_offset(eps) :
                    for disc in self._iterator_discriminant(eps, offset) :
                        maass_coeffs[(disc, eps, theta)] = \
                             sum( a**(weight-1) *
                                  vvform[self._reduce_vector_valued_index((theta[0]//a, theta[1]//a))][vv_expfactor * disc // a**2]
                                  for a in divisors(eps))
        else :
            ## The factor for the exponent of the basis of vector valued forms
            ## and the factor D in the formula for the discriminant are combined
            ## here 
            vv_expfactor = (2 * vv_expfactor) // (- self.__D) 

            if self.__D // 4 % 2 == 0 :
                for eps in self._iterator_content() : 
                    for (theta, offset) in self._semireduced_vector_valued_indices_with_discriminant_offset(eps) :
                        for disc in self._iter_discriminant(eps, offset) :
                            maass_coeffs[(disc, eps, theta)] = \
                                 sum( a**(weight-1) * (1 if (theta[0] + theta[1] - 1) % 4 == 0 else -1) *
                                      vvform[self._reduce_vector_valued_index((theta[0]//a, theta[1]//a))][vv_expfactor * disc // a**2]
                                      for a in divisors(eps))
            else :
                for eps in self._iterator_content() : 
                    for (theta, offset) in self._semireduced_vector_valued_indices_with_discriminant_offset(eps) :
                        for disc in self._iter_discriminant(eps, offset) :
                            maass_coeffs[(disc, eps, theta)] = \
                                 sum( a**(weight-1) * (1 if (theta[1] - 1) % 4 == 0 else -1) *
                                      vvform[self._reduce_vector_valued_index((theta[0]//a, theta[1]//a))][vv_expfactor * disc // a**2]
                                      for a in divisors(eps) )
        lift_coeffs = dict()
        ## TODO: Check whether this is correct. Add the character as an argument.
        for ((a,b1,b2,c), eps, disc) in self.precision().iter_positive_forms_for_character_with_content_and_discriminant(for_character = with_character) :
            (theta1, theta2) = self._reduce_vector_valued_index((b1/eps, b2/eps))
            theta = (eps * theta1, eps * theta2)
            try:
                lift_coeffs[(a,b1,b2,c)] = maass_coeffs[(disc, eps, theta)]
#.........这里部分代码省略.........

def enumerate_totallyreal_fields_all(n, B, verbose=0, return_seqs=False,
                                     return_pari_objects=True):
    r"""
    Enumerates *all* totally real fields of degree ``n`` with discriminant
    at most ``B``, primitive or otherwise.

    INPUT:

    - ``n`` -- integer, the degree
    - ``B`` -- integer, the discriminant bound
    - ``verbose`` -- boolean or nonnegative integer or string (default: 0)
      give a verbose description of the computations being performed. If
      ``verbose`` is set to ``2`` or more then it outputs some extra
      information. If ``verbose`` is a string then it outputs to a file
      specified by ``verbose``
    - ``return_seqs`` -- (boolean, default False) If ``True``, then return
      the polynomials as sequences (for easier exporting to a file). This
      also returns a list of four numbers, as explained in the OUTPUT
      section below.
    - ``return_pari_objects`` -- (boolean, default: True) if both
      ``return_seqs`` and ``return_pari_objects`` are ``False`` then it
      returns the elements as Sage objects; otherwise it returns pari
      objects.

    EXAMPLES::

        sage: enumerate_totallyreal_fields_all(4, 2000)
        [[725, x^4 - x^3 - 3*x^2 + x + 1],
        [1125, x^4 - x^3 - 4*x^2 + 4*x + 1],
        [1600, x^4 - 6*x^2 + 4],
        [1957, x^4 - 4*x^2 - x + 1],
        [2000, x^4 - 5*x^2 + 5]]
        sage: enumerate_totallyreal_fields_all(1, 10)
        [[1, x - 1]]

    TESTS:

    Each of the outputs must be elements of Sage if ``return_pari_objects``
    is set to ``False``::

        sage: enumerate_totallyreal_fields_all(2, 10)
        [[5, x^2 - x - 1], [8, x^2 - 2]]
        sage: enumerate_totallyreal_fields_all(2, 10)[0][1].parent()
        Interface to the PARI C library
        sage: enumerate_totallyreal_fields_all(2, 10, return_pari_objects=False)[0][1].parent()
        Univariate Polynomial Ring in x over Rational Field


    In practice most of these will be found by
    :func:`~sage.rings.number_field.totallyreal.enumerate_totallyreal_fields_prim`,
    which is guaranteed to return all primitive fields but often returns
    many non-primitive ones as well. For instance, only one of the five
    fields in the example above is primitive, but
    :func:`~sage.rings.number_field.totallyreal.enumerate_totallyreal_fields_prim`
    finds four out of the five (the exception being `x^4 - 6x^2 + 4`).

    The following was fixed in :trac:`13101`::

        sage: enumerate_totallyreal_fields_all(8, 10^6)  # long time (about 2 s)
        []
    """

    S = []
    counts = [0,0,0,0]
    if len(divisors(n)) > 4:
        raise ValueError("Only implemented for n = p*q with p,q prime")
    for d in divisors(n):
        if d > 1 and d < n:
            Sds = enumerate_totallyreal_fields_prim(d, int(math.floor((1.*B)**(1.*d/n))), verbose=verbose)
            for i in range(len(Sds)):
                if verbose:
                    print "="*80
                    print "Taking F =", Sds[i][1]
                F = NumberField(ZZx(Sds[i][1]), 't')
                T = enumerate_totallyreal_fields_rel(F, n/d, B, verbose=verbose, return_seqs=return_seqs)
                if return_seqs:
                    for i in range(4):
                        counts[i] += T[0][i]
                    S += [[t[0],pari(t[1]).Polrev()] for t in T[1]]
                else:
                    S += [[t[0],t[1]] for t in T]
                j = i+1
                for E in enumerate_totallyreal_fields_prim(n/d, int(math.floor((1.*B)**(1./d)/(1.*Sds[i][0])**(n*1./d**2)))):
                    for EF in F.composite_fields(NumberField(ZZx(E[1]), 'u')):
                        if EF.degree() == n and EF.disc() <= B:
                            S.append([EF.disc(), pari(EF.absolute_polynomial())])
    S += enumerate_totallyreal_fields_prim(n, B, verbose=verbose)
    S.sort(cmp=lambda x, y: cmp(x[0], y[0]) or cmp(x[1], y[1]))
    weed_fields(S)

    # Output.
    if verbose:
        saveout = sys.stdout
        if isinstance(verbose, str):
            fsock = open(verbose, 'w')
            sys.stdout = fsock
        # Else, print to screen
        print "="*80
        print "Polynomials tested: {}".format(counts[0])
        print ( "Polynomials with discriminant with large enough square"
#.........这里部分代码省略.........