def perpendicular_bisector(self): #UHP
        r"""
        Return the perpendicular bisector of the hyperbolic geodesic ``self``
        if that geodesic has finite length.

        EXAMPLES::

            sage: UHP = HyperbolicPlane().UHP()
            sage: g = UHP.random_geodesic()
            sage: h = g.perpendicular_bisector()
            sage: c = lambda x: x.coordinates()
            sage: bool(c(g.intersection(h)[0]) - c(g.midpoint()) < 10**-9)
            True

        Infinite geodesics cannot be bisected::

            sage: UHP.get_geodesic(0, 1).perpendicular_bisector()
            Traceback (most recent call last):
            ...
            ValueError: the length must be finite
        """
        if self.length() == infinity:
            raise ValueError("the length must be finite")
        start = self._start.coordinates()
        d = self._model._dist_points(start, self._end.coordinates()) / 2
        S = self.complete()._to_std_geod(start)
        T1 = matrix([[exp(d/2), 0], [0, exp(-d/2)]])
        s2 = sqrt(2) * 0.5
        T2 = matrix([[s2, -s2], [s2, s2]])
        isom_mtrx = S.inverse() * (T1 * T2) * S # We need to clean this matrix up.
        if (isom_mtrx - isom_mtrx.conjugate()).norm() < 5*EPSILON: # Imaginary part is small.
            isom_mtrx = (isom_mtrx + isom_mtrx.conjugate()) / 2 # Set it to its real part.
        H = self._model.get_isometry(isom_mtrx)
        return self._model.get_geodesic(H(self._start), H(self._end))

def walsh_matrix(m0):
    """
    This is the generator matrix of a Walsh code. The matrix of
    codewords correspond to a Hadamard matrix.

    EXAMPLES::

        sage: walsh_matrix(2)
        [0 0 1 1]
        [0 1 0 1]
        sage: walsh_matrix(3)
        [0 0 0 0 1 1 1 1]
        [0 0 1 1 0 0 1 1]
        [0 1 0 1 0 1 0 1]
        sage: C = LinearCode(walsh_matrix(4)); C
        [16, 4] linear code over GF(2)
        sage: C.spectrum()
        [1, 0, 0, 0, 0, 0, 0, 0, 15, 0, 0, 0, 0, 0, 0, 0, 0]

    This last code has minimum distance 8.

    REFERENCES:

    - http://en.wikipedia.org/wiki/Hadamard_matrix
    """
    m = int(m0)
    if m == 1:
        return matrix(GF(2), 1, 2, [ 0, 1])
    if m > 1:
        row2 = [x.list() for x in walsh_matrix(m-1).augment(walsh_matrix(m-1)).rows()]
        return matrix(GF(2), m, 2**m, [[0]*2**(m-1) + [1]*2**(m-1)] + row2)
    raise ValueError("%s must be an integer > 0."%m0)

def matrix_multiplicative_order(m):
    r"""
    Return the order of the 2x2 matrix ``m``.
    """
    if m.is_one():
        return Integer(1)
    elif m.det() != 1 and m.det() != -1:
        return Infinity

    # now we compute the potentially preserved quadratic form
    # i.e. looking for A such that m^t A m = A
    m00 = m[0,0]
    m01 = m[0,1]
    m10 = m[1,0]
    m11 = m[1,1]
    M = matrix(m.base_ring(),
        [[m00**2, m00*m10, m10**2],
         [m00*m01, m00*m11, m10*m11],
         [m01**2, m01*m11, m11**2]])

    # might there be several solutions ? (other than scaling)... should not
    try:
        v = (M-identity_matrix(3)).solve_right()
    except ValueError: # no solution
        return False

    raise NotImplementedError("your matrix is conjugate to an orthogonal matrix but the angle might not be rational.. to be terminated.")

    # then we conjugate and check if the angles are rational
    # we need to take a square root of a symmetric matrix... this is not implemented!
    A = matrix(m.base_ring(), [[v[0],v[1]],[v[1],v[2]]])

    def parity_check_matrix(self):
        r"""
        Returns a parity check matrix of ``self``.

        This matrix is computed directly from :func:`original_code`.

        EXAMPLES::

            sage: set_random_seed(42)
            sage: C = codes.RandomLinearCode(9, 5, GF(7))
            sage: Ce = codes.ExtendedCode(C)
            sage: Ce.parity_check_matrix()
            [1 1 1 1 1 1 1 1 1 1]
            [1 0 0 0 2 1 6 6 4 0]
            [0 1 0 0 6 1 6 1 0 0]
            [0 0 1 0 3 2 6 2 1 0]
            [0 0 0 1 4 5 4 3 5 0]
        """
        F = self.base_ring()
        zero = F.zero()
        one = F.one()
        H = self.original_code().parity_check_matrix()
        nr, nc = H.nrows(), H.ncols()
        Hlist = H.list()
        v = matrix(F, nr + 1, 1, [one] + [zero] * nr)
        return matrix(F, nr + 1, nc, [one] * nc + Hlist).augment(v)

    def read_matrix(self, filename):
        r"""
        Read a matrix in 4ti2 format from the file ``filename`` in
        directory ``directory()``.

        INPUT:

        - ``filename`` - The name of the file to read from.

        OUTPUT:
        The data from the file as a matrix over `\ZZ`.

        EXAMPLES::

            sage: from sage.interfaces.four_ti_2 import four_ti_2
            sage: four_ti_2.write_matrix([[1,2,3],[3,4,6]], "test_file")
            sage: four_ti_2.read_matrix("test_file")
            [1 2 3]
            [3 4 6]
        """
        from sage.matrix.constructor import matrix
        try:
            f = open(os.path.join(self.directory(), filename))
            lines = f.readlines()
            f.close()
        except IOError:
            return matrix(ZZ, 0, 0)

        nrows, ncols = map(ZZ, lines.pop(0).strip().split())
        return matrix(ZZ, nrows, ncols,
                      [map(ZZ, line.strip().split()) for line in lines
                       if line.strip() != ""])

    def __init__(self, R, elements):
        """
        Initialize ``self``.

        EXAMPLES::

            sage: R.<x,y,z> = QQ[]
            sage: K = KoszulComplex(R, [x,y])
            sage: TestSuite(K).run()
        """
        # Generate the differentials
        self._elements = elements
        n = len(elements)
        I = range(n)
        diff = {}
        zero = R.zero()
        for i in I:
            M = matrix(R, binomial(n,i), binomial(n,i+1), zero)
            j = 0
            for comb in itertools.combinations(I, i+1):
                for k,val in enumerate(comb):
                    r = rank(comb[:k] + comb[k+1:], n, False)
                    M[r,j] = (-1)**k * elements[val]
                j += 1
            M.set_immutable()
            diff[i+1] = M
        diff[0] = matrix(R, 0, 1, zero)
        diff[0].set_immutable()
        diff[n+1] = matrix(R, 1, 0, zero)
        diff[n+1].set_immutable()
        ChainComplex_class.__init__(self, ZZ, ZZ(-1), R, diff)

    def random_isometry(self, preserve_orientation=True, **kwargs):
        r"""
        Return a random isometry in the Upper Half Plane model.

        INPUT:

        - ``preserve_orientation`` -- if ``True`` return an
          orientation-preserving isometry

        OUTPUT:

        - a hyperbolic isometry

        EXAMPLES::

            sage: A = HyperbolicPlane().UHP().random_isometry()
            sage: B = HyperbolicPlane().UHP().random_isometry(preserve_orientation=False)
            sage: B.preserves_orientation()
            False
        """
        [a,b,c,d] = [RR.random_element() for k in range(4)]
        while abs(a*d - b*c) < EPSILON:
            [a,b,c,d] = [RR.random_element() for k in range(4)]
        M = matrix(RDF, 2,[a,b,c,d])
        M = M / (M.det()).abs().sqrt()
        if M.det() > 0:
            if not preserve_orientation:
                M = M * matrix(2,[0,1,1,0])
        elif preserve_orientation:
            M = M * matrix(2,[0,1,1,0])
        return self._Isometry(self, M, check=False)

    def __repr__(self):
        r"""
        Return string representation.

        OUTPUT:

        String.

        EXAMPLES::

            sage: from sage.geometry.polyhedron.double_description import \
            ....:     DoubleDescriptionPair, StandardAlgorithm
            sage: A = matrix(QQ, [(1,0,1), (0,1,1), (-1,-1,1)])
            sage: DD = StandardAlgorithm(A).run()
            sage: DD.__repr__()
            'Double description pair (A, R) defined by\n    [ 1  0  1]
             [ 2/3 -1/3 -1/3]\nA = [ 0  1  1],   R = [-1/3  2/3 -1/3]\n
             [-1 -1  1]        [ 1/3  1/3  1/3]'
        """
        from sage.typeset.ascii_art import ascii_art
        from sage.matrix.constructor import matrix
        s = ascii_art('Double description pair (A, R) defined by')
        A = ascii_art(matrix(self.A))
        A._baseline = (len(self.A) / 2)
        A = ascii_art('A = ') + A
        R = ascii_art(matrix(self.R).transpose())
        if len(self.R) > 0:
            R._baseline = (len(self.R[0]) / 2)
        else:
            R._baseline = 0
        R = ascii_art('R = ') + R
        return str(s * (A + ascii_art(',   ') + R))

    def parity_check_matrix(self):
        r"""
        Returns a parity check matrix of ``self``.

        This matrix is computed directly from :func:`original_code`.

        EXAMPLES::

            sage: C = LinearCode(matrix(GF(2),[[1,0,0,1,1],\
                                               [0,1,0,1,0],\
                                               [0,0,1,1,1]]))
            sage: C.parity_check_matrix()
            [1 0 1 0 1]
            [0 1 0 1 1]
            sage: Ce = codes.ExtendedCode(C)
            sage: Ce.parity_check_matrix()
            [1 1 1 1 1 1]
            [1 0 1 0 1 0]
            [0 1 0 1 1 0]
        """
        F = self.base_ring()
        zero = F.zero()
        one = F.one()
        H = self.original_code().parity_check_matrix()
        nr, nc = H.nrows(), H.ncols()
        Hlist = H.list()
        v = matrix(F, nr + 1, 1, [one] + [zero] * nr)
        M = matrix(F, nr + 1, nc, [one] * nc + Hlist).augment(v)
        M.set_immutable()
        return M

    def __call__(self, n, modulus=0):

        """
        Give the nth term of a binary recurrence sequence, possibly mod some modulus.

        INPUT:

        - ``n`` -- an integer (the index of the term in the binary recurrence sequence)

        - ``modulus`` -- a natural number (optional --  default value is 0)

        OUTPUT:

        - An integer (the nth term of the binary recurrence sequence modulo ``modulus``)

        EXAMPLES::

            sage: R = BinaryRecurrenceSequence(3,3,2,1)
            sage: R(2)
            9
            sage: R(101)
            16158686318788579168659644539538474790082623100896663971001
            sage: R(101,12)
            9
            sage: R(101)%12
            9

        """
        R = Integers(modulus)
        F = matrix(R, [[0,1],[self.c,self.b]])            # F*[u_{n}, u_{n+1}]^T = [u_{n+1}, u_{n+2}]^T (T indicates transpose).
        v = matrix(R, [[self.u0],[self.u1]])
        return list(F**n*v)[0][0]

    def reflection_involution(self):
        r"""
        Return the isometry of the involution fixing the geodesic ``self``.

        EXAMPLES::

            sage: UHP = HyperbolicPlane().UHP()
            sage: g1 = UHP.get_geodesic(0, 1)
            sage: g1.reflection_involution()
            Isometry in UHP
            [ 1  0]
            [ 2 -1]
            sage: UHP.get_geodesic(I, 2*I).reflection_involution()
            Isometry in UHP
            [ 1  0]
            [ 0 -1]
        """
        x, y = [real(k.coordinates()) for k in self.ideal_endpoints()]
        if x == infinity:
            M = matrix([[1, -2*y], [0, -1]])
        elif y == infinity:
            M = matrix([[1, -2*x], [0, -1]])
        else:
            M = matrix([[(x+y)/(y-x), -2*x*y/(y-x)], [2/(y-x), -(x+y)/(y-x)]])
        return self._model.get_isometry(M)

    def symmetric_matrix(self):
        r"""
        The symmetric matrix `M` such that `(x y z) M (x y z)^t`
        is the defining equation of ``self``.

        EXAMPLES ::

            sage: R.<x, y, z> = QQ[]
            sage: C = Conic(x^2 + x*y/2 + y^2 + z^2)
            sage: C.symmetric_matrix()
            [  1 1/4   0]
            [1/4   1   0]
            [  0   0   1]

            sage: C = Conic(x^2 + 2*x*y + y^2 + 3*x*z + z^2)
            sage: v = vector([x, y, z])
            sage: v * C.symmetric_matrix() * v
            x^2 + 2*x*y + y^2 + 3*x*z + z^2
        """
        [a,b,c,d,e,f] = self.coefficients()
        if self.base_ring().characteristic() == 2:
            if b == 0 and c == 0 and e == 0:
                return matrix([[a,0,0],[0,d,0],[0,0,f]])
            raise ValueError, "The conic self (= %s) has no symmetric matrix " \
                              "because the base field has characteristic 2" % \
                              self
        from sage.matrix.constructor import matrix
        return matrix([[  a , b/2, c/2 ],
                       [ b/2,  d , e/2 ],
                       [ c/2, e/2,  f  ]])

    def bigraphical(self, G, A=None, K=QQ, names=None):
        r"""
        Return a bigraphical hyperplane arrangement.

        INPUT:

        - ``G`` -- graph

        - ``A`` -- list, matrix, dictionary (default: ``None``
          gives semiorder), or the string 'generic'

        - ``K`` -- field (default: `\QQ`)

        - ``names`` -- tuple of strings or ``None`` (default); the
          variable names for the ambient space

        OUTPUT:

        The hyperplane arrangement with hyperplanes `x_i - x_j =
        A[i,j]` and `x_j - x_i = A[j,i]` for each edge `v_i, v_j` of
        ``G``.  The indices `i,j` are the indices of elements of
        ``G.vertices()``.

        EXAMPLES::

            sage: G = graphs.CycleGraph(4)
            sage: G.edges()
            [(0, 1, None), (0, 3, None), (1, 2, None), (2, 3, None)]
            sage: G.edges(labels=False)
            [(0, 1), (0, 3), (1, 2), (2, 3)]
            sage: A = {0:{1:1, 3:2}, 1:{0:3, 2:0}, 2:{1:2, 3:1}, 3:{2:0, 0:2}}
            sage: HA = hyperplane_arrangements.bigraphical(G, A)
            sage: HA.n_regions()
            63
            sage: hyperplane_arrangements.bigraphical(G, 'generic').n_regions()
            65
            sage: hyperplane_arrangements.bigraphical(G).n_regions()
            59

        REFERENCES:

        ..  [BigraphicalArrangements] S. Hopkins, D. Perkinson.
            "Bigraphical Arrangements".
            :arxiv:`1212.4398`
        """
        n = G.num_verts()
        if A is None:  # default to G-semiorder arrangement
            A = matrix(K, n, lambda i, j: 1)
        elif A == 'generic':
            A = random_matrix(ZZ, n, x=10000)
            A = matrix(K, A)
        H = make_parent(K, n, names)
        x = H.gens()
        hyperplanes = []
        for e in G.edges():
            i = G.vertices().index(e[0])
            j = G.vertices().index(e[1])
            hyperplanes.append( x[i] - x[j] - A[i][j])
            hyperplanes.append(-x[i] + x[j] - A[j][i])
        return H(*hyperplanes)

    def orthonormal_1(dim_n=5):
        """
        A matrix of rational approximations to orthonormal vectors to
        ``(1,...,1)``.

        INPUT:

        - ``dim_n`` - the dimension of the vectors

        OUTPUT:

        A matrix over ``QQ`` whose rows are close to an orthonormal
        basis to the subspace normal to ``(1,...,1)``.

        EXAMPLES::

            sage: from sage.geometry.polyhedron.library import Polytopes
            sage: m = Polytopes.orthonormal_1(5)
            sage: m
            [ 70711/100000   -7071/10000             0             0             0]
            [    1633/4000     1633/4000 -81649/100000             0             0]
            [   7217/25000    7217/25000    7217/25000  -43301/50000             0]
            [ 22361/100000  22361/100000  22361/100000  22361/100000  -44721/50000]
        """
        pb = []
        for i in range(0,dim_n-1):
            pb.append([1.0/(i+1)]*(i+1) + [-1] + [0]*(dim_n-i-2))
        m = matrix(RDF,pb)
        new_m = []
        for i in range(0,dim_n-1):
            new_m.append([RDF(100000*q/norm(m[i])).ceil()/100000 for q in m[i]])
        return matrix(QQ,new_m)

def find_kadziela_matrices(M,T):
    '''
    The matrix M describes the relation between periods (A,B,D)^t
    and the periods (A0,B0)^t, where (A,B,D) are the periods of
    the Teitelbaum periods, and (A0,B0) are the Darmon ones.
           (A,B,D)^t = M * (A0,B0)^t
    The matrix T describes the action of Hecke on homology.
    That is, the first column of T describes the image of T
    on the first basis vector.

    The output are matrices X and Y such that
         X * matrix(2,2,[A,B,B,D]) = matrix(2,2,[A0,B0,C0,D0]) * Y

    '''
    a, b, c, d, e, f = M.list()
    x, y, z, t = T.list()
    #     1,  2,  3,  4,  5,  6,  7,  8
    r1 = [a,  c,  0,  0, -1,  0,  0,  0]
    r2 = [b,  d,  0,  0,  0,  0, -1,  0]
    r3 = [c,  e,  0,  0,  0, -1,  0,  0]
    r4 = [d,  f,  0,  0,  0,  0,  0, -1]
    r5 = [0,  0,  a,  c,  0,  0, -1,  0]
    r6 = [0,  0,y*b,y*d, -z,  0,x-t,  0]
    r7 = [0,  0,  c,  e,  0,  0,  0, -1]
    r8 = [0,  0,y*d,y*f,  0, -z,  0,x-t]
    AA = matrix(ZZ,8,8,[r1,r2,r3,r4,r5,r6,r7,r8])
    if AA.rank() == 8:
        raise ValueError('Not isogenous')
    r = AA.right_kernel().matrix().rows()[0].list()
    X = matrix(ZZ,2,2,r[:4])
    Y = matrix(ZZ,2,2,r[4:])
    return X, Y

    def __init__(self, n):
        r"""
        Hecke triangle group (2, n, infinity).
        Namely the von Dyck group corresponding to the triangle group
        with angles (pi/2, pi/n, 0).

        INPUT:

        - ``n``   - ``infinity`` or an integer greater or equal to ``3``.

        OUTPUT:

        The Hecke triangle group for the given parameter ``n``.

        EXAMPLES::

            sage: G = HeckeTriangleGroup(12)
            sage: G
            Hecke triangle group for n = 12
            sage: G.category()
            Category of groups
        """

        self._n = n
        self._T = matrix(AA, [[1,self.lam()],[0,1]])
        self._S = matrix(AA, [[0,-1],[1,0]])

        FinitelyGeneratedMatrixGroup_generic.__init__(self, ZZ(2), AA, [self._S, self._T])

 def Reverse(self):
     A1 = matrix(3, [1,1,1, 0,1,0, 0,0,1])
     A2 = matrix(3, [1,0,0, 1,1,1, 0,0,1])
     A3 = matrix(3, [1,0,0, 0,1,0, 1,1,1])
     R = matrix(3, [0,1,1, 1,0,1, 1,1,0])
     gens = {1:A1, 2:A2, 3:A3, 4:R}
     return MatrixCocycle(gens)

    def Cassaigne_accelerated(self, order=3):
        r"""
        EXAMPLES::

            sage: from slabbe.matrix_cocycle import cocycles
            sage: c = cocycles.Cassaigne_accelerated(order=3)
            sage: c
            Cocycle with 6 gens over Language of finite words over
            alphabet ['11', '121', '12^{2}1', '212', '21^{2}2', '22']

        """
        C1 = matrix(3, [1,1,0, 0,0,1, 0,1,0])
        C2 = matrix(3, [0,1,0, 1,0,0, 0,1,1])
        C = {1:C1, 2:C2}
        gens = {}
        for i in range(order):
            for (a,b) in [(1,2), (2,1)]:
                if i == 0:
                    code = '{}{}'.format(a,a)
                elif i == 1:
                    code = '{}{}{}'.format(a,b,a)
                else:
                    code = '{}{}^{{{}}}{}'.format(a,b,i,a)
                gens[code] = C[a]*C[b]**i*C[a]
        return MatrixCocycle(gens)

    def image_mod_n(self):
        r"""
        Return the image of this group in `SL(2, \ZZ / N\ZZ)`.

        EXAMPLE::

            sage: Gamma0(3).image_mod_n()
            Matrix group over Ring of integers modulo 3 with 2 generators (
            [2 0]  [1 1]
            [0 2], [0 1]
            )

        TEST::

            sage: for n in [2..20]:
            ...     for g in Gamma0(n).gamma_h_subgroups():
            ...       G = g.image_mod_n()
            ...       assert G.order() == Gamma(n).index() / g.index()
        """
        N = self.level()
        if N == 1:
            raise NotImplementedError("Matrix groups over ring of integers modulo 1 not implemented")
        gens = [matrix(Zmod(N), 2, 2, [x, 0, 0, Zmod(N)(1)/x]) for x in self._generators_for_H()]
        gens += [matrix(Zmod(N),2,[1,1,0,1])]
        return MatrixGroup(gens)

    def kernel_vector(self, way='LLL', verbose=False):
        r"""
        todo: clean this

        EXAMPLES::

            sage: from slabbe import ChristoffelGraph
            sage: C = ChristoffelGraph((2,5,7))
            sage: C.kernel_vector()
            [(-1, -1, 1), (3, -4, 0)]

        """
        from sage.arith.misc import gcd
        if way == 'vect_gcd':
            a,b,c = self._v
            gcd_ac = gcd(a,c)
            gcd_bc = gcd(b,c)
            U = ua,ub,uc = vector((c,0,-a)) / gcd(a,c)
            V = va,vb,vc = vector((0,c,-b)) / gcd(b,c)
            rows = U,V
        elif way == 'echelon':
            a,b,c = self._v
            m = matrix(ZZ, 4, [1,1,1,c,0,-a,0,c,-b,b,-a,0])
            me = m.echelon_form()
            if verbose:
                print(me)
            rows = me[1],me[2]
        elif way == 'LLL':
            dim = self.dimension()
            if dim == 3:
                a,b,c = self._v
                M = matrix(ZZ, 4, [1,1,1,c,0,-a,0,c,-b,b,-a,0])
            elif dim == 4:
                a,b,c,d = self._v
                M = matrix(ZZ, 7, (1,1,1,1,b,-a,0,0,c,0,-a,0,0,c,-b,0,d,0,0,-a,0,d,0,-b,0,0,d,-c))
            else:
                raise ValueError("dimension (=%s) must be 3 or 4" % dim)
            rows = M.LLL().rows()
            VS = rows[0].parent()
            zero = VS(0)
            un = VS((1,)*dim)
            assert zero in rows, "(0,0,0) not in LLL result"
            assert un in rows, "(1,1,1) not in LLL result"
            while zero in rows: rows.remove(zero)
            while un in rows: rows.remove(un)
        elif way == 'vect':
            a,b,c = self._v
            U = ua,ub,uc = vector((c,0,-a))
            V = va,vb,vc = vector((0,c,-b))
            rows = U,V
        else:
            raise ValueError("unknown way")
        R = matrix(rows)
        if sum(map(abs, R.minors(dim-1))) != sum(map(abs,self._v)):
            print(R)
            print(R.minors(dim-1))
            print(sum(map(abs, R.minors(dim))))
            print(sum(map(abs,self._v)))
            raise Exception("The result (=%s) is false " % rows)
        return rows