def are_equivalent(self, x, y, trans = False):
        r""" 
        Test whether or not cusps x and y are equivalent modulo self.  If self
        has a reduce_cusp() method, use that; otherwise do a slow explicit
        test. 

        If trans = False, returns True or False. If trans = True, then return
        either False or an element of self mapping x onto y.

        EXAMPLE::

            sage: Gamma0(7).are_equivalent(Cusp(1/3), Cusp(0), trans=True)
            [  3  -1]
            [-14   5]
            sage: Gamma0(7).are_equivalent(Cusp(1/3), Cusp(1/7))
            False
        """
        x = Cusp(x)
        y = Cusp(y)
        if not trans:
            try:
                xr = self.reduce_cusp(x)
                yr = self.reduce_cusp(y)
                if xr != yr:
                    return False
                if xr == yr:
                    return True
            except NotImplementedError:
                pass

        from all import SL2Z 
    
        vx = lift_to_sl2z(x.numerator(),x.denominator(), 0)
        dx = SL2Z([vx[2], -vx[0], vx[3], -vx[1]])
        vy = lift_to_sl2z(y.numerator(),y.denominator(), 0)
        dy = SL2Z([vy[2], -vy[0], vy[3], -vy[1]])
    
        for i in xrange(self.index()):
            # Note that the width of any cusp is bounded above by the index of self.
            # If self is congruence, then the level of self is a much better bound, but
            # this method is written to work with non-congruence subgroups as well,
            if dy * SL2Z([1,i,0,1])*(~dx) in self:
                if trans:
                    return dy * SL2Z([1,i,0,1]) * ~dx
                else:
                    return True
            elif (self.is_odd() and dy * SL2Z([-1,-i,0,-1]) * ~dx in self):
                if trans:
                    return dy * SL2Z([-1,-i,0,-1]) * ~dx
                else:
                    return True
        return False

    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 coset_reps(self):
        r"""
        Return representatives for the right cosets of this congruence
        subgroup in `{\rm SL}_2(\ZZ)` as a generator object.

        Use ``list(self.coset_reps())`` to obtain coset reps as a
        list.

        EXAMPLES::

            sage: list(Gamma0(5).coset_reps())
            [
            [1 0]  [ 0 -1]  [1 0]  [ 0 -1]  [ 0 -1]  [ 0 -1]
            [0 1], [ 1  0], [1 1], [ 1  2], [ 1  3], [ 1  4]
            ]
            sage: list(Gamma0(4).coset_reps())
            [
            [1 0]  [ 0 -1]  [1 0]  [ 0 -1]  [ 0 -1]  [1 0]
            [0 1], [ 1  0], [1 1], [ 1  2], [ 1  3], [2 1]
            ]
            sage: list(Gamma0(1).coset_reps())
            [
            [1 0]
            [0 1]
            ]
        """
        from all import SL2Z
        N = self.level()
        if N == 1: # P1List isn't very happy working modulo 1
            yield SL2Z([1,0,0,1])
        else:
            for z in sage.modular.modsym.p1list.P1List(N):
                yield SL2Z(lift_to_sl2z(z[0], z[1], N))

    def cusp_data(self, c):
        r"""
        Return a triple (g, w, t) where g is an element of self generating the
        stabiliser of the given cusp, w is the width of the cusp, and t is 1 if
        the cusp is regular and -1 if not.

        EXAMPLES::

            sage: Gamma1(4).cusp_data(Cusps(1/2))
            (
            [ 1 -1]
            [ 4 -3], 1, -1
            )
        """
        c = Cusp(c)

        # first find an element of SL2Z sending infinity to the given cusp
        w = lift_to_sl2z(c.denominator(), c.numerator(), 0)
        g = SL2Z([w[3], w[1], w[2],w[0]])

        for d in xrange(1,1+self.index()):
            if g * SL2Z([1,d,0,1]) * (~g) in self:
                return (g * SL2Z([1,d,0,1]) * (~g), d, 1)
            elif g * SL2Z([-1,-d,0,-1]) * (~g) in self:
                return (g * SL2Z([-1,-d,0,-1]) * (~g), d, -1)
        raise ArithmeticError("Can't get here!")

def lift_to_gamma1(g, m, n):
    r"""
    If ``g = [a,b,c,d]`` is a list of integers defining a `2 \times 2` matrix
    whose determinant is `1 \pmod m`, return a list of integers giving the
    entries of a matrix which is congruent to `g \pmod m` and to
    `\begin{pmatrix} 1 & * \\ 0 & 1 \end{pmatrix} \pmod n`. Here `m` and `n`
    must be coprime.

    Here `m` and `n` should be coprime positive integers. Either of `m` and `n`
    can be `1`. If `n = 1`, this still makes perfect sense; this is what is
    called by the function :func:`~lift_matrix_to_sl2z`. If `m = 1` this is a
    rather silly question, so we adopt the convention of always returning the
    identity matrix.

    The result is always a list of Sage integers (unlike ``lift_to_sl2z``,
    which tends to return Python ints).

    EXAMPLE::

        sage: from sage.modular.local_comp.liftings import lift_to_gamma1
        sage: A = matrix(ZZ, 2, lift_to_gamma1([10, 11, 3, 11], 19, 5)); A
        [371  68]
        [ 60  11]
        sage: A.det() == 1
        True
        sage: A.change_ring(Zmod(19))
        [10 11]
        [ 3 11]
        sage: A.change_ring(Zmod(5))
        [1 3]
        [0 1]
        sage: m = list(SL2Z.random_element())
        sage: n = lift_to_gamma1(m, 11, 17)
        sage: assert matrix(Zmod(11), 2, n) == matrix(Zmod(11),2,m)
        sage: assert matrix(Zmod(17), 2, [n[0], 0, n[2], n[3]]) == 1
        sage: type(lift_to_gamma1([10,11,3,11],19,5)[0])
        <type 'sage.rings.integer.Integer'>

    Tests with `m = 1` and with `n = 1`::

        sage: lift_to_gamma1([1,1,0,1], 5, 1)
        [1, 1, 0, 1]
        sage: lift_to_gamma1([2,3,11,22], 1, 5)
        [1, 0, 0, 1]
    """
    if m == 1:
        return [ZZ(1),ZZ(0),ZZ(0),ZZ(1)]
    a,b,c,d = [ZZ(x) for x in g]
    if not (a*d - b*c) % m == 1:
        raise ValueError( "Determinant is {0} mod {1}, should be 1".format((a*d - b*c) % m, m) )
    c2 = crt(c, 0, m, n)
    d2 = crt(d, 1, m, n)
    a3,b3,c3,d3 = map(ZZ, lift_to_sl2z(c2,d2,m*n))
    r = (a3*b - b3*a) % m
    return [a3 + r*c3, b3 + r*d3, c3, d3]

    def gamma0_coset_reps(self):
        r"""
        Return a set of coset representatives for self \\ Gamma0(N), where N is
        the level of self.

        EXAMPLE::

            sage: GammaH(108, [1,-1]).gamma0_coset_reps()
            [[1 0] [0 1], [-43 -45] [108 113], [ 31  33] [108 115], [-49 -54]
            [108 119], [ 25  28] [108 121], [-19 -22] [108 125], [-17 -20] [108
            127], [ 47  57] [108 131], [ 13  16] [108 133], [ 41  52] [108
            137], [  7   9] [108 139], [-37 -49] [108 143], [-35 -47] [108
            145], [ 29  40] [108 149], [ -5  -7] [108 151], [ 23  33] [108
            155], [-11 -16] [108 157], [ 53  79] [108 161]]
        """
        from all import Gamma0
        N = self.level()
        G = Gamma0(N)
        s = []
        return [G(lift_to_sl2z(0, d.lift(), N)) for d in _GammaH_coset_helper(N, self._list_of_elements_in_H())]