# 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).