def new_local_splitting_big(q, l, prec):
    #######################################################
    #Same as _local_splitting from BTQuotient sage package
    #Except without enumerating entire graph
    #######################################################
    #INPUT: q = Discriminant of our quaternion algebra
    #       prec = l-adic precision 
    #       prec should be 2*(distance from starting vertex)
    #OUTPUT: representations of i,j,k as l-adic 2x2 matrices
    ########################################################
    A = QuaternionAlgebra(-1,-q)
    ZZl = Zp(l,prec)
    a = ZZl(-1)
    b = ZZl(-q)
    M = MatrixSpace(ZZl, 2)
    if a.is_square():
        alpha = a.sqrt()
        II = M([alpha, 0, 2*alpha, -alpha])
        JJ = M([b, -b, b-1, -b])
    else:
        II = M([0,a,1,0])
        z = 0
        JJ = 0
        while(JJ == 0):
            c = a*z*z + b
            if c.is_square():
                x = c.sqrt()
                JJ = M([x,-a*z,z,-x])
            else:
                z += 1
    KK = II*JJ
    return II, JJ, KK


def new_local_splitting_map_big(q,prec):
    #######################################################
    #Same as _local_splitting_map from BTQuotient sage package
    #Except without enumerating entire graph
    #######################################################
    #INPUT: q = Discriminant of our quaternion algebra
    #       prec = l-adic precision 
    #        prec should be 2*(distance from starting vertex)
    #OUTPUT: a map that takes elements of the quaternion algebra
    #        and outputs them as l-adic 2x2 matrices
    ####################################################### 
    I,J,K = new_local_splitting_big(q,l,prec)
    def phi(n):
        R = I.parent()
        v = n.coefficient_tuple()
        return R(v[0] + I*v[1]+ J*v[2] + K*v[3])
    return phi

def new_Iotainv(Obasis, q, l, prec):
    #######################################################
    #Same as Iotainv from BTQuotient sage package,
    #Except without enumerating entire graph,
    #######################################################
    #INPUT: basis of the order we fixed
    #       q = Discriminant of our quaternion algebra
    #       prec = l-adic precision\n",
    #       prec should be 2*(distance from starting vertex)
    #OUTPUT: Matrix associated to the inverse of embedding iota: O --> M2(Q_l)
    ######################################################
    B = Obasis
    phi = new_local_splitting_map_big(q,prec)
    lN = l**prec
    Iotamod = Matrix(Zmod(l ** prec), 4, 4,
                          [phi(B[kk])[ii, jj] for ii in range(2)
                           for jj in range(2) for kk in range(4)])
    Iotainv_lift = Iotamod.inverse().lift()
    Mat_44 = MatrixSpace(ZZ, 4, 4)
    Iotainv = Mat_44([Iotainv_lift[ii, jj] % lN for ii in range(4) for jj in range(4)])
    return Iotainv

def new_target_QF(q,l,dir,dis,repeat):
    prec = 2*dis+1

    X1 = Matrix(ZZ,2,[1,0,0,0])
    X2 = Matrix(ZZ,2,[0,1,0,0])
    X3 = Matrix(ZZ,2,[0,0,1,0])
    X4 = Matrix(ZZ,2,[0,0,0,1])
    X = [X1,X2,X3,X4]

    v0 = Matrix(ZZ,2,[1,0,0,1])
    v0adj = v0.adjugate()
    #choosing a direction - v1
    #creating the matrix that defines the whole isogeny
    v1 = dir^repeat
    
    # These matrices live in M2(Q_l)
    vecM = [v1*X[ii]*v0adj for ii in range(4)] 
    
    A.<i,j,k> = QuaternionAlgebra(q)
    v = [1] + A.gens()
    
    # Choose a maximal order O (it will be associated to v0).
    O = A.maximal_order()
    Obasis = O.basis()
    # Or we can also fix the basis we want
    Obasis = [A(1),i,(1+k)/2,(i+j)/2]     # O=Z[lambda1, lambda2, lambda3, lambda4]
    
    # iota: O --> M2(Q_l) This depends on the basis of O.
    # iotainv: M2(Q_l) --> O
    emb = new_Iotainv(Obasis,q,l,prec)
    R = MatrixSpace(ZZ,4)
    
    # embed the basis vecM to a basis in B
    Maux = emb * R([[vecM[ii][jj, kk] for ii in range(4)]  #defining an element of the resulting endomophism ring ai + bj + ck + d
                    for jj in range(2) for kk in range(2)])
    Maug = R(Maux).augment(R(l^dis))
    M = Maug.transpose()
    E = M.echelon_form().submatrix(0,0,4,4)
    Et = E.transpose() # columns encode the Eichler order associated to (v0,v1), wrt the basis lambda1,...,lambda4
    
    # Change of basis to get a basis for Et wrt 1,i,j,k
    
    C = Matrix(QQ, 4, 4, [[x[0],x[1],x[2],x[3]] for x in Obasis]).transpose()
    O_eichler = C*Et
    aux  = O_eichler.transpose()
    O_eichler_basis = [aux[ii][0] + aux[ii][1]*A.0 + aux[ii][2]*A.1 + aux[ii][3]*A.2 for ii in range(4)]
    
    #QFmat is the quadratic form matrix for the norm form of the starting order
    B = [A(Obasis[tt]) for tt in range(4)]
    Mat_44 = MatrixSpace(ZZ, 4, 4)
    OQuadForm = QuadraticForm(Mat_44([(B[ii] * B[jj].conjugate()).reduced_trace() for ii in range(4) for jj in range(4)]))
    QFmat = OQuadForm.matrix()
    
    #Quadratic form matrix for the Eichler order Et
    newQFmat = E*QFmat*Et #quadratic form matrix of the new order
    
    # Write the norm form
    Pol.<a,b,c,d> = PolynomialRing(QQ,4)
    poli = 0
    for ii in range(4):
        for jj in range(4):
            poli = poli + 1/2*newQFmat[ii][jj]*Pol.gen(ii)*Pol.gen(jj)
            
    return O_eichler_basis, newQFmat, poli