#Computes Eisenstein series and Poincare square series for the dual Weil representation

#most of this algorithm consists of implementing of the p-adic generating function of Cowan, Katz, White for calculating certain Igusa zeta functions

@cached_function
def local_normal_form_with_change_vars(Qmatrix,p):
    #basically copied from the built-in function local_normal_form(), but it also computes a change-of-variables to the local normal form. The input is a Gram matrix "Qmatrix" and an odd prime "p". The output is a diagonal quadratic form Q_Jordan and a change-of-basis matrix M.
    Q = QuadraticForm(Qmatrix)
    I = range(Q.dim())
    assert p != 2
    M = identity_matrix(QQ,Q.dim())
    Q_Jordan = DiagonalQuadraticForm(ZZ,[])
    while Q.dim() > 0:
        n = Q.dim()
        (min_i, min_j) = Q.find_entry_with_minimal_scale_at_prime(p)
        if min_i == min_j:
            min_val = valuation(2 * Q[min_i, min_j], p)
        else:
            min_val = valuation(Q[min_i, min_j], p)
        if (min_i == min_j):
            Q.swap_variables(0, min_i, in_place = True)
            M.swap_rows(I[0],I[min_i])
        else:
            Q.swap_variables(0, min_i, in_place = True)
            Q.swap_variables(1, min_j, in_place = True)
            M.swap_rows(I[0],I[min_i])
            M.swap_rows(I[1],I[min_j])
            Q.add_symmetric(1, 0, 1, in_place = True)
            M.add_multiple_of_row(I[0],I[1],1)
        min_scale = p ** min_val
        a = 2 * Q[0,0]
        for j in range(1, n):
            b = Q[0, j]
            g = GCD(a, b)
            Q.multiply_variable(ZZ(a/g), j, in_place = True)
            Q.add_symmetric(ZZ(-b/g), j, 0, in_place = True)
            M.rescale_row(I[j],a/g)
            M.add_multiple_of_row(I[j],I[0],-b/g)
        Q_Jordan = Q_Jordan + Q.extract_variables(range(1))
        I.remove(I[0])
        Q = Q.extract_variables(range(1, n))
    return [Q_Jordan,M]


def IsospectralNormalForm(Q,L,p): #Takes a quadratic form "Q" and linear form "L" and outputs an isospectral normal form over the odd prime "p". The quadratic part of the result is represented as a list "quads" = [a_1,...,a_n] (corresponding to the quadratic form a_1 x_1^2 + ... + a_n x_n^2); the linear part "linsgcd" is a number b (corresponding to the linear form y->b*y) and the constant term of the result is "const".
    [D,P] = local_normal_form_with_change_vars(Q.matrix(),p)
    L = L*P.transpose()
    lins = []
    quads = []
    consts = []
    for j in range(Q.dim()):
        a = D[j,j]
        b = L[0,j] / 2
        if b.valuation(p) < a.valuation(p):
            lins.append(b)
        else:
            quads.append(a)
            consts.append(-b^2 / (4*a))
    const = sum(consts)
    if len(lins) > 0:
        linsgcd = gcd(lins)
        return [quads,linsgcd,const]
    else:
        return [quads,0,const]


def IARD(a,r,d,p): #Computes the helper function I_a(r,d) with respect to the odd prime p, as in Proposition 20 of jacobi_eisenstein.pdf
    assert p != 2
    R.<t> = PolynomialRing(QQ)
    if r % 2 == 0:
        kron = kronecker_symbol((-1)^(r/2) * d,p)
        iard = 1 - p^(-r/2) * kron
        if a % p == 0:
            iard = iard * (1 + p^(-r/2) * kron * t) *(1 - 1/p) / (1 - t/p)
        else:
            iard = iard * ((1 - 1/p) / (1 - t/p) + p^(-r/2) * kron)
    else:
        iard = 1
        if a % p == 0:
            iard = iard * (1 - t*p^(-r)) * (1 - p^(-1)) / (1 - p^(-1) * t)
        else:
            kron = kronecker_symbol(a*d*(-1)^((1+r)/2),p)
            iard = iard * (1 + p^(-(1+r)/2)*kron*t)*(1 - p^(-1)) / (1 - p^(-1) * t) - p^(-r) - p^(-(1+r)/2) * kron
    return iard

def IgusaZeta(Q,L,c,p): #Computes the Igusa (local) zeta function of the polynomial Q(x) + L(x) + c over the *odd* prime p. The output is a power series in the variable "t" (interpreted as t = p^(-s)).
    [wrongquads,linsgcd,const] = IsospectralNormalForm(Q,L,p)
    c = const + c/2
    u = denominator(c/2)
    c = c*u
    linsgcd = linsgcd*u
    quads = []
    for a in wrongquads:
        quads.append(a*u)
    R.<t> = PolynomialRing(QQ)
    Z = 0*t
    dold = []
    rold = []
    for j in range(len(quads)):
        a = quads[j]
        pi = a.valuation(p)
        while pi > len(dold) - 1:
            dold.append(1)
            rold.append(0)
        dold[pi] = dold[pi]* ZZ(a / p^pi)
        rold[pi] = rold[pi] + 1
    d = []
    r = []
    qpower = []
    if linsgcd == 0:
        if c != 0:
            maxrange = max(len(dold),c.valuation(p)+1)
        else:
            maxrange = len(dold)
    else:
        if c != 0:
            maxrange = max(len(dold),c.valuation(p)+1,linsgcd.valuation(p))
        else:
            maxrange = max(len(dold),linsgcd.valuation(p))
    for k in range(maxrange+1):
        if k >= len(dold):
            dold.append(1)
            rold.append(0)
        l = k
        while l <= maxrange:
            if l > len(d) - 1:
                d.append(1)
                r.append(0)
                qpower.append(0)
            if (l % 2 == k % 2):
                d[l] = d[l] * dold[k]
                r[l] = r[l] + rold[k]
            if l != k:
                qpower[l] = qpower[l] + r[k]
            l = l + 1
    w = len(dold) - 1
    if linsgcd == 0 and c == 0:
            R = sum(rold)
            Z = (t^(w-1)*IARD(0,r[w-1],d[w-1],p) / p^(qpower[w-1]) + t^w *IARD(0,r[w],d[w],p) / p^(qpower[w]) )/ (1 - t^2 / p^R)
            for i in range(w-1):
                Z = Z + IARD(0,r[i],d[i],p)*t^i / p^(qpower[i])
            return Z
    elif linsgcd.valuation(p) > c.valuation(p):
        kappa = c.valuation(p)
        Z = t^kappa / p^(qpower[kappa + 1])
        for j in range(kappa + 1):
            Z = Z + t^j * IARD(ZZ(c/p^j),r[j],d[j],p) / p^(qpower[j])
        return Z
    elif linsgcd != 0:
        lamda = linsgcd.valuation(p)
        Z = (t^lamda / p^(qpower[lamda])) * ((1 - 1/p) / (1 - t/p))
        for i in range(lamda):
            Z = Z + t^i * IARD(0,r[i],d[i],p) / p^(qpower[i])
        return Z
    
    

@cached_function
def TwoadicClassify(Qmatrix): #Classifies unimodular quadratic forms over Z2. Every form is equivalent to a direct sum of at most two squares ax^2, and at most one "elliptic plane" 2(x^2 + xy + y^2), and any number of hyperbolic planes 2xy. "Qmatrix" is a Gram matrix of the input quadratic form. This doesn't verify whether Qmatrix is unimodular.
    J = twonf_with_change_vars_correct(Qmatrix)
    J = J[0].matrix()
    Squares = []
    Ell = 0
    Hyp = 0
    k = 0
    while k < J.nrows():
        if k < J.nrows() - 1:
            if J[k,k+1] != 0:
                if J[k,k]*J[k+1,k+1] % 8 == 0:
                    Hyp = Hyp + 1
                else:
                    Ell = Ell + 1
                k = k + 2
            else:
                Squares.append((J[k,k] / 2) % 8)
                Squares.append((J[k+1,k+1] / 2) % 8)
                k = k + 2
        else:
            Squares.append((J[k,k] / 2) % 8)
            k = k + 1
    while len(Squares) > 2:
        [a,b,c] = sorted(Squares[0:3])
        Squares = Squares[3:len(Squares)]
        if [a,b,c] in [[1,3,5],[1,1,7],[5,5,7]]:
            Squares.append(1)
            Hyp = Hyp + 1
        elif [a,b,c] in [[1,3,7],[3,3,5],[5,7,7]]:
            Squares.append(3)
            Hyp = Hyp + 1
        elif [a,b,c] in [[1,5,7],[3,5,5],[1,1,3]]:
            Squares.append(5)
            Hyp = Hyp + 1
        elif [a,b,c] in [[3,5,7],[1,7,7],[1,3,3]]:
            Squares.append(7)
            Hyp = Hyp + 1
        elif [a,b,c] in [[3,3,3],[3,7,7]]:
            Squares.append(1)
            Ell = Ell + 1
        elif [a,b,c] in [[1,1,1],[1,5,5]]:
            Squares.append(3)
            Ell = Ell + 1
        elif [a,b,c] in [[7,7,7],[3,3,7]]:
            Squares.append(5)
            Ell = Ell + 1
        elif [a,b,c] in [[5,5,5],[1,1,5]]:
            Squares.append(7)
            Ell = Ell + 1
    r = floor(Ell/2)
    Hyp = Hyp + 2*r
    Ell = Ell % 2
    return [Squares,Ell,Hyp]
    
def IaQQQrank(Q0): #Helper function
    return len(Q0[0]) + 2*Q0[1] + 2*Q0[2]
    
def Ig(a,u,j): #This is Ig(z^(a + 2^b Z_2 + 2^v Z_2) in the notation of [CKW] section 3.4
    R.<t> = PolynomialRing(QQ)
    if u == Infinity:
        if a.valuation(2) >= j:
            return t^j / (2-t)
        else:
            return t^(a.valuation(2))
    else:
        return Ig(a,Infinity,min(u,j))
    
def hatHQ(a,u,Q0): #the helper function H1
    R.<t> = PolynomialRing(QQ)
    r = IaQQQrank(Q0)
    if len(Q0[0]) == 0:
        return (1 - 2^(-r))*Ig(a,u,1)
    else:
        return Ig(a,u,0) - 2^(-r)*Ig(a,u,1)
    
def tildeHQ(a,u,Q0): #the helper function H2
    R.<t> = PolynomialRing(QQ)
    eps = (-1)^Q0[1]
    r = IaQQQrank(Q0)
    if len(Q0[0]) == 0:
        return (1 - eps*2^(-r/2))*(Ig(a,u,1) + eps*2^(-r/2)*Ig(a,u,2))
    elif len(Q0[0]) == 1:
        b = Q0[0][0]
        return (1 - eps*2^(-(r-1)/2))*Ig(a,u,0) - 2^(-r) * Ig(a,u,2) + eps*2^(-(r+1)/2) * (Ig(a,u,2) + Ig(a+b,u,2))
    else:
        b = Q0[0][0]
        c = Q0[0][1]
        if (b+c) % 4 == 0:
            return Ig(a,u,0) - eps*2^(-r/2) * Ig(a,u,1) + (eps*2^(-r/2) - 2^(-r)) * Ig(a,u,2)
        else:
            return (1 - eps*2^(-(r-2)/2))*Ig(a,u,0) + eps*2^(-r/2)*Ig(a,u,1) + eps*2^(-r/2)*Ig(a+b,u,2) - 2^(-r)*Ig(a,u,2)
    
def tildeHQdiff(a,u,Q0): #the helper function H3
    R.<t> = PolynomialRing(QQ)
    eps = (-1)^Q0[1]
    r = IaQQQrank(Q0)
    if len(Q0[0]) == 0:
        return t - t
    if len(Q0[0]) == 1:
        b = Q0[0][0]
        return 2^(-r)*(Ig(a+b,u,3) - Ig(a+b,u,2))
    if len(Q0[0]) == 2:
        b = Q0[0][0]
        c = Q0[0][1]
        if (b+c) % 4 == 0:
            sig = (-1)^((b+c)/4)
            return sig * 2^(-r) * (Ig(a,u,3) - Ig(a,u,2))
        else:
            return -2^(1-r) * Ig(a,u,1) + 2^(-r) * (Ig(a,u,2) + Ig(a+b+c,u,3))
    

def IaQQQ(a,u,U0,U1,U2): #Computes the helper functions I_a(U_0,U_1,U_2) as in Appendix B
    Q0 = TwoadicClassify(U0.matrix()) #basically Q0 = U0. This doesn't seem to be broken so I won't fix it.
    Q1 = TwoadicClassify(U1.matrix())
    Q2 = TwoadicClassify(U2.matrix())
    R.<t> = PolynomialRing(QQ)
    if len(Q1[0]) > 0 and len(Q2[0]) > 0:
        return hatHQ(a,u,Q0)
    elif len(Q2[0]) > 0:
        return tildeHQ(a,u,Q0)
    eps1 = (-1)^Q1[1]
    r1 = IaQQQrank(Q1)
    if len(Q1[0]) == 0:
        return tildeHQ(a,u,Q0) + eps1*2^(-r1/2) * tildeHQdiff(a,u,Q0)
    elif len(Q1[0]) == 1:
        c = Q1[0][0]
        return hatHQ(a,u,Q0) + eps1*2^(-(r1+1)/2) * (tildeHQdiff(a,u,Q0) + tildeHQdiff(a+2*c,u,Q0))
    elif len(Q1[0]) == 2:
        c = Q1[0][0]
        d = Q1[0][1]
        if (c + d) % 4 == 0:
            return hatHQ(a,u,Q0) + eps1 * 2^(-r1/2) * tildeHQdiff(a,u,Q0)
        else:
            return hatHQ(a,u,Q0) + eps1 * 2^(-r1/2) * tildeHQdiff(a+2*c,u,Q0)

#@cached_function
def twonf_with_change_vars_correct(Qmatrix): #basically copied from the built-in function local_normal_form(), specialized to p=2, but it also computes a change-of-variables to the local normal form. it skips the "Cassels' proof" step. this is accounted for later.
    Q = QuadraticForm(Qmatrix)
    #optimistic name for luck. this function has been hard to debug
    I = range(Q.dim())
    P = identity_matrix(QQ,Q.dim())
    Q_Jordan = DiagonalQuadraticForm(ZZ,[])
    while Q.dim() > 0:
        n = Q.dim()
        (min_i,min_j) = Q.find_entry_with_minimal_scale_at_prime(2)
        if min_i == min_j:
            min_val = valuation(2 * Q[min_i,min_j],2)
        else:
            min_val = valuation(Q[min_i, min_j],2)
        if (min_i == min_j):
            block_size = 1
            Q.swap_variables(0,min_i,in_place = True)
            P.swap_rows(I[0],I[min_i])
        else:
            Q.swap_variables(0, min_i, in_place = True)
            Q.swap_variables(1, min_j, in_place = True)
            P.swap_rows(I[0],I[min_i])
            P.swap_rows(I[1],I[min_j])
            block_size = 2
        min_scale = 2^(min_val)
        if (block_size == 1):
            a = 2 * Q[0,0]
            for j in range(block_size, n):
                b = Q[0, j]
                g = GCD(a, b)
                Q.multiply_variable(ZZ(a/g), j, in_place = True)
                Q.add_symmetric(ZZ(-b/g), j, 0, in_place = True)
                P.rescale_row(I[j],a/g)
                P.add_multiple_of_row(I[j],I[0],-b/g)
            NewQ = QuadraticForm(matrix([a]))
        else:
            a1 = 2*Q[0,0]
            a2 = Q[0,1]
            b2 = 2*Q[1,1]
            big_det = (a1*b2 - a2^2)
            small_det = big_det / (min_scale * min_scale)
            for j in range(block_size,n):
                a = Q[0,j]
                b = Q[1,j]
                Q.multiply_variable(big_det,j,in_place = True)
                Q.add_symmetric(ZZ(-(a*b2 - b*a2)),j,0,in_place = True)
                Q.add_symmetric(ZZ(-(-a*a2 + b*a1)),j,1,in_place = True)
                Q.divide_variable(ZZ(min_scale^2),j,in_place = True)
                P.rescale_row(I[j],big_det)
                P.add_multiple_of_row(I[j],I[0],-a*b2+b*a2)
                P.add_multiple_of_row(I[j],I[1],a*a2-b*a1)
                P.rescale_row(I[j],min_scale^(-2))
        Q_Jordan = Q_Jordan + Q.extract_variables(range(block_size))
        for j in range(block_size):
            I.remove(I[0])
        Q = Q.extract_variables(range(block_size, n))
    return [Q_Jordan,P]

@cached_function
def TwoadicIsospectralNormalForm(Q,L): #reduces the polynomial Q(x) + L(x) to an isospectral normal form. In the output NewQ is a quadratic form in Jordan normal form and "linsgcd" is a constant (which should be understood as a linear form in a variable which does not appear in NewQ), and "const" is an additional constant term.
    [D,P] = twonf_with_change_vars_correct(Q.matrix())
    L = L*P.transpose()
    lins = []
    NewQ = DiagonalQuadraticForm(ZZ,[])
    consts = []
    j = 0
    while j < Q.dim():
        if j < Q.dim() - 1 and D[j,j+1] != 0:
            l = L[0,j]
            m = L[0,j+1]
            A = 2*D[j,j]
            B = 2*D[j,j+1]
            C = 2*D[j+1,j+1]
            assert min(A.valuation(2),C.valuation(2)) >= B.valuation(2)
            if min(l.valuation(2),m.valuation(2)) >= B.valuation(2):
                NewQ = NewQ + D.extract_variables([j,j+1])
                U = matrix([[2*A,B],[B,2*C]])
                v = U.inverse()*vector([-l,-m])
                consts.append(A*v[0]^2 + B*v[0]*v[1] + C*v[1]^2 + l*v[0] + m*v[1]) #complete the square
            else:
                lins.append(l)
                lins.append(m)
            j = j + 2
        else:
            a = 2*D[j,j]
            b = L[0,j]
            if b.valuation(2) < a.valuation(2):
                lins.append(b)
            elif b.valuation(2) == a.valuation(2):
                lins.append(2*b)
            else:
                NewQ = NewQ + D.extract_variables([j])
                consts.append(-b^2 / (4*a))
            j = j + 1
    const = sum(consts)
    if len(lins) > 0:
        linsgcd = gcd(lins)
        return [NewQ,linsgcd,const]
    else:
        return [NewQ,0,const]

@cached_function
def twoadic_jordan_blocks_correct(Qmatrix): #decompose a Gram matrix which is in Jordan normal form into a list of unimodular Jordan blocks and valuations. The suffix "_correct" suggests that I found this difficult but I no longer remember why.
    Q = QuadraticForm(Qmatrix)
    Jordans = []
    Jordanlist = []
    valuations = []
    i = 0
    n = Q.dim()
    while (i < n):
        if (i == n-1) or (Q[i,i+1] == 0):
            block_size = 1
        else:
            block_size = 2
        if block_size == 1:
            scale = valuation(Q[i,i], 2) + 1
            if len(valuations) > 0 and scale == valuations[len(valuations) -1]:
                e = valuations.index(scale)
                Jordans[e] = Jordans[e] + QuadraticForm(ZZ,Q.extract_variables([i]).matrix() / 2^(scale - 1))
            else:
                valuations.append(scale)
                Jordans.append(QuadraticForm(ZZ,Q.extract_variables([i]).matrix() / 2^(scale - 1)))
            i = i + 1
        else:
            scale = valuation(Q[i,i+1], 2)
            if len(valuations) > 0 and scale == valuations[len(valuations) -1]:
                e = valuations.index(scale)
                Jordans[e] = Jordans[e] + QuadraticForm(ZZ,Q.extract_variables([i,i+1]).matrix() / 2^(scale))
            else:
                valuations.append(scale)
                Jordans.append(QuadraticForm(ZZ,Q.extract_variables([i,i+1]).matrix() / 2^(scale)))
            i = i + 2
    for i in range(len(valuations)):
        Jordanlist.append([valuations[i],Jordans[i]])
    return Jordanlist


def TwoadicIgusaZeta(Q,L,c): #as in the above function IgusaZeta(Q,L,c,p), but now p=2.
    bigr = Q.dim()
    [NewQ,linsgcd,const] = TwoadicIsospectralNormalForm(Q,L)
    const = const + c
    R.<t> = PolynomialRing(QQ)
    Z = 0*t
    jordanblocks = twoadic_jordan_blocks_correct(NewQ.matrix())
    Qs = []
    vb = linsgcd.valuation(2)
    vc = const.valuation(2)
    if len(jordanblocks) > 0:
        maxrange = jordanblocks[len(jordanblocks)-1][0] + 1
    else:
        maxrange = 2
    if linsgcd != 0:
        maxrange = max(maxrange,vb+1)
    if const != 0:
        maxrange = max(maxrange,vc+3)
    maxrange = max(maxrange,2)
    jordans = [DiagonalQuadraticForm(ZZ,[])] * maxrange
    for k in range(len(jordanblocks)):
        jordans[jordanblocks[k][0]] = jordans[jordanblocks[k][0]] + jordanblocks[k][1]
    rs = []
    Zeroform = DiagonalQuadraticForm(ZZ,[])
    Sqform = DiagonalQuadraticForm(ZZ,[1])
    qpowers = []
    for k in range(maxrange):
        Qs.append(DiagonalQuadraticForm(ZZ,[]))
        rs.append(0)
        l = k
        while l >= 0:
            Qs[k] = Qs[k] + jordans[l]
            rs[k] = rs[k] + jordans[l].dim()
            l = l - 2
        if k == 0:
            qpowers.append(1)
        else:
            qpowers.append(qpowers[k-1] * 2^(rs[k-1]))
    if vb == Infinity and vc == Infinity:
        r = sum(rs)
        w = max(len(Qs) - 1,1)
        Z = (t^(w-1) * IaQQQ(0,Infinity,Qs[w-1],Qs[w],Zeroform) * qpowers[w-1]^(-1) + t^w * IaQQQ(0,Infinity,Qs[w],Qs[w-1],Zeroform) * qpowers[w]^(-1)) / (1 - t^2 * 2^(-bigr))
        for i in range(0,w-1):
            Z = Z + t^i * IaQQQ(0,Infinity,Qs[i],Qs[i+1],jordans[i+2]) / qpowers[i]
        return Z
    elif vb < Infinity and vb <= vc:
        Z = t^vb / (qpowers[vb] * (2 - t))
        for i in range(vb - 2):
            Z = Z + t^i * IaQQQ(0,Infinity,Qs[i],Qs[i+1],jordans[i+2]) / qpowers[i]
        if vb >= 2:
            Z = Z + t^(vb-2) * IaQQQ(0,2,Qs[vb-2],Qs[vb-1],Sqform) / qpowers[vb-2]
        if vb >= 1:
            Z = Z + t^(vb-1) * IaQQQ(0,1,Qs[vb-1],Sqform,Sqform) / qpowers[vb-1]
        return Z
    elif vb < Infinity and vc < vb and vb <= vc + 2:
        Z = t^vc / qpowers[vc+1]
        for i in range(vb - 2):
            Z = Z + t^i * IaQQQ(const/2^i,Infinity,Qs[i],Qs[i+1],jordans[i+2]) / qpowers[i]
        if vb >= 2:
            Z = Z + t^(vb-2) * IaQQQ(const/2^(vb-2),2,Qs[vb-2],Qs[vb-1],Sqform) / qpowers[vb-2]
        if vb >= 1 and vb <= vc + 1:
            Z = Z + t^(vb-1) * IaQQQ(const/2^(vb-1),1,Qs[vb-1],Sqform,Sqform) / qpowers[vb-1]
        return Z
    else:
        Z = t^vc / qpowers[vc+1]
        for i in range(vc+1):
            Z = Z + t^i * IaQQQ(const/2^i,Infinity,Qs[i],Qs[i+1],jordans[i+2]) / qpowers[i]
        return Z


def Lvalue(g,n,S,p,k): #the Euler factors in the Eisenstein series
    R.<t> = PolynomialRing(QQ)
    Q = QuadraticForm(S)
    d = Q.dim()
    L = 2*matrix(g*S)
    c = 2*n + g*S*g
    if p == 2:
        IgZ = (1 - TwoadicIgusaZeta(Q,L,c)) * (1 - t/2) / (1-t)
        return IgZ(2^(1 + d/2 - k))
    else:
        IgZ = (1 - t*IgusaZeta(Q,L,c,p)) * (1 - t/p) / (1 - t)
        return IgZ(p^(1 + d/2 - k))

def LvalueNopole(g,S): #weight k=1 weirdness
    R.<t> = PolynomialRing(QQ)
    Q = QuadraticForm(S)
    d = Q.dim()
    L = 2*matrix(g*S)
    c = g*S*g
    IgZ = (1 - TwoadicIgusaZeta(Q,L,c)) * (1 - t/2) * (1 - 2^(-d/2)*t) / (1 - t)
    return IgZ(2^(d/2))
    
def JacobiLValue(g,b,m,r,n,S,p,k): #the Euler factors in the Jacobi Eisenstein series
    R.<t> = PolynomialRing(QQ)
    Srows = S.nrows()
    Snew = matrix(Srows+1)
    tildeb = S*b
    shiftm = m + b*S*b/2
    for i in range(Srows):
        for j in range(i+1):
            Snew[i,j] = S[i,j]
            Snew[j,i] = S[j,i]
        Snew[i,Srows] = tildeb[i]
        Snew[Srows,i] = tildeb[i]
    Snew[Srows,Srows] = 2*shiftm
    Q = QuadraticForm(Snew)
    d = Q.dim()
    c = 2*n + g*S*g
    Lv = vector(list(g*S) + [r + g*S*b])
    L = -2*matrix(Lv)
    if p == 2:
        IgZ = (1 - TwoadicIgusaZeta(Q,L,c)) * (1 - t/2) / (1 - t)
        return IgZ(2^(2 + (d-1)/2 - k))
    else:
        IgZ = (1 - t*IgusaZeta(Q,L,c,p)) * (1 - t/p) / (1-t)
        return IgZ(p^(2+ (d-1)/2 - k))
    
def JacobiLValueNopole(g,b,m,r,n,p,S): #weight k=3/2
    R.<t> = PolynomialRing(QQ)
    Srows = S.nrows()
    Snew = matrix(Srows+1)
    tildeb = S*b
    shiftm = m + b*S*b/2
    for i in range(Srows):
        for j in range(i+1):
            Snew[i,j] = S[i,j]
            Snew[j,i] = S[j,i]
        Snew[i,Srows] = tildeb[i]
        Snew[Srows,i] = tildeb[i]
    Snew[Srows,Srows] = 2*shiftm
    Q = QuadraticForm(Snew)
    d = Q.dim()
    c = 2*n + g*S*g
    Lv = vector(list(g*S) + [r + g*S*b])
    L = -2*matrix(Lv)
    if p == 2:
        IgZ = (1 - TwoadicIgusaZeta(Q,L,c)) * (1 - t/2) * (1 - 2^(-d/2)*t) / (1 - t)
        return IgZ(2^(d/2))
    else:
        IgZ = (1 - t*IgusaZeta(Q,L,c,p)) * (1 - t/p) * (1 - p^(-d/2)*t) / (1 - t)
        return IgZ(p^(d/2))

@cached_function
def quadratic_L_function__correct(k,D): #quadratic_L_function__exact() makes D squarefree. This function undoes that.
    d = squarefree_part(D)
    val=quadratic_L_function__exact(k,d)
    if d % 4 == 3:
        d = d * 4
    for pe in factor(D):
        p = pe[0]
        if d % p != 0:
            val = val * (1 - kronecker_symbol(d,p)*p^(-k))
    return val

@cached_function
def EisensteinCoefficient(g,n,k,S): #coefficient c(n,g) of Bruinier's Eisenstein series of weight k and Gram matrix S. WARNING: when k = 3/2 or 2 the result is generally not a modular form.
    if not (n + g*S*g / 2).is_integer():
        raise ValueError('bad value of n')
    if k < 1 or not (2*k + QuadraticForm(S).signature()) % 4 == 0:
        raise ValueError('bad value of k')
    if n < 0:
        return 0
    if n == 0:
        if denominator(g) == 1:
            if k != 1:
                return 1
            if k == 1:
                return 1 - WeightOneZeroVal(g,S)
        else:
            return 0
    dets = S.determinant()
    Q = QuadraticForm(S)
    firstfac = 2^k * pi^k * n^(k-1) * (-1)^((2*k + Q.signature())/ 4) / (sqrt(abs(dets))*gamma(k)) #unnecessary
    e = S.nrows()
    dgamma = denominator(g)
    if e % 2 == 1:
        n0 = 1
        for pe in factor(n):
            p = pe[0]
            if 2*dets % p == 0:
                n0 = n0 * p^(pe[1])
            else:
                n0 = n0 * p^(pe[1] % 2)
        bigD = 2*(-1)^((e+1)/2)*n0*dets*dgamma^2
        secondfac = quadratic_L_function__correct(ZZ(k-1/2),bigD) / zeta(2*k - 1)
        thirdfac = 1
        for pe in factor(2*n*dgamma^2*abs(dets)):
            p = pe[0]
            Eulerfactor = Lvalue(g,n,S,p,k)
            thirdfac = thirdfac * (1 - kronecker_symbol(bigD,p)*p^(1/2 - k))*Eulerfactor / (1 - p^(1 - 2*k))
        return QQ((firstfac * secondfac * thirdfac).canonicalize_radical())
    else:
        D = (-1)^(e/2) * dets
        secondfac = 1 / quadratic_L_function__correct(k,D)
        thirdfac = 1
        for pe in factor(2*n*dgamma^2 * abs(dets)):
            p = pe[0]
            Eulerfactor = Lvalue(g,n,S,p,k)
            thirdfac = thirdfac * Eulerfactor / (1 - kronecker_symbol(D,p)*p^(-k))
        return QQ((firstfac * secondfac * thirdfac).canonicalize_radical())

    
def WeightOneZeroVal(oldg,oldS): #corrected constant term in Eisenstein series of weight one
    if oldS.nrows() == 2: #dumb hack
        S = matrix([[oldS[0,0],oldS[0,1],0,0],[oldS[1,0],oldS[1,1],0,0],[0,0,0,1],[0,0,1,0]])
        g = vector([oldg[0],oldg[1],0,0])
    else:
        S = oldS
        g = oldg
    Q = QuadraticForm(S)
    if not frac(g*S*g/2) == 0:
        return 0
    s = pi*(-1)^((2 - Q.signature())  /4) / sqrt(abs(S.determinant()))
    D = -abs(S.determinant())
    s = s / quadratic_L_function__correct(1,D)
    problem = []
    for pe in factor(D):
        p = pe[0]
        try:
            Lval = Lvalue(g,0,S,p,1)
        except ZeroDivisionError:
            problem.append(p)
            Lval = LvalueNopole(g,S)
        s = s * Lval
    modD = D
    for p in problem:
        modD = modD / p^(D.valuation(p))
    s = s * quadratic_L_function__correct(0,modD)
    return s

def DiscriminantGroup(S): #representatives of discriminant group of S
    from itertools import product
    L = []
    List = []
    [D,U,V] = S.smith_form()
    for k in range(0,D.nrows()):
        L.append(list(vector(range(D[k,k])) / D[k,k]))
    for r in product(*L):
        g = V*vector(r)
        for k in range(len(g)):
            g[k] = frac(g[k])
        List.append(g)
    return List

@cached_function
def ReducedDiscriminantGroup(S): #discriminant group of S with only one representative of {+/- g} for each g in the discriminant group
    G = DiscriminantGroup(S)
    L = []
    for g in G:
        u = []
        for k in range(len(g)):
            u.append(frac(-g[k]))
        u = vector(u)
        if not u in L:
            L.append(g)
    return L

def Jacobi(g,b,m,n,r,k,S): #the coefficient (n,r,g) of the Jacobi Eisenstein series for the b-Schroedinger representation of index m, weight k, Gram matrix S. does not check input validity so be careful
    if 4*m*n - r^2 < 0:
        return 0
    if 4*m*n == r^2:
        lam = r / (2*m)
        if lam.is_integer():
            if denominator(g - lam*b) == 1:
                return 1
        return 0
    Q = QuadraticForm(S)
    sig = 2*k + Q.signature()
    dets = S.determinant()
    firstfac = (-1)^(sig/4) * pi^(k - 1/2) * (4*m*n - r^2)^(k - 3/2) / (2^(k-2) * m^(k-1) * gamma(k - 1/2) * sqrt(abs(dets)))
    dbeta = denominator(b)
    dgamma = denominator(g)
    tilden = dbeta^2 * dgamma^2 * (n - r^2 / (4*m))
    e = S.nrows()
    badprimes = [2]
    for pe in factor(tilden):
        badprimes.append(pe[0])
    for pe in factor(m):
        badprimes.append(pe[0])
    for pe in factor(abs(dets)):
        badprimes.append(pe[0])
    badprimes = list(set(badprimes))
    if e % 2 == 0:
        bigD = m*dbeta^2 * (-1)^(e/2 + 1) * tilden * dets * (4*m)^2
        for p in badprimes:
            bigD = bigD * p^2
        secondfac = quadratic_L_function__correct(k-1,bigD) / zeta(2*k - 2)
        thirdfac = 1
        for p in badprimes:
            Eulerfactor = JacobiLValue(g,b,m,r,n,S,p,k)
            thirdfac = thirdfac * Eulerfactor * (1 - p^(1-k)*kronecker_symbol(bigD,p)) / (1 - p^(2 - 2*k))
        return QQ((firstfac * secondfac * thirdfac).canonicalize_radical())
    else:
        bigD = 2*dbeta^2 * (-1)^((e+1)/2) * dets*m
        for p in badprimes:
            bigD = bigD * p^2
        secondfac = 1 / quadratic_L_function__correct(ZZ(k-1/2),bigD)
        thirdfac = 1
        for p in badprimes:
            Eulerfactor = JacobiLValue(g,b,m,r,n,S,p,k)
            thirdfac = thirdfac * Eulerfactor / (1 - p^(1/2 - k) * kronecker_symbol(bigD,p))
        return QQ((firstfac * secondfac * thirdfac).canonicalize_radical())

@cached_function
def PSS(g,b,m,n,k,S): #the coefficient (n,g) of the Poincare square series of index (m,b), weight k, Gram matrix S.
    if not (m + b*S*b / 2).is_integer():
        raise ValueError('bad value of m')
    if not (n + g*S*g / 2).is_integer():
        raise ValueError('bad value of n')
    if k < 3/2 or not (2*k + QuadraticForm(S).signature()) % 4 == 0:
        raise ValueError('bad value of k')
    if n < 0:
        return 0
    assert m >= 0
    s = 0
    for shiftr in range(floor(g*S*b - 2*sqrt(m*n)),floor(2*sqrt(m*n) + g*S*b)+1):
        r = shiftr - g*S*b
        s = s + Jacobi(g,b,m,n,r,k,S)
    if k == 5/2 and n > 0:
        s = s + weight5halfPSSfix(g,b,m,n,S)
    if k == 2 and n > 0: #result is NOT always a modular form
        s = s + weight2PSSfix(g,b,m,n,S)
    if k == 3/2:
        s = s + weight3halfPSSfix(g,b,m,n,S)
    return s

def weight5halfPSSfix(g,b,m,n,S): #add a theta derivative
    dbeta = denominator(b)
    dgamma = denominator(g)
    s = 0
    r = sqrt(4*m*n)
    shiftr = r + g*S*b
    if shiftr.is_integer():
        dets = S.determinant()
        Q = QuadraticForm(S)
        epsilon = 24*n*(-1)^((5 + Q.signature()) / 4) / sqrt(2*m*abs(dets))
        badprimes = [2]
        for pe in factor(m):
            badprimes.append(pe[0])
        for pe in factor(abs(dets)):
            badprimes.append(pe[0])
        badprimes = list(set(badprimes))
        e = S.nrows()
        bigD = 2*m * (-1)^((e+1)/2)* dets
        if bigD.is_square():
            for p in badprimes:
                EulerFactor = JacobiLValue(g,b,m,r,n,S,p,5/2)
                epsilon = epsilon * EulerFactor * (1 - p^(-1)) / (1 - p^(-2))
            s = s + epsilon
    shiftrdown = r - g*S*b
    if shiftrdown.is_integer():
        dets = S.determinant()
        Q = QuadraticForm(S)
        epsilon = 24*n*(-1)^((5 + Q.signature()) / 4) / sqrt(2*m*abs(dets))
        badprimes = [2]
        for pe in factor(m):
            badprimes.append(pe[0])
        for pe in factor(abs(dets)):
            badprimes.append(pe[0])
        badprimes = list(set(badprimes))
        e = S.nrows()
        bigD = 2*m * (-1)^((e+1)/2)* dets
        if bigD.is_square():
            for p in badprimes:
                EulerFactor = JacobiLValue(g,b,m,-r,n,S,p,5/2)
                epsilon = epsilon * EulerFactor * (1 - p^(-1)) / (1 - p^(-2))
            s = s + epsilon
    return s

def helper_B(x):
    return x - (floor(x) - floor(-x)) / 2

def DimensionFormula(k,S): #returns a list [dim M, dim S], where dim M is the dimension of modular forms of weight k and dim S is the dimension of cusp forms.
    r = QuadraticForm(S).signature()
    assert (2*k + r) % 4 == 0
    if k < 0:
        return 0
    if k <= 2: #slow
        L1 = []
        L2 = []
        Spans = CuspSpan(k+12,S)
        for g in ReducedDiscriminantGroup(S):
            offset = frac(g*S*g/2)
            n = 1 - offset
            U = []
            V = []
            for mb in Spans:
                m = mb[0]
                b = vector(mb[1])
                E = EisensteinCoefficient(g,n,k+12,S)
                P = PSS(g,b,m,n,k+12,S)
                U.append(E-P)
                if offset != 0:
                    V.append(E-P)
            L1.append(U)
            if offset != 0:
                L2.append(V)
        return [len(Spans)-matrix(L2).rank(), len(Spans)-matrix(L1).rank()]
    G = DiscriminantGroup(S)
    dets = abs(S.determinant())
    size = 0
    gausssum1 = 0
    gausssum2 = 0
    gausssum3 = 0
    alpha4 = 0
    B1 = 0
    B2 = 0
    for g in G:
        d = denominator(g)
        gSg = g*S*g
        B1 = B1 + helper_B(gSg / 2)
        if 2 % d == 0:
            size = size + 1
            B2 = B2 + helper_B(gSg / 2)
        else:
            size = size + 1/2
        if (gSg / 2).is_integer():
            if 2 % d == 0:
                alpha4 = alpha4 + 1
            else:
                alpha4 = alpha4 + 1/2
        gausssum1 = gausssum1 + exp(pi*I*gSg)
        gausssum2 = gausssum2 + exp(2*pi*I*gSg)
        gausssum3 = gausssum3 + exp(-3*pi*I*gSg)
    alphaT = (alpha4 + B1 + B2) / 2
    moddim = size*(k - 1)/12 - (exp(pi*I*(4*k + 3*r - 10) / 12) * (gausssum1 + gausssum3)).real() / (3*sqrt(3*dets)) + exp(pi*I*(2*k + r) / 4) * gausssum2.real() / (4*sqrt(dets)) + alphaT
    return [round(moddim.n()),round((moddim-alpha4).n())]

def CuspSpan(k,S): #span cusp forms in weight at least 5/2. Output is a list of [m,b] where PSS(...,b,m,...) - E(...) is a basis. 
    if k < 5/2 or not (2*k + QuadraticForm(S).signature()) % 4 == 0: #todo: add weight 3/2 and weight 2 using Ehlen and Skoruppa's algorithm
        raise ValueError('bad value of k')
    [dimMk,dimSk] = DimensionFormula(k,S)
    if dimSk == 0:
        return []
    sturmbound = floor(k / 12) + 1
    L = []
    G = ReducedDiscriminantGroup(S)
    notdone = True
    m = 1
    bs = []
    oldrank = 0
    while notdone:
        for b in G:
            U = []
            m_offset = frac(b*S*b/2)
            for g in G:
                offset = frac(g*S*g/2)
                for n in range(1,sturmbound+1):
                    U.append(EisensteinCoefficient(g,n-offset,k,S) - PSS(g,b,m-m_offset,n-offset,k,S))
            newL = list(L)
            newL.append(U)
            newrank = matrix(newL).rank()
            if newrank > oldrank:
                L = newL
                oldrank = newrank
                bs.append([m-m_offset,b])
            if oldrank == dimSk:
                notdone = False
                break
        m = m + 1
    return bs

@cached_function
def FundaSol(D): #fundamental solution to pell equation x^2 - Dy^2 = 1
    K.<alpha> = QuadraticField(D)
    G = K.unit_group()
    u = G.1
    u = K(u)
    if u < 0:
        u = -u
    if u < 1:
        u = u^(-1)
    if u.norm() == -1:
        u = u^2
    j = 1
    newu = u
    inZsqrtD = False
    while not inZsqrtD:
        x = realpart(D,newu)
        y = imagpart(D,newu)
        if x.is_integer() and y.is_integer():
            inZsqrtD = True
            break
        newu = newu*u
    return newu


def ismultiple(D,mu,x): #is x a multiple of some element of mu
    K.<alpha> = QuadraticField(D)
    f = D / D.squarefree_part()
    for y in mu:
        z = x/y
        zbar = z.norm() / z
        if realpart(D,z).is_integer() and imagpart(D,z).is_integer():
            return True
    return False

def mus(D,C): #generate representatives of orbits of solutions to a pell-type equation
    K.<alpha> = QuadraticField(D)
    R.<x> = PolynomialRing(K)
    RQ.<y> = PolynomialRing(QQ)
    u = FundaSol(D)
    bound = floor(2*sqrt(C*abs(u)).n()) + 1
    mu = []
    if is_square(C):
        mu.append(sqrt(C))
    for lam in range(bound):
        p = x^2 + lam*x + C
        if not p.is_irreducible() and RQ(p).is_irreducible():
            mu.append(p.roots()[0][0])
            mu.append(p.roots()[1][0])
    return mu


def minimus(D,C): #generate a minimal set of minimal representatives of orbits of solutions to a pell-type equation
    minimuus = []
    muus = mus(D,C)
    for muu in muus:
        if muu != 0:
            if muu < 0:
                muu = -muu
            eps0 = FundaSol(D)
            notdone = True
            muconj = C / muu
            if muu < muconj:
                [muu,muconj] = [muconj,muu]
            while notdone:
                if muu < muconj * eps0^2:
                    notdone = False
                    break
                muu = muu*eps0^(-1)
                muconj = muconj*eps0
        if not ismultiple(D,minimuus,muu) and not ismultiple(D,minimuus,muconj):
            minimuus.append(muu)
    return minimuus

def realpart(D,x):
    K.<alpha> = QuadraticField(D)
    return K(x).trace()

def imagpart(D,x):
    K.<alpha> = QuadraticField(D)
    return QQ((2*x - realpart(D,x)) / alpha)

def congmu(D,muu,modulus,remainder): #find a minimal solution to Pell-type equation in the orbit of muu congruent to remainder mod modulus
    im_modulus = modulus * D / D.squarefree_part()
    eps0 = FundaSol(D)
    newmuu = eps0^(-1) * muu
    L = []
    while True:
        newmuu = eps0*newmuu
        R_part = realpart(D.squarefree_part(),newmuu)
        I_part = imagpart(D.squarefree_part(),newmuu)
        if R_part % modulus == remainder:
            return newmuu
        RI_parts = [R_part % modulus, I_part % im_modulus]
        if RI_parts in L:
            return 0
        L.append(RI_parts)
        
def congfundasol(D,muu,modulus): #congruent fundamental solution for modulus and orbit of muu
    if muu == 0:
        return 0
    eps0 = FundaSol(D)
    R_part = realpart(D,muu)
    I_part = imagpart(D,muu)
    eps = eps0
    newmuu = muu
    while True:
        newmuu = newmuu*eps0
        newR_part = realpart(D,newmuu)
        newI_part = imagpart(D,newmuu)
        if (R_part - newR_part) % modulus == 0:
            return eps
        eps = eps*eps0

def corrector(D,C,muu,eps0): #the correction term associated to the class "muu" in the weight two formula for PSS
    if muu == 0:
        return 0
    K.<alpha> = QuadraticField(D)
    muconj = C / muu
    c = (eps0 + eps0^(-1))/2
    d = (eps0 - c) / alpha
    e = (muu + muconj) / 2
    f = (muu - e) / alpha
    s = 4*D*(f + 2*d*e / (eps0 - 1).norm())
    if ismultiple(D,[muu],muconj) or muu.is_integer():
        s = s/2
    return s

def weight2PSSfix(g,b,m,n,S): #difference between weight 2 PSS coefficient and naive sum
    s = 0
    if n <= 0:
        return 0
    dbeta = denominator(b)
    dgamma = denominator(g)
    D = abs(S.determinant())
    sig = QuadraticForm(S).signature()
    if not is_square(D):
        modulus = dbeta*dgamma
        remainder1 = (-dbeta*dgamma*(g*S*b)) % modulus
        remainder2 = (dbeta*dgamma*(g*S*b)) % modulus
        C = dbeta^2 * dgamma^2 * m*n
        muus = minimus(D,C)
        for muu in muus:
            congmuu1 = congmu(D,muu,modulus,remainder1)
            congmuu2 = congmu(D,muu,modulus,remainder2)
            if congmuu1 != 0 or congmuu2 != 0:
                eps1 = congfundasol(D,congmuu1,modulus)
                eps2 = congfundasol(D,congmuu2,modulus)
                if congmuu1 != 0:
                    r = realpart(D,congmuu1) / (denominator(b)*denominator(g))
                if congmuu2 != 0:
                    r = -realpart(D,congmuu2) / (denominator(b)*denominator(g))
                badprimes = [2]
                if not r^2 - 4*m*n == 0:
                    for pe in factor(r^2 / m - 4*n):
                        badprimes.append(pe[0])
                for pe in factor(m):
                    badprimes.append(pe[0])
                for pe in factor(D):
                    badprimes.append(pe[0])
                badprimes = list(set(badprimes))
                J = QQ(1)
                for p in badprimes:
                    J = J*JacobiLValue(g,b,m,r,n,S,p,2) / (1 + p^(-1))
                corr = (corrector(D,C,congmuu1,eps1) + corrector(D,C,congmuu2,eps2)) / 2
                s = s + J*corr
        s = -6*(-1)^((4+sig)/4)*s / (m*D*denominator(b)*denominator(g))
        return QQ(s)
    else:
        for tilder in range(-2*dbeta*m*(n+1)-1,2*dbeta*m*(n+1)+1):
            r = tilder - frac(g*S*b)
            tilden = dbeta^2 * dgamma^2 * (n - r^2 / (4*m))
            badprimes = [2]
            if tilden != 0:
                for pe in factor(tilden):
                    badprimes.append(pe[0])
            for pe in factor(m):
                badprimes.append(pe[0])
            for pe in factor(D):
                badprimes.append(pe[0])
            fac = 1
            badprimes = list(set(badprimes))
            for p in badprimes:
                fac = fac * (1 - p^(-1)) / (1 - p^(-2))
                fac = fac * JacobiLValue(g,b,m,r,n,S,p,2)
            if is_square(D* (r^2 - 4*m*n)):
                if r^2 - 4*m*n != 0:
                    s = s - 6*(-1)^(sig/4) * fac * (abs(r) - sqrt(r^2 - 4*m*n))
                else:
                    s = s - 3*(-1)^(sig/4) * fac * (abs(r) - sqrt(r^2 - 4*m*n))
        return QQ(s / (m*sqrt(D)))
    
def weight3halfPSSfix(g,b,m,n,S): #add theta
    positive_r = sqrt(4*m*n)
    negative_r = -sqrt(4*m*n)
    dgamma = denominator(g)
    dbeta = denominator(b)
    e = S.nrows()
    s = 0
    if (positive_r + frac(g*S*b)).is_integer():
        dets = S.determinant()
        Q = QuadraticForm(S)
        epsilon = (-1)^((1 - Q.signature()) / 4) / sqrt(2*m*abs(dets))
        dbeta = denominator(b)
        dgamma = denominator(g)
        badprimes = [2]
        bigD = 2*dbeta^2 * m * (-1)^((e+1)/2)* dets
        for pe in factor(m):
            badprimes.append(pe[0])
        for pe in factor(abs(dets)):
            badprimes.append(pe[0])
        badprimes = list(set(badprimes))
        problem = []
        for p in badprimes:
            try:
                Lval = JacobiLValue(g,b,m,positive_r,n,S,p,3/2)
            except ZeroDivisionError:
                problem.append(p)
                Lval = JacobiLValueNopole(g,b,m,positive_r,n,p,S)
            epsilon = epsilon * Lval
        modD = bigD
        for p in problem:
            modD = modD / p^(bigD.valuation(p))
        epsilon = epsilon * quadratic_L_function__correct(0,modD) / quadratic_L_function__correct(1,bigD)
        s = s + epsilon*pi
    if (negative_r + frac(g*S*b)).is_integer() and n != 0:
        dets = S.determinant()
        Q = QuadraticForm(S)
        epsilon = (-1)^((1 - Q.signature()) / 4) / sqrt(2*m*abs(dets))
        dbeta = denominator(b)
        dgamma = denominator(g)
        badprimes = [2]
        bigD = 2*dbeta^2 * m * (-1)^((e+1)/2)* dets
        for pe in factor(m):
            badprimes.append(pe[0])
        for pe in factor(abs(dets)):
            badprimes.append(pe[0])
        badprimes = list(set(badprimes))
        problem = []
        for p in badprimes:
            try:
                Lval = JacobiLValue(g,b,m,negative_r,n,S,p,3/2)
            except ZeroDivisionError:
                problem.append(p)
                Lval = JacobiLValueNopole(g,b,m,negative_r,n,p,S)
            epsilon = epsilon * Lval
        modD = bigD
        for p in problem:
            modD = modD / p^(bigD.valuation(p))
        epsilon = epsilon * quadratic_L_function__correct(0,modD) / quadratic_L_function__correct(1,bigD)
        s = s + epsilon*pi
    return -s

def EisensteinBound(S): #k = e/2 + 2. S should be NEGATIVE definite. Returns a number t such that the absolute values of coefficients of Eisenstein series are bounded below by t*n^(k-1)
    dets = S.determinant()
    e = S.nrows()
    k = e/2 + 2
    if e % 2 == 0:
        D = (-1)^(e/2) * dets
        k = ZZ(k)
        t = (2*pi)^k * (2 - (1 - 2^(1-k))*zeta(k-1)) / (sqrt(abs(dets)) * gamma(k) * quadratic_L_function__correct(k,D))
        for pe in factor(2*dets):
            p = pe[0]
            t = t * (1 - p^(-1))
        return t.n()
    if e % 2 == 1:
        t = (2*pi)^k / (5 * (1 - 2^(-3-e)) * gamma(e/2 + 2) * zeta(e/2 + 3/2) * sqrt(abs(dets)))
        for pe in factor(dets):
            p = pe[0]
            t = t * (1 - p^(-1)) / (1 - p^(-3 - e))
        return t.n()

def BorcherdsPrincipalPart(wt,S): #assume S is negative definite!! Tries to find Borcherds products of weight "wt" and returns their principal parts
    from itertools import product
    e = S.nrows()
    k = e/2 + 2
    if e % 2 == 0:
        k = ZZ(k)
    Spans = CuspSpan(k,S)
    vs = []
    for i in range(len(Spans)):
        vs.append([0])
    gn = []
    r = [-wt]
    bound = (2*wt/EisensteinBound(S))^(1 / (e/2 + 1))
    for g in ReducedDiscriminantGroup(S):
        multiplyby2 = False
        for l in range(len(g)):
            if not (2*g[l]).is_integer():
                multiplyby2 = True
        offset = frac(g*S*g/2)
        for n in range(ceil(offset),floor(bound + offset)+1):
            if n - offset > 0:
                E = EisensteinCoefficient(g,n-offset,k,S)
                u = -E/2
                if multiplyby2:
                    u = u * 2
                if True:
                    r.append(u)
                    gn.append([g,n-offset])
                    for j in range(len(Spans)):
                        mb = Spans[j]
                        m = mb[0]
                        b = vector(mb[1])
                        coeff = E - PSS(g,b,m,n-offset,k,S)
                        if multiplyby2:
                            coeff = coeff*2
                        vs[j].append(coeff)
    positive = []
    for i in range(len(r)-1): #positivity condition
        v = [0]*len(r)
        v[i+1] = 1
        gni = gn[i]
        g = gni[0]
        n = gni[1]
        for l in range(1,floor(bound / n) + 1):
            newg1 = g*l
            newg2 = g*l
            for g_index in range(e):
                newg1[g_index] = frac(newg1[g_index])
                newg2[g_index] = frac(-newg2[g_index])
            try:
                new_index = gn.index([newg1,n*l^2])
            except ValueError:
                try:
                    new_index = gn.index([newg2,n*l^2])
                except ValueError:
                    new_index = -1
            if new_index >= 0:
                v[new_index + 1] = 1
        positive.append(v)
    vs.append(r)
    points = Polyhedron(ieqs = positive, eqns = vs).integral_points()
    Output = []
    for point in points:
        Output.append(principal_part(point,gn,e,wt))
    return Output

def principal_part(beta,gn,e,wt): #formats principal parts as an easier-to-read string
    e0 = '(' + (''.join(str(i) + ',' for i in list(vector([0]*e))))[:-1] + ')'
    output = str(2*wt) + "*e_" + e0
    for j in range(len(beta)):
        u1 = '(' + (''.join(str(i) + ',' for i in list(gn[j][0])))[:-1] + ')'
        if beta[j] != 0:
            if beta[j] != 1:
                output = output + " + " + str(beta[j]) + "*q^(" + str(-gn[j][1]) + ")*e_" + u1
            if beta[j] == 1:
                output = output + " + q^(" + str(-gn[j][1]) + ")*e_" + u1
        multiplyby2 = False
        g = vector([QQ(0)]*e)
        for i in range(e):
            g[i] = frac(-gn[j][0][i])
            if g[i] != gn[j][0][i]:
                multiplyby2 = True
        u2 = '(' + (''.join(str(i) + ',' for i in list(g)))[:-1] + ')'
        if multiplyby2 and beta[j] != 0:
            if beta[j] != 1:
                output = output + " + " + str(beta[j]) + "*q^(" + str(-gn[j][1]) + ")*e_" + u2
            if beta[j] == 1:
                output = output + " + q^(" + str(-gn[j][1]) + ")*e_" + u2
    return output

def parseInput(inputs,S): #parses the "easy-to-read" string back into a list. This could be done better.
    inputsplit = inputs.split()
    constant_term = inputsplit.pop(0)
    const = sage_eval(constant_term[:constant_term.find('*')])
    RDS = ReducedDiscriminantGroup(S)
    R.<q> = PowerSeriesRing(QQ)
    U = [R(0)]*(len(RDS))
    U[0] = R(const)
    order = 0
    for i in range(len(inputsplit)):
        if i % 2 == 1:
            term = inputsplit[i]
            N = term.find('e')
            if term[0] == 'q':
                exponent = sage_eval(term[3:N-2])
                constg = 1
            else:
                h = 0
                while term[h] != 'q':
                    h = h + 1
                constg = sage_eval(term[:h-1])
                exponent = sage_eval(term[h+3:N-2])
            if S.nrows() > 1:
                g = vector(sage_eval(term[N+2:len(term)]))
            if S.nrows() == 1:
                g = vector([sage_eval(term[N+2:len(term)])])
            offset = frac(g*S*g/2)
            exponent = ZZ(exponent - offset)
            order = min(order,exponent)
            try:
                indexg = RDS.index(g)
                U[indexg] = U[indexg] + constg*q^exponent
            except ValueError:
                pass
    U.append(-order)
    return U


def WeylVector(inputs,S): #assume S is negative definite! Corrected formula. Outputs the Weyl vector associated to an "inputs" (a string in the format of the result of a BorcherdsPrincipalPart) associated to the Gram matrix S. 
    from itertools import product
    min_eig = -max(S.eigenvalues()).n()
    e = S.nrows() + 2
    v = vector(QQ,[0]*e)
    R.<q> = PowerSeriesRing(QQ)
    RDS = ReducedDiscriminantGroup(S)
    U = parseInput(inputs,S)
    order = U.pop()
    F = R(0)
    coffsum = U[0][0] / 24
    vecsum = vector(QQ,[0]*(e-2))
    for i in range(len(U)):
        g = RDS[i]
        g2 = vector(QQ,[0]*(e-2))
        for j in range(e-2):
            g2[j] = frac(1 - g[j])
        multby2 = False
        for j in range(e-2):
            if not (2*g[j]).is_integer():
                multby2 = True
        theta = R(0)
        L = []
        for j in range(e-2):
            L.append(range(-floor(sqrt(order/min_eig))-3,ceil(sqrt(order/min_eig))+3))
        for x in product(*L):
            u = vector(x) + g
            u2 = vector(x) + g2
            theta = theta + q^(ZZ(-u*S*u/2 + frac(u*S*u/2)))
            if multby2:
                theta = theta + q^(ZZ(-u2*S*u2/2 + frac(u2*S*u2/2)))
            positive = True
            for j in range(len(u)):
                if (S*u)[j] < 0:
                    break
                if (S*u)[j] > 0:
                    positive = False
                    break
            if u != vector([0]*S.nrows()):
                coffsum = coffsum + U[i][u*S*u/2 - frac(u*S*u/2)] / 24
                if positive: #seems to be correct now
                    vecsum = vecsum + U[i][u*S*u/2 - frac(u*S*u/2)] * u / 2
            positive = True
            for j in range(len(u)):
                if (S*u2)[j] < 0:
                    break
                if (S*u2)[j] > 0:
                    positive = False
                    break
            if u2 != u and u2 != vector([0]*S.nrows()):
                coffsum = coffsum + U[i][u2*S*u2 / 2 - frac(u2*S*u2/2)] / 24
                if positive:
                    vecsum = vecsum + U[i][u2*S*u2/2 - frac(u2*S*u2/2)] * u2 / 2
        theta = theta + O(q^(order+1))
        F = F + U[i]*theta
    return vector([coffsum]+list(vecsum)+[-(F*eisenstein_series_qexp(2,order+1))[0]])

def InputFunction(inputs,prec,S): #the input function into Borcherds lift up to precision O(q^prec). "Inputs" should be formatted as in the output of BorcherdsPrincipalPart. S should be a *negative-definite* Gram matrix. 
    R.<q> = PowerSeriesRing(QQ)
    RDS = ReducedDiscriminantGroup(S)
    U = parseInput(inputs,S)
    V = [R(0)]*len(U)
    mult = ZZ(U.pop())
    weight = 12*mult -S.nrows() / 2
    if S.nrows() % 2 == 0:
        weight = ZZ(weight)
    for i in range(len(U)):
        V[i] = delta_qexp(prec)^mult * U[i]
    Spans = CuspSpan(weight,-S)
    dimn = ZZ(mult * len(RDS))
    M = matrix(QQ,dimn,len(Spans)+1)
    v = vector([0]*dimn)
    for i in range(len(RDS)):
        for j in range(mult):
            indexval = i*mult + j
            g = RDS[i]
            offset = frac(g*(-S)*g/2)
            if offset == 0:
                offset = 1
            n = 1 + j - offset
            E = EisensteinCoefficient(g,n,weight,-S)
            M[indexval,0] = E
            for k in range(len(Spans)):
                mb = Spans[k]
                m = mb[0]
                b = vector(mb[1])
                M[indexval,k+1] = PSS(g,b,m,n,weight,-S)
            v[indexval] = (V[i])[j]
    combo = M.solve_right(v)
    outputs = []
    for g in RDS:
        f = R(0)
        offset = frac(g*(-S)*g/2)
        for n in range(prec + mult):
            s = combo[0]*EisensteinCoefficient(g,n-offset,weight,-S)
            for i in range(len(Spans)):
                mb = Spans[i]
                m = mb[0]
                b = vector(mb[1])
                s = s + combo[i+1]*PSS(g,b,m,n-offset,weight,-S)
            f = f + s*q^n
        f = f + O(q^(prec + mult))
        outputs.append([g,-offset,f / delta_qexp(prec+mult)^mult])
    return outputs