// Magma code and data related to "Modular curves of prime power level with infinitely many rational points",
// by Andrew V. Sutherland and David Zywina.

load "RationalFunctions.m";
load "ModularFunctions.m";


CongSubgroup := recformat<N:RngIntElt, label:MonStgElt, gens:SeqEnum, H:GrpMat, index:RngIntElt,  cusps,widths:SeqEnum, hauptmodul, h:RngSerPuisElt, J:FldFunRatUElt, valueatcusps, sup :SeqEnum>;

/* Record for a congruence subgroup Gamma of genus 0.
   N:              level of Gamma
   label:          label of Gamma given by Cummin and Pauli (describes Gamma up to conjugacy in PGL_2(Z).
   H:              image of Gamma in SL_2(Z/NZ).
   gens:           generators for subgroup H of SL_2(Z/NZ).
   index:          index of Gamma in SL_2(Z).
   cusps:          sequence of matrices A in SL_2(Z) such that A*infty represent the cusps of X_Gamma; 
                   the first cusp will be the one at infinity.
   widths:         width of the cusps.

   hauptmodul:     An abstract description of a hauptmodul for Gamma.  It is in one of the following forms:
                   1) A sequence of records of type SiegelPower; represents the product of those Siegel powers.
                   2) A sequence of sequence as in (1); represents the sum of products of Siegel powers.

   h:              q-expansion for the hauptmodul.   
   J:              rational function such that J(h) is the modular j-invariant.
   valueatcusps:  value of the hauptmodul at each of the cusps
   
   sup:            in a list of congruence subgroups will give label of congruence subgroups Gamma' containing Gamma.
*/
   

            
/*
The following is a complete list of the congruence subgroups of genus 0 with prime power level, up to conjugacy in PGL_2(Z), indexed by their 
Cummin-Pauli name given.  Our specific choices are different than theirs; we have chosen so that there are lots 
of inclusions between the groups.   See http://www.uncg.edu/mat/faculty/pauli/congruence/congruence.html

WARNING: we have commented out the level 11 entry.
*/

CPlist:=AssociativeArray();     

CPlist["1A0"]:= rec<CongSubgroup | N:=1,  label:="1A0", gens:=[ ]>;
CPlist["2A0"]:= rec<CongSubgroup | N:=2,  label:="2A0", gens:=[ [1,1,1,0] ]>;
CPlist["2B0"]:= rec<CongSubgroup | N:=2,  label:="2B0", gens:=[ [1,0,1,1] ]>;
CPlist["2C0"]:= rec<CongSubgroup | N:=2,  label:="2C0", gens:=[ [1,0,0,1] ]>;
CPlist["3A0"]:= rec<CongSubgroup | N:=3,  label:="3A0", gens:=[ [2,0,0,2], [0,2,1,0], [1,2,2,2] ]>;
CPlist["3B0"]:= rec<CongSubgroup | N:=3,  label:="3B0", gens:=[ [1,0,1,1], [2,0,0,2] ]>;
CPlist["3C0"]:= rec<CongSubgroup | N:=3,  label:="3C0", gens:=[ [2,0,0,2], [0,2,1,0] ]>;
CPlist["3D0"]:= rec<CongSubgroup | N:=3,  label:="3D0", gens:=[ [2,0,0,2] ]>;
CPlist["4A0"]:= rec<CongSubgroup | N:=4,  label:="4A0", gens:=[ [3,0,0,3], [2,1,1,1], [1,2,1,3] ]>;
CPlist["4B0"]:= rec<CongSubgroup | N:=4,  label:="4B0", gens:=[ [1,0,1,1], [3,0,0,3] ]>;
CPlist["4C0"]:= rec<CongSubgroup | N:=4,  label:="4C0", gens:=[ [3,2,1,1], [1,0,2,1], [3,0,2,3] ]>;
CPlist["4D0"]:= rec<CongSubgroup | N:=4,  label:="4D0", gens:=[ [2,3,3,3], [3,1,1,2], [1,3,3,2], [3,0,0,3], [2,1,1,1] ]>;
CPlist["4E0"]:= rec<CongSubgroup | N:=4,  label:="4E0", gens:=[ [1,0,2,1], [3,0,0,3] ]>;
CPlist["4F0"]:= rec<CongSubgroup | N:=4,  label:="4F0", gens:=[ [7,0,0,7], [5,0,0,5], [3,6,7,1], [5,0,4,5], [1,6,1,7] ]>;
CPlist["4G0"]:= rec<CongSubgroup | N:=4,  label:="4G0", gens:=[ [3,0,0,3] ]>;
CPlist["5A0"]:= rec<CongSubgroup | N:=5,  label:="5A0", gens:=[ [2,0,0,3], [0,-1,1,0], [1,2,1,3] ]>;
CPlist["5B0"]:= rec<CongSubgroup | N:=5,  label:="5B0", gens:=[ [1,0,1,1], [2,0,0,3] ]>;
CPlist["5C0"]:= rec<CongSubgroup | N:=5,  label:="5C0", gens:=[ [2,0,0,3], [3,2,4,3] ]>;
CPlist["5D0"]:= rec<CongSubgroup | N:=5,  label:="5D0", gens:=[ [1,0,1,1], [4,0,0,4] ]>;
CPlist["5E0"]:= rec<CongSubgroup | N:=5,  label:="5E0", gens:=[ [2,0,0,3], [0,-1,1,0] ]>;
CPlist["5F0"]:= rec<CongSubgroup | N:=5,  label:="5F0", gens:=[ [3,2,4,3] ]>;
CPlist["5G0"]:= rec<CongSubgroup | N:=5,  label:="5G0", gens:=[ [2,0,0,3] ]>;
CPlist["5H0"]:= rec<CongSubgroup | N:=5,  label:="5H0", gens:=[ [4,0,0,4] ]>;
CPlist["7A0"]:= rec<CongSubgroup | N:=7,  label:="7A0", gens:=[ [0,2,3,0 ], [2,6,5,5 ], [4,6,0,2 ], [6,0,0,6 ], [6,2,6,1] ]>;
CPlist["7B0"]:= rec<CongSubgroup | N:=7,  label:="7B0", gens:=[ [3,0,2,5 ], [5,4,3,4 ], [6,0,0,6] ]>;
CPlist["7C0"]:= rec<CongSubgroup | N:=7,  label:="7C0", gens:=[ [6,3,4,1], [5,6,0,3], [6,0,0,6], [1,5,1,6] ]>;
CPlist["7D0"]:= rec<CongSubgroup | N:=7,  label:="7D0", gens:=[ [0,2,3,0 ], [1,4,3,6 ], [2,6,5,5 ], [6,0,0,6] ]>;
CPlist["7E0"]:= rec<CongSubgroup | N:=7,  label:="7E0", gens:=[ [6,0,0,6 ], [6,5,2,3] ]>;
CPlist["7F0"]:= rec<CongSubgroup | N:=7,  label:="7F0", gens:=[ [0,2,3,0 ], [4,5,3,4 ], [6,0,0,6] ]>;
CPlist["7G0"]:= rec<CongSubgroup | N:=7,  label:="7G0", gens:=[ [6,3,4,1], [6,0,0,6], [1,5,1,6] ]>;
CPlist["8A0"]:= rec<CongSubgroup | N:=8,  label:="8A0", gens:=[ [1,4,4,1], [7,0,0,7], [7,2,3,1], [5,2,1,7], [5,0,4,5], [5,4,0,5], [1,7,3,6] ]>;
CPlist["8B0"]:= rec<CongSubgroup | N:=8,  label:="8B0", gens:=[ [7,0,0,7], [5,0,0,5], [5,2,7,3], [5,0,4,5], [1,6,1,7] ]>;
CPlist["8C0"]:= rec<CongSubgroup | N:=8,  label:="8C0", gens:=[ [1,0,1,1], [7,0,0,7], [5,0,0,5] ]>;
CPlist["8D0"]:= rec<CongSubgroup | N:=8,  label:="8D0", gens:=[ [7,0,0,7], [5,0,0,5], [5,0,4,5], [1,6,1,7], [3,0,2,3] ]>;
CPlist["8E0"]:= rec<CongSubgroup | N:=8,  label:="8E0", gens:=[ [2,1,1,5], [7,0,0,7], [5,0,4,5], [5,4,0,5], [1,7,3,6] ]>;
CPlist["8F0"]:= rec<CongSubgroup | N:=8,  label:="8F0", gens:=[ [6,1,5,1], [7,0,0,7], [5,0,0,5], [7,7,2,1] ]>;
CPlist["8G0"]:= rec<CongSubgroup | N:=8,  label:="8G0", gens:=[ [1,0,2,1], [7,0,0,7], [5,0,0,5] ]>;
CPlist["8H0"]:= rec<CongSubgroup | N:=8,  label:="8H0", gens:=[ [7,0,0,7], [5,0,0,5], [5,0,4,5], [1,6,1,7] ]>;
CPlist["8I0"]:= rec<CongSubgroup | N:=8,  label:="8I0", gens:=[ [1,0,1,1], [7,0,0,7] ]>;
CPlist["8J0"]:= rec<CongSubgroup | N:=8,  label:="8J0", gens:=[ [7,0,2,7], [3,4,6,3], [1,0,4,1], [1,0,6,1] ]>;
CPlist["8K0"]:= rec<CongSubgroup | N:=8,  label:="8K0", gens:=[ [5,6,5,3], [7,0,0,7], [5,0,4,5], [1,4,0,1] ]>;
CPlist["8L0"]:= rec<CongSubgroup | N:=8,  label:="8L0", gens:=[ [7,4,6,7], [5,6,5,3], [7,0,0,7], [1,0,4,1] ]>;
CPlist["8M0"]:= rec<CongSubgroup | N:=8,  label:="8M0", gens:=[ [7,0,0,7], [7,2,3,1], [1,7,3,6] ]>;
CPlist["8N0"]:= rec<CongSubgroup | N:=8,  label:="8N0", gens:=[ [5,4,4,5], [7,0,0,7], [7,0,4,7], [3,4,4,3], [3,4,0,3], [1,0,4,1], [5,4,0,5] ]>;
CPlist["8O0"]:= rec<CongSubgroup | N:=8,  label:="8O0", gens:=[ [1,0,2,1], [7,0,0,7] ]>;
CPlist["8P0"]:= rec<CongSubgroup | N:=8,  label:="8P0", gens:=[ [5,6,5,3], [7,0,0,7], [1,0,4,1] ]>;
CPlist["9A0"]:= rec<CongSubgroup | N:=9,  label:="9A0", gens:=[ [0,8,1,0 ], [1,1,7,8 ], [2,6,6,5 ], [7,0,0,4 ], [7,6,6,4] ]>;
CPlist["9B0"]:= rec<CongSubgroup | N:=9,  label:="9B0", gens:=[ [1,0,1,1], [2,0,0,5] ]>;
CPlist["9C0"]:= rec<CongSubgroup | N:=9,  label:="9C0", gens:=[ [4,0,3,7], [7,0,3,4], [8,0,0,8], [7,6,7,1] ]>;
CPlist["9D0"]:= rec<CongSubgroup | N:=9,  label:="9D0", gens:=[ [0,8,1,0 ], [4,6,6,7 ], [5,3,3,2 ], [7,0,0,4] ]>;
CPlist["9E0"]:= rec<CongSubgroup | N:=9,  label:="9E0", gens:=[ [0,8,1,0 ], [4,3,6,7 ], [7,0,0,4 ], [8,0,0,8] ]>;
CPlist["9F0"]:= rec<CongSubgroup | N:=9,  label:="9F0", gens:=[ [0,8,1,0], [0,1,8,0], [1,2,3,7], [4,1,1,5], [8,0,0,8] ]>;
CPlist["9G0"]:= rec<CongSubgroup | N:=9,  label:="9G0", gens:=[ [0,8,1,0 ], [1,1,7,8 ], [4,6,6,7 ], [5,3,3,2] ]>;
CPlist["9H0"]:= rec<CongSubgroup | N:=9,  label:="9H0", gens:=[ [4,0,0,7 ], [4,6,3,7 ], [8,0,0,8] ]>;
CPlist["9I0"]:= rec<CongSubgroup | N:=9,  label:="9I0", gens:=[ [1,0,1,1], [8,0,0,8] ]>;
CPlist["9J0"]:= rec<CongSubgroup | N:=9,  label:="9J0", gens:=[ [1,0,3,1], [8,0,0,8], [7,6,7,1] ]>;
//CPlist["11A0"]:= rec<CongSubgroup | N:=11,  label:="11A0", gens:=[ [8,6,2,3 ], [9,2,4,1 ], [10,0,0,10] ]>;
CPlist["13A0"]:= rec<CongSubgroup | N:=13,  label:="13A0", gens:=[ [1,11,3,8 ], [7,4,7,6 ], [9,7,2,6 ], [12,0,0,12] ]>;
CPlist["13B0"]:= rec<CongSubgroup | N:=13,  label:="13B0", gens:=[ [0,4,3,2 ], [4,11,0,10 ], [5,3,7,7] ]>;
CPlist["13C0"]:= rec<CongSubgroup | N:=13,  label:="13C0", gens:=[ [0,4,3,2 ], [6,9,6,7 ], [7,7,2,4] ]>;
CPlist["16A0"]:= rec<CongSubgroup | N:=16,  label:="16A0", gens:=[ [11,4,0,3], [13,4,0,5], [0,11,13,12], [13,8,12,5], [9,8,0,9], [9,0,8,9], [1,11,7,14] ]>;
CPlist["16B0"]:= rec<CongSubgroup | N:=16,  label:="16B0", gens:=[ [11,4,2,11], [9,0,0,9], [13,14,13,3], [15,0,0,15], [5,0,4,13], [1,8,4,1], [1,0,8,1] ]>;
CPlist["16C0"]:= rec<CongSubgroup | N:=16,  label:="16C0", gens:=[ [1,0,1,1], [15,0,0,15], [5,0,0,13] ]>;
CPlist["16D0"]:= rec<CongSubgroup | N:=16,  label:="16D0", gens:=[ [1,0,1,1], [17,0,0,17], [9,0,0,9], [11,8,1,11], [-1,0,0,-1] ]>;
CPlist["16E0"]:= rec<CongSubgroup | N:=16,  label:="16E0", gens:=[ [1,0,12,1], [9,0,14,9], [9,0,0,9], [5,6,13,3], [15,0,0,15], [9,0,8,9], [5,8,8,13] ]>;
CPlist["16F0"]:= rec<CongSubgroup | N:=16,  label:="16F0", gens:=[ [1,11,7,14], [11,4,0,3], [13,4,0,5] ]>;
CPlist["16G0"]:= rec<CongSubgroup | N:=16,  label:="16G0", gens:=[ [9,0,14,9], [9,0,0,9], [11,8,8,3], [-1,0,0,-1] ]>;
CPlist["16H0"]:= rec<CongSubgroup | N:=16,  label:="16H0", gens:=[ [1,0,1,1], [17,0,0,17], [9,0,0,9], [-1,0,0,-1] ]>;
CPlist["25A0"]:= rec<CongSubgroup | N:=25,  label:="25A0", gens:=[ [1,0,1,1], [2,0,0,13] ]>;
CPlist["25B0"]:= rec<CongSubgroup | N:=25,  label:="25B0", gens:=[ [1,0,1,1], [4,0,0,19] ]>;
CPlist["27A0"]:= rec<CongSubgroup | N:=27,  label:="27A0", gens:=[ [26,0,0,26], [7,18,26,13], [13,9,2,16], [19,0,0,10], [19,0,3,10], [1,0,9,1] ]>;
CPlist["32A0"]:= rec<CongSubgroup | N:=32,  label:="32A0", gens:=[ [1,0,26,1], [1,0,8,1], [15,0,10,15], [21,16,10,29], [17,0,24,17], [25,0,26,9], [17,0,29,17], [17,0,16,17], [1,0,20,1], [21,16,28,29], [1,0,16,1], [25,0,4,9 ] ]>;


// Computes the group H, the index of Gamma in SL(2,Z) and the cusp data of X_Gamma for each congruence subgroup Gamma.

for k in Keys(CPlist) do
	Gamma:=CPlist[k];
	if k eq "1A0" then
		Gamma`index:=1;  Gamma`cusps:=[[1,0,0,1]];  Gamma`widths:=[1];
	else
		N:=Gamma`N;
		SL2:=SL(2,Integers(N));
		Gamma`H:=sub< SL2 | Gamma`gens cat [[-1,0,0,-1]]>;
		Gamma`index:= Index(SL2,Gamma`H);
     		Gamma`cusps, Gamma`widths := CuspData(Gamma`H);
	end if;
	CPlist[k]:=Gamma;
end for;


// For a congruence subgroup Gamma, determines the congruence subgroups Gamma' in the list that contain, and do not equal, Gamma.   
// We then order them by decreasing index in SL(2,Z).

for k  in Keys(CPlist) do
    Gamma:=CPlist[k]; 
    if Gamma`label eq "1A0" then
       CPlist[k]`sup:=[]; continue k;
    end if; 
    N:=Gamma`N; H:=Gamma`H; 
    SL2:=SL(2,Integers(N));   

    Gamma`sup:=["1A0"];
    for k_ in Keys(CPlist) do
        if k_ eq "1A0" then continue k_; end if;
        if k_ eq k then continue k_; end if;
 
        Gamma_:=CPlist[k_];
        N_:=Gamma_`N; H_:=Gamma_`H;

        if N mod N_ eq 0 then
           SL2_:=SL(2,Integers(N_));
           red:=hom<SL2->SL2_|[SL2_!SL2.i: i in [1..#Generators(SL2)]]>;
           if H subset (H_ @@ red) then
              Gamma`sup:=Gamma`sup cat [Gamma_`label];
           end if;
        end if;
    end for;
    CPlist[k]:=Gamma;
end for;

// In our list of congruence subgroups, we determine which are subgroups of each other.

for k in Keys(CPlist) do
    sup:=CPlist[k]`sup; 
    ind:=[CPlist[k_]`index: k_ in sup];
    ParallelSort(~ind, ~sup); 
    Reverse(~sup);
    CPlist[k]`sup:=sup;
end for;



// Compute a hauptmodul in terms of Siegel functions for each congruence subgroup in our list (except SL(2,Z)).

for k in Keys(CPlist) diff {"1A0"} do    
    CPlist[k]`hauptmodul:=FindHauptmodul(CPlist[k]`H);
end for;


// Compute the values of the hauptmodul of X_Gamma at each of its cusps when X_Gamma has at least 2 cusps.
// (In this case, our constructed hauptmodul is expressed as a product of Siegel functions)

for k in Keys(CPlist) do
    Gamma:=CPlist[k];  
    
    if #Gamma`cusps ge 2 then
       N:=Gamma`N;  K:=CyclotomicField(N); 
       Gamma`valueatcusps:=[];
       for A in Gamma`cusps do 
           w:=LeadingTermSiegel( ActOnSiegel(Gamma`hauptmodul,A) );
           if w[2] lt 0 then
              Gamma`valueatcusps:=Gamma`valueatcusps cat [[1,0]];
           elif w[2] gt 0 then
              Gamma`valueatcusps:=Gamma`valueatcusps cat [[0,1]];
           else
              Gamma`valueatcusps:=Gamma`valueatcusps cat [[K!w[1],1]];
           end if;
       end for;
    else
       Gamma`valueatcusps:=[[1,0]];
    end if;

    CPlist[k]:=Gamma;
end for;


// Compute q-expansion of the hauptmoduls; we compute enough terms for our computations.

R<q>:=PuiseuxSeriesRing(Rationals());
CPlist["1A0"]`h:=jInvariant(q);

for k in Keys(CPlist) diff {"1A0"} do  
    Gamma:=CPlist[k];  N:=Gamma`N;  
    R<q>:=PuiseuxSeriesRing(CyclotomicField(N));
    Gamma`h:=R!SiegelExpansion(Gamma`hauptmodul,7);  
    CPlist[k]:=Gamma;
end for;


// For the computed hauptmodul h of X_Gamma, computes the rational function J
// such that J(h) is the modular j-invariant. 

keys:=[k: k in Keys(CPlist)];
ind:=[CPlist[k]`index:k in keys];
ParallelSort(~ind, ~keys);

for k in keys do
    Gamma:=CPlist[k];  N:=Gamma`N;
    F<t>:=FunctionField(CyclotomicField(N));
    if k eq "1A0" then
       CPlist[k]`J:=t; continue k; 
    end if;
    
    Gamma_:=CPlist[ Gamma`sup[1] ];
    f:=F!FindRelation(Gamma`h,Gamma_`h, Gamma`index div Gamma_`index);

    CPlist[k]`J:=Evaluate(F!Gamma_`J, f);
end for;



