<jupyter_start><jupyter_text>This sage worksheet accompanies my paper  (joint with <NAME> and <NAME>)  ['The Manin constant and the modular degree'](https://arxiv.org/abs/1911.09446)   Lattices of modular formsHere we check that for $N=40$ and $N=80$, the lattices of modular forms $S_2(\Gamma_0(N),\mathbb Z)$ and $H^0(X_0(N),\Omega)$ are not equal. If they were we would have been able to extend Proposition 6.10 and show that $X_0(N)$ has rational singularities for all levels $N$ with $\mathrm{val}_2(N)=3$ and $\mathrm{val}_2(N)=4$.  Requirements: SageMath and the function get_expansions from the python package [products-of-eisenstein-series](https://github.com/michaelneururer/products-of-eisenstein-series) (joint work with <NAME>) to calculate the Fourier expansions.<jupyter_code>def hecke_basis_matrix(N,k=2):
    r"""
    INPUT: 
            - int N, level
            - int k, weight
    OUTPUT: Matrix that contains the Fourier expansions of the basis of normalised Hecke eigenforms of S_2(Gamma_0(N))
    WARNING: This is only implemented if all Hecke eigenforms have rational Fourier coefficients.
    """
    sturm = sturm_bound(N,k)
    Pow.<q> = PowerSeriesRing(QQ)
    res = []
    for d in N.divisors():
        try:
            eigenspace = Newforms(d,k)
        except ValueError:
            raise NotImplementedError('This command only works if all Hecke eigenforms of level N are rational')
        for n in eigenspace:
            for m in ZZ(N/d).divisors():
                res.append(map(lambda x: m^(k/2)*x,(n.qexp(sturm+1)(q^m)+O(q^(sturm+1))).padded_list()))
    return Matrix(res)

def cusp_form_basis_matrix(N,k=2):
    r"""
    INPUT: 
            - int N, level
            - int k, weight
    OUTPUT: Matrix that contains the Fourier expansiosn of the basis of cusp forms with integral Fourier coefficients produced by Sage
    """
    sturm = sturm_bound(N,k)
    return Matrix([c.coefficients(range(sturm+1)) for c in CuspForms(N,k).basis()])<jupyter_output><empty_output><jupyter_text>The case N = 40<jupyter_code>H_40 = hecke_basis_matrix(40); 
print 'H_40 = '
print H_40
e1,e2,e3 = CuspForms(40).basis();
print 'e1 = ', e1, ', e2 = ', e2, ', e3 = ', e3
M_40 = cusp_form_basis_matrix(40); 
print 'M_40 = '
print M_40
B_40 = H_40.solve_left(M_40); B_40
B_40*H_40 == M_40<jupyter_output><empty_output><jupyter_text>The matrix $H_{40}$ contains as rows the Fourier expansions of $g$, the newform of level $20$, $g_2 = 2g(q^2)$, and $f$, the newform of level $40$. The matrix $M_{40}$ contains the Fourier expansions of a basis $e_1 = q+O(q^4), e_2 = q^2 + O(q^4), e_3 = q^3 + O(q^4)$ of $S_2(\Gamma_0(40),\mathbb{Z})$.From the matrix $B_{40}$, that solves $B_{40}\cdot H_{40} = M_{40}$ we can read off that $e_1 = f, e_2 = g_2$, and $e_3 = \frac{f-g}{2}$. Let $L = S_2(\Gamma_0(40),\mathbb{Z})$ and $$L' = \mathbb{Z}g \oplus \mathbb{Z}g_2 \oplus f\subseteq H^0(X_0(40),\Omega)\subset S_2(\Gamma_0(40),\mathbb{Z}).$$In order to show $L' = H^0(X_0(40),\Omega)$ it suffices to show $e_3\notin H^0(X_0(40),\Omega)$. We do this by applying Proposition 5.13 and that $v_2(e_3|_{1/2})=-2$.<jupyter_code>count = 1
valuations_12 = []
for e in [e1,e2,e3]:
    print 'Calculating Fourier expansions for e{}'.format(count)
    print 'Expansion at cusp 1/2:'
    s2 = get_expansion(e, mat = Matrix(2,2,[1,2,2,5]),prec=8);
    print s2
    K2 = s2.series.base_ring() #Number field
    P2 = K2.prime_above(2)
    print('Denominators Cusp 1/2')
    print [valuation(s2[i],P2)/P2.ramification_index() for i in range(10)]
    print 'Minimal valuation', min([valuation(s2[i],P2)/P2.ramification_index() for i in range(10)])
    valuations_12.append(min([valuation(s2[i],P2)/P2.ramification_index() for i in range(10)]))
    count+=1
print 'All valuations at 1/2'
print valuations_12<jupyter_output>Calculating Fourier expansions for e1
Expansion at cusp 1/2:
1/10*zeta80^20*q10 + 1/10*zeta80^20*q10^5 + 2/5*zeta80^20*q10^7 - 3/10*zeta80^20*q10^9 - 2/5*zeta80^20*q10^11 - 1/5*zeta80^20*q10^13 + 1/5*zeta80^20*q10^17 - 2/5*zeta80^20*q10^19 - 2/5*zeta80^20*q10^23 + 1/10*zeta80^20*q10^25 - 1/5*zeta80^20*q10^29 + 4/5*zeta80^20*q10^31 + 2/5*zeta80^20*q10^35 + 3/5*zeta80^20*q10^37 - 3/5*zeta80^20*q10^41 + 4/5*zeta80^20*q10^43 - 3/10*zeta80^20*q10^45 - 2/5*zeta80^20*q10^47 + 9/10*zeta80^20*q10^49 + 3/5*zeta80^20*q10^53 - 2/5*zeta80^20*q10^55 + 2/5*zeta80^20*q10^59 - 1/5*zeta80^20*q10^61 - 6/5*zeta80^20*q10^63 - 1/5*zeta80^20*q10^65 - 4/5*zeta80^20*q10^67 - 3/5*zeta80^20*q10^73 - 8/5*zeta80^20*q10^77 + O(q10^80)

Denominators Cusp 1/2
[+Infinity, -1, +Infinity, +Infinity, +Infinity, -1, +Infinity, 1, +Infinity, -1]
Minimal valuation -1
Calculating Fourier expansions for e2
Expansion at cusp 1/2:
-1/20*zeta80^20*q10 - 1/10*zeta80^20*q10^3 + 1/20*zeta80^20*q10^5 + 1/10*zeta80^20*q10^7 - 1/20*ze[...]<jupyter_text>We check that $e_2+e_3$ is in $H^0(X_0(40),\Omega)$<jupyter_code>s0 = get_expansion(e2+e3, mat = Matrix(2,2,[0,1,-1,0])); print s0
K0 = s0.series.base_ring() #Number field
P0 = K0.prime_above(2)
print 'Denominator at 0'
print min([valuation(s0[i],P0)/P0.ramification_index() for i in range(10)])
s2 = get_expansion(e2+e3, mat = Matrix(2,2,[1,0,2,1])); print s2
K2 = s2.series.base_ring() #Number field
P2 = K2.prime_above(2)
print 'Denominator at 1/2'
print min([valuation(s2[i],P2)/P2.ramification_index() for i in range(10)])
s4 = get_expansion(e2+e3, mat = Matrix(2,2,[1,0,4,1])); print s4
K4 = s4.series.base_ring() #Number field
P4 = K4.prime_above(2)
print 'Denominator at 1/4'
print min([valuation(s4[i],P4)/P4.ramification_index() for i in range(10)])
s8 = get_expansion(e2+e3, mat = Matrix(2,2,[5,3,8,5])); print s8
K8 = s8.series.base_ring() #Number field
P8 = K8.prime_above(2)
print 'Denominator at 1/8'
print min([valuation(s8[i],P8)/P8.ramification_index() for i in range(10)])<jupyter_output>-1/40*q40 + 1/40*q40^2 + 1/40*q40^3 - 1/20*q40^6 + 1/40*q40^7 + 1/40*q40^9 - 1/40*q40^10 - 1/20*q40^11 + 1/20*q40^14 - 1/40*q40^15 + 1/20*q40^17 + 1/40*q40^18 + 1/20*q40^21 - 1/8*q40^23 - 1/40*q40^25 + 1/20*q40^26 - 1/20*q40^27 - 1/20*q40^29 + 1/20*q40^30 + 3/20*q40^31 - 3/20*q40^34 + 3/40*q40^35 - 1/10*q40^37 - 1/10*q40^38 + 1/20*q40^39 - 1/10*q40^42 + 9/40*q40^43 + 1/20*q40^45 + 3/20*q40^46 + 1/40*q40^47 - 3/40*q40^49 + 1/40*q40^50 - 3/20*q40^51 + 1/10*q40^54 - 1/20*q40^55 - 1/10*q40^57 + 3/20*q40^58 - 1/10*q40^59 - 1/10*q40^62 - 7/40*q40^63 + 1/20*q40^65 - 1/8*q40^67 + 3/20*q40^69 - 1/20*q40^70 + 3/20*q40^71 + 1/20*q40^73 + 1/20*q40^74 + 1/40*q40^75 + 1/5*q40^77 - 1/10*q40^78 - 1/10*q40^79 + O(q40^80)

Denominator at 0
-3
-1/10*zeta80^8*q10^2 + 1/10*zeta80^12*q10^3 + 1/10*zeta80^20*q10^5 + 1/5*zeta80^24*q10^6 - 3/10*zeta80^28*q10^7 + (-1/5*zeta80^28 + 1/5*zeta80^20 - 1/5*zeta80^12 + 1/5*zeta80^4)*q10^9 - 1/10*q10^10 - 1/5*zeta80^4*q10^11 + 1/5*zeta80^12*q10^13 + 1/5*zeta80^16*q10^14[...]<jupyter_text>The case N = 80We repeat the same calculations as before<jupyter_code>H_80 = hecke_basis_matrix(80); 
print 'H_80 = '
print H_80
f1,f2,f3,f4,f5,f6,f7 = CuspForms(80, prec = 10).basis();
print 'f1 = ', f1, ', f2 = ', f2, ', f3 = ', f3
print 'f4 = ', f4, ', f5 = ', f5, ', f6 = ', f6, ', f7 = ', f7
M_80 = cusp_form_basis_matrix(80); 
print 'M_80 = '
print M_80
B_80 = H_80.solve_left(M_80); 
print 'B_80 = '
print B_80<jupyter_output>H_80 = 
[ 0  1  0 -2  0 -1  0  2  0  1  0  0  0  2  0  2  0 -6  0 -4  0 -4  0  6  0]
[ 0  0  2  0  0  0 -4  0  0  0 -2  0  0  0  4  0  0  0  2  0  0  0  0  0  0]
[ 0  0  0  0  4  0  0  0  0  0  0  0 -8  0  0  0  0  0  0  0 -4  0  0  0  0]
[ 0  1  0  0  0  1  0 -4  0 -3  0  4  0 -2  0  0  0  2  0  4  0  0  0  4  0]
[ 0  0  2  0  0  0  0  0  0  0  2  0  0  0 -8  0  0  0 -6  0  0  0  8  0  0]
[ 0  1  0  0  0  1  0  4  0 -3  0 -4  0 -2  0  0  0  2  0 -4  0  0  0 -4  0]
[ 0  1  0  2  0 -1  0 -2  0  1  0  0  0  2  0 -2  0 -6  0  4  0 -4  0 -6  0]
f1 =  q - q^9 + O(q^10) , f2 =  q^2 + O(q^10) , f3 =  q^3 + O(q^10)
f4 =  q^4 + O(q^10) , f5 =  q^5 - 2*q^9 + O(q^10) , f6 =  q^6 + O(q^10) , f7 =  q^7 + O(q^10)
M_80 = 
[ 0  1  0  0  0  0  0  0  0 -1  0  0  0  0  0  0  0 -2  0  0  0 -2  0  0  0]
[ 0  0  1  0  0  0  0  0  0  0  1  0  0  0 -4  0  0  0 -3  0  0  0  4  0  0]
[ 0  0  0  1  0  0  0  0  0  0  0 -1  0  0  0 -1  0  0  0  1  0  0  0 -4  0]
[ 0  0  0  0  1  0  0  0  0  0  0  0 -2  0  0  0  0 [...]<jupyter_text>None of the forms $f_i$ are in $H^0(X_0(80),\Omega)$. But $2f_1,2f_2,8f_3,4f_4,2f_5,4f_6,8f_7$ are. As can be seen from the calculation below (the last few lines after 'SUMMARY' suffice)<jupyter_code>count = 1
valuations_0 = []
valuations_12 = []
valuations_14 = []
valuations_18 = []
for f in [f1,f2,f3,f4,f5,f6,f7]:
    print 'Calculating Fourier expansions for f{}'.format(count)
    print 'Expansion at cusp 0:'
    s0 = get_expansion(f, mat = Matrix(2,2,[0,-1,1,0]));
    print s0
    K0 = s0.series.base_ring() #Number field
    P0 = K0.prime_above(2)
    print('Denominators Cusp 0')
    print [valuation(s0[i],P0)/P0.ramification_index() for i in range(10)]
    print 'Minimal valuation', min([valuation(s0[i],P0)/P0.ramification_index() for i in range(10)])
    valuations_0.append(min([valuation(s0[i],P0)/P0.ramification_index() for i in range(10)]))
    print 'Expansion at cusp 1/2:'
    s2 = get_expansion(f, mat = Matrix(2,2,[1,2,2,5]));
    print s2
    K2 = s2.series.base_ring() #Number field
    P2 = K2.prime_above(2)
    print('Denominators Cusp 1/2')
    print [valuation(s2[i],P2)/P2.ramification_index() for i in range(10)]
    print 'Minimal valuation', min([valuation(s2[i],P2)/P2.ramification_index() for i in range(10)])
    valuations_12.append(min([valuation(s2[i],P2)/P2.ramification_index() for i in range(10)]))
    s4 = get_expansion(f, mat = Matrix(2,2,[1,1,4,5]));
    print 'Expansion at cusp 1/4:'
    print s4
    K4 = s4.series.base_ring() #Number field
    P4 = K4.prime_above(2)
    print('Denominators Cusp 1/4')
    print [valuation(s4[i],P4)/P4.ramification_index() for i in range(10)]
    print 'Minimal valuation', min([valuation(s4[i],P4)/P4.ramification_index() for i in range(10)])
    valuations_14.append(min([valuation(s4[i],P4)/P4.ramification_index() for i in range(10)]))
    s8 = get_expansion(f, mat = Matrix(2,2,[1,0,8,1]));
    print s8
    K8 = s8.series.base_ring() #Number field
    P8 = K8.prime_above(2)
    print('Denominators Cusp 1/8')
    print [valuation(s8[i],P8)/P8.ramification_index() for i in range(10)]
    print 'Minimal valuation', min([valuation(s8[i],P8)/P8.ramification_index() for i in range(10)])
    valuations_18.append(min([valuation(s8[i],P8)/P8.ramification_index() for i in range(10)]))
    count+=1
print 'SUMMARY:'
print 'All valuations at 0'
print valuations_0
print 'All valuations at 1/2'
print valuations_12
print 'All valuations at 1/4'
print valuations_14
print 'All valuations at 1/8'
print valuations_18<jupyter_output>Calculating Fourier expansions for f1
Expansion at cusp 0:
-1/160*q80 - 1/160*q80^2 - 1/160*q80^3 - 1/80*q80^4 - 1/160*q80^7 + 1/160*q80^9 - 1/160*q80^10 + 1/80*q80^11 + 1/40*q80^12 + 1/40*q80^14 + 1/160*q80^15 + 1/80*q80^17 + 3/160*q80^18 + 1/80*q80^20 + 1/80*q80^21 - 1/40*q80^22 + 1/32*q80^23 - 1/160*q80^25 + 1/80*q80^26 + 1/80*q80^27 - 1/40*q80^28 - 1/80*q80^29 - 3/80*q80^31 - 1/80*q80^34 - 3/160*q80^35 - 1/80*q80^36 - 1/40*q80^37 - 1/40*q80^38 - 1/80*q80^39 - 9/160*q80^43 + 1/80*q80^45 - 1/40*q80^46 - 1/160*q80^47 - 3/160*q80^49 - 1/160*q80^50 + 3/80*q80^51 - 1/40*q80^52 + 1/80*q80^55 - 1/40*q80^57 + 1/80*q80^58 + 1/40*q80^59 - 1/40*q80^60 + 1/20*q80^62 + 7/160*q80^63 + 1/80*q80^65 + 1/32*q80^67 + 3/40*q80^68 + 3/80*q80^69 + 1/40*q80^70 - 3/80*q80^71 + 1/80*q80^73 - 3/80*q80^74 - 1/160*q80^75 + 1/20*q80^76 + 1/20*q80^77 + 1/40*q80^79 + 1/160*q80^81 + 3/80*q80^82 - 1/32*q80^83 + 1/20*q80^84 - 1/40*q80^85 + 1/20*q80^86 - 3/80*q80^87 + 3/80*q80^89 + 3/160*q80^90 + 3/80*q80^91 - 3/40*q[...]<jupyter_text>We show that $f_1 + f_2 + f_5$ is in $H^0(X_0(80),\Omega)$. The same calculation can be done for $f = 2(f_3+f_7)$ to show that this is also in $H^0(X_0(80),\Omega)$.<jupyter_code>f = f1+f2+f5
print 'Expansion at cusp 0:'
s0 = get_expansion(f, mat = Matrix(2,2,[0,-1,1,0]));
print s0
K0 = s0.series.base_ring() #Number field
P0 = K0.prime_above(2)
print('Denominators Cusp 0')
print [valuation(s0[i],P0)/P0.ramification_index() for i in range(10)]
print 'Minimal valuation', min([valuation(s0[i],P0)/P0.ramification_index() for i in range(10)])
valuations_0.append(min([valuation(s0[i],P0)/P0.ramification_index() for i in range(10)]))
s2 = get_expansion(f, mat = Matrix(2,2,[1,2,2,5]));
print s2
K2 = s2.series.base_ring() #Number field
P2 = K2.prime_above(2)
print('Denominators Cusp 1/2')
print [valuation(s2[i],P2)/P2.ramification_index() for i in range(10)]
print 'Minimal valuation', min([valuation(s2[i],P2)/P2.ramification_index() for i in range(10)])
s4 = get_expansion(f, mat = Matrix(2,2,[1,1,4,5]));
print 'Expansion at cusp 1/4:'
print s4
K4 = s4.series.base_ring() #Number field
P4 = K4.prime_above(2)
print('Denominators Cusp 1/4')
print [valuation(s4[i],P4)/P4.ramification_index() for i in range(10)]
print 'Minimal valuation', min([valuation(s4[i],P4)/P4.ramification_index() for i in range(10)])
s8 = get_expansion(f, mat = Matrix(2,2,[1,0,8,1]));
print 'Expansion at cusp 1/8:'
print s8
K8 = s8.series.base_ring() #Number field
P8 = K8.prime_above(2)
print('Denominators Cusp 1/8')
print [valuation(s8[i],P8)/P8.ramification_index() for i in range(10)]
print 'Minimal valuation', min([valuation(s8[i],P8)/P8.ramification_index() for i in range(10)])
s16 = get_expansion(f, mat = Matrix(2,2,[1,0,16,1]));
print 'Expansion at cusp 1/16:'
print s16
K16 = s16.series.base_ring() #Number field
P16 = K16.prime_above(2)
print('Denominators Cusp 1/16')
print [valuation(s16[i],P16)/P16.ramification_index() for i in range(10)]
print 'Minimal valuation', min([valuation(s16[i],P16)/P16.ramification_index() for i in range(10)])<jupyter_output>Expansion at cusp 0:
-1/80*q80 - 1/80*q80^2 - 1/80*q80^5 + 3/80*q80^9 - 1/80*q80^10 + 1/40*q80^13 + 1/20*q80^14 - 1/40*q80^17 + 3/80*q80^18 - 1/20*q80^22 - 1/80*q80^25 + 1/40*q80^26 + 1/40*q80^29 - 1/40*q80^34 - 3/40*q80^37 - 1/20*q80^38 + 3/40*q80^41 + 3/80*q80^45 - 1/20*q80^46 - 9/80*q80^49 - 1/80*q80^50 - 3/40*q80^53 + 1/40*q80^58 + 1/40*q80^61 + 1/10*q80^62 + 1/40*q80^65 + 1/20*q80^70 + 3/40*q80^73 - 3/40*q80^74 + 1/5*q80^77 - 9/80*q80^81 + 3/40*q80^82 - 1/40*q80^85 + 1/10*q80^86 + 3/40*q80^89 + 3/80*q80^90 - 1/20*q80^94 + 7/40*q80^97 - 9/80*q80^98 - 3/40*q80^101 - 3/40*q80^106 - 7/40*q80^109 - 1/20*q80^110 - 9/40*q80^113 - 3/40*q80^117 + 1/20*q80^118 - 1/16*q80^121 + 1/40*q80^122 - 1/80*q80^125 - 3/20*q80^126 + 1/40*q80^130 + 1/5*q80^133 - 1/10*q80^134 - 1/8*q80^137 + 1/40*q80^145 + 3/40*q80^146 + 1/8*q80^149 + 3/40*q80^153 + 1/5*q80^154 + 1/40*q80^157 + O(q80^160)

Denominators Cusp 0
[+Infinity, -4, -4, +Infinity, +Infinity, -4, +Infinity, +Infinity, +Infinity, -4]
Minimal valua[...]<jupyter_text>Checking that $S_2(\Gamma_0(80),\mathbb{Z})$ is not equal to $H^0(X_0(80,\Omega)$Here we calculate the index of $L'$, the lattice generated by Hecke eigenforms in $L$, the lattice $S_2(\Gamma_0(80),\mathbb{Z})$. It is 2048. We also calculate the index of the lattice $$ which is contained in $H^0(X_0(80,\Omega)$ by previous calculations. This index is $1024$, showing that $L'$ is strictly smaller than $H^0(X_0(80,\Omega)$.<jupyter_code>G = identity_matrix(7)
basis = map(lambda x: list(x),(B_80^(-1)).rows())
hecke_lattice = IntegralLattice(G, basis)
hecke_lattice
sqrt(-hecke_lattice.discriminant())
sublattice1_basis= [[2,0,0,0,0,0,0],[0,2,0,0,0,0,0],[0,0,8,0,0,0,0],[0,0,0,4,0,0,0],[0,0,0,0,2,0,0],[0,0,0,0,0,4,0],[0,0,0,0,0,0,8]]
sublattice1 = IntegralLattice(G, sublattice1_basis) #A small lattice that is contained in H^0(X_0(80),Omega)
sublattice2_basis = [[2,0,0,0,0,0,0],[1,1,0,0,1,0,0],[0,0,8,0,0,0,0],[0,0,0,4,0,0,0],[0,0,0,0,2,0,0],[0,0,0,0,0,4,0],[0,0,2,0,0,0,2]]
sublattice2 = IntegralLattice(G, sublattice2_basis) #A slightly bigger lattice
sqrt(-sublattice_of_K2.discriminant())<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18968417118295647, 'fraction_numerical': 0.17735002803214353, 'mean_word_length': 3.8945121951219512, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 32, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '493310', 'license': 'MIT', 'n_tokens_mistral': 8891, 'n_tokens_neox': 7582, 'n_words': 1615, 'path': 'manin_constant/lattices_of_modular_forms.ipynb', 'repo_name': 'michaelneururer/notebooks-research', 'size': '16082'}