// Magma verification script for "Modular curves of prime power level with infinitely many rational points",
// by Andrew V. Sutherland and David Zywina.

// This source file verifies the various computational assertions made in the paper.

// Compute hauptmoduls of various congruence subgroups; this data is in an associative array named "CPlist".
print " "; "Computing hauptmoduls (this could take several minutes).";
load "ComputeHauptmoduls.m";


function  IsFixedByMatrix(Gamma,psi,A :prec:=3)  
/*
        Input: A matrix A in GL(2,Z/NZ), a rational function psi in K_N(t) of degree 1, and a 
        congruence subgroup Gamma (as produced by ComputeHauptmoduls.txt) of genus 0. 

        The quantity prec is used to indicated how many terms of the hauptmodul to compute; the function
        will call itself with a larger value of prec if needed.

        Output: Return true when the modular function psi(h) is fixed by the left action of A; false otherwise.              
*/
     
    N:=#BaseRing(Parent(A));
    K:=CyclotomicField(N);
    F<t>:=FunctionField(K);
    R<q>:=PuiseuxSeriesRing(K);

    d:=Integers()!Determinant(A);
    num:=Numerator(psi); den:=Denominator(psi);
    psi_d:=( Conjugate(Coefficient(num,1),d) *t + Conjugate(Coefficient(num,0),d) );
    psi_d:=psi_d/( Conjugate(Coefficient(den,1),d) *t + Conjugate(Coefficient(den,0),d) );             

    h :=Gamma`hauptmodul;
    hA:=ActOnSiegel(h, Transpose(A));

    if #Gamma`cusps ge 3 then
       for i in [1..Minimum(#Gamma`cusps,3)] do

           b:=LeadingTermSiegel(ActOnSiegel(h ,Gamma`cusps[i]));   C1:=b[1]; v1:=b[2];
           b:=LeadingTermSiegel(ActOnSiegel(hA,Gamma`cusps[i]));   C2:=b[1]; v2:=b[2]; 
           if v1 gt 0 then P1:=[0,1]; elif v1 lt 0 then P1:=[1,0]; else P1:=[C1,1]; end if;
           if v2 gt 0 then P2:=[0,1]; elif v2 lt 0 then P2:=[1,0]; else P2:=[C2,1]; end if;
           // P1 and P2 in P^1(K) are the values of h and hA, respectively, at the i-th cusp

           if EvaluateRationalFunction(psi,P1) ne EvaluateRationalFunction(psi_d,P2) then
              return false;
           end if;    
       end for;
       return true;
    end if;

  
    f:=Evaluate(F!psi,R!SiegelExpansion(h,prec));
    fA:=R!Evaluate(psi_d,R!SiegelExpansion(hA, prec));  // transpose needed.
    if AbsolutePrecision(fA-f) le 2*Gamma`widths[1]/Gamma`N then
       return IsFixedByMatrix(Gamma,psi,A:prec:=prec+1);  // try again by computing more terms of the q-expansion
    end if;
    if Valuation(fA-f) gt 2*Gamma`widths[1]/Gamma`N then 
       return true; 
    end if;                              
           
    return false;
end function;


print " "; "Verifying that the pairs (G,J) from the paper work, i.e., X_G is isomorphic to P^1  and the morphism to the j-line is given by the rational function J(t) in Q(t).   (This could take a while).";


levels:=Sort([N: N in {gl2tab[k]`gl2level : k in Keys(gl2tab)} | N gt 1]);
// We check by increasing level N.

for N in levels do
for k in Sort([k: k in Keys(gl2tab) | gl2tab[k]`gl2level eq N]) do

    r:=gl2tab[k];  
    print r`label;
    
    // The record r gives a group G and a rational function J.
    GL2:=GL(2,Integers(N));
    SL2:=SL(2,Integers(N)); 
    G:=sub<GL2|r`gens>;

    K:=CyclotomicField(N);
    F<t>:=FunctionField(K);
    J :=F!eval(r`jmap);
 
    /* We let Gamma be the genus 0 congruence subgroup associated to G; it consists of A in SL(2,Z) 
       for which A^t modulo N lies in G. We allow ourselves to replace G by a conjugate to accomplish this.
       Associated to Gamma we have a precomputed hauptmodul h and a rational function J_ in K(t) such that J_(h)=j.
    */
    
    Gamma:=CPlist[ r`sl2label ];
    J_:=F!Gamma`J;

    SL2_:=SL(2,Integers(Gamma`N));
    red:=hom<SL2->SL2_|[SL2_!SL2.i: i in [1..#Generators(SL2)]]>;
    H:= sub<SL2_|{Transpose(A):A in Gamma`H}> @@ red;
 
    _,B:=IsConjugate(GL2,G meet SL2,H);
    G:=sub<GL2|{B^(-1)*A*B : A in Generators(G)}>;
    assert G meet SL2 eq H;

    // We now compute the set Psi of degree 1 rational functions psi(t) in K(t) for which J_(t)=J(psi(t)).
    if #Gamma`cusps ge 3 then     
       xi:=Gamma`valueatcusps[3];
       triple:=[[1,0],[0,1],xi];
       b, Psi:= SameRationalParametrization(F!Gamma`J,J,K : findall:=true, triple:=triple); 
    else    
       b, Psi:= SameRationalParametrization(F!Gamma`J,J,K : findall:=true);
    end if;
    assert b;    

    
    // We finally verify that the function field of X_G is Q(psi(h)) for some psi in Psi.

    found:=false;
    for psi in Psi do
        U,f:=UnitGroup(Integers(N));
        S:={Random({g: g in G| Determinant(g) eq f(u)}) : u in Generators(U)};
        // S and H generate G.

        if &and[IsFixedByMatrix(Gamma,psi,A): A in S] then
           found:=true;
           break psi;
        end if;
    end for;
    if not found then
       "Invalid pair (G,J)"; assert false;
    end if;


end for;
end for;

g0Verify();
    
load "g1groups.m";
g1Verify();

print "Done.";
