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 _cycle_type(self, s):
        """
        EXAMPLES::

            sage: cis = species.PartitionSpecies().cycle_index_series()
            sage: [cis._cycle_type(p) for p in Partitions(3)]
            [[3, 1, 1], [2, 1, 1, 1], [1, 1, 1, 1, 1]]
            sage: cis = species.PermutationSpecies().cycle_index_series()
            sage: [cis._cycle_type(p) for p in Partitions(3)]
            [[3, 1, 1, 1], [2, 2, 1, 1], [1, 1, 1, 1, 1, 1]]
            sage: cis = species.SetSpecies().cycle_index_series()
            sage: [cis._cycle_type(p) for p in Partitions(3)]
            [[1], [1], [1]]
        """
        if s == []:
            return self._card(0)
        res = []
        for k in range(1, self._upper_bound_for_longest_cycle(s)+1):
            e = 0
            for d in divisors(k):
                m = moebius(d)
                if m == 0:
                    continue
                u = s.power(k/d)
                e += m*self.count(u)
            res.extend([k]*int(e/k))
        res.reverse()
        return Partition(res)

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 _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.linear_combination((self._p([d]),moebius(k//d))
                                    for d in divisors(k))

    def possible_orders(self):
        """
        Return the possible orders of this torsion subgroup, computed from
        a known divisor and multiple of the order.

        EXAMPLES::

            sage: J0(11).rational_torsion_subgroup().possible_orders()
            [5]
            sage: J0(33).rational_torsion_subgroup().possible_orders()
            [100, 200]

        Note that this function has not been implemented for `J_1(N)`,
        though it should be reasonably easy to do so soon (see Conrad,
        Edixhoven, Stein)::

            sage: J1(13).rational_torsion_subgroup().possible_orders()
            Traceback (most recent call last):
            ...
            NotImplementedError: torsion multiple only implemented for Gamma0
        """
        try:
            return self._possible_orders
        except AttributeError:
            pass
        u = self.multiple_of_order()
        l = self.divisor_of_order()
        assert u % l == 0
        O = [l * d for d in divisors(u//l)]
        self._possible_orders = O
        return O

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.arith.all 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 arith_prod_coeff(n):
            if n == 0:
                res = p.zero()
            else:
                index_set = ((d, n // d) for d in divisors(n))
                res = sum(arith_prod_sf(self.coefficient(i), g.coefficient(j)) for i,j in index_set)

            # Build a list which has res in the `n`th slot and 0's before and after
            # to feed to sum_generator
            res_in_seq = [p.zero()]*n + [res, p.zero()]

            return self.parent(res_in_seq)

    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],
             8: [1],
             10: [1, 49],
             12: [1],
             15: [1, 49],
             16: [1],
             20: [1, 49],
             24: [1, 191],
             30: [1, 49, 137, 233],
             40: [1, 7, 49, 103],
             48: [1, 191],
             60: [1, 49, 137, 233],
             80: [1, 7, 49, 103],
             120: [1, 7, 49, 103, 137, 191, 233, 239],
             240: [1, 7, 49, 103, 137, 191, 233, 239]}
            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 cardinality(self):
        """
        TESTS::

            sage: [ LyndonWords(3,i).cardinality() for i in range(1, 11) ]
            [3, 3, 8, 18, 48, 116, 312, 810, 2184, 5880]
        """
        if self._k == 0:
            return Integer(1)
        else:
            s = Integer(0)
            for d in divisors(self._k):
                s += moebius(d)*(self._n**(self._k/d))
        return s/self._k

    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]
        """
        from six.moves import range
        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 = {d: sum(bv[i] * long_etas[i].r(d)
                           for i in range(r.nrows()))
                    for d in divisors(N)}
            short_etas.append(self(dict))
        return short_etas

    def eval_at_permutation_roots_on_generators(self, k, rho):
        r"""
        Evaluate `p_k` at eigenvalues of permutation matrix.

        This function evaluates a symmetric function ``p([k])``
        at the eigenvalues of a permutation matrix with cycle
        structure ``\rho``.

        This function evaluates a `p_k` at the roots of unity

        .. MATH::

            \Xi_{\rho_1},\Xi_{\rho_2},\ldots,\Xi_{\rho_\ell}

        where

        .. MATH::

            \Xi_{m} = 1,\zeta_m,\zeta_m^2,\ldots,\zeta_m^{m-1}

        and `\zeta_m` is an `m` root of unity.
        This is characterized by `p_k[ A , B ] = p_k[A] + p_k[B]` and
        `p_k[ \Xi_m ] = 0` unless `m` divides `k` and `p_{rm}[\Xi_m]=m`.

        INPUT:

        - ``k`` -- a non-negative integer
        - ``rho`` -- a partition or a list of non-negative integers

        OUTPUT:

        - an element of the base ring

        EXAMPLES::

            sage: p = SymmetricFunctions(QQ).p()
            sage: p.eval_at_permutation_roots_on_generators(3, [6])
            0
            sage: p.eval_at_permutation_roots_on_generators(3, [3])
            3
            sage: p.eval_at_permutation_roots_on_generators(3, [1])
            1
            sage: p.eval_at_permutation_roots_on_generators(3, [3,3])
            6
            sage: p.eval_at_permutation_roots_on_generators(3, [1,1,1,1,1])
            5
        """
        return self.base_ring().sum(d*list(rho).count(d) for d in divisors(k))

    def pAdicEisensteinSeries(self, ring, prec=20):
        r"""
        Calculate the q-expansion of the p-adic Eisenstein series of given
        weight-character, normalised so the constant term is 1.

        EXAMPLE::

            sage: kappa = pAdicWeightSpace(3)(3, DirichletGroup(3,QQ).0)
            sage: kappa.pAdicEisensteinSeries(QQ[['q']], 20)
            1 - 9*q + 27*q^2 - 9*q^3 - 117*q^4 + 216*q^5 + 27*q^6 - 450*q^7 + 459*q^8 - 9*q^9 - 648*q^10 + 1080*q^11 - 117*q^12 - 1530*q^13 + 1350*q^14 + 216*q^15 - 1845*q^16 + 2592*q^17 + 27*q^18 - 3258*q^19 + O(q^20)
        """
        if not self.is_even():
            raise ValueError("Eisenstein series not defined for odd weight-characters")
        q = ring.gen()
        s = ring(1) + 2*self.one_over_Lvalue() * sum([sum([self(d)/d for d in divisors(n)]) * q**n for n in xrange(1, prec)])
        return s.add_bigoh(prec)

    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 ZZ.zero()

        n = sum(le)

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

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.arith.all 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 = next(it)
                while not Dp.is_coprime(N) or not (Dp*D).is_principal():
                    Dp = next(it)
                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 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)
        doctest:...: DeprecationWarning: lift2smallest_field2 will be removed in a future release of Sage. Consider using sage.coding.code_constructions._lift2smallest_field instead, though this is private and may be removed in the future without deprecation warning. If you care about this functionality being in Sage, consider opening a Trac ticket for promoting the function to public.
        See http://trac.sagemath.org/21165 for details.
        (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
    """
    deprecation(
        21165,
        "lift2smallest_field2 will be removed in a future release of Sage. Consider using sage.coding.code_constructions._lift2smallest_field instead, though this is private and may be removed in the future without deprecation warning. If you care about this functionality being in Sage, consider opening a Trac ticket for promoting the function to public.",
    )
    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 _cis_iterator(self, base_ring):
        r"""
        The cycle index series of the species of cyclic permutations is
        given by

        .. MATH::

             -\sum_{k=1}^\infty \phi(k)/k * log(1 - x_k)


        which is equal to

        .. MATH::

             \sum_{n=1}^\infty \frac{1}{n} * \sum_{k|n} \phi(k) * x_k^{n/k}

        .

        EXAMPLES::

            sage: P = species.CycleSpecies()
            sage: cis = P.cycle_index_series()
            sage: cis.coefficients(7)
            [0,
             p[1],
             1/2*p[1, 1] + 1/2*p[2],
             1/3*p[1, 1, 1] + 2/3*p[3],
             1/4*p[1, 1, 1, 1] + 1/4*p[2, 2] + 1/2*p[4],
             1/5*p[1, 1, 1, 1, 1] + 4/5*p[5],
             1/6*p[1, 1, 1, 1, 1, 1] + 1/6*p[2, 2, 2] + 1/3*p[3, 3] + 1/3*p[6]]
        """
        from sage.combinat.sf.sf import SymmetricFunctions
        p = SymmetricFunctions(base_ring).power()

        zero = base_ring(0)

        yield zero
        for n in _integers_from(1):
            res = zero
            for k in divisors(n):
                res += euler_phi(k)*p([k])**(n//k)
            res /= n
            yield self._weight*res

    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