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 _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 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 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 1
        else:
            s = 0
            for d in divisors(self.k):
                s += moebius(d)*(self.n**(self.k/d))
        return s/self.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 __iter__(self):
        r"""
        Iterate over ``self``.

        EXAMPLES::

            sage: PTC = PrimarySimilarityClassTypes(2)
            sage: PTC.cardinality()
            3
        """
        n = self._n
        if self._min[0].divides(n):
            for par in Partitions(ZZ(n / self._min[0]), starting=self._min[1]):
                yield self.element_class(self, self._min[0], par)
        for d in filter(lambda d: d > self._min[0], divisors(n)):
            for par in Partitions(ZZ(n / d)):
                yield self.element_class(self, d, par)

def primitives(n, invertible=False, q=None):
    """
    Return the number of similarity classes of simple matrices
    of order n with entries in a finite field of order ``q``.
    This is the same as the number of irreducible polynomials
    of degree d.

    If ``invertible`` is ``True``, then only the number of
    similarity classes of invertible matrices is returned.

    .. NOTE::

        All primitive classes are invertible unless ``n`` is `1`.

    INPUT:

    - ``n`` -- a positive integer

    - ``invertible`` -- boolean; if set, only number of non-zero classes is returned

    - ``q`` -- an integer or an indeterminate

    OUTPUT:

    - a rational function of the variable ``q``

    EXAMPLES::

        sage: from sage.combinat.similarity_class_type import primitives
        sage: primitives(1)
        q
        sage: primitives(1, invertible = True)
        q - 1
        sage: primitives(4)
        1/4*q^4 - 1/4*q^2
        sage: primitives(4, invertible = True)
        1/4*q^4 - 1/4*q^2
    """
    if q is None:
        q = QQ["q"].gen()
    p = sum([moebius(n / d) * q ** d for d in divisors(n)]) / n
    if invertible and n == 1:
        return p - 1
    else:
        return p

    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 cardinality(self):
        """
        Returns the number of Lyndon words with the evaluation e.

        EXAMPLES::

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

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

        ::

            sage: comps = [[],[2,2],[3,2,7],[4,2]]+Compositions(4).list()
            sage: lws = [ LyndonWords(comp) for comp in comps]
            sage: all( [ lw.cardinality() == len(lw.list()) for lw in lws] )
            True
        """
        evaluation = self.e
        le = __builtin__.list(evaluation)
        if len(evaluation) == 0:
            return 0

        n = sum(evaluation)

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

def __find_eisen_chars_gamma1(N, k):
    """
    Find all triples `(\psi_1, \psi_2, t)` that give rise to an Eisenstein series of weight `k` on
    `\Gamma_1(N)`.

    EXAMPLES::

        sage: pars = sage.modular.modform.eis_series.__find_eisen_chars_gamma1(12, 4)
        sage: [(x[0].values_on_gens(), x[1].values_on_gens(), x[2]) for x in pars]
        [((1, 1), (1, 1), 1),
        ((1, 1), (1, 1), 2),
        ((1, 1), (1, 1), 3),
        ((1, 1), (1, 1), 4),
        ((1, 1), (1, 1), 6),
        ((1, 1), (1, 1), 12),
        ((1, 1), (-1, -1), 1),
        ((-1, -1), (1, 1), 1),
        ((-1, 1), (1, -1), 1),
        ((1, -1), (-1, 1), 1)]

        sage: pars =  sage.modular.modform.eis_series.__find_eisen_chars_gamma1(12, 5)
        sage: [(x[0].values_on_gens(), x[1].values_on_gens(), x[2]) for x in pars]
        [((1, 1), (-1, 1), 1),
        ((1, 1), (-1, 1), 3),
        ((-1, 1), (1, 1), 1),
        ((-1, 1), (1, 1), 3),
        ((1, 1), (1, -1), 1),
        ((1, 1), (1, -1), 2),
        ((1, 1), (1, -1), 4),
        ((1, -1), (1, 1), 1),
        ((1, -1), (1, 1), 2),
        ((1, -1), (1, 1), 4)]
    """
    pairs = []
    s = (-1)**k
    G = dirichlet.DirichletGroup(N)
    E = list(G)
    parity = [c(-1) for c in E]
    for i in range(len(E)):
        for j in range(i,len(E)):
            if parity[i]*parity[j] == s and N % (E[i].conductor()*E[j].conductor()) == 0:
                chi, psi = __common_minimal_basering(E[i], E[j])
                if k != 1:
                    pairs.append((chi, psi))
                    if i!=j: pairs.append((psi,chi))
                else:
                    # if weight is 1 then (chi, psi) and (chi, psi) are the
                    # same form
                    if psi.is_trivial() and not chi.is_trivial():
                        # need to put the trivial character first to get the L-value right
                        pairs.append((psi, chi))
                    else:
                        pairs.append((chi, psi))
        #end fors
    #end if

    triples = []
    D = divisors(N)
    for chi, psi in pairs:
        c_chi = chi.conductor()
        c_psi = psi.conductor()
        D = divisors(N/(c_chi * c_psi))
        if (k==2 and chi.is_trivial() and psi.is_trivial()):
            D.remove(1)
        chi, psi = __common_minimal_basering(chi, psi)
        for t in D:
            triples.append((chi, psi, t))
    return triples

def __find_eisen_chars(character, k):
    """
    Find all triples `(\psi_1, \psi_2, t)` that give rise to an Eisenstein series of the given weight and character.

    EXAMPLES::

        sage: sage.modular.modform.eis_series.__find_eisen_chars(DirichletGroup(36).0, 4)
        []

        sage: pars =  sage.modular.modform.eis_series.__find_eisen_chars(DirichletGroup(36).0, 5)
        sage: [(x[0].values_on_gens(), x[1].values_on_gens(), x[2]) for x in pars]
        [((1, 1), (-1, 1), 1),
        ((1, 1), (-1, 1), 3),
        ((1, 1), (-1, 1), 9),
        ((1, -1), (-1, -1), 1),
        ((-1, 1), (1, 1), 1),
        ((-1, 1), (1, 1), 3),
        ((-1, 1), (1, 1), 9),
        ((-1, -1), (1, -1), 1)]
    """
    N = character.modulus()
    if character.is_trivial():
        if k%2 != 0:
            return []
        char_inv = ~character
        V = [(character, char_inv, t) for t in divisors(N) if t>1]
        if k != 2:
            V.insert(0,(character, char_inv, 1))
        if is_squarefree(N):
            return V
        # Now include all pairs (chi,chi^(-1)) such that cond(chi)^2 divides N:
        # TODO: Optimize -- this is presumably way too hard work below.
        G = dirichlet.DirichletGroup(N)
        for chi in G:
            if not chi.is_trivial():
                f = chi.conductor()
                if N % (f**2) == 0:
                    chi = chi.minimize_base_ring()
                    chi_inv = ~chi
                    for t in divisors(N//(f**2)):
                        V.insert(0, (chi, chi_inv, t))
        return V


    eps = character
    if eps(-1) != (-1)**k:
        return []
    eps = eps.maximize_base_ring()
    G = eps.parent()

    # Find all pairs chi, psi such that:
    #
    #  (1) cond(chi)*cond(psi) divides the level, and
    #
    #  (2) chi*psi == eps, where eps is the nebentypus character of self.
    #
    # See [Miyake, Modular Forms] Lemma 7.1.1.

    K = G.base_ring()
    C = {}

    t0 = misc.cputime()

    for e in G:
        m = Integer(e.conductor())
        if m in C:
            C[m].append(e)
        else:
            C[m] = [e]

    misc.verbose("Enumeration with conductors.",t0)

    params = []
    for L in divisors(N):
        misc.verbose("divisor %s"%L)
        if L not in C:
            continue
        GL = C[L]
        for R in divisors(N/L):
            if R not in C:
                continue
            GR = C[R]
            for chi in GL:
                for psi in GR:
                    if chi*psi == eps:
                        chi0, psi0 = __common_minimal_basering(chi, psi)
                        for t in divisors(N//(R*L)):
                            if k != 1 or ((psi0, chi0, t) not in params):
                                params.append( (chi0,psi0,t) )
    return params

    def dimension_new_cusp_forms(self, k=2, eps=None, p=0, algorithm="CohenOesterle"):
        r"""
        Dimension of the new subspace (or `p`-new subspace) of cusp forms of
        weight `k` and character `\varepsilon`.

        INPUT:

        - ``k`` - an integer (default: 2)

        - ``eps`` - a Dirichlet character

        -  ``p`` - a prime (default: 0); just the `p`-new subspace if given

        - ``algorithm`` - either "CohenOesterle" (the default) or "Quer". This
          specifies the method to use in the case of nontrivial character:
          either the Cohen--Oesterle formula as described in Stein's book, or
          by Moebius inversion using the subgroups GammaH (a method due to
          Jordi Quer).

        EXAMPLES::

            sage: G = DirichletGroup(9)
            sage: eps = G.0^3
            sage: eps.conductor()
            3
            sage: [Gamma1(9).dimension_new_cusp_forms(k, eps) for k in [2..10]]
            [0, 0, 0, 2, 0, 2, 0, 2, 0]
            sage: [Gamma1(9).dimension_cusp_forms(k, eps) for k in [2..10]]
            [0, 0, 0, 2, 0, 4, 0, 6, 0]
            sage: [Gamma1(9).dimension_new_cusp_forms(k, eps, 3) for k in [2..10]]
            [0, 0, 0, 2, 0, 2, 0, 2, 0]

        Double check using modular symbols (independent calculation)::

            sage: [ModularSymbols(eps,k,sign=1).cuspidal_subspace().new_subspace().dimension()  for k in [2..10]]
            [0, 0, 0, 2, 0, 2, 0, 2, 0]
            sage: [ModularSymbols(eps,k,sign=1).cuspidal_subspace().new_subspace(3).dimension()  for k in [2..10]]
            [0, 0, 0, 2, 0, 2, 0, 2, 0]

        Another example at level 33::

            sage: G = DirichletGroup(33)
            sage: eps = G.1
            sage: eps.conductor()
            11
            sage: [Gamma1(33).dimension_new_cusp_forms(k, G.1) for k in [2..4]]
            [0, 4, 0]
            sage: [Gamma1(33).dimension_new_cusp_forms(k, G.1, algorithm="Quer") for k in [2..4]]
            [0, 4, 0]
            sage: [Gamma1(33).dimension_new_cusp_forms(k, G.1^2) for k in [2..4]]
            [2, 0, 6]
            sage: [Gamma1(33).dimension_new_cusp_forms(k, G.1^2, p=3) for k in [2..4]]
            [2, 0, 6]

        """

        if eps == None:
            return GammaH_class.dimension_new_cusp_forms(self, k, p)

        N = self.level()
        eps = DirichletGroup(N)(eps)

        from all import Gamma0

        if eps.is_trivial():
            return Gamma0(N).dimension_new_cusp_forms(k, p)

        from congroup_gammaH import mumu

        if p == 0 or N%p != 0 or eps.conductor().valuation(p) == N.valuation(p):
            D = [eps.conductor()*d for d in divisors(N//eps.conductor())]
            return sum([Gamma1_constructor(M).dimension_cusp_forms(k, eps.restrict(M), algorithm)*mumu(N//M) for M in D])
        eps_p = eps.restrict(N//p)
        old = Gamma1_constructor(N//p).dimension_cusp_forms(k, eps_p, algorithm)
        return self.dimension_cusp_forms(k, eps, algorithm) - 2*old

def hecke_operator_on_qexp(f, n, k, eps=None, prec=None, check=True, _return_list=False):
    r"""
    Given the `q`-expansion `f` of a modular form with character
    `\varepsilon`, this function computes the image of `f` under the
    Hecke operator `T_{n,k}` of weight `k`.

    EXAMPLES::

        sage: M = ModularForms(1,12)
        sage: hecke_operator_on_qexp(M.basis()[0], 3, 12)
        252*q - 6048*q^2 + 63504*q^3 - 370944*q^4 + O(q^5)
        sage: hecke_operator_on_qexp(M.basis()[0], 1, 12, prec=7)
        q - 24*q^2 + 252*q^3 - 1472*q^4 + 4830*q^5 - 6048*q^6 + O(q^7)
        sage: hecke_operator_on_qexp(M.basis()[0], 1, 12)
        q - 24*q^2 + 252*q^3 - 1472*q^4 + 4830*q^5 - 6048*q^6 - 16744*q^7 + 84480*q^8 - 113643*q^9 - 115920*q^10 + 534612*q^11 - 370944*q^12 - 577738*q^13 + O(q^14)

        sage: M.prec(20)
        20
        sage: hecke_operator_on_qexp(M.basis()[0], 3, 12)
        252*q - 6048*q^2 + 63504*q^3 - 370944*q^4 + 1217160*q^5 - 1524096*q^6 + O(q^7)
        sage: hecke_operator_on_qexp(M.basis()[0], 1, 12)
        q - 24*q^2 + 252*q^3 - 1472*q^4 + 4830*q^5 - 6048*q^6 - 16744*q^7 + 84480*q^8 - 113643*q^9 - 115920*q^10 + 534612*q^11 - 370944*q^12 - 577738*q^13 + 401856*q^14 + 1217160*q^15 + 987136*q^16 - 6905934*q^17 + 2727432*q^18 + 10661420*q^19 - 7109760*q^20 + O(q^21)

        sage: (hecke_operator_on_qexp(M.basis()[0], 1, 12)*252).add_bigoh(7)
        252*q - 6048*q^2 + 63504*q^3 - 370944*q^4 + 1217160*q^5 - 1524096*q^6 + O(q^7)

        sage: hecke_operator_on_qexp(M.basis()[0], 6, 12)
        -6048*q + 145152*q^2 - 1524096*q^3 + O(q^4)

    An example on a formal power series::

        sage: R.<q> = QQ[[]]
        sage: f = q + q^2 + q^3 + q^7 + O(q^8)
        sage: hecke_operator_on_qexp(f, 3, 12)
        q + O(q^3)
        sage: hecke_operator_on_qexp(delta_qexp(24), 3, 12).prec()
        8
        sage: hecke_operator_on_qexp(delta_qexp(25), 3, 12).prec()
        9

    An example of computing `T_{p,k}` in characteristic `p`::

        sage: p = 199
        sage: fp = delta_qexp(prec=p^2+1, K=GF(p))
        sage: tfp = hecke_operator_on_qexp(fp, p, 12)
        sage: tfp == fp[p] * fp
        True
        sage: tf = hecke_operator_on_qexp(delta_qexp(prec=p^2+1), p, 12).change_ring(GF(p))
        sage: tfp == tf
        True
    """
    if eps is None:
        # Need to have base_ring=ZZ to work over finite fields, since
        # ZZ can coerce to GF(p), but QQ can't.
        eps = DirichletGroup(1, base_ring=ZZ)[0]
    if check:
        if not (is_PowerSeries(f) or is_ModularFormElement(f)):
            raise TypeError("f (=%s) must be a power series or modular form" % f)
        if not is_DirichletCharacter(eps):
            raise TypeError("eps (=%s) must be a Dirichlet character" % eps)
        k = Integer(k)
        n = Integer(n)
    v = []

    if prec is None:
        if is_ModularFormElement(f):
            # always want at least three coefficients, but not too many, unless
            # requested
            pr = max(f.prec(), f.parent().prec(), (n + 1) * 3)
            pr = min(pr, 100 * (n + 1))
            prec = pr // n + 1
        else:
            prec = (f.prec() / ZZ(n)).ceil()
            if prec == Infinity:
                prec = f.parent().default_prec() // n + 1

    if f.prec() < prec:
        f._compute_q_expansion(prec)

    p = Integer(f.base_ring().characteristic())
    if k != 1 and p.is_prime() and n.is_power_of(p):
        # if computing T_{p^a} in characteristic p, use the simpler (and faster)
        # formula
        v = [f[m * n] for m in range(prec)]
    else:
        l = k - 1
        for m in range(prec):
            am = sum([eps(d) * d ** l * f[m * n // (d * d)] for d in divisors(gcd(n, m)) if (m * n) % (d * d) == 0])
            v.append(am)
    if _return_list:
        return v
    if is_ModularFormElement(f):
        R = f.parent()._q_expansion_ring()
    else:
        R = f.parent()
    return R(v, prec)

def _cl_term(n, R = RationalField()):
    """
    Compute the order-n term of the cycle index series of the virtual species `\Omega`,
    the compositional inverse of the species `E^{+}` of nonempty sets.

    EXAMPLES::

        sage: from sage.combinat.species.generating_series import _cl_term
        sage: [_cl_term(i) for i in range(4)]
        [0, p[1], -1/2*p[1, 1] - 1/2*p[2], 1/3*p[1, 1, 1] - 1/3*p[3]]
    """

    n = Integer(n) #check that n is an integer

    p = SymmetricFunctions(R).power()

    res = p.zero()
    if n == 1:
        res = p([1])
    elif n > 1:
        res = 1/n * ((-1)**(n-1) * p([1])**n - sum(d * p([Integer(n/d)]).plethysm(_cl_term(d, R)) for d in divisors(n)[:-1]))

    return res