def _find_cusps(self):
        r"""
        Calculate the reduced representatives of the equivalence classes of
        cusps for this group. Adapted from code by Ron Evans.

        EXAMPLE::

            sage: Gamma(8).cusps() # indirect doctest
            [0, 1/4, 1/3, 3/8, 1/2, 2/3, 3/4, 1, 4/3, 3/2, 5/3, 2, 7/3, 5/2, 8/3, 3, 7/2, 11/3, 4, 14/3, 5, 6, 7, Infinity]
        """
        n = self.level()
        C=[QQ(x) for x in xrange(n)]

        n0=n//2
        n1=(n+1)//2

        for r in xrange(1, n1):
            if r > 1 and gcd(r,n)==1:
                C.append(ZZ(r)/ZZ(n))
            if n0==n/2 and gcd(r,n0)==1:
                C.append(ZZ(r)/ZZ(n0))

        for s in xrange(2,n1):
            for r in xrange(1, 1+n):
                if GCD_list([s,r,n])==1:
                    # GCD_list is ~40x faster than gcd, since gcd wastes loads
                    # of time initialising a Sequence type.
                    u,v = _lift_pair(r,s,n)
                    C.append(ZZ(u)/ZZ(v))

        return [Cusp(x) for x in sorted(C)] + [Cusp(1,0)]

def _lift_pair(U,V,N):
    r"""
    Utility function. Given integers ``U, V, N``, with `N \ge 1` and `{\rm
    gcd}(U, V, N) = 1`, return a pair `(u, v)` congruent to `(U, V) \bmod N`,
    such that `{\rm gcd}(u,v) = 1`, `u, v \ge 0`, `v` is as small as possible,
    and `u` is as small as possible for that `v`.

    *Warning*: As this function is for internal use, it does not do a
    preliminary sanity check on its input, for efficiency. It will recover
    reasonably gracefully if ``(U, V, N)`` are not coprime, but only after
    wasting quite a lot of cycles!

    EXAMPLE::

        sage: from sage.modular.arithgroup.congroup_gamma import _lift_pair
        sage: _lift_pair(2,4,7)
        (9, 4)
        sage: _lift_pair(2,4,8) # don't do this
        Traceback (most recent call last):
        ...
        ValueError: (U, V, N) must be coprime
    """
    u = U % N
    v = V % N
    if v == 0:
        if u == 1:
            return (1,0)
        else:
            v = N
    while gcd(u, v) > 1:
        u = u+N
        if u > N*v: raise ValueError, "(U, V, N) must be coprime"
    return (u, v)

    def primitive(self, signed=True):
        """
        Return hyperplane defined by primitive equation.

        INPUT:

        - ``signed`` -- boolean (optional, default: ``True``); whether
          to preserve the overall sign

        OUTPUT:

        Hyperplane whose linear expression has common factors and
        denominators cleared. That is, the same hyperplane (with the
        same sign) but defined by a rescaled equation. Note that
        different linear expressions must define different hyperplanes
        as comparison is used in caching.

        If ``signed``, the overall rescaling is by a positive constant
        only.

        EXAMPLES::

            sage: H.<x,y> = HyperplaneArrangements(QQ)
            sage: h = -1/3*x + 1/2*y - 1;  h
            Hyperplane -1/3*x + 1/2*y - 1
            sage: h.primitive()
            Hyperplane -2*x + 3*y - 6
            sage: h == h.primitive()
            False
            sage: (4*x + 8).primitive()
            Hyperplane x + 0*y + 2

            sage: (4*x - y - 8).primitive(signed=True)   # default
            Hyperplane 4*x - y - 8
            sage: (4*x - y - 8).primitive(signed=False)
            Hyperplane -4*x + y + 8
        """
        from sage.rings.all import lcm, gcd
        coeffs = self.coefficients()
        try:
            d = lcm([x.denom() for x in coeffs])
            n = gcd([x.numer() for x in coeffs])
        except AttributeError:
            return self
        if not signed:
            for x in coeffs:
                if x > 0:
                    break
                if x < 0: 
                    d = -d
                    break
        parent = self.parent()
        d = parent.base_ring()(d)
        n = parent.base_ring()(n)
        if n == 0:
            n = parent.base_ring().one()
        return parent(self * d / n)

    def _normalize_coefficients(self, c):
        """
        If our coefficient ring is the field of fractions over a univariate
        polynomial ring over the rationals, then we should clear both the
        numerator and denominator of the denominators of their
        coefficients.
        
        EXAMPLES::
        
            sage: P = JackPolynomialsP(QQ)
            sage: t = P.base_ring().gen()
            sage: a = 2/(1/2*t+1/2)
            sage: P._normalize_coefficients(a)
            4/(t + 1)
            sage: a = 1/(1/3+1/6*t)
            sage: P._normalize_coefficients(a)
            6/(t + 2)
            sage: a = 24/(4*t^2 + 12*t + 8)
            sage: P._normalize_coefficients(a)
            6/(t^2 + 3*t + 2)
        """
        BR = self.base_ring()
        if is_FractionField(BR) and BR.base_ring() == QQ:
            denom = c.denominator()
            numer = c.numerator()

            #Clear the denominators
            a = lcm([i.denominator() for i in denom.coeffs()])
            b = lcm([i.denominator() for i in numer.coeffs()])
            l = Integer(a).lcm(Integer(b))
            denom *= l
            numer *= l

            #Divide through by the gcd of the numerators
            a = gcd([i.numerator() for i in denom.coeffs()])
            b = gcd([i.numerator() for i in numer.coeffs()])
            l = Integer(a).gcd(Integer(b))

            denom = denom / l
            numer = numer / l
            
            return c.parent()(numer, denom)
        else:
            return c

    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 WP(self, *q, **kw):
        # Specific keyword arguments instead of **kw would be preferable,
        # later versions of Python might support specific (optional) keyword
        # arguments after *q.
        r"""
        Construct weighted projective `n`-space over a field.

        INPUT:

        - ``q`` -- a sequence of positive integers relatively prime to
          one another. The weights ``q`` can be given either as a list
          or tuple, or as positional arguments.

        Two keyword arguments:

        - ``K`` -- a field (default: `\QQ`).
        - ``names`` -- string or list (tuple) of strings (default 'z+'). See
          :func:`~sage.schemes.toric.variety.normalize_names` for
          acceptable formats.

        OUTPUT:

        - A :class:`toric variety
          <sage.schemes.toric.variety.ToricVariety_field>`.
          If `q=(q_0,\dots,q_n)`, then the output is the weighted projective
          space `\mathbb{P}(q_0,\dots,q_n)` over `K`. ``names`` are the names
          of the generators of the homogeneous coordinate ring.

        EXAMPLES:

        A hyperelliptic curve `C` of genus 2 as a subscheme of the weighted
        projective plane `\mathbb{P}(1,3,1)`::

            sage: X = toric_varieties.WP([1,3,1], names='x y z')
            sage: X.inject_variables()
            Defining x, y, z
            sage: g = y^2-(x^6-z^6)
            sage: C = X.subscheme([g]); C
            Closed subscheme of 2-d toric variety covered by 3 affine patches defined by:
              -x^6 + z^6 + y^2
        """
        if len(q) == 1:
            # tuples and lists of weights are acceptable input
            if isinstance(q[0], (list, tuple)):
                q = q[0]
        q = list(q)
        m = len(q)
        # allow case q=[1]? (not allowed presently)
        if m < 2:
            raise ValueError("more than one weight must be provided (got %s)" % q)
        for i in range(m):
            try:
                q[i] = ZZ(q[i])
            except (TypeError):
                raise TypeError("the weights (=%s) must be integers" % q)
            if q[i] <= 0:
                raise ValueError("the weights (=%s) must be positive integers" % q)
        if not gcd(q) == 1:
            raise ValueError("the weights (=%s) must be relatively prime" % q)

        # set default values for K and names
        K = QQ
        names = "z+"
        for key in kw:
            if key == "K":
                K = kw["K"]
                if K not in _Fields:
                    raise TypeError("K (=%r) must be a field" % K)
            elif key == "names":
                names = kw["names"]
                names = normalize_names(names, m, DEFAULT_PREFIX)
            else:
                raise TypeError("got an unexpected keyword argument %r" % key)

        L = ToricLattice(m)
        L_sub = L.submodule([L(q)])
        Q = L / L_sub
        rays = []
        cones = []
        w = range(m)
        L_basis = L.basis()
        for i in w:
            b = L_basis[i]
            v = Q.coordinate_vector(Q(b))
            rays = rays + [v]
            w_c = w[:i] + w[i + 1 :]
            cones = cones + [tuple(w_c)]
        fan = Fan(cones, rays)
        return ToricVariety(fan, coordinate_names=names, base_field=K)

    def _compute_lattice(self, rational_only=False, rational_subgroup=False):
        r"""
        Return a list of vectors that define elements of the rational
        homology that generate this finite subgroup.
        
        INPUT:
        
        
        -  ``rational_only`` - bool (default: False); if
           ``True``, only use rational cusps.
        
        
        OUTPUT:
        
        
        -  ``list`` - list of vectors
        
        
        EXAMPLES::
        
            sage: J = J0(37)
            sage: C = sage.modular.abvar.cuspidal_subgroup.CuspidalSubgroup(J)
            sage: C._compute_lattice()
            Free module of degree 4 and rank 4 over Integer Ring
            Echelon basis matrix:
            [  1   0   0   0]
            [  0   1   0   0]
            [  0   0   1   0]
            [  0   0   0 1/3]
            sage: J = J0(43)
            sage: C = sage.modular.abvar.cuspidal_subgroup.CuspidalSubgroup(J)
            sage: C._compute_lattice()
            Free module of degree 6 and rank 6 over Integer Ring
            Echelon basis matrix:
            [  1   0   0   0   0   0]
            [  0 1/7   0 6/7   0 5/7]
            [  0   0   1   0   0   0]
            [  0   0   0   1   0   0]
            [  0   0   0   0   1   0]
            [  0   0   0   0   0   1]
            sage: J = J0(22)
            sage: C = sage.modular.abvar.cuspidal_subgroup.CuspidalSubgroup(J)
            sage: C._compute_lattice()
            Free module of degree 4 and rank 4 over Integer Ring
            Echelon basis matrix:
            [1/5 1/5 4/5   0]
            [  0   1   0   0]
            [  0   0   1   0]
            [  0   0   0 1/5]
            sage: J = J1(13)
            sage: C = sage.modular.abvar.cuspidal_subgroup.CuspidalSubgroup(J)
            sage: C._compute_lattice()
            Free module of degree 4 and rank 4 over Integer Ring
            Echelon basis matrix:
            [ 1/19     0     0  9/19]
            [    0  1/19  1/19 18/19]
            [    0     0     1     0]
            [    0     0     0     1]
        
        We compute with and without the optional
        ``rational_only`` option.
        
        ::
        
            sage: J = J0(27); G = sage.modular.abvar.cuspidal_subgroup.CuspidalSubgroup(J)
            sage: G._compute_lattice()
            Free module of degree 2 and rank 2 over Integer Ring
            Echelon basis matrix:
            [1/3   0]
            [  0 1/3]
            sage: G._compute_lattice(rational_only=True)
            Free module of degree 2 and rank 2 over Integer Ring
            Echelon basis matrix:
            [1/3   0]
            [  0   1]
        """
        A = self.abelian_variety()
        Cusp = A.modular_symbols()
        Amb  = Cusp.ambient_module()
        Eis  = Amb.eisenstein_submodule()

        C = Amb.cusps()
        N = Amb.level()
        
        if rational_subgroup:
            # QQ-rational subgroup of cuspidal subgroup
            assert A.is_ambient()
            Q = Cusp.abvarquo_rational_cuspidal_subgroup()
            return Q.V()
        
        if rational_only:
            # subgroup generated by differences of rational cusps
            if not is_Gamma0(A.group()):
                raise NotImplementedError, 'computation of rational cusps only implemented in Gamma0 case.'
            if not N.is_squarefree():
                data = [n for n in range(2,N) if gcd(n,N) == 1]
                C = [c for c in C if is_rational_cusp_gamma0(c, N, data)]

        v = [Amb([infinity, alpha]).element() for alpha in C]
        cusp_matrix = matrix(QQ, len(v), Amb.dimension(), v)
#.........这里部分代码省略.........

    def _subpoly_parallel_facets(self):
        """
        Generator for all lattice sub-polyhedra with parallel facets.

        In a sub-polyhedron `Y\subset X` not all edges of `Y` need to
        be parallel to `X`. This method iterates over all
        sub-polyhedra where they are parallel, up to an overall
        translation of the sub-polyhedron. Degenerate sub-polyhedra of
        dimension strictly smaller are included.

        OUTPUT:

        A generator yielding `\ZZ`-polyhedra. By construction, each
        facet of the returned polyhedron is parallel to one of the
        facets of ``self``.

        EXAMPLES::

            sage: X = Polyhedron(vertices=[(0,0), (0,1), (1,0), (1,1)])
            sage: X._subpoly_parallel_facets()
            <generator object _subpoly_parallel_facets at 0x...>
            sage: for p in X._subpoly_parallel_facets():
            ...       print p.Vrepresentation()
            (A vertex at (0, 0),)
            (A vertex at (0, -1), A vertex at (0, 0))
            (A vertex at (-1, 0), A vertex at (0, 0))
            (A vertex at (-1, -1), A vertex at (-1, 0), A vertex at (0, -1), A vertex at (0, 0))

        TESTS::

            sage: X = Polyhedron(vertices=[(0,), (3,)])
            sage: [ p.vertices() for p in X._subpoly_parallel_facets() ]
            [(A vertex at (0),),
             (A vertex at (-1), A vertex at (0)),
             (A vertex at (-2), A vertex at (0)),
             (A vertex at (-3), A vertex at (0))]
            sage: list( Polyhedron(vertices=[[0,0]])._subpoly_parallel_facets() )
            [A 0-dimensional polyhedron in ZZ^2 defined as the convex hull of 1 vertex]
            sage: list( Polyhedron()._subpoly_parallel_facets() )
            [The empty polyhedron in ZZ^0]
        """
        if self.dim()>2 or not self.is_compact():
            raise NotImplementedError('only implemented for bounded polygons')
        from sage.geometry.polyhedron.plot import cyclic_sort_vertices_2d
        vertices = cyclic_sort_vertices_2d(self.vertices())
        n = len(vertices)
        if n==1:  # single point
            yield self
            return
        edge_vectors = []
        for i in range(0,n):
            v = vertices[(i+1) % n].vector() - vertices[i].vector()
            d = gcd(list(v))
            v_prim = (v/d).change_ring(ZZ)
            edge_vectors.append([ v_prim*i for i in range(d+1) ])
        origin = self.ambient_space().zero()
        parent = self.parent()
        from sage.combinat.cartesian_product import CartesianProduct
        for edges in CartesianProduct(*edge_vectors):
            v = []
            point = origin
            for e in edges:
                point += e
                v.append(point)
            if point!=origin:   # does not close up, not a subpolygon
                continue
            yield parent([v, [], []], None)

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 dimension__jacobi_scalar_f(k, m, f) :
    if moebius(f) != (-1)**k :
        return 0
    
    ## We use chapter 6 of Skoruppa's thesis
    ts = filter(lambda t: gcd(t, m // t) == 1, m.divisors())
    
    ## Eisenstein part
    eis_dimension = 0
    
    for t in ts :
        eis_dimension +=   moebius(gcd(m // t, f)) \
                        * (t // t.squarefree_part()).isqrt() \
                        * (2 if (m // t) % 4 == 0 else 1) 
    eis_dimension = eis_dimension // len(ts)
    
    if k == 2 and f == 1 :
        eis_dimension -= len( (m // m.squarefree_part()).isqrt().divisors() ) 
    
    ## Cuspidal part
    cusp_dimension = 0
    
    tmp = ZZ(0)
    for t in ts :
        tmp += moebius(gcd(m // t, f)) * t
    tmp = tmp / len(ts)
    cusp_dimension += tmp * (2 * k - 3) / ZZ(12)
    print "1: ", cusp_dimension
    
    if m % 2 == 0 :
        tmp = ZZ(0)
        for t in ts :
            tmp += moebius(gcd(m // t, f)) * kronecker_symbol(-4, t)
        tmp = tmp / len(ts)
        
        cusp_dimension += 1/ZZ(2) * kronecker_symbol(8, 2 * k - 1) * tmp
        print "2: ", 1/ZZ(2) * kronecker_symbol(8, 2 * k - 1) * tmp
        
    tmp = ZZ(0)
    for t in ts :
        tmp += moebius(gcd(m // t, f)) * kronecker_symbol(t, 3)
    tmp = tmp / len(ts)
    if m % 3 != 0 :
        cusp_dimension += 1 / ZZ(3) * kronecker_symbol(k, 3) * tmp
        print ": ", 1 / ZZ(3) * kronecker_symbol(k, 3) * tmp
    elif k % 3 == 0 :
        cusp_dimension += 2 / ZZ(3) * (-1)**k * tmp
        print "3: ", 2 / ZZ(3) * (-1)**k * tmp
    else :
        cusp_dimension += 1 / ZZ(3) * (kronecker_symbol(k, 3) + (-1)**(k - 1)) * tmp
        print "3: ", 1 / ZZ(3) * (kronecker_symbol(k, 3) + (-1)**(k - 1)) * tmp
    
    tmp = ZZ(0)
    for t in ts :
        tmp +=   moebius(gcd(m // t, f)) \
               * (t // t.squarefree_part()).isqrt() \
               * (2 if (m // t) % 4 == 0 else 1)
    tmp = tmp / len(ts)
    cusp_dimension -= 1 / ZZ(2) * tmp
    print "4: ", -1 / ZZ(2) * tmp
    
    tmp = ZZ(0)
    for t in ts :
        tmp +=   moebius(gcd(m // t, f)) \
               * sum(   (( len(BinaryQF_reduced_representatives(-d, True))
                           if d not in [3, 4] else ( 1 / ZZ(3) if d == 3 else 1 / ZZ(2) ))
                         if d % 4 == 0 or d % 4 == 3 else 0 )
                      * kronecker_symbol(-d, m // t)
                      * ( 1 if (m // t) % 2 != 0 else
                          ( 4 if (m // t) % 4 == 0 else 2 * kronecker_symbol(-d, 2) ))
                      for d in (4 * m).divisors() )
    tmp = tmp / len(ts)
    cusp_dimension -= 1 / ZZ(2) * tmp
    print "5: ", -1 / ZZ(2) * tmp
    
    if k == 2 :
        cusp_dimension += len( (m // f // (m // f).squarefree_part()).isqrt().divisors() )
    
    return eis_dimension + cusp_dimension

def charpolys(v, B, filename=None):
    """
    Compute characteristic polynomials of T_P for primes P with norm
    <= B coprime to the level, for all spaces of Hilbert modular forms
    for all the levels in v.

    INPUT:

    - `v` -- list of positive integers
    - `B` -- positive integer
    - ``filename`` -- optional string; if given, output is also written
      to that file (in addition to stdout).

    OUTPUT:

    - outputs a table with rows corresponding to the ideals
      of Q(sqrt(5)) with norm in v, and optionally creates a file

    EXAMPLES::

        sage: from sage.modular.hilbert.sqrt5_tables import charpolys
        sage: out = charpolys([1..20], 10)
        4 2 ... [(5,x-6),(3,x-10)]
        5 -2*a+1 ... [(2,x-5),(3,x-10)]
        9 3 ... [(2,x-5),(5,x-6)]
        11 -3*a+1 ... [(2,x-5),(5,x-6),(3,x-10)]
        11 -3*a+2 ... [(2,x-5),(5,x-6),(3,x-10)]
        16 4 ... [(5,x-6),(3,x-10)]
        19 -4*a+1 ... [(2,x-5),(5,x-6),(3,x-10)]
        19 -4*a+3 ... [(2,x-5),(5,x-6),(3,x-10)]
        20 -4*a+2 ... [(3,x-10)]

        sage: out = charpolys([20, 11], 10)
        20 -4*a+2 ... [(3,x-10)]
        11 -3*a+1 ... [(2,x-5),(5,x-6),(3,x-10)]
        11 -3*a+2 ... [(2,x-5),(5,x-6),(3,x-10)]

    Test writing to a file::

        sage: if os.path.exists('tmp_table.txt'): os.unlink('tmp_table.txt')
        sage: out = charpolys([20, 11], 10, 'tmp_table.txt')
        20 -4*a+2 ... [(3,x-10)]
        11 -3*a+1 ... [(2,x-5),(5,x-6),(3,x-10)]
        11 -3*a+2 ... [(2,x-5),(5,x-6),(3,x-10)]
        sage: r = open('tmp_table.txt').read()
        sage: 'x-10' in r
        True
        sage: r.count('\n')
        3
        sage: os.unlink('tmp_table.txt')    
    """
    if len(v) == 0:
        return ""
    out = ""
    F = open(filename, "a") if filename else None
    P = [p for p in ideals_of_bounded_norm(B) if p.is_prime()]
    for N in ideals_of_norm(v):
        t = cputime()
        H = IcosiansModP1ModN(N)
        T = [
            (p.smallest_integer(), H.hecke_matrix(p).fcp()) for p in P if gcd(Integer(p.norm()), Integer(N.norm())) == 1
        ]
        tm = "%.2f" % cputime(t)
        s = "%s %s %s %s" % (N.norm(), no_space(reduced_gen(N)), tm, no_space(T))
        print s
        out += s + "\n"
        if F:
            F.write(s + "\n")
            F.flush()
    return out

def saturation(A, proof=True, p=0, max_dets=5):
    """
    Compute a saturation matrix of A.

    INPUT:

    - A     -- a matrix over ZZ
    - proof -- bool (default: True)
    - p     -- int (default: 0); if not 0 only guarantees that output is
      p-saturated
    - max_dets -- int (default: 4) max number of dets of submatrices to
      compute.

    OUTPUT:

    matrix -- saturation of the matrix A.

    EXAMPLES::

        sage: from sage.matrix.matrix_integer_dense_saturation import saturation
        sage: A = matrix(ZZ, 2, 2, [3,2,3,4]); B = matrix(ZZ, 2,3,[1,2,3,4,5,6]); C = A*B
        sage: C
        [11 16 21]
        [19 26 33]
        sage: C.index_in_saturation()
        18
        sage: S = saturation(C); S
        [11 16 21]
        [-2 -3 -4]
        sage: S.index_in_saturation()
        1
        sage: saturation(C, proof=False)
        [11 16 21]
        [-2 -3 -4]
        sage: saturation(C, p=2)
        [11 16 21]
        [-2 -3 -4]
        sage: saturation(C, p=2, max_dets=1)
        [11 16 21]
        [-2 -3 -4]
    """
    # Find a submatrix of full rank and instead saturate that matrix.
    r = A.rank()
    if A.is_square() and r == A.nrows():
        return identity_matrix(ZZ, r)
    if A.nrows() > r:
        P = []
        while len(P) < r:
            P = matrix_integer_dense_hnf.probable_pivot_rows(A)
        A = A.matrix_from_rows(P)

    # Factor out all common factors from all rows, just in case.
    A = copy(A)
    A._factor_out_common_factors_from_each_row()

    if A.nrows() <= 1:
        return A

    A, zero_cols = A._delete_zero_columns()

    if max_dets > 0:
        # Take the GCD of at most num_dets randomly chosen determinants.
        nr = A.nrows()
        nc = A.ncols()
        d = 0
        trials = min(binomial(nc, nr), max_dets)
        already_tried = []
        while len(already_tried) < trials:
            v = random_sublist_of_size(nc, nr)
            tm = verbose("saturation -- checking det condition on submatrix")
            d = gcd(d, A.matrix_from_columns(v).determinant(proof=proof))
            verbose("saturation -- got det down to %s" % d, tm)
            if gcd(d, p) == 1:
                return A._insert_zero_columns(zero_cols)
            already_tried.append(v)

        if gcd(d, p) == 1:
            # already p-saturated
            return A._insert_zero_columns(zero_cols)

        # Factor and p-saturate at each p.
        # This is not a good algorithm, because all the HNF's in it are really slow!
        #
        # tm = verbose('factoring gcd %s of determinants'%d)
        # limit = 2**31-1
        # F = d.factor(limit = limit)
        # D = [p for p, e in F if p <= limit]
        # B = [n for n, e in F if n > limit]  # all big factors -- there will only be at most one
        # assert len(B) <= 1
        # C = B[0]
        # for p in D:
        #    A = p_saturation(A, p=p, proof=proof)

    # This is a really simple but powerful algorithm.
    # FACT: If A is a matrix of full rank, then hnf(transpose(A))^(-1)*A is a saturation of A.
    # To make this practical we use solve_system_with_difficult_last_row, since the
    # last column of HNF's are typically the only really big ones.
    B = A.transpose().hermite_form(include_zero_rows=False, proof=proof)
    B = B.transpose()

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

    def multiple_of_order(self, maxp=None):
        """
        Return a multiple of the order of this torsion group.

        The multiple is computed using characteristic polynomials of Hecke
        operators of odd index not dividing the level.

        INPUT:


        -  ``maxp`` - (default: None) If maxp is None (the
           default), return gcd of best bound computed so far with bound
           obtained by computing GCD's of orders modulo p until this gcd
           stabilizes for 3 successive primes. If maxp is given, just use all
           primes up to and including maxp.


        EXAMPLES::

            sage: J = J0(11)
            sage: G = J.rational_torsion_subgroup()
            sage: G.multiple_of_order(11)
            5
            sage: J = J0(389)
            sage: G = J.rational_torsion_subgroup(); G
            Torsion subgroup of Abelian variety J0(389) of dimension 32
            sage: G.multiple_of_order()
            97
            sage: [G.multiple_of_order(p) for p in prime_range(3,11)]
            [92645296242160800, 7275, 291]
            sage: [G.multiple_of_order(p) for p in prime_range(3,13)]
            [92645296242160800, 7275, 291, 97]
            sage: [G.multiple_of_order(p) for p in prime_range(3,19)]
            [92645296242160800, 7275, 291, 97, 97, 97]

        ::

            sage: J = J0(33) * J0(11) ; J.rational_torsion_subgroup().order()
            Traceback (most recent call last):
            ...
            NotImplementedError: torsion multiple only implemented for Gamma0

        The next example illustrates calling this function with a larger
        input and how the result may be cached when maxp is None::

            sage: T = J0(43)[1].rational_torsion_subgroup()
            sage: T.multiple_of_order()
            14
            sage: T.multiple_of_order(50)
            7
            sage: T.multiple_of_order()
            7
        """
        if maxp is None:
            try:
                return self.__multiple_of_order
            except AttributeError:
                pass
        bnd = ZZ(0)
        A = self.abelian_variety()
        if A.dimension() == 0:
            T = ZZ(1)
            self.__multiple_of_order = T
            return T
        N = A.level()
        if not (len(A.groups()) == 1 and is_Gamma0(A.groups()[0])):
            # to generalize to this case, you'll need to
            # (1) define a charpoly_of_frob function:
            #       this is tricky because I don't know a simple
            #       way to do this for Gamma1 and GammaH.  Will
            #       probably have to compute explicit matrix for
            #       <p> operator (add to modular symbols code),
            #       then compute some charpoly involving
            #       that directly...
            # (2) use (1) -- see my MAGMA code.
            raise NotImplementedError("torsion multiple only implemented for Gamma0")
        cnt = 0
        if maxp is None:
            X = Primes()
        else:
            X = prime_range(maxp+1)
        for p in X:
            if (2*N) % p == 0:
                continue

            f = A.hecke_polynomial(p)
            b = ZZ(f(p+1))

            if bnd == 0:
                bnd = b
            else:
                bnd_last = bnd
                bnd = ZZ(gcd(bnd, b))
                if bnd == bnd_last:
                    cnt += 1
                else:
                    cnt = 0
                if maxp is None and cnt >= 2:
                    break

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