def iter_positive_forms_with_content(self) :
        if self.__disc is infinity :
            raise ValueError, "infinity is not a true filter index"
        
        
        if self.__reduced :        
            for a in xrange(1,isqrt(self.__disc // 3) + 1) :
                for b in xrange(a+1) :
                    g = gcd(a, b)
                    for c in xrange(a, (b**2 + (self.__disc - 1))//(4*a) + 1) :
                        yield (a,b,c), gcd(g,c)
        else :
            maxtrace = floor(5*self.__disc / 15 + sqrt(self.__disc)/2)
            for a in xrange(1, maxtrace + 1) :
                for c in xrange(1, maxtrace - a + 1) :
                    g = gcd(a,c)
                    
                    Bu = isqrt(4*a*c - 1)

                    di = 4*a*c - self.__disc
                    if di >= 0 :
                        Bl = isqrt(di) + 1 
                    else :
                        Bl = 0
                    
                    for b in xrange(-Bu, -Bl + 1) :
                        yield (a,b,c), gcd(g,b)
                    for b in xrange(Bl, Bu + 1) :
                        yield (a,b,c), gcd(g,b)
        #! if self.__reduced

        raise StopIteration

    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
        `\mathbf{P}^1(\QQ)`.  These are returned in a reduced
        form; see self.reduce_cusp for the definition of reduced.

        ALGORITHM:
            Lemma 3.2 in Cremona's 1997 book shows that for the action
            of Gamma1(N) on "signed projective space"
            `\Q^2 / (\Q_{\geq 0}^+)`, we have `u_1/v_1 \sim u_2 / v_2`
            if and only if `v_1 = v_2 \bmod N` and `u_1 = u_2 \bmod
            gcd(v_1, N)`. It follows that every orbit has a
            representative `u/v` with `v \le N` and `0 \le u \le
            gcd(v, N)`.  We iterate through all pairs `(u,v)`
            satisfying this.

            Having found a set containing at least one of every
            equivalence class modulo Gamma1(N), we can be sure of
            picking up every class modulo GammaH(N) since this
            contains Gamma1(N); and the reduce_cusp call does the
            checking to make sure we don't get any duplicates.

        EXAMPLES::

            sage: Gamma1(5)._find_cusps()
            [0, 2/5, 1/2, Infinity]
            sage: Gamma1(35)._find_cusps()
            [0, 2/35, 1/17, 1/16, 1/15, 1/14, 1/13, 1/12, 3/35, 1/11, 1/10, 1/9, 4/35, 1/8, 2/15, 1/7, 1/6, 6/35, 1/5, 3/14, 8/35, 1/4, 9/35, 4/15, 2/7, 3/10, 11/35, 1/3, 12/35, 5/14, 13/35, 2/5, 3/7, 16/35, 17/35, 1/2, 8/15, 4/7, 3/5, 9/14, 7/10, 5/7, 11/14, 4/5, 6/7, 9/10, 13/14, Infinity]
            sage: Gamma1(24)._find_cusps() == Gamma1(24).cusps(algorithm='modsym')
            True
            sage: GammaH(24, [13,17])._find_cusps() == GammaH(24,[13,17]).cusps(algorithm='modsym')
            True
        """

        s = []
        hashes = []
        N = self.level()

        for d in xrange(1, 1+N):
            w = N.gcd(d)
            M = int(w) if w > 1 else 2
            for a in xrange(1,M):
                if gcd(a, w) != 1:
                    continue
                while gcd(a, d) != 1:
                    a += w
                c = self.reduce_cusp(Cusp(a,d))
                h = hash(c)
                if not h in hashes:
                    hashes.append(h)
                    s.append(c)
        return sorted(s)

def num_cusps_of_width(N, d):
    r"""
    Return the number of cusps on `X_0(N)` of width d.

    INPUT:


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

    -  ``d`` - (integer): an integer dividing N, the cusp
       width


    EXAMPLES::

        sage: [num_cusps_of_width(18,d) for d in divisors(18)]
        [1, 1, 2, 2, 1, 1]
    """
    try:
        N = ZZ(N)
        d = ZZ(d)
        assert N>0
        assert d>0
        assert ((N % d) == 0)
    except TypeError:
        raise TypeError, "N and d must be integers"
    except AssertionError:
        raise AssertionError, "N and d must be positive integers with d|N"

    return euler_phi(gcd(d, N//d))

def _normalize_H(H, level):
    """
    Normalize representatives for a given subgroup H of the units
    modulo level.

    NOTE: This function does *not* make any attempt to find a minimal
    set of generators for H. It simply normalizes the inputs for use
    in hashing.

    EXAMPLES::

        sage: sage.modular.arithgroup.congroup_gammaH._normalize_H([23], 10)
        [3]
        sage: sage.modular.arithgroup.congroup_gammaH._normalize_H([1,5], 7)
        [5]
        sage: sage.modular.arithgroup.congroup_gammaH._normalize_H([4,18], 14)
        Traceback (most recent call last):
        ...
        ArithmeticError: The generators [4, 18] must be units modulo 14
        sage: sage.modular.arithgroup.congroup_gammaH._normalize_H([3,17], 14)
        [3]
        sage: sage.modular.arithgroup.congroup_gammaH._normalize_H([-1,7,9], 10)
        [7, 9]
    """
    H = [ZZ(h) for h in H]
    for h in H:
        if gcd(h, level) > 1:
            raise ArithmeticError, 'The generators %s must be units modulo %s'%(H, level)
    H = list(set([h%level for h in H]))
    H.sort()
    if 1 in H:
        H.remove(1)
    return H

def _list_subgroup(N, gens):
    r"""
    Given an integer ``N`` and a list of integers ``gens``, return a list of
    the elements of the subgroup of `(\ZZ / N\ZZ)^\times` generated by the
    elements of ``gens``.

    EXAMPLE::

        sage: sage.modular.arithgroup.congroup_gammaH._list_subgroup(11, [3])
        [1, 3, 4, 5, 9]
    """
    H = set([1])
    N = int(N)
    for g in gens:
        if gcd(g, N) != 1:
            raise ValueError, "gen (=%s) is not in (Z/%sZ)^*"%(g,N)
        gk = int(g) % N
        sbgrp = [gk]
        while not (gk in H):
            gk = (gk * g)%N
            sbgrp.append(gk)
        H = set([(x*h)%N for x in sbgrp for h in H])
    H = list(H)
    H.sort()
    return H

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 order_at_cusp(self, cusp):
        r"""
        Return the order of vanishing of self at the given cusp.

        INPUT:


        -  ``cusp`` -  a CuspFamily object


        OUTPUT:


        - an integer


        EXAMPLES::

            sage: e = EtaProduct(2, {2:24, 1:-24})
            sage: e.order_at_cusp(CuspFamily(2, 1)) # cusp at infinity
            1
            sage: e.order_at_cusp(CuspFamily(2, 2)) # cusp 0
            -1
        """
        if not isinstance(cusp, CuspFamily):
            raise TypeError("Argument (=%s) should be a CuspFamily" % cusp)
        if cusp.level() != self.level():
            raise ValueError("Cusp not on right curve!")
        return 1/ZZ(24)/gcd(cusp.width(), self.level()//cusp.width()) * sum( [ell*self.r(ell)/cusp.width() * (gcd(cusp.width(), self.level()//ell))**2  for ell in self._keys] )

def num_cusps_of_width(N, d):
    r"""
    Return the number of cusps on `X_0(N)` of width d.

    INPUT:


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

    -  ``d`` - (integer): an integer dividing N, the cusp
       width


    EXAMPLES::

        sage: [num_cusps_of_width(18,d) for d in divisors(18)]
        [1, 1, 2, 2, 1, 1]
        sage: num_cusps_of_width(4,8)
        Traceback (most recent call last):
        ...
        ValueError: N and d must be positive integers with d|N
    """
    N = ZZ(N)
    d = ZZ(d)
    if N <= 0 or d <= 0 or (N % d) != 0:
        raise ValueError("N and d must be positive integers with d|N")

    return euler_phi(gcd(d, N//d))

    def _test__jacobi_predicted_taylor_coefficients(fs, q_precision) :
        r"""
        Given a list of power series, which are the corrected Taylor coefficients
        of a Jacobi form, return the renormalized uncorrected ones, assuming that
        all but one `f` vanish.
        
        INPUT:
        
        - ``fs`` -- A list of power series.
        
        - ``q_precision`` -- An integer.
        
        OUPUT:
        
        - A list of power series.
        
        TESTS:
        
        See jacobi_form_by_taylor_expansion.
        """
        from sage.rings.arith import gcd
        
        R = PowerSeriesRing(ZZ, 'q'); q = R.gen(0)
        
        diff = lambda f: f.derivative().shift(1)
        normalize = lambda f: f / gcd(f.list()) if f != 0 else f
        diffnorm = lambda f,l: normalize(reduce(lambda a, g: g(a), l*[diff], f))

        taylor_coefficients = list()
        allf = R(0)
        for f in fs :
            allf = f(q_precision) + diffnorm(allf, 1)            
            taylor_coefficients.append(allf)

        return taylor_coefficients

    def random_element(self, bound=100, *args, **kwds):
        r"""
        Return a random element of `{\rm SL}_2(\ZZ)` with entries whose
        absolute value is strictly less than bound (default 100).
        Additional arguments and keywords are passed to the random_element
        method of ZZ.

        (Algorithm: Generate a random pair of integers at most bound. If they
        are not coprime, throw them away and start again. If they are, find an
        element of `{\rm SL}_2(\ZZ)` whose bottom row is that, and
        left-multiply it by `\begin{pmatrix} 1 & w \\ 0 & 1\end{pmatrix}` for
        an integer `w` randomly chosen from a small enough range that the
        answer still has entries at most bound.)

        It is, unfortunately, not true that all elements of SL2Z with entries <
        bound appear with equal probability; those with larger bottom rows are
        favoured, because there are fewer valid possibilities for w.

        EXAMPLES::

            sage: SL2Z.random_element()
            [60 13]
            [83 18]
            sage: SL2Z.random_element(5)
            [-1  3]
            [ 1 -4]

        Passes extra positional or keyword arguments through::

            sage: SL2Z.random_element(5, distribution='1/n')
            [ 1 -4]
            [ 0  1]
        """
        if bound <= 1: raise ValueError("bound must be greater than 1")
        c = ZZ.random_element(1-bound, bound, *args, **kwds)
        d = ZZ.random_element(1-bound, bound, *args, **kwds)
        if gcd(c,d) != 1: # try again
            return self.random_element(bound, *args, **kwds)
        else:
            a,b,c,d = lift_to_sl2z(c,d,0)
            whi = bound
            wlo = bound
            if c > 0:
                whi = min(whi, ((bound - a)/ZZ(c)).ceil())
                wlo = min(wlo, ((bound + a)/ZZ(c)).ceil())
            elif c < 0:
                whi = min(whi, ((bound + a)/ZZ(-c)).ceil())
                wlo = min(wlo, ((bound - a)/ZZ(-c)).ceil())

            if d > 0:
                whi = min(whi, ((bound - b)/ZZ(d)).ceil())
                wlo = min(wlo, ((bound + b)/ZZ(d)).ceil())
            elif d < 0:
                whi = min(whi, ((bound + b)/ZZ(-d)).ceil())
                wlo = min(wlo, ((bound - b)/ZZ(-d)).ceil())

            w = ZZ.random_element(1-wlo, whi, *args, **kwds)
            a += c*w
            b += d*w
            return self([a,b,c,d])

 def eval(self, expansion, weight = None) :
     if weight is None :
         try :
             weight = expansion.weight()
         except AttributeError :
             raise ValueError, "weight must be defined for the Hecke action"
     
     precision = expansion.precision()
     if precision.is_infinite() :
         precision = expansion._bounding_precision()
     else :
         precision = precision._hecke_operator(self.__l)
     characters = expansion.non_zero_components()
     expansion_level = precision.level()
     
     if gcd(expansion_level, self.__l) != 1 :
         raise ValueError, "Level of expansion and of Hecke operator must be coprime."
     
     hecke_expansion = dict()
     for ch in characters :
         hecke_expansion[ch] = dict( (k, self.hecke_coeff(expansion, ch, k, weight, expansion_level))
                                     for k in precision )
     
     result = expansion.parent()._element_constructor_(hecke_expansion)
     result._set_precision(expansion.precision()._hecke_operator(self.__l))
     
     return result

    def get_representatives( self, t, N) :
        r"""
        A helper function used in hecke_coeff that computes the right
        coset representatives of $\Gamma^0(t) \cap \Gamma_0(N) \ Gamma_0(N)$ where
        $\Gamma^0(t)$ is the subgroup of $SL(2,Z)$ where the upper right hand
        corner is divisible by $t$.

        NOTE
            We use the bijection $\Gamma^0(t)\SL(2,Z) \rightarrow P^1(\Z/t\Z)$
            given by $A \mapsto [1:0]A$.
        """
        if t == 1 : return [(1,0,0,1)]
        
        rep_list = []
        
        for (x,y) in P1List(t):
            ## we know that (N, x, y) = 1
            N1 = gcd(N, x)
            if N1 != 1 :
                x = x + t * N // N1

            ## we calculate a pair c,d satisfying a minimality condition
            ## to make later multiplications cheaper
            (_, d, c) = Integer(x)._xgcd(Integer(N * y), minimal=True)
            
            #print (x, y, -N * c, d)
            rep_list.append((x, y, -N * c, d))
                
        return rep_list

def circle_drops(A,B):
    # Drops going around the unit circle for those A and B.
    # See http://user.math.uzh.ch/dehaye/thesis_students/Nicolas_Wider-Integrality_of_factorial_ratios.pdf
    # for longer description (not original, better references exist)    
    marks = lcm(lcm(A),lcm(B))    
    tmp = [0 for i in range(marks)]
    for a in A:
#        print tmp
        for i in range(a):
            if gcd(i, a) == 1:
                tmp[i*marks/a] -= 1
    for b in B:
#        print tmp
        for i in range(b):
            if gcd(i, b) == 1:
                tmp[i*marks/b] += 1
#    print tmp
    return [sum(tmp[:j]) for j in range(marks)]

def blift(LF,Li,p,S=None):
    r"""
    Search for a solution to the given list of inequalities. If found, lift the solution to
    an appropriate valuation. See Lemma 3.3.6 in [Molnar, M.Sc. thesis]

    INPUT:

    - ``LF`` -- a list of integer polynomials in one variable (the normalized coefficients)

    - ``Li`` -- an integer, the bound on coefficients

    - ``p`` -- a prime

    OUTPUT:

    - Boolean -- whether or not the lift is successful

    - integer -- the lift

    EXAMPLES::

        sage: R.<b> = PolynomialRing(QQ)
        sage: from sage.schemes.projective.endPN_minimal_model import blift
        sage: blift([8*b^3 + 12*b^2 + 6*b + 1, 48*b^2 + 483*b + 117, 72*b + 1341, -24*b^2 + 411*b + 99, -144*b + 1233, -216*b], 2, 3)
        (True, 4)
    """

    P = LF[0].parent()
    #Determine which inequalities are trivial, and scale the rest, so that we only lift
    #as many times as needed.
    keepScaledIneqs = [scale(P(coeff),Li,p) for coeff in LF if coeff != 0]
    keptVals = [i[2] for i in keepScaledIneqs if i[0]]
    if keptVals != []:
        #Determine the valuation to lift until.
        liftval = max(keptVals)
    else:
        #All inequalities are satisfied.
        return True,1
    if S is None:
        S = PolynomialRing(Zmod(p),'b')
    keptScaledIneqs = [S(i[1]) for i in keepScaledIneqs if i[0]]
    #We need a solution for each polynomial on the left hand side of the inequalities,
    #so we need only find a solution for their gcd.
    g = gcd(keptScaledIneqs)
    rts = g.roots(multiplicities=False)
    for r in rts:
        #Recursively try to lift each root
        r_initial=QQ(r)
        newInput = P([r_initial,p])
        LG = [F(newInput) for F in LF]
        lift,lifted = blift(LG,Li,p,S=S)
        if lift:
            #Lift successful.
            return True,r_initial+ p*lifted
    #Lift non successful.
    return False,0

 def __contains__(self, f) :
     if self.__disc is infinity :
         return True
     
     (a, b, c) = f
     disc = 4*a*c - b**2
     if disc == 0 :
         return gcd([a, b, c]) < self._indefinite_content_bound()
     else :
         return disc < self.__disc

 def arith_prod_of_partitions(l1, l2):
     # Given two partitions `l_1` and `l_2`, we construct a new partition `l_1 \\boxtimes l_2` by
     # the following procedure: each pair of parts `a \\in l_1` and `b \\in l_2` contributes
     # `\\gcd (a, b)`` parts of size `\\lcm (a, b)` to `l_1 \\boxtimes l_2`. If `l_1` and `l_2`
     # are partitions of integers `n` and `m`, respectively, then `l_1 \\boxtimes l_2` is a
     # partition of `nm`.
     term_iterable = chain.from_iterable( repeat(lcm(pair), times=gcd(pair)) for pair in product(l1, l2) )
     term_list = sorted(term_iterable, reverse=True)
     res = Partition(term_list)
     return res

    def _new_group_from_level(self, level):
        r"""
        Return a new group of the same type (Gamma0, Gamma1, or
        GammaH) as self of the given level. In the case that self is of type
        GammaH, we take the largest H inside `(\ZZ/ \text{level}\ZZ)^\times`
        which maps to H, namely its inverse image under the natural reduction
        map.

        EXAMPLES::

            sage: G = Gamma0(20)
            sage: G._new_group_from_level(4)
            Congruence Subgroup Gamma0(4)
            sage: G._new_group_from_level(40)
            Congruence Subgroup Gamma0(40)

            sage: G = Gamma1(10)
            sage: G._new_group_from_level(6)
            Traceback (most recent call last):
            ...
            ValueError: one level must divide the other

            sage: G = GammaH(50,[7]); G
            Congruence Subgroup Gamma_H(50) with H generated by [7]
            sage: G._new_group_from_level(25)
            Congruence Subgroup Gamma_H(25) with H generated by [7]
            sage: G._new_group_from_level(10)
            Congruence Subgroup Gamma0(10)
            sage: G._new_group_from_level(100)
            Congruence Subgroup Gamma_H(100) with H generated by [7, 57]
        """
        from congroup_gamma0 import is_Gamma0
        from congroup_gamma1 import is_Gamma1
        from congroup_gammaH import is_GammaH
        from all import Gamma0, Gamma1, GammaH

        N = self.level()
        if (level % N) and (N % level):
            raise ValueError, "one level must divide the other"
        if is_Gamma0(self):
            return Gamma0(level)
        elif is_Gamma1(self):
            return Gamma1(level)
        elif is_GammaH(self):
            H = self._generators_for_H()
            if level > N:
                d = level // N
                diffs = [N * i for i in range(d)]
                newH = [h + diff for h in H for diff in diffs]
                return GammaH(level, [x for x in newH if gcd(level, x) == 1])
            else:
                return GammaH(level, [h % level for h in H])
        else:
            raise NotImplementedError

    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 __init__(self, parent, polys, check=True):
        """
        The Python constructor.

        See :class:`SchemeMorphism_polynomial` for details.

        EXAMPLES::

            sage: A2.<x,y> = AffineSpace(QQ,2)
            sage: H = A2.Hom(A2)
            sage: H([x-y, x*y])
            Scheme endomorphism of Affine Space of dimension 2 over Rational Field
              Defn: Defined on coordinates by sending (x, y) to
                    (x - y, x*y)
        """
        if check:
            if not isinstance(polys, (list, tuple)):
                raise TypeError, "polys (=%s) must be a list or tuple"%polys
            source_ring = parent.domain().coordinate_ring()
            target = parent.codomain().ambient_space()
            if len(polys) != target.ngens():
                raise ValueError, "there must be %s polynomials"%target.ngens()
            try:
                polys = [source_ring(poly) for poly in polys]
            except TypeError:
                raise TypeError, "polys (=%s) must be elements of %s"%(polys,source_ring)
            from sage.rings.quotient_ring import QuotientRing_generic
            if isinstance(source_ring, QuotientRing_generic):
                lift_polys = [f.lift() for f in polys]
            else:
                lift_polys = polys
            from sage.schemes.generic.projective_space import is_ProjectiveSpace
            if is_ProjectiveSpace(target):
                # if the codomain is a subscheme of projective space,
                # then we need to make sure that polys have no common
                # zeros 
                if isinstance(source_ring, QuotientRing_generic):
                    # if the coordinate ring of the domain is a
                    # quotient by an ideal, we need to check that the
                    # gcd of polys and the generators of the ideal is 1
                    gcd_polys = lift_polys + list(source_ring.defining_ideal().gens())
                else:
                    # if the domain is affine space, we just need to
                    # check the gcd of polys
                    gcd_polys = polys
                from sage.rings.arith import gcd
                if gcd(gcd_polys) != 1:
                    raise ValueError, "polys (=%s) must not have common factors"%polys
            polys = Sequence(lift_polys)
            polys.set_immutable()
            # Todo: check that map is well defined (how?)
        self.__polys = polys
        SchemeMorphism.__init__(self, parent)

    def eval_twisted_symbol_on_Da(self, a): # rename! should this be in modsym?
        """
        Returns `\Phi_{\chi}(\{a/p}-{\infty})` where `Phi` is the OMS and
        `\chi` is a the quadratic character corresponding to self

        INPUT:
            - ``a`` -- integer in range(p)

        OUTPUT:

        The distribution `\Phi_{\chi}(\{a/p\}-\{\infty\})`.

        EXAMPLES:
            sage: from sage.modular.pollack_stevens.space import ps_modsym_from_elliptic_curve
            sage: E = EllipticCurve('57a')
            sage: p = 5
            sage: prec = 4
            sage: phi = ps_modsym_from_elliptic_curve(E)
            sage: ap = phi.Tq_eigenvalue(p,prec)
            sage: Phi = phi.p_stabilize_and_lift(p,ap = ap, M = prec, algorithm='stevens')
            sage: L = pAdicLseries(Phi)
            sage: L.eval_twisted_symbol_on_Da(1)
            (2 + 2*5 + 2*5^2 + 2*5^3 + O(5^4), 2 + 3*5 + 2*5^2 + O(5^3), 4*5 + O(5^2), 3 + O(5))

            sage: from sage.modular.pollack_stevens.space import ps_modsym_from_elliptic_curve
            sage: E = EllipticCurve('40a4')
            sage: p = 7
            sage: prec = 4
            sage: phi = ps_modsym_from_elliptic_curve(E)
            sage: ap = phi.Tq_eigenvalue(p,prec)
            sage: Phi = phi.p_stabilize_and_lift(p,ap = ap, M = prec, algorithm='stevens')
            sage: L = pAdicLseries(Phi)
            sage: L.eval_twisted_symbol_on_Da(1)
            (4 + 6*7 + 3*7^2 + O(7^4), 2 + 7 + O(7^3), 4 + 6*7 + O(7^2), 6 + O(7))

        """
        symb = self.symb()
        p = symb.parent().prime()
        S0p = Sigma0(p)
        Dists = symb.parent().coefficient_module()
        M = Dists.precision_cap()
        p = Dists.prime()
        twisted_dist = Dists.zero_element()
        m_map = symb._map
        D = self._quadratic_twist
        for b in range(1, abs(D) + 1):
            if gcd(b, D) == 1:
                M1 = S0p([1, (b / abs(D)) % p**M, 0, 1])
                new_dist = m_map(M1 * M2Z([a, 1, p, 0]))*M1
                new_dist = new_dist.scale(kronecker(D, b)).normalize()
                twisted_dist = twisted_dist + new_dist
                #ans = ans + self.eval(M1 * M2Z[a, 1, p, 0])._right_action(M1)._lmul_(kronecker(D, b)).normalize()
        return twisted_dist.normalize()