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

// This file contains various utility functions for computing invariants of subgroups of GL(Z/NZ)

// Given a subgroup H of GL(2,Z/nZ), computes the index of det(H) in (Z/nZ)*
function GL2DetIndex(H)
    M,pi:=MultiplicativeGroup(BaseRing(H));
    return Index(M,sub<M|[Inverse(pi)(Determinant(h)):h in Generators(H)]>);
end function;

// Compute the abelian group invariants of a Z/nZ-submodule of Z/nZ x Z/nZ (which Magma apparently doesn't now is an abelian group)
function GL2ModuleInvariants(V)
    if Dimension(V) eq 0 then return []; end if;
    if Dimension(V) eq 1 then return [#V]; end if;
    assert Dimension(V) le 2;
    r1:=#sub<V|V.1>;  r2:=#sub<V|V.2>;
    return [GCD(r1,r2),LCM(r1,r2)];
end function;

// Given invariants A and B of two Z/nZ-submodules of Z/nZ x Z/nZ, return true if A is isomorphic to a submodule of B
function GL2IsSubModule(A,B)
    i:=#B-#A;
    if i lt 0 then return false; end if;
    for j in [1..#A] do if not IsDivisibleBy(B[i+j],A[j]) then return false; end if; end for;
    return true;
end function;

// Given a subgroup of GL(2,Z/nZ), computes the invariants of the sub-module of Z/nZ x Z/nZ fixed by the left action of G (returns a list [], [a], or [a,b] with a|b|n)
function GL2FixModule(H)
    V := Eigenspace(Identity(H),1);
    for h in Generators(H) do V:= V meet Eigenspace(Transpose(h),1); end for;   // take transpose so that we  are using left actions (Magma default is right action!)
    return GL2ModuleInvariants(V);
end function;

// Given a subgroup H of GL(2,Z/nZ), returns true if H contains an element corresponding to complex conjugation
// This is tested by checking det = -1, tr = 0, and fixing a maximal cyclic subgroup (note: this is faster than checking for an element conjugate to [1,0,0,-1] or [1,1,0,-1])
function GL2ContainsCC(H)
    return #[h:h in H|Determinant(h) eq -1 and Trace(h) eq 0 and GL2IsSubModule([#BaseRing(H)],GL2FixModule(sub<H|h>))] gt 0;
end function;

// Given a subgroup G of GL(2,Z/nZ), return nu2(Gamma), where Gamma is the PSL(2,Z) congruence subgroup Gamma corresponding to G (which depends only on the intersection of G with SL(2,/Z/nZ) mod +/-1)
// That is, we take the inverse image of G meet SL(2,Z/nZ) in SL(2,Z)/{-1} and count the number of elliptic points of order 2, i.e. Gamma-inequivalent SL(2,Z)-translates of i that have non-trivial stabilizers in Gamma
// This amounts to counting right cosets fixed by the stablizer of i, which is generated by e2=[0,1,-1,0]
function GL2Nu2(G)
    SL2 := SL(2,BaseRing(G));
    H := sub<SL2|G meet SL2,[-1,0,0,-1]>;
    R:=RightTransversal(SL2,H);
    e2:=SL2![0,1,-1,0];
    return #{g:g in R|g*e2*g^-1 in H};
end function;
    
// Given a subgroup G of GL(2,Z/nZ), return nu3(Gamma), where Gamma is the PSL(2,Z) congruence subgroup Gamma corresponding to G (which depends only on the intersection of G with SL(2,/Z/nZ) mod +/-1)
// That is, we take the inverse image of G meet SL(2,Z/nZ) in SL(2,Z)/{-1} and count the number of elliptic points of order 3, i.e. Gamma-inequivalent SL(2,Z)-translates of rho that have non-trivial stabilizers in Gamma
// This amounts to counting right cosets fixed by the stablizer of z_3, which is generated by e2=[0,1,-1,0]
function GL2Nu3(G)
    SL2 := SL(2,BaseRing(G));
    H := sub<SL2|G meet SL2,[-1,0,0,-1]>;
    R:=RightTransversal(SL2,H);
    e3:=SL2![0,1,-1,-1];
    return #{g:g in R|g*e3*g^-1 in H};
end function;

// Given a subgroup G of GL(2,Z/nZ), return the number of cusp orbits for Gamma, where Gamma is the PSL(2,Z) congruence subgroup Gamma corresponding to G (which depends only on the intersection of G with SL(2,/Z/nZ) mod +/-1)
// That is, we take the inverse image of G meet SL(2,Z/nZ) in SL(2,Z)/{-1} and count the number of orbits of infty, equivalently, the number of orbits of the right-coset space of Gamma under the right-action of T=[1,1,0,1]
function GL2CuspCount(G)
    SL2 := SL(2,BaseRing(G));
    H := sub<SL2|G meet SL2,[-1,0,0,-1]>;
    R:=[r:r in RightTransversal(SL2,H)];
    T:=SL2![1,1,0,1];
    c:=0;
    while #R ne 0 do
        // pick a right coset rep
        r := R[#R];  R:=Prune(R); c+:=1;
        // compute its orbit under right-multiplication by T, removing each coset rep in the orbit from R as we go
        while true do
            r := r*T;
            S:=[i:i in [1..#R]|R[i]*r^-1 in H]; // figure out which coset in our list is equal to the coset represented by r
            if #S eq 0 then break; end if;     // if we can't find it, its because we already removed it (i.e we cycled back to our starting point)
            assert #S eq 1;
            R:=Remove(R,S[1]);
        end while;
    end while;
    return c;
end function;

// Given a subgroup G of GL(2,Z/nZ), return the genus of the modular curve X_G', where G'=<G,-I>
function GL2Genus(G)
    SL2:=SL(2,BaseRing(G));
    H:=sub<SL2|G meet SL2,[-1,0,0,-1]>;
    g := 1 + Index(SL2,H)/12 - GL2Nu2(H)/4 - GL2Nu3(H)/3 - GL2CuspCount(H)/2;
    return Integers()!g;
end function;
