I have been working on this program, trying to implement the algorithm for level N on page 9 of this paper, but I am running into trouble, and I can't figure out why. When I run this program, if I set m, which is basically the precision, to 6, it runs fine. However, if I increase the precision to m=12, I get an error right at the very end, when trying to find the matrix A. It says that there is no solution. However, according to the mathematical proof in the paper, this should not be the case, so I don't know what is wrong.

I have spent hours troubleshooting this, and cannot find any errors. Can anyone help with this? Can anyone tell where my implementation differs from the algorithm I am trying to implement?

I'm not sure if this is the kind of question that should be asked on this, but I just can't fix this and need some advice, and don't know where else to ask.

My code is as follows:

given N>=2, p>= 5 prime not dividing N, integer k and positive integer m:

p=5
k=10
m=12
N=2

Compute the k0 and j

((k0,j),)=Polyhedron(ieqs=[[0,1,0], [p-2,-1,0]], eqns=[[-k, 1,p-1]]).integral_points()

Compute n

n=floor(((p+1)/(p-1))*(m-1))

for i=0,...,n, compute di, the dimension of the space of classical modular forms of level 1 and weight k0+i(p-1)

d=list('d_%d' % s for s in range(0,n))
for i in range(0,n):
    d[i]=dimension_modular_forms(Gamma1(N),k0+i*(p-1)) #not sure if I should
    #use dimension_cusp_forms instead

Compute the m_i:

M=list(var('m_%d' % s) for s in range(0,n))
M[0]=d[0]
for i in range(1,n):
    M[i]=d[i]-d[i-1]

Compute lo (l in the paper):

lo=sum(M[i] for i in range(0,n))

Compute the working precision mp (m' in the paper):

mp=m+floor(n/(p+1))

compute the sturm bound l (l' in the paper)

l=ModularForms(Gamma1(N),k0+n*(p-1)).sturm_bound()

STEP 2

For each 0<=i<=n, denote by M_k0+i(p-1) a row reduced basis of q-expansions in Z[[q]]/(p^mp,q^lp) of the space of classical modular forms of weight k0+i(p-1) and level N

R.<q> = PowerSeriesRing(ZZ,default_prec=2*(p*l))
Poly.<x>=PolynomialRing(ZZ)

First, define a reduction modulo (p^mp,q^(l*p)):

def ordred2(a,p):
    b=Poly(a)
    b+= x^(2*(l*p)) #this is just to make absolutely sure we have enough coefficients to avoid issues with list index
    ordred2_b=b.coefficients(sparse=False)
    ordred2_coeff=list('ordred2_coeff_%d' % d for d in range(0,(l*p)))
    for i in range(0,(l*p)):
        ordred2_coeff[i]=Mod(ordred2_b[i],p^mp)
    return R(ordred2_coeff)

Now, define a function that returns the row-reduced basis of q-expansions:

def Mrr(k0,i):
    return ModularForms(Gamma1(N),k0+i*(p-1)).q_echelon_basis(prec=(l*p))

Now calculate all the bases, reduce modulo (p^mp,q^(l*p)), and store them in a list:

ListBasis=list('M_%d' % s for s in range(0,n))
for i in range(0,n):
    if k0+i*(p-1) == 0:
        Mrri=list('Mrri_%d' % s for s in range(0,1))
        Mrri[0]=1
        ListBasis[i]=Mrri
    else:
        Mrri=list('Mrri_%d' % s for s in range(0,len(Mrr(k0,i))))
        for y in range(0,len(Mrr(k0,i))):
            Mrri[y]=ordred2(Mrr(k0,i)[y],p) #the ordred2 is to do the reduction modulo (p^mp,q^(l*p))
        ListBasis[i]=Mrri

Now define W(i). If there is a mistake, it is probably in this part. The way it is defined in the paper is: for all i, set W(i)=[]. Then, for each 1<=w<=d[i], if the degree of the lowest term of the w-th element in M_(k0+i(p-1)) does not occur as the degree of the lowest term of ANY element in M_(k0+(i-1)(p-1))

recall that we calculated these M earlier, as elements of ListBasis. M_(k0+i(p-1))=ListBasis[i]

basically we add it when the degree of the lowest term is "new"

def W(i):
    F=list('f_%d' % s for s in range(0,M[i]))
    t=0
    if i==0:
        F=ListBasis[0] #since there is nothing to compare it to, but clearly the degree of the lowest term is different
                       #than N/A, so it should be added to W(0)
    else:
        while t<M[i]: #M[i] is the difference in degree between two successive spaces of modular forms, so it will be
                      #the number of new degrees of lowest terms that appear. I use < because indexing starts at 0
            for w in range(0,d[i]):
                if all(R(ListBasis[i][w]).polynomial().ord() != R(ListBasis[i-1][y]).polynomial().ord() for y in range(0,d[i-1])):
                    #this checks the degrees of the lowest terms
                    F[t]=ListBasis[i][w]
                    t+=1
                else:
                    asdfggfsbgb=1 #placeholder for it to do nothing, just having Else gives me a syntax error
    return F

STEP 3

compute the q-expansion in z[[q]]/(p^m',q^lp) of the eisenstein series E_p-1(q)

eisen=ordred2(eisenstein_series_qexp(p-1,p*l,normalization='constant'),p)

Now define the elements e_i,s from the paper, and do so by defining them in a double list:

Matrixe=[[1 for s in range(0,M[i])] for i in range(0,n)]

for i in range(0,n):
        if M[i]==0:
            adsaffdsa=3 #placeholder for it to do nothing, just having Else gives me a syntax error
        else:
            for s in range (0,M[i]):
                Matrixe[i][s]=ordred2((ordred2(p^floor(i/(p+1))*(eisen)^(-i),p))*(ordred2(W(i)[s],p)),p)

Step 4

Define eisenq = E_p-1(q^p)

eisenq=ordred2(eisenstein_series_qexp(p-1,l*p,normalization='constant').subs({q:q^p}),p)

Now, use this to define G(q) as in the paper:

def G(q):
    return eisen/eisenq

Now, define the u_i,s as in the paper:

Matrixu=[[1 for s in range(0,M[i])] for i in range(0,n)]

for i in range(0,n):
        if M[i]==0:
            adsaffdsa=3 #placeholder for it to do nothing, just having Else gives me a syntax error
        else:
            for s in range (0,M[i]):
                Matrixu[i][s]=ordred2((G(q)^j)*Matrixe[i][s],p)

Step 5

calculate the t_i,s as in the paper:

Matrixt=[[1 for s in range(0,M[i])] for i in range(0,n)]
for i in range(0,n):
        if M[i]==0:
            adsaffdsa=3 #placeholder for it to do nothing, just having Else gives me a syntax error
        else:
            for s in range (0,M[i]):
                abcd=(R(Matrixu[i][s])+q^(l*p+1)).polynomial().coefficients(sparse=False)
                List_t=list('M_%d' % s for s in range(0,l))
                for y in range(0,l):
                    List_t[y]=abcd[p*y] #this whole section is applying the operator U_p to the elements of Matrixu
                Matrixt[i][s]=R(List_t)

Step 6

Define the matrix T

First, define the base ring:

S=Integers(p^mp)

Now initialize the matrix:

T=matrix(ring=S,nrows=lo,ncols=l)

now define a new function to reduce modulo q^l. This will return a list of coefficients, rather than a polynomial.

def ordred3(a,p):
    b=Poly(a)
    b+= x^(2*(l*p))
    ordred3_b=b.coefficients(sparse=False)
    ordred3_coeff=list('ordred3_coeff_%d' % s for s in range(0,(l)))
    for i in range(0,(l)):
        ordred3_coeff[i]=Mod(ordred3_b[i],p^mp)
    return ordred3_coeff

Now calculate the entries of T: Each row of T correcponds to a non-zero element of Matrixt, and the l entries in that row correspond to the coefficients

jv=0 #this will keep track of the row of the matrix T
i=0 #this will keep track of the i in t_i,s
s=0 #this will keep track of the s in t_i,s
#i could theoretically be increased past the size of Matrixt, but this will not happen because the math is such that
#jv will reach lo and the loop will terminate either before i reaches its limit or at the same time it reaches its limit
while jv<lo:
        if M[i]==0:
            i += 1 
        else:
            while s<M[i]:
                for k in range(0,l): # this will keep track of the column, and thus also the degree of the corresponding
                                     # coefficient
                    T[jv,k]=ordred3(Matrixt[i][s],p)[k]
                s+=1
                jv+=1
            else:
                s=0
                i+=1

Define the matrix E, basically the same method as for T

E=matrix(ring=S,nrows=lo,ncols=l)

j=0
i=0
s=0
while j<lo:
    if M[i]==0:
        i += 1
    else:
        while s<M[i]:
            for k in range(0,l):
                E[j,k]=ordred3(Matrixe[i][s],p)[k]
            s+=1
            j+=1

        else:
            s=0
            i+=1

Solve T=AE for A: If this has no solution, multiply by P and try again.

try:
    E.solve_left(T)
    A=E.solve_left(T)
except:
    E.solve_left(p*T)
    A=E.solve_left(p*T)

Define one last function for reduction:

def ordred4(b,p):
    ordred4_b=Poly(b).coefficients(sparse=False)
    ordred4_coeff=list('ordred4_coeff_%d' % s for s in range(0,len(ordred4_b)))
    for i in range(0,len(ordred4_b)):
        ordred4_coeff[i]=Mod(ordred4_b[i],p^m)
    return Poly(ordred4_coeff)

Calculate the determinant of Id-A*x, where x is just a variable, Modulo p^m

ordred4((identity_matrix(S,lo)-A*x).determinant(),p)