def _split_hyperbolic(L) :
            cur_cor = 2
            Lcor = L + cur_cor * identity_matrix(L.nrows())
            while not is_positive_definite(Lcor) :
                cur_cor += 2
                Lcor = L + cur_cor * identity_matrix(L.nrows())

            a = FreeModule(ZZ, L.nrows()).gen(0)
            if a * L * a >= 0 :
                Lcor_length_inc = max(3, a * L * a)
                cur_Lcor_length = Lcor_length_inc
                while True :
                    short_vectors = flatten(QuadraticForm(Lcor).short_vector_list_up_to_length( a * Lcor * a )[cur_Lcor_length - Lcor_length_inc: cur_Lcor_length], max_level = 1)
                    for a in short_vectors :
                        if a * L * a < 0 :
                            break
                    else :
                        continue
                    break
            n = -a * L * a // 2

            short_vectors = E8.short_vector_list_up_to_length(n + 1)[-1]

            for v in short_vectors :
                for w in short_vectors :
                    if v * E8_gram * w == 2 * n - 1 :
                        LE8_mat = L.block_sum(E8_gram)
                        v_form = vector( list(a) + list(v) ) * LE8_mat
                        w_form = vector( list(a) + list(w) ) * LE8_mat
                        Lred_basis = matrix(ZZ, [v_form, w_form]).right_kernel().basis_matrix().transpose()
                        Lred_basis = matrix(ZZ, Lred_basis)

                        return Lred_basis.transpose() * LE8_mat * Lred_basis

def _explode_embedding_list(v0,M,emblist,power = 1):
    p = v0.codomain().base_ring().prime()
    list_embeddings = []
    for tau0,gtau_orig in emblist:
        gtau = gtau_orig**power
        verbose('gtau = %s'%gtau)
        ## First method
        for u1 in is_represented_by_unit(M,ZZ(gtau[0,0]),p):
            u_M1 = matrix(QQ,2,2,[u1**-1,0,0,u1])
            gtau1 = u_M1 * gtau
            tau01 = tau0 / (u1**2)
            if gtau1[0,0] % M == 1:
                list_embeddings.append((tau01,gtau1,1))
            elif gtau1[0,0] % M == -1:
                list_embeddings.append((tau01,-gtau1,1))
        ## Second method
        if M > 1:
            a_inv = ZZ((1/Zmod(M)(gtau[0,0])).lift())
            for u2 in is_represented_by_unit(M,a_inv,p):
                u_M2 = matrix(QQ,2,2,[u2,0,0,u2**-1])
                gtau2 = u_M2 * gtau
                tau02 = u2**2 * tau0 #act_flt(u_M2,tau0)
                if gtau2[0,0] % M == 1:
                    list_embeddings.append((tau02,gtau2,1))
                elif gtau1[0,0] % M == -1:
                    list_embeddings.append((tau02,-gtau2,1))
    verbose('Found %s embeddings...'%len(list_embeddings))
    return list_embeddings

def matrix_of_independent_rows(field, rows, width):

    M = matrix(field, rows, sparse=True)
    N = matrix(field, 0, width, sparse=True)
    NE = copy(N)

    for i in range(M.nrows()):
        NE2 = NE.stack(M[i, :])
        NE2.echelonize()
        if not NE2[-1, :].is_zero():
            NE = NE2
            N = N.stack(M[i, :])

    return N

    def weil_representation(self) :
        r"""
        OUTPUT:
        
        - A pair of matrices corresponding to T and S.
        """
        disc_bilinear = lambda a, b: (self._dual_basis * vector(QQ, a.lift())) * self._L * (self._dual_basis * vector(QQ, b.lift()))
        disc_quadratic = lambda a: disc_bilinear(a, a) / ZZ(2)
        
        zeta_order = ZZ(lcm([8, 12, prod(self.invariants())] + map(lambda ed: 2 * ed, self.invariants())))
        K = CyclotomicField(zeta_order); zeta = K.gen()

        R = PolynomialRing(K, 'x'); x = R.gen()
#        sqrt2s = (x**2 - 2).factor()
#        if sqrt2s[0][0][0].complex_embedding().real() > 0 :        
#            sqrt2  = sqrt2s[0][0][0]
#        else : 
#            sqrt2  = sqrt2s[0][1]
        Ldet_rts = (x**2 - prod(self.invariants())).factor()
        if Ldet_rts[0][0][0].complex_embedding().real() > 0 :
            Ldet_rt  = Ldet_rts[0][0][0] 
        else :
            Ldet_rt  = Ldet_rts[0][0][0]
                
        Tmat  = diagonal_matrix( K, [zeta**(zeta_order*disc_quadratic(a)) for a in self] )
        Smat = zeta**(zeta_order / 8 * self._L.nrows()) / Ldet_rt  \
               * matrix( K,  [ [ zeta**ZZ(-zeta_order * disc_bilinear(gamma,delta))
                                 for delta in self ]
                               for gamma in self ])
        
        return (Tmat, Smat)

 def __init__(self, G, V, trivial_action = False):
     self._group = G
     self._coeffmodule = V
     self._trivial_action = trivial_action
     self._gen_pows = []
     self._gen_pows_neg = []
     if trivial_action:
         self._acting_matrix = lambda x, y: matrix(V.base_ring(),V.dimension(),V.dimension(),1)
         gens_local = [ (None, None) for g in G.gens() ]
     else:
         def acting_matrix(x,y):
             try:
                 return V.acting_matrix(x,y)
             except:
                 return V.acting_matrix(G.embed(x.quaternion_rep,V.base_ring().precision_cap()), y)
         self._acting_matrix = acting_matrix
         gens_local = [ (g, g**-1) for g in G.gens() ]
     onemat = G(1)
     try:
         dim = V.dimension()
     except AttributeError:
         dim = len(V.basis())
     one = Matrix(V.base_ring(),dim,dim,1)
     for g, ginv in gens_local:
         A = self._acting_matrix(g, dim)
         self._gen_pows.append([one, A])
         Ainv = self._acting_matrix(ginv, dim)
         self._gen_pows_neg.append([one, Ainv])
     Parent.__init__(self)
     return

def _hecke_operator_on_basis(B, V, n, k, eps):
    """
    Does the work for hecke_operator_on_basis once the input
    is normalized.

    EXAMPLES::

        sage: hecke_operator_on_basis(ModularForms(1,16).q_expansion_basis(30), 3, 16) # indirect doctest
        [   -3348        0]
        [       0 14348908]

    The following used to cause a segfault due to accidentally
    transposed second and third argument (#2107)::

        sage: B = victor_miller_basis(100,30)
        sage: t2 = hecke_operator_on_basis(B, 100, 2)
        Traceback (most recent call last):
        ...
        ValueError: The given basis vectors must be linearly independent.
    """
    prec = V.degree()
    TB = [hecke_operator_on_qexp(f, n, k, eps, prec, check=False, _return_list=True)
                for f in B]
    TB = [V.coordinate_vector(w) for w in TB]
    return matrix(V.base_ring(), len(B), len(B), TB, sparse=False)

 def space(self):
     r'''
     Calculates the homology space as a Z-module.
     '''
     verb = get_verbose()
     set_verbose(0)
     V = self.coefficient_module()
     R = V.base_ring()
     Vdim = V.dimension()
     G = self.group()
     gens = G.gens()
     ambient = R**(Vdim * len(gens))
     if self.trivial_action():
         cycles = ambient
     else:
         # Now find the subspace of cycles
         A = Matrix(R, Vdim, 0)
         for g in gens:
             for v in V.gens():
                 A = A.augment(matrix(R,Vdim,1,list(vector(g**-1 * v - v))))
         K = A.right_kernel_matrix()
         cycles = ambient.submodule([ambient(list(o)) for o in K.rows()])
     boundaries = []
     for r in G.get_relation_words():
         grad = self.twisted_fox_gradient(G(r).word_rep)
         for v in V.gens():
             boundaries.append(cycles(ambient(sum([list(a * vector(v)) for a in grad],[]))))
     boundaries = cycles.submodule(boundaries)
     ans = cycles.quotient(boundaries)
     set_verbose(verb)
     return ans

    def _matrix_(self, R):
        r"""
        Return Sage matrix from this matlab element.

        EXAMPLES::

            sage: A = matlab('[1,2;3,4]')       # optional - matlab
            sage: matrix(ZZ, A)                 # optional - matlab
            [1 2]
            [3 4]
            sage: A = matlab('[1,2;3,4.5]')     # optional - matlab
            sage: matrix(RR, A)                 # optional - matlab
            [1.00000000000000 2.00000000000000]
            [3.00000000000000 4.50000000000000]

            sage: a = matlab('eye(50)')         # optional - matlab
            sage: matrix(RR, a)                 # optional - matlab
            50 x 50 dense matrix over Real Field with 53 bits of precision

        """
        from sage.matrix.all import matrix
        matlab = self.parent()
        entries = matlab.strip_answer(matlab.eval("mat2str({0})".format(self.name())))
        entries = entries.strip()[1:-1].replace(';', ' ')
        entries = map(R, entries.split(' '))
        nrows, ncols = map(int, str(self.size()).strip().split())
        m = matrix(R, nrows, ncols, entries)
        return m

 def element_of_norm(self,N,use_magma = False,return_all = False, radius = None, max_elements = None, local_condition = None): # in nonsplitcartan
     try:
         if return_all:
             return [self._element_of_norm[N]]
         else:
             return self._element_of_norm[N]
     except (AttributeError,KeyError):
         pass
     if not hasattr(self,'_element_of_norm'):
         self._element_of_norm  = dict([])
     eps = self.eps
     q = self.q
     M = self.level
     llinv = (self.GFq(N)**-1).lift()
     if M != 1:
         while llinv % M != 1:
             llinv += q
     found = False
     for a,b in product(range(q*M),repeat = 2):
         if (a**2*llinv - b**2*eps*llinv - 1) % (q*M) == 0 and (b*eps) % M == 0:
             verbose('Found a=%s, b=%s'%(a,b))
             found = True
             break
     assert found
     m0 = matrix(ZZ,2,2,[a,b*llinv,b*eps,a*llinv])
     a,b,c,d = lift(m0,q*M).list()
     candidate = self.B([a,N*b,c,N*d])
     assert self._is_in_order(candidate)
     assert candidate.determinant() == N
     set_immutable(candidate)
     self._element_of_norm[N] = candidate
     if return_all:
         return [candidate]
     else:
         return candidate

    def hecke_matrix(self, l, use_magma = True, g0 = None, with_torsion = False): # l can be oo
        verb = get_verbose()
        set_verbose(0)
        if with_torsion:
            dim = len(self.gens())
            gens = self.gens()
        else:
            dim = self.rank()
            gens = self.free_gens()
        R = self.coefficient_module().base_ring()
        M = matrix(R,dim,dim,0)
        coeff_dim = self.coefficient_module().dimension()

        for j,cycle in enumerate(gens):
            # Construct column j of the matrix
            new_col = vector(self.apply_hecke_operator(cycle, l, use_magma = use_magma, g0 = g0))
            if with_torsion:
                M.set_column(j,list(new_col))
            else:
                try:
                    invs = self.space().invariants()
                    M.set_column(j,[o for o,a in zip(new_col,invs) if a == 0])
                except AttributeError:
                    M.set_column(j,list(new_col))
        set_verbose(verb)
        return M

 def transition_matrix(self, basis, n):
     """
     Returns the transitions matrix between self and basis for the
     homogenous component of degree n.
     
     EXAMPLES::
     
         sage: HLP = HallLittlewoodP(QQ)
         sage: s   = SFASchur(HLP.base_ring())
         sage: HLP.transition_matrix(s, 4)
         [             1             -t              0            t^2           -t^3]
         [             0              1             -t             -t      t^3 + t^2]
         [             0              0              1             -t            t^3]
         [             0              0              0              1 -t^3 - t^2 - t]
         [             0              0              0              0              1]
         sage: HLQ = HallLittlewoodQ(QQ)
         sage: HLQ.transition_matrix(s,3)
         [                        -t + 1                        t^2 - t                     -t^3 + t^2]
         [                             0                  t^2 - 2*t + 1           -t^4 + t^3 + t^2 - t]
         [                             0                              0 -t^6 + t^5 + t^4 - t^2 - t + 1]
         sage: HLQp = HallLittlewoodQp(QQ)
         sage: HLQp.transition_matrix(s,3)
         [      1       0       0]
         [      t       1       0]
         [    t^3 t^2 + t       1]
     """
     P = sage.combinat.partition.Partitions_n(n)
     Plist = P.list()
     m = []
     for row_part in Plist:
         z = basis(self(row_part))
         m.append( map( lambda col_part: z.coefficient(col_part), Plist ) )
     return matrix(m)

    def Cube_deformation(self,k, names=None):
        r""" 
        Construct, for each `k\in\ZZ_{\geq 0}`, a toric variety with
        `\ZZ_k`-torsion in the Chow group.

        The fans of this sequence of toric varieties all equal the
        face fan of a unit cube topologically, but the
        ``(1,1,1)``-vertex is moved to ``(1,1,2k+1)``. This example
        was studied in [FS]_.

        INPUT: 

        - ``k`` -- integer. The case ``k=0`` is the same as
          :meth:`Cube_face_fan`.
        
        - ``names`` -- string. Names for the homogeneous
          coordinates. See
          :func:`~sage.schemes.toric.variety.normalize_names`
          for acceptable formats.

        OUTPUT:

        A :class:`toric variety
        <sage.schemes.toric.variety.ToricVariety_field>`
        `X_k`. Its Chow group is `A_1(X_k)=\ZZ_k`.

        EXAMPLES::

            sage: X_2 = toric_varieties.Cube_deformation(2)
            sage: X_2
            3-d toric variety covered by 6 affine patches
            sage: X_2.fan().ray_matrix()
            [ 1  1 -1 -1 -1 -1  1  1]
            [ 1 -1  1 -1 -1  1 -1  1]
            [ 5  1  1  1 -1 -1 -1 -1]
            sage: X_2.gens()
            (z0, z1, z2, z3, z4, z5, z6, z7)

        REFERENCES:

        ..  [FS]
            William Fulton, Bernd Sturmfels, "Intersection Theory on
            Toric Varieties", http://arxiv.org/abs/alg-geom/9403002
        """
        # We are going to eventually switch off consistency checks, so we need
        # to be sure that the input is acceptable.
        try:
            k = ZZ(k)   # make sure that we got a "mathematical" integer
        except TypeError:
            raise TypeError("cube deformations X_k are defined only for "
                            "non-negative integer k!\nGot: %s" % k)
        if k < 0:
            raise ValueError("cube deformations X_k are defined only for "
                             "non-negative k!\nGot: %s" % k)
        rays = lambda kappa: matrix([[ 1, 1, 2*kappa+1],[ 1,-1, 1],[-1, 1, 1],[-1,-1, 1],
                                       [-1,-1,-1],[-1, 1,-1],[ 1,-1,-1],[ 1, 1,-1]])
        cones = [[0,1,2,3],[4,5,6,7],[0,1,7,6],[4,5,3,2],[0,2,5,7],[4,6,1,3]]
        fan = Fan(cones, rays(k))
        return ToricVariety(fan, coordinate_names=names)

 def set_wp(self, wp):
     epsinv = matrix(QQ,2,2,[0,-1,self.p,0])**-1
     set_immutable(wp)
     assert is_in_Gamma0loc(self.embed(wp,20) * epsinv, det_condition = False)
     assert all((self.is_in_Gpn_order(wp**-1 * g * wp) for g in self.Gpn_Obasis()))
     assert self.is_in_Gpn_order(wp)
     self._wp = wp
     return self._wp

def normalize_square_matrices(matrices):
    """
    Find a common space for all matrices.

    OUTPUT:

    A list of matrices, all elements of the same matrix space.

    EXAMPLES::

        sage: from sage.groups.matrix_gps.finitely_generated import normalize_square_matrices
        sage: m1 = [[1,2],[3,4]]
        sage: m2 = [2, 3, 4, 5]
        sage: m3 = matrix(QQ, [[1/2,1/3],[1/4,1/5]])
        sage: m4 = MatrixGroup(m3).gen(0)
        sage: normalize_square_matrices([m1, m2, m3, m4])
        [
        [1 2]  [2 3]  [1/2 1/3]  [1/2 1/3]
        [3 4], [4 5], [1/4 1/5], [1/4 1/5]
        ]
    """
    deg = []
    gens = []
    for m in matrices:
        if is_MatrixGroupElement(m):
            deg.append(m.parent().degree())
            gens.append(m.matrix())
            continue
        if is_Matrix(m):
            if not m.is_square():
                raise TypeError('matrix must be square')
            deg.append(m.ncols())
            gens.append(m)
            continue
        try:
            m = list(m)
        except TypeError:
            gens.append(m)
            continue
        if isinstance(m[0], (list, tuple)):
            m = [list(_) for _ in m]
            degree = ZZ(len(m))
        else:
            degree, rem = ZZ(len(m)).sqrtrem()
            if rem!=0:
                raise ValueError('list of plain numbers must have square integer length')
        deg.append(degree)
        gens.append(matrix(degree, degree, m))
    deg = set(deg)
    if len(set(deg)) != 1:
        raise ValueError('not all matrices have the same size')
    gens = Sequence(gens, immutable=True)
    MS = gens.universe()
    if not is_MatrixSpace(MS):
        raise TypeError('all generators must be matrices')
    if MS.nrows() != MS.ncols():
        raise ValueError('matrices must be square')
    return gens

    def _compute_padic_splitting(self,prec):
        verbose('Entering compute_padic_splitting')
        prime = self.p
        if self.seed is not None:
            self.magma.eval('SetSeed(%s)'%self.seed)
        R = Qp(prime,prec+10) #Zmod(prime**prec) #
        B_magma = self.Gn._B_magma
        a,b = self.Gn.B.invariants()
        if self._matrix_group:
            self._II = matrix(R,2,2,[1,0,0,-1])
            self._JJ = matrix(R,2,2,[0,1,1,0])
            goodroot = self.F.gen().minpoly().change_ring(R).roots()[0][0]
            self._F_to_local = self.F.hom([goodroot])
        else:
            verbose('Calling magma pMatrixRing')
            if self.F == QQ:
                _,f = self.magma.pMatrixRing(self.Gn._O_magma,prime*self.Gn._O_magma.BaseRing(),Precision = 20,nvals = 2)
                self._F_to_local = QQ.hom([R(1)])
            else:
                _,f = self.magma.pMatrixRing(self.Gn._O_magma,sage_F_ideal_to_magma(self.Gn._F_magma,self.ideal_p),Precision = 20,nvals = 2)
                try:
                    self._goodroot = R(f.Image(B_magma(B_magma.BaseRing().gen(1))).Vector()[1]._sage_())
                except SyntaxError:
                    raise SyntaxError("Magma has trouble finding local splitting")
                self._F_to_local = None
                for o,_ in self.F.gen().minpoly().change_ring(R).roots():
                    if (o - self._goodroot).valuation() > 5:
                        self._F_to_local = self.F.hom([o])
                        break
                assert self._F_to_local is not None
            verbose('Initializing II,JJ,KK')
            v = f.Image(B_magma.gen(1)).Vector()
            self._II = matrix(R,2,2,[v[i+1]._sage_() for i in xrange(4)])
            v = f.Image(B_magma.gen(2)).Vector()
            self._JJ = matrix(R,2,2,[v[i+1]._sage_() for i in xrange(4)])
            v = f.Image(B_magma.gen(3)).Vector()
            self._KK = matrix(R,2,2,[v[i+1]._sage_() for i in xrange(4)])
            self._II , self._JJ = lift_padic_splitting(self._F_to_local(a),self._F_to_local(b),self._II,self._JJ,prime,prec)
        self.Gn._F_to_local = self._F_to_local
        if not self.use_shapiro():
            self.Gpn._F_to_local = self._F_to_local

        self._KK = self._II * self._JJ
        self._prec = prec
        return self._II, self._JJ, self._KK

    def positive_quadratic_form(self) :
        r"""
        Find the Gram matrix of a positive definite quadratic form
        which is stabily equivalent to `L`.
        """
        E8_gram = matrix(ZZ, 8,
                [2, -1, 0, 0,  0, 0, 0, 0,
                -1, 2, -1, 0,  0, 0, 0, 0,
                0, -1, 2, -1,  0, 0, 0, -1,
                0, 0, -1, 2,  -1, 0, 0, 0,
                0, 0, 0, -1,  2, -1, 0, 0,
                0, 0, 0, 0,  -1, 2, -1, 0,
                0, 0, 0, 0,  0, -1, 2, 0,
                0, 0, -1, 0,  0, 0, 0, 2])
        E8 = QuadraticForm(E8_gram)
        L = self._L
        
        ## This is a workaround since GP crashes
        is_positive_definite = lambda L: all( L[:n,:n].det() > 0 for n in range(1, L.nrows() + 1) )
        def _split_hyperbolic(L) :
            cur_cor = 2
            Lcor = L + cur_cor * identity_matrix(L.nrows())
            while not is_positive_definite(Lcor) :
                cur_cor += 2
                Lcor = L + cur_cor * identity_matrix(L.nrows())

            a = FreeModule(ZZ, L.nrows()).gen(0)
            if a * L * a >= 0 :
                Lcor_length_inc = max(3, a * L * a)
                cur_Lcor_length = Lcor_length_inc
                while True :
                    short_vectors = flatten(QuadraticForm(Lcor).short_vector_list_up_to_length( a * Lcor * a )[cur_Lcor_length - Lcor_length_inc: cur_Lcor_length], max_level = 1)
                    for a in short_vectors :
                        if a * L * a < 0 :
                            break
                    else :
                        continue
                    break
            n = -a * L * a // 2

            short_vectors = E8.short_vector_list_up_to_length(n + 1)[-1]

            for v in short_vectors :
                for w in short_vectors :
                    if v * E8_gram * w == 2 * n - 1 :
                        LE8_mat = L.block_sum(E8_gram)
                        v_form = vector( list(a) + list(v) ) * LE8_mat
                        w_form = vector( list(a) + list(w) ) * LE8_mat
                        Lred_basis = matrix(ZZ, [v_form, w_form]).right_kernel().basis_matrix().transpose()
                        Lred_basis = matrix(ZZ, Lred_basis)

                        return Lred_basis.transpose() * LE8_mat * Lred_basis

        while not is_positive_definite(L) :
            L = _split_hyperbolic(L)
        
        return L

    def padic_regulator(self,prec=20):
        r"""
        Computes the canonical `p`-adic regulator on the extended Mordell-Weil group as in [MTT]
        (with the correction of [Wer] and sign convention in [SW].)
        The `p`-adic Birch and Swinnerton-Dyer conjecture
        predicts that this value appears in the formula for the leading term of the
        `p`-adic L-function.

        INPUT:

        - ``prec`` - the `p`-adic precision, default is 20.

        REFERENCES:

        - [MTT] B. Mazur, J. Tate, and J. Teitelbaum,
          On `p`-adic analogues of the conjectures of Birch and
          Swinnerton-Dyer, Inventiones mathematicae 84, (1986), 1-48.

        - [Wer] Annette Werner, Local heights on abelian varieties and rigid analytic unifomization,
          Doc. Math. 3 (1998), 301-319.

        - [SW] William Stein and Christian Wuthrich, Computations About Tate-Shafarevich Groups
          using Iwasawa theory, preprint 2009.

        EXAMPLES::

            sage: eq = EllipticCurve('130a1').tate_curve(5)
            sage: eq.padic_regulator()
            2*5^-1 + 1 + 2*5 + 2*5^2 + 3*5^3 + 3*5^6 + 5^7 + 3*5^9 + 3*5^10 + 3*5^12 + 4*5^13 + 3*5^15 + 2*5^16 + 3*5^18 + 4*5^19 + O(5^20)

        """
        prec = prec + 4

        K = Qp(self._p, prec=prec)
        rank = self._E.rank()
        if rank == 0:
            return K(1)

        if not self.is_split():
            raise NotImplementedError("The p-adic regulator is not implemented for non-split multiplicative reduction.")


        basis = self._E.gens()
        M = matrix.matrix(K, rank, rank, 0)

        height =   self.padic_height(prec= prec)
        point_height = [height(P) for P in basis]
        for i in range(rank):
            for j in range(i+1, rank):
                M[i, j] = M[j, i] = (- point_height[i] - point_height[j] + height(basis[i] + basis[j]))/2
        for i in range(rank):
            M[i,i] = point_height[i]

        return M.determinant()

 def generate_wp_candidates(self, p, ideal_p, **kwargs):
     eps = self.eps
     q = self.q
     for a,b in product(range(q),repeat = 2):
         if (p**2*a**2 - b**2*eps - p) % (q) == 0:
             verbose('Found a=%s, b=%s'%(a,b))
             break
     c = (self.GFq(p)**-1 * b * eps).lift()
     d = a
     a,b,c,d = lift(matrix(ZZ,2,2,[p*a,b,c,d]),q).list()
     Fp = FiniteField(p)
     if c % p == 0:
         c += a*q
         d += b*q
     assert c % p != 0
     r = (Fp(-a)*Fp(c*q)**-1).lift()
     a += q*c*r
     b += q*d*r
     ans = matrix(ZZ,2,2,[a,b,p*c,p*d])
     ans.set_immutable()
     yield ans

def matrix_generating_function(m):
    r"""
    returns the generating function of each scalar as a matrix
    """
    num_var = 10
    dimx = m.nrows()
    dimy = m.ncols()
    d = dict()
    for x in range(dimx):
        for y in range(dimy):
            d[(x,y)]=m[x,y].get_generating_function()
    return matrix(PolynomialRing(ZZ,['x'+str(i) for i in range(num_var)]),d)

def lift(A, N):
    r"""
    Lift a matrix A from SL_m(Z/NZ) to SL_m(Z).

    Follows Shimura, Lemma 1.38, p21.

    """
    assert A.is_square()
    assert A.determinant() != 0
    m = A.nrows()
    if m == 1:
        return identity_matrix(1)

    D, U, V = A.smith_form()
    if U.determinant() == -1 :
        U = matrix(2,2,[-1,0,0,1])* U
    if V.determinant() == -1 :
        V = V *matrix(2,2,[-1,0,0,1])
    D = U*A*V
    assert U.determinant() == 1
    assert V.determinant() == 1
    a = [ D[i, i] for i in range(m) ]
    b = prod(a[1:])
    W = identity_matrix(m)
    W[0, 0] = b
    W[1, 0] = b-1
    W[0, 1] = 1
    X = identity_matrix(m)
    X[0, 1] = -a[1]
    Ap = D.parent()(D)
    Ap[0, 0] = 1
    Ap[1, 0] = 1-a[0]
    Ap[1, 1] *= a[0]
    assert (W*U*A*V*X).change_ring(Zmod(N)) == Ap.change_ring(Zmod(N))
    Cp = diagonal_matrix(a[1:])
    Cp[0, 0] *= a[0]
    C = lift(Cp, N)
    Cpp = block_diagonal_matrix(identity_matrix(1), C)
    Cpp[1, 0] = 1-a[0]
    return (~U * ~W * Cpp * ~X * ~V).change_ring(ZZ)