# This program has been on my website since July 2013, but the proof
# behind the program hasn't been written down until now (2020). In
# a joint work with Hanson Smith we give two proofs for the algorithm.
#
# The preprint is at: https://arxiv.org/abs/2004.13644
# The numbers e_i(N) and f_i(N) at the end of page 3 in the preprint are
# computed by these two programs:
#
# e_i(N) = IG(i,N)
# f_i(N) = DegreeCusp(i,N)
#
IG := proc(i,N)
if i = 2 and N = 4 then
1
else
N/igcd(i,N)
fi
end:
DegreeCusp := proc(i,N)
local d;
options remember;
d := numtheory['phi']( igcd(i,N) );
if i=0 or i=N/2 then
ceil(d/2)
else
d
fi
end:
# MinFormula(k, t) returns the piecewise-linear function v_k(t)
# from the main theorem, computed using Remark 4.3.
#
MinFormula := proc(N,t) local i;
if N < 2 then
error "wrong input"
elif N = 2 then
4 * t - 1
elif N = 3 then
9 * min(t, 1/3) - 8*t
else
add(N * Phi(igcd(i,N)) * (min(t, i/N) - 4*(i/N)*(1-i/N)*t), i=1..iquo(N-1, 2))
fi
end:
Phi := proc(d)
numtheory['phi'](d) / d
end:
# Returns the divisor of modular unit F_k in X1(N), i.e. it returns the valuation
# of F_k in X1(N) at cusp C_i for each i from 0 .. floor(N/2).
#
# The reason for the _bc notation is because before this program, I computed
# divisors in another way that was much slower, so to save time I put these
# divisors into tables (the files "cusp_divisors" and "cusp_divisors_large").
#
# But the basis of modular units (the file "gonality" lists that basis) that
# was used in the table "cusp_divisors" is different from the basis used here.
# Hence the subscript _bc to indicate that basis.
#
Divisor_F_bc := proc(N,k)
local t, i, f;
options remember;
f := MinFormula(k, t);
[seq(eval(f, t=i/N) * IG(i,N), i=0..floor(N/2))]
end:
# Note: the piecewise linear function f = v_k(t) = MinFormula(k, t)
# is {unweighted} and to get the actual valuation at cusp C_i you
# still have to multiply that by e_i(N) = IG(i,N), see the main theorem.
# Basis of divisors of modular units. The Z-SPAN is the same as that in the
# table "cusp_divisors" but the basis differs. For table "cusp_divisors", the
# basis is given in the file "gonality", whereas the basis used here is the
# F_2, F_3,... in the paper.
#
L_bc := proc(N)
local k;
[seq(Divisor_F_bc(N,k), k=2..floor(N/2)+1)]
end:
# Returns the structure of the cuspidal part of the class group:
#
CuspidalClassGroupStructure := proc(N)
local M,i;
M := LinearAlgebra['SmithForm'](Matrix(L_bc(N)));
subs(1=NULL,[seq(M[i,i], i=1..floor(N/2))])
end:
# Maarten Derickx implemented a Sage version of this program as well as numerous
# other useful programs such as the CuspidalClassGroupStructure for X1(m, m*n).