 /*
    Magma code related to the paper :
        "Computing images of Galois representations attached to elliptic curves", Forum of Mathematics, Sigma 4 (2016) e4 (79 pages), http://dx.doi.org/10.1017/fms.2015.33
    Copyright (c) 2015 by Andrew V. Sutherland
    You are free to use and/or modify this code under GPL version 2 or any later version (but please be sure to cite the paper)
    
    This file contains functions for identifying, constructing, and enumerating subgroups of GL2(Z/pZ).
    It includes an implementation of Algorithms 1 and 2 from the paper, counting formulas derived in Section 3, and functions for determining the subgroup labels described in section 7.
    
    The main functions of note are
    
        GL2Counts(p) -- returns the number of subgroup conjugacy classes in GL2(Z/pZ), implements the formuals derived in Section 3 of the paper
        GL2Subgroups(p) -- returns a complete list of the subgroups of GL2(Z/pZ) up to conjugacy, using Algorithm 2 of the paper
        GL2SubgroupLabel(H) -- returns a text label that uniquely identifies the conjugacy class of H of GL2(Z/pZ), as described in section 7 of the paper
        GL2SubgroupFromLabel(s) -- given a text label s identifying a subgroup of GL2(Z/pZ) constructs a subgroup with generators as decribed in Section 7 of the paper
        GL2LocalConjugate(H) -- given a subgroup of GL2(Z/pZ) or its label, returns the label of the unique non-conjugate locally conjugate group, if any, or the empty string
        GL2Twists(H) -- given a subgroup of GL2(Z/pZ) or its label, returns a set containing the label of G:=<H,-1> and the index 2 subgroups of G that do not contain -1 (this set contains the labels of all subgroups that arise for twists of E with G_E(p)=H (including H itself) and will have cardinality at most 3, see Section 5.4 of the paper)
*/

// code to count subgroups of GL_2(Z/pZ) up to conjugacy
alpha:=function(n) f:=Factorization(n);  return &*[&+[(2*e[2]-2*i+1)*e[1]^i: i in [0..e[2]]]:e in f]; end function;
beta:=function(n) f:=Factorization(n); e2:=Valuation(n,2); if e2 eq 0 then b2:=1; else b2:=(2*e2^2-2*e2+3); end if; return b2*&*([(e[2]+1)^2:e in f|e[1] ne 2] cat [1]); end function;
tau:=function(n) if Denominator(n) ne 1 then return 0; end if; return #Divisors(Integers()!n); end function;
    
GL2Counts :=function(p)
    assert IsPrime(p);
    if p eq 2 then return 4; end if;
    Sn := tau(p-1);
    Bn := alpha(p-1);
    Zn := tau(p-1);
    sCn := ExactQuotient(alpha(p-1)+beta(p-1),2);
    nsCn:=tau(p^2-1);
    sDn := beta(p-1)-tau(p-1)+tau((p-1)/2)^2-tau(ExactQuotient(p-1,2))-1/2*( 1+KroneckerSymbol(-1,p))*tau((p-1)/4);
    nsDn := tau(p^2-1)-tau(p-1)+tau((p^2-1)/4)-tau((p-1)/2);
    DDn:=tau(ExactQuotient(p-1,2^Valuation(p-1,2)));
    if p le 3 then EA4n:= 0; else EA4n := tau((p-1)/2) + 1/2*(1+KroneckerSymbol(-3,p))*tau((p-1)/6); end if;
    if p le 3 then ES4n:= 0; else ES4n := (1-1/4*(1-KroneckerSymbol(2,p))*(1-KroneckerSymbol(-2,p)))*tau((p-1)/2) + 1/2*(1+KroneckerSymbol(-1,p))*tau((p-1)/4); end if;
    if p le 5 then EA5n:= 0; else EA5n := 1/2*(1+KroneckerSymbol(5,p))*tau((p-1)/2); end if;
    tot := Sn + Bn + sCn + nsCn - Zn + sDn + nsDn - DDn+ EA4n + ES4n + EA5n;
    //print p,"&",Sn,"&",Bn,"&",sCn,"&",nsCn,"&",Zn,"&",sDn,"&",nsDn,"&",DDn,"&",EA4n,"&",ES4n,"&",EA5n,"&",tot,"\\\\";
    return tot;
end function;

dets:=function(H) return IsTrivial(H) select 1 else LCM([Order(Determinant(h)):h in Generators(H)]); end function;
chi:=function(g) return KroneckerSymbol(Integers()!(Trace(g)^2-4*Determinant(g)),Characteristic(BaseRing(g))); end function;

// return array whose nth entry is the least positive integer with order n mod p (requires n|(p-1))
MinReps :=function(p)
    assert IsPrime(p);
    r := PrimitiveRoot(p);
    a:=Integers(p)!r;
    A:=[p:i in [1..p-1]];
    for i in [1..p-1] do
        o:=ExactQuotient(p-1,GCD(p-1,i));
        A[o]:=Min(A[o],Integers()!a);
        a:=r*a;
    end for;
    return A;
end function;

Scalars:=function(p)
    return Center(GL(2,p));
end function;
    
IsScalar := function(H)
    for g in Generators(H) do if g[1][2] ne 0 or g[2][1] ne 0 or g[1][1] ne g[2][2] then return false; end if; end for;
    return true;
end function;

SplitCartan:=function(p)
    local r;
    r:=PrimitiveRoot(p);
    return sub<GL(2,p)|[1,0,0,r],[r,0,0,1]>;
end function;

Borel:=function(p)
    return sub<GL(2,p)|SplitCartan(p),[1,1,0,1]>;
end function;

NonSplitCartan:=function(p)
    if p eq 2 then return sub<GL(2,p)|[1,1,1,0]>; end if;
    r:=PrimitiveRoot(p);
    G:=GL(2,p);
    while true do
        x:=Random(GF(p));  y:=Random(GF(p));
        if y eq 0 then continue; end if;
        H:=sub<G|[x,r*y,y,x]>;
        if #H eq p^2-1 then break; end if;
    end while;
    return H;
end function;

NormalizerSplitCartan:=function(p)
    return sub<GL(2,p)|SplitCartan(p),[0,1,1,0]>;
end function;

NormalizerNonSplitCartan:=function(p)
    return sub<GL(2,p)|NonSplitCartan(p),[1,0,0,-1]>;
end function;

IsDiagonal:=function(a)
    return a[1][2] eq 0 and a[2][1] eq 0;
end function;

// given diagonalizabe H, returns conjugate of H that is diagonal with minimal integers a and b such that H~<[a,0,0,1/a],[b,0,0,r/b]> where r is a minimal generator for det(H)
SplitCartanConjugate := function(H:M:=[])
    local p,H0,a,H1,h,g,b;
    assert IsAbelian(H);
    assert {chi(h):h in Generators(H)} subset {0,1};
    p:=Characteristic(BaseRing(H));
    if p eq 2 then assert #H eq 1; return H; end if;
    H0:=H meet SL(2,p);
    if #M ne p-1 then M:=MinReps(p); end if;
    a:= M[#H0];                                         // least a such that <[a,0,0,1/a]> is conjugate to H meet SL(2,p)
    r:=M[dets(H)];                                      // least r that generates det(H)
    H1,pi:=quo<H|H0>;
    h:=Inverse(pi)([h:h in H1|Order(h) eq #H1][1]);
    g:=[g:g in sub<H|h>|Determinant(g) eq r][1];
    b:= Min([Min([e[1]:e in Eigenvalues(h0*g)]):h0 in H0]);
    return sub<GL(2,p)|[a,0,0,1/a],[b,0,0,r/b]>,Integers()!a,Integers()!b;
end function;

// given diagonal H, returns minimal integers a and b such that H=<[a,0,0,1/a], [b,0,0,r/b]>, where r is minimal generator for det(H)
DiagonalSubgroupGenerators:=function(H:M:=[])
    for h in Generators(H) do assert IsDiagonal(h); end for;
    p:=#BaseRing(H);
    H0:=H meet SL(2,p);
    if #M ne p-1 then M:=MinReps(p); end if;
    a:= M[#H0];                                         // least a such that <[a,0,0,1/a]> is conjugate to H meet SL(2,p)
    r:=M[dets(H)];                                      // least r that generates det(H)
    H1,pi:=quo<H|H0>;
    h:=Inverse(pi)([h:h in H1|Order(h) eq #H1][1]);
    g:=[g:g in sub<H|h>|Determinant(g) eq r][1];
    b := Min([(h0*g)[1][1]:h0 in H0]);
    return sub<GL(2,p)|[a,0,0,1/a],[b,0,0,r/b]>,Integers()!a,Integers()!b;
end function;

// given H contained in some Borel subgroup, returns conjugate of H in upper-triangular Borel
BorelConjugate := function(H:M:=[])
    local p,s,n,a,b,G,H1;
    p:=#BaseRing(H);
    if not IsDivisibleBy(Order(H),p) then return SplitCartanConjugate(H:M:=M); end if;
    // Get an element of order p
    while true do h:=Random(H); n:=Order(h); if IsDivisibleBy(n,p) then h:=h^ExactQuotient(n,p); break; end if; end while;
    G:=GL(2,p);
    b,A:=IsConjugate(G,h,G![1,1,0,1]);
    assert b;
    BH:=Conjugate(H,A);
    D:=sub<G|[G![h[1][1],0,0,h[2][2]]:h in Generators(BH)]>;
    D,a,b:=DiagonalSubgroupGenerators(D:M:=M); 
    HH:=sub<G|D,G![1,1,0,1]>;
    assert BH eq HH;
    return HH,a,b;
end function;

// given H contained in a non-split Cartan returns conjugate subgroup of standard NonSplitCartan generated by [a,b*s,b,a] where s is the minimal primitive root mod p and (a,b) is lexicograpahically minimal
NonSplitCartanConjugate := function(H:M:=[])
    local p,n,F,u,m,r,a,s,d,b,G,GL2,ZL2;
    p:=#BaseRing(H);
    assert p gt 2;
    b,g:=IsCyclic(H);
    assert b and chi(g) le 0;
    GL2:=GL(2,p);
    ZL2:=Scalars(p);
    HZ:=sub<H|g^ExactQuotient(Order(g),GCD(Order(g),p-1))>;
    if #M ne p-1 then M:=MinReps(p); end if;
    if IsScalar(H) then
        z:=M[#HZ];
        G:=sub<H|[z,0,0,z]>;
        assert G eq H;
        return G,z,0;
    end if;
    n:=ExactQuotient(#H,#HZ);
    F:=GF(p^2);
    u:=Integers(p)!(2+Trace(RootOfUnity(n,F)));
    m:=dets(H);
    r:=Integers(p)!M[m];
    maxdets := [r^i: i in [1..m]|GCD(m,i) eq 1];
    s:=Integers(p)!M[p-1];
    if u eq 0 then
        a:=0;
        B:=[Sqrt(-d/s):d in maxdets];
        B:= B cat [-b:b in B];
        I:=[i:i in [1..#B]|sub<ZL2|(GL2![0,B[i]*s,B[i],0])^2> eq HZ];
        b:=Min([B[i]:i in I]);
    else
        A:=[Sqrt(d*u)/2:d in maxdets];
        A:=A cat [-a:a in A];
        B:=[Sqrt((a^2-4*a^2/u)/s):a in A];
        I:=[i:i in [1..#A]|sub<ZL2|(GL2![A[i],B[i]*s,B[i],A[i]])^n> eq HZ];
        a:=Min([A[i]:i in I]);
        b:=Sqrt((a^2-4*a^2/u)/s);
        b:=Min([b,-b]);
    end if;
    G:=sub<GL2|[a,b*s,b,a]>;
    assert IsConjugate(GL2,G,H);
    return G,Integers()!a,Integers()!b;
end function;
    
// Given H with dihedral projective image conjugate to a subgroup of the normalizer of a split Cartan, returns conjugate subgroup in the normalizer of the standard SplitCartan
// If H meet SL(2,p) is in the Cartan, returns minimal a and b such that H=<[a,0,0,1/a],[0,b,-r/b,0]> where r is minimal generator of det(H)
// Otherwise returns minimal a,b,c such that H=<[a,0,0,1/a],[0,b,-1/b,0],[0,c,-r/c,0]>
NormalizerSplitCartanConjugate := function(H:M:=[])
    local p,GL2,SL2,ZL2,G,G0,H1,B,S,a,b,c,d,e,g;
    p:=#BaseRing(H);
    GL2:=GL(2,p);
    SL2:=SL(2,p);
    HZ:=H meet Scalars(p);
    // pick an index 2 abelian subgroup of H; if there is more than one choice (as permittedy by Lemma 3.16), opt for the cyclic choice.
    S:=[G`subgroup:G in MaximalSubgroups(H:IsAbelian:=true,IndexEqual:=2)|{chi(g):g in Generators(G`subgroup)} subset {0,1}];
    assert #S gt 0;
    i:=0; for j in [1..#S] do if IsCyclic(S[j]) then i:=j; break; end if; end for;
    G0:= i ne 0 select S[i] else S[1];
    if #M ne p-1 then M:=MinReps(p); end if;
    G,a,b:=SplitCartanConjugate(G0:M:=M);    // conjugate it into the SplitCartan
    B:=GL2![b,0,0,M[dets(G)]/b];                    // save generator of det(G) for use below
    // Now use Corollary 3.17 to distinguish the two subgroups of C^+ that intersect C in G, one will have det(G)=-det(H-G0) and the other will have det(G) and -det(H-G0) disjoint
    g:=GL2![0,1,1,0];
    repeat x:=Random(H); until not x in G0;
    if IsDivisibleBy(dets(G),Order(-Determinant(x)))  then
        e:=0; H1:=sub<GL2|G,g>;
    else
        d:=GL2![1,0,0,M[p-1]];  e:=ExactQuotient(p-1,#HZ);  g:=g*d^e; H1:=sub<GL2|G,g>;
    end if;
//    assert IsConjugate(GL2,H,H1);
    r:=GF(p)!M[dets(H)];
    S:=[h:h in Generators(H1 meet SL2)|not h in G];
    g:=[B^i*g:i in [1..dets(G)]|Determinant(B^i*g) eq r][1];
    c:=Min([(h*g)[1][2]:h in sub<GL2|[a,0,0,1/a]>]);
    if #S eq 0 then
        G:=sub<GL2|[a,0,0,1/a],[0,c,-r/c,0]>;
        assert IsConjugate(GL2,G,H);
        return G,Integers()!a,Integers()!c,0;
    else
        b:=Min([(S[1]*h)[1][2]:h in sub<GL2|[a,0,0,1/a]>]);
        G:=sub<GL2|[a,0,0,1/a],[0,b,-1/b,0],[0,c,-r/c,0]>;
        assert IsConjugate(GL2,G,H);
        return G,Integers()!a,Integers()!b,Integers()!c;
    end if;
end function;
    
// Given H with dihedral projective image conjugate to a subgroup of the normalizer of a non-split Cartan, returns conjugate subgroup in the normalizer of the standard NonSplitCartan
// If H meet SL(2,p) is in the Cartan, returns minimal a and b such that H=<[a,b*s,b,a],[1,0,0,-1]*d^e> where s is minimal generator of det(H) and d is a generator for the NonSplitCartan
NormalizerNonSplitCartanConjugate := function(H:M:=[])
    local p,GL2,G0,G,a,b,c,d,e,g;
    p:=#BaseRing(H);
    GL2:=GL(2,p);
    G0:=[G`subgroup:G in MaximalSubgroups(H:IsAbelian:=true,IndexEqual:=2)|{chi(g):g in Generators(G`subgroup)} subset {0,-1}][1];   // take any maximal abelian index 2 subgroup
    if #M ne p-1 then M:=MinReps(p); end if;
    G,a,b:=NonSplitCartanConjugate(G0:M:=M);    // conjugate it into the NonSplitCartan
    HZ:=G meet Scalars(p);
    // Now use Corollary 3.17 to distinguish the two subgroups of C^+ that intersect C in G, one will have det(G)=-det(H-G0) and the other will have det(G) and -det(H-G0) disjoint
    g:=GL2![1,0,0,-1];
    repeat x:=Random(H); until not x in G0;
    if IsDivisibleBy(dets(G),Order(-Determinant(x)))  then
        e:=0; H1:=sub<GL2|G,g>;
    else
        d:=NonSplitCartan(p).1;  e:=ExactQuotient(p-1,#HZ);  H1:=sub<GL2|G,g*d^e>;
    end if;
//    assert IsConjugate(GL2,H,H1);
    return H1,Integers()!a,Integers()!b,Integers()!e;
end function;

GroupLetters:=["G","B","Cs","Cn","Ns","Nn","A4","S4","A5"];

// id format is [p,d,n,a,b,c] where p is characteristic, d is (p-1)/#det(G), n+1 is an index into GroupLetters, and a,b,c are nonnegative integers as defined in the paper

GroupLabelFromId:=function(id)
    assert #id ge 3;
    assert IsDivisibleBy(id[1]-1,id[2]);
    // special case for S4, needed for backward compatibiliy, suppress the .2 (which is the typical case and only case with full det)
    if id[3] eq 7 and id[4] eq 2 then
        if id[2] eq 1 then return Sprintf("%o%o",id[1],GroupLetters[8]); else return Sprintf("%o%o[%o]",id[1],GroupLetters[8],id[2]); end if;
    end if;
    if id[2] eq 1 then
        if #id eq 3 then return Sprintf ("%o%o",id[1],GroupLetters[id[3]+1]); end if;
        if #id eq 4 then return Sprintf ("%o%o.%o",id[1],GroupLetters[id[3]+1],id[4]); end if;
        if #id eq 5 then return Sprintf ("%o%o.%o.%o",id[1],GroupLetters[id[3]+1],id[4],id[5]); end if;
        if #id eq 6 then return Sprintf ("%o%o.%o.%o.%o",id[1],GroupLetters[id[3]+1],id[4],id[5],id[6]); end if;
    else
        if #id eq 3 then return Sprintf ("%o%o[%o]",id[1],GroupLetters[id[3]+1],id[2]); end if;
        if #id eq 4 then return Sprintf ("%o%o.%o[%o]",id[1],GroupLetters[id[3]+1],id[4],id[2]); end if;
        if #id eq 5 then return Sprintf ("%o%o.%o.%o[%o]",id[1],GroupLetters[id[3]+1],id[4],id[5],id[2]); end if;
        if #id eq 6 then return Sprintf ("%o%o.%o.%o.%o[%o]",id[1],GroupLetters[id[3]+1],id[4],id[5],id[6],id[2]); end if;
    end if;
    assert false;
end function;

GroupIdFromLabel:=function(label)
    local j,p,B,C,d,A,id;
    for i:=1 to #GroupLetters do
        j:=Position(label,GroupLetters[i]);
        if j eq 0 then continue; end if;
        p:=StringToInteger(Substring(label,1,j-1));
        B:=Split(label,"[");
        if #B gt 1 then
            C:=Split(B[2],"]")[1];
            d:=StringToInteger(C);
        else
            d:=1;
        end if;
        A:=Split(B[1],".");
        id:=[p,d,i-1] cat [StringToInteger(A[k]):k in [2..#A]];
        return id;
    end for;
    print "Invalid subgroup label", label;
    assert false;
end function;

GL2SubgroupID := function(H:PH:=[],M:=[])
    local p,GL2,SL2,ZL2,PG,pi,d,G,a,b;
    p := #BaseRing(H);
    if p eq 2 then
        if #H eq 1 then return [p,1,2]; end if;
        if #H eq 2 then return [p,1,1]; end if;
        if #H eq 3 then return [p,1,3]; end if;
        assert #H eq 6;
        return [p,1,0];
    end if;
    assert IsPrime(2);
    H:=ChangeRing(H,GF(p));
    GL2:=GL(2,p);
    SL2:=SL(2,p);
    d := ExactQuotient(p-1,ExactQuotient(#H, #(H meet SL2)));
    if SL2 subset H then return [p,d,0]; end if;
    HZ:=H meet Scalars(p);
    if #M ne p-1 then M:=MinReps(p); end if;
    if #H mod p eq 0 then
        if d*#H eq p*(p-1)^2 then return [p,d,1]; end if;
        G,a,b:=BorelConjugate(H);
        return [p,d,1,a,b];
    end if;
    if ExactQuotient(#H,#HZ) le 60 then
        PH:=quo<H|H meet HZ>;
        if IsIsomorphic(PH,AlternatingGroup(4)) then
            if dets(H) eq dets(HZ) then
                return [p,d,6,1];
            else
                assert 3*dets(HZ) eq dets(H);
                assert p mod 3 eq 1;
                return [p,d,6,3];
            end if;
        end if;
        if IsIsomorphic(PH,SymmetricGroup(4)) then
            if dets(H) eq 2*dets(HZ) then
                return [p,d,7,2];
            else
                assert d gt 1;
                assert dets(HZ) eq dets(H);
                return [p,d,7,1];
            end if;
        end if;
        if IsIsomorphic(PH,AlternatingGroup(5)) then
            assert p^2 mod 5 eq 1;
            assert dets(H) eq dets(HZ);
            return [p,d,8,1];
        end if;
    end if;
    if #PH eq 0 then PH:=quo<H|H meet HZ>; end if;
    if IsCyclic(PH) then
        if {chi(h):h in Generators(H)} subset {0,1} then
            if d*#H eq (p-1)^2 then return [p,d,2]; end if;
            G,a,b:=SplitCartanConjugate(H);
            return [p,d,2,a,b];
        else
            assert {chi(h): h in Generators(H)} subset {0,-1};
            if d*#H eq p^2-1 then return [p,d,3]; end if;
            H1,a,b:=NonSplitCartanConjugate(H);
            return [p,d,3,a,b];
        end if;
    end if;
    S:= [G`subgroup:G in MaximalSubgroups(H:IsAbelian:=true,IndexEqual:=2)|{chi(g):g in Generators(G`subgroup)}subset {0,1}];
    if #S ge 1 then
        H1,a,b,c:=NormalizerSplitCartanConjugate(H);
        if d*#H1 eq 2*(p-1)^2 then return [p,d,4]; end if;
        if c  gt 0 then return [p,d,4,a,b,c]; else return [p,d,4,a,b]; end if;
    end if;
    H1,a,b,c:=NormalizerNonSplitCartanConjugate(H);
    if d*#H1 eq 2*(p^2-1) then return [p,d,5]; end if;
    if c gt 0 then return [p,d,5,a,b,c]; else return [p,d,5,a,b]; end if;
end function;

GL2SubgroupLabel := function(H:PH:=[],M:=[])
    return GroupLabelFromId(GL2SubgroupID(H:PH:=PH,M:=M));
end function;
    
MatFrob := function(A);
    p:=Characteristic(BaseRing(A));
    return Parent(A)![A[1][1]^p,A[1][2]^p,A[2][1]^p,A[2][2]^p];
end function;

// Given a subgroup G of GL(2,p^2) that is conjugate to a subgroup H of GL(2,p), returns such an H, using the algorithm in Glasby and Howlett (Writing representations over minimal fields, Comm. Alg. 25 (1997), 1703-1711).
ConjugateToRationalSubgroup := function(G)
    local F,p,r,x,C,mu,R,v,X,A,H;
    F:=BaseRing(G);
    if Degree(F) eq 1 then return G; end if;
    assert Degree(F) eq 2;
    p:=Characteristic(F);
    r := [];
    for g in Generators(G) do
        r:=Append(r,[g[1][1]-g[1][1]^p,-g[2][1]^p,g[1][2],0]);
        r:=Append(r,[-g[1][2]^p,g[1][1]-g[2][2]^p,0,g[1][2]]);
        r:=Append(r,[g[2][1],0,g[2][2]-g[1][1]^p,-g[2][1]^p]);
        r:=Append(r,[0,g[2][1],-g[1][2]^p,g[2][2]-g[2][2]^p]);
    end for;
    while true do
        x:=Random(NullspaceOfTranspose(Matrix(r)));
        C:=MatrixRing(F,2)![x[i]:i in [1..Degree(x)]];
        if IsInvertible(C) then C:=GL(2,F)!C; break; end if;
    end while;
    for g in Generators(G) do if C^-1*g*C ne MatFrob(g) then print C, g; assert false; end if; end for;
    mu := C*MatFrob(C);
    assert mu[1][1] eq mu[2][2] and mu[1][2] eq 0 and mu[2][1] eq 0;
    mu := GF(p)!mu[1][1];
    b,v:=NormEquation(F,mu);
    C:=GL(2,F)![1/v,0,0,1/v]*C;
    assert C*MatFrob(C) eq Identity(G);
    while true do
        X:=Random(MatrixRing(F,2));
        A:=X+C*MatFrob(X);
        if not IsInvertible(A) then continue; end if;
        A:=GL(2,F)!A;
        H:=Conjugate(G,A);
        for h in Generators(H) do assert MatFrob(h) eq h; end for;
        return sub<GL(2,p)|Generators(H)>;
    end while;
end function;

//. Implements Algorithm 1 of the paper based on Thm 5.5 of Flannery-O'Brien (Linear groups of small degree over finite fields, Internat. J. Algebra Comput.  15 (2005), 467--502),
// d indicates the index of det(H) in (Z/pZ)* and a = [det(H):det(Z(H))] (slightly different inputs than in the paper)
ConstructA4Subgroup := function(p,d,a:verify:=false)
    local F,G,w,s,t,z,H,v,Z,PG,pi;
    assert p ge 5;  // we don't consider H with pi(H)=A4 an exceptional subgroup of GL(2,3) because in this case H is equal to SL(2,3)
    assert a eq 1 or ((p mod 3 eq 1) and a eq 3);
    assert IsDivisibleBy(p-1,a*d);
    if p mod 4 eq 1 then
        F:=GF(p);
    else
        F:=GF(p^2);
    end if;
    G:=GL(2,F);
    w:=RootOfUnity(4,F);
    s:=G![(w-1)/2,(w-1)/2,(w+1)/2,-(w+1)/2];
    t:=G![w,0,0,-w];
    z:=F!PrimitiveRoot(p);
    if d mod 2 eq 0 then e:=ExactQuotient(d,2); else e:=d; end if;   // forces -1 into Z(H)
    if a eq 1 then
        H:=ConjugateToRationalSubgroup(sub<G|s,t,[z^e,0,0,z^e]>);
    else
        v:=z^e;
        v := G![v,0,0,v];
        H:=ConjugateToRationalSubgroup(sub<G|v*s,t>);
    end if;
    // the isomorphism test below is slow so only do it when asked to
    if verify then
        Z:=Center(GL(2,p));
        PG,pi:=quo<GL(2,p)|Z>;
        assert IsIsomorphic(pi(H),AlternatingGroup(4)) and ExactQuotient(p-1,dets(H)) eq d and a*dets(H meet Z) eq dets(H);
    end if;
    return H;
end function;

//. Implements Thm 5.8 of Flannery-O'Brien
// d indicates the index of det(H) in (Z/pZ)* and a = [det(H):det(Z(H))] (slightly different inputs than in the paper)
ConstructS4Subgroup := function(p,d,a:verify:=false)
    local F,G,w,s,c,t,z,zzmH,u,uz,Z,PG,pi;
    assert p ge 5;  // we don't consider H with pi(H)=A4 an exceptional subgroup of GL(2,3) because in this case H is equal to GL(2,3)
    assert IsDivisibleBy(p-1,a*d);
    if a eq 1 then assert p mod 8 eq 1 or p mod 8 eq 7; else assert a eq 2; end if;
    if p mod 8 eq 1 or p mod 8 eq 7 then assert d mod 2 eq 0; end if;
    if p mod 8 eq 5 then assert ExactQuotient(p-1,d) mod 4 eq 0; end if;
    F:=GF(p^2);
    G:=GL(2,F);
    w:=RootOfUnity(4,F);
    s:=G![(w-1)/2,(w-1)/2,(w+1)/2,-(w+1)/2];
    c:=Sqrt(F!2);
    t:=G![(1+w)/c,0,0,(1-w)/c];
    z:=F!PrimitiveRoot(p);
    if d mod 2 eq 0 then e:=ExactQuotient(d,2); else e:=d; end if;   // forces -1 into Z(H)
    zz:= G![z^e,0,0,z^e];
    if a eq 1 then
        H:=ConjugateToRationalSubgroup(sub<G|s,t,zz>);
    else
        if p mod 8 eq 1 then zz := zz^2; e:=d; end if;
        u := Sqrt(z^e);
        uz:= G![u,0,0,u];
        H:=ConjugateToRationalSubgroup(sub<G|s,uz*t,zz>);
    end if;
    // the isomorphism test below is slow so only do it when asked to
    if verify then
        Z:=Center(GL(2,p));
        PG,pi:=quo<GL(2,p)|Z>;
        assert IsIsomorphic(pi(H),SymmetricGroup(4)) and ExactQuotient(p-1,dets(H)) eq d and a*dets(H meet Z) eq dets(H);
    end if;
    return H;
end function;

//. Implements Thm 5.11 of Flannery-O'Brien
// d indicates the index of det(H) in (Z/pZ)* and a = [det(H):det(Z(H))], which must be 1 (slightly different inputs than in the paper)
ConstructA5Subgroup := function(p,d,a:verify:=false)
    local e,s,t,z,w,F,G,H,Z;
    assert IsDivisibleBy(p-1,a*d);
    assert p^2 mod 5 eq 1 and d mod 2 eq 0 and a eq 1;
    F:=GF(p^2);
    G:=GL(2,F);
    w:=RootOfUnity(4,F);
    b := Sqrt(F!5);
    s:=G![(w-1)/2,(w-1)/2,(w+1)/2,-(w+1)/2];
    t:=G![w,0,0,-w];
    v:=G![w/2,(1-b)/4-w*(1+b)/4,(-1+b)/4-w*(1+b)/4,-w/2];
    z:=F!PrimitiveRoot(p);
    if d mod 2 eq 0 then e:=ExactQuotient(d,2); else e:=d; end if;   // forces -1 into Z(H)
     H:=ConjugateToRationalSubgroup(sub<G|s,t,v,[z^e,0,0,z^e]>);
    // the isomorphism test below is slow so only do it when asked to
    if verify then
        Z:=Center(GL(2,p));
        PG,pi:=quo<GL(2,p)|Z>;
        assert IsIsomorphic(pi(H),AlternatingGroup(5)) and ExactQuotient(p-1,dets(H)) eq d and a*dets(H meet Z) eq dets(H);
    end if;
    return H;
end function;

// returns the GL2 subgroup corresponding to the specified label -- for the sake of speed this function does not verify that the label is the unique label for the group (as specified in Section 7 of the paper), you can apply GL2SubgroupLabel to the result to check this
GL2SubgroupFromLabel := function(id:M:=[])
    local p,x,y,s,a,b,c,G;
    if Type(id) eq MonStgElt then id := GroupIdFromLabel(id); end if;
    p:=id[1]; d:=id[2];
    if #M ne p-1 then M:=MinReps(p); end if;
    r:=Integers(p)!M[ExactQuotient(p-1,d)];
    e:=Integers(p)!M[p-1];
    G:=GL(2,p);
    if id[3] eq 0 then return sub<G|[1,1,0,1],[1,0,1,1],[1,0,0,r]>; end if;
    if id[3] eq 1 then
        if #id eq 3 then a:= e; b:=1; else assert #id eq 5; a:=id[4]; b:=id[5]; end if;
        return sub<G|[a,0,0,1/a],[b,0,0,r/b],[1,1,0,1]>;
    end if;
    if id[3] eq 2 then
        if #id eq 3 then a:= e; b:=1; else assert #id eq 5; a:=id[4]; b:=id[5]; end if;
        return sub<G|[a,0,0,1/a],[b,0,0,r/b]>;
    end if;
    if id[3] eq 3 then
        if p eq 2 then assert #id eq 3; return sub<G|[0,1,1,1]>; end if;
        if #id eq 3 then
            return sub<G|NonSplitCartan(p).1^d>;
        else
            assert #id eq 5;
            return sub<G|[id[4],id[5]*e,id[5],id[4]]>;
        end if;
    end if;
    if id[3] eq 4 then
        if #id eq 3 then return  sub<G|[e,0,0,1/e],[1,0,0,r/1],[0,-1,1,0]>; end if;
        assert #id ge 5;
        if #id eq 5 then return sub<G|[id[4],0,0,1/id[4]],[0,id[5],-r/id[5],0]>; end if;
        return sub<G|[id[4],0,0,1/id[4]],[0,id[5],-1/id[5],0],[0,id[6],-r/id[6],0]>;
    end if;
    if id[3] eq 5 then 
        if #id eq 3 then 
            if d eq 1 then
                return NormalizerNonSplitCartan(p);
            else
                h:=NonSplitCartan(p).1;  n:=Log(Determinant(h),GF(p)!-1);
                return  sub<G|h^d,h^n*G![1,0,0,-1]>;
            end if;
        end if;
        assert #id ge 5;
        if #id eq 5 then return sub<G|[id[4],e*id[5],id[5],id[4]],[1,0,0,-1]>; end if;
        d:=NonSplitCartan(p).1;
        return sub<G|[id[4],e*id[5],id[5],id[4]],G![1,0,0,-1]*d^id[6]>;
    end if;
    if id[3] eq 6 then return ConstructA4Subgroup(p,d,id[4]); end if;
    if id[3] eq 7 then if #id lt 4 then c:=2; else c:=id[4]; end if; return ConstructS4Subgroup(p,d,c); end if;
    if id[3] eq 8 then return ConstructA5Subgroup(p,d,1); end if;
    assert false;
end function;

// returns list of subgroups containing SL(2,p)
SL2Subgroups := function(p)
    local r;
    r := GF(p)!PrimitiveRoot(p);
    return [sub<GL(2,p)|[1,1,0,1],[1,0,1,1],[1,0,0,r^a]>:a in Divisors(p-1)];
end function;

// returns list subgroups of the Borel not in the split Cartan, using Lemma 5
BSubgroups := function(p)
    local r,S;
    r := GF(p)!PrimitiveRoot(p);
    S:=[];
    for a,b in Divisors(p-1) do
        S:= S cat [sub<GL(2,p)|[r^a,0,0,1/r^a],[r^(i*ExactQuotient(a,GCD(a,b))),0,0,r^ExactQuotient(p-1,b)/r^(i*ExactQuotient(a,GCD(a,b)))],[1,1,0,1]>: i in [0..GCD(a,b)-1]];
    end for;
    return S;
end function;
    
// function to filter conjugate subgroups of the split Cartan
IsMinCsConjugate := function(n,a,b,i)
    local c,d,j;
    try
        c := ExactQuotient(a,GCD(a,b));
        d := ExactQuotient(n,b);
        j := Solution(Integers(a)!c, Integers(a)!(d-i*c));
        return j ge i;
    catch e
        return true;
    end try;
end function;

// returns list of non-conjugate subgroups of split Cartan, including scalar subgroups
CsSubgroups := function(p)
    local r,S;
    r := GF(p)!PrimitiveRoot(p);
    S:=[];
    for a,b in Divisors(p-1) do
        S:= S cat [sub<GL(2,p)|[r^a,0,0,1/r^a],[r^(i*ExactQuotient(a,GCD(a,b))),0,0,r^ExactQuotient(p-1,b)/r^(i*ExactQuotient(a,GCD(a,b)))]>: i in [0..GCD(a,b)-1] | IsMinCsConjugate(p-1,a,b,i)];
    end for;
    return S;
end function;

// returns list of non-scalar subgroups of non-split Cartan
CnSubgroups := function(p)
    local g,G;
    G:=GL(2,p);   
    if p eq 2 then
        g:=G![1,1,1,0];
    else
        g:=NonSplitCartan(p).1;
    end if;
    return [sub<G|g^n>:n in Divisors(p^2-1)|not IsDivisibleBy(n,p+1)];
end function;

IsNormalCsSubgroup := function(n,a,b,i)
        return IsDivisibleBy(ExactQuotient(2*a*i,GCD(a,b))-ExactQuotient(n,b),a);
end function;

NormalCsSubgroups := function(p)
    local r,S;
    r := GF(p)!PrimitiveRoot(p);
    S:=[];
    for a,b in Divisors(p-1) do
        S:= S cat [sub<GL(2,p)|[r^a,0,0,1/r^a],[r^(i*ExactQuotient(a,GCD(a,b))),0,0,r^ExactQuotient(p-1,b)/r^(i*ExactQuotient(a,GCD(a,b)))]>: i in [0..GCD(a,b)-1] | IsNormalCsSubgroup(p-1,a,b,i)];
    end for;
    return S;
end function;

// Check for duplicate dihedral subgroups of C_s+ (which are necessarily isomorphic to the Klein 4-group) , using Lemma 18 (keep the one where H is cyclic)
IsKleinDuplicateNCs := function(H,nZ)
    local p,g;
    p:=#BaseRing(H);
    g:=GL(2,BaseRing(H))![1,0,0,-1];
    return ExactQuotient(p-1,nZ) mod 2 eq 0 and 2*nZ eq #H and g in H;
end function;

// returns a list of subgroups of the normalizer of the split Cartan with dihedral projective image
NCsSubgroups := function(p)
    local G,Z,S,SZ,gamma,delta,nI,T1,T2;
    G:=GL(2,p);
    Z:=Scalars(p);
    S:=[H:H in NormalCsSubgroups(p)|not IsScalar(H)];
/*
    This code works but isn't significantly faster than letting Magma compute the intersection of H and Z, so we  take the simpler approach
    ZZ:=Integers();
    g:=G![1,0,0,PrimitiveRoot(p)];
    T:=[0:i in [1..p]];
    for e in [1..p-1] do T[1+ZZ!(Trace(a)^2/Determinant(a)) where a:=g^e] := ExactQuotient(p-1,GCD(p-1,e)); end for;
    time SZ:= [sub<Z|#Generators(H) eq 1 select Z.1^ExactQuotient((p-1)*T[1+ZZ!(Trace(H.1)^2/Determinant(H.1))],#H) else Z.1^ExactQuotient((p-1)*LCM(T[1+ZZ!(Trace(H.1)^2/Determinant(H.1))],T[1+ZZ!(Trace(H.2)^2/Determinant(H.2))]),#H)>:H in S];
*/
    SZ:=[H meet Z:H in S];
    gamma:=G![0,1,1,0];
    delta:=G![1,0,0,PrimitiveRoot(p)];
    nI:=G![-1,0,0,-1];
    T1 := (p mod 4 ne 1) select [sub<G|H,gamma>:H in S] else [sub<G|S[i],gamma>:i in [1..#S]| not IsKleinDuplicateNCs(S[i],#SZ[i])];// duplicates can only occur in the det(H) = - det(G-H) case
    T2:=[sub<G|S[i],gamma*delta^ExactQuotient(p-1,#SZ[i])>:i in [1..#S] | nI in S[i] and dets(S[i]) eq dets(SZ[i])];
    return T1 cat T2;
end function;

// returns a list of subgroups of the normalizer of the nonsplit Cartan with dihedral proejctive image
NCnSubgroups := function(p)
    local G,Z,S,SZ,gamma,delta,nI,T1,T2,NCs;
    G:=GL(2,p);
    Z:=Scalars(p);
    S:=[H:H in CnSubgroups(p)];
    SZ:=[sub<H|H.1^ExactQuotient(Order(H.1),GCD(Order(H.1),p-1))>:H in S];  // fast computation of scalar intersections -- this is much faster than letting Magma compute H meet Z
    gamma:=G![1,0,0,-1];
    delta:=NonSplitCartan(p).1;
    nI:=G![-1,0,0,-1];
    T1:=[sub<G|H,gamma>:H in S];
    T2:=[sub<G|S[i],gamma*delta^ExactQuotient(p-1,#SZ[i])>:i in [1..#S] | nI in S[i] and dets(S[i]) eq dets(SZ[i])];
    NCs:=NormalizerSplitCartan(p);
    return [G:G in T1 cat T2|not G subset NCs];
end function;
   

//. Implements Thm 5.5 of Flannery-O'Brien (Linear groups of small degree over finite fields, Internat. J. Algebra Comput.  15 (2005), 467--502),
A4Subgroups := function(p)
    local F,GG,w,s,t,z,S,v;
    if p lt 5 then return[]; end if;    // don't count A4 in GL(2,3) containing SL(2,3)
    if p mod 4 eq 1 then
        F:=GF(p);
    else
        F:=GF(p^2);
    end if;
    GG:=GL(2,F);
    w:=RootOfUnity(4,F);
    s:=GG![(w-1)/2,(w-1)/2,(w+1)/2,-(w+1)/2];
    t:=GG![w,0,0,-w];
    z:=F!PrimitiveRoot(p);
    S:=[];
    for d in Divisors(p-1) do
        if ExactQuotient(p-1,d) mod 2 eq 1 then continue; end if;
        S:=Append(S,ConjugateToRationalSubgroup(sub<GG|s,t,[z^d,0,0,z^d]>));
        if d mod 3 eq 0 then
            v:=z^ExactQuotient(d,3);
            v := GG![v,0,0,v];
            S:=Append(S,ConjugateToRationalSubgroup(sub<GG|v*s,t>));
        end if;
    end for;
    return S;
end function;

//. Implements Thm 5.8 of Flannery-O'Brien
S4Subgroups := function(p)
    local F,GG,w,s,a,t,z,S,zz,u,uz;
    if p lt 5 then return[]; end if;    // don't count S4 == GL(2,3)
    F:=GF(p^2);
    GG:=GL(2,F);
    w:=RootOfUnity(4,F);
    s:=GG![(w-1)/2,(w-1)/2,(w+1)/2,-(w+1)/2];
    a:=Sqrt(F!2);
    t:=GG![(1+w)/a,0,0,(1-w)/a];
    z:=F!PrimitiveRoot(p);
    S:=[];
    for d in Divisors(p-1) do
        if ExactQuotient(p-1,d) mod 2 eq 1 then continue; end if;
        zz:= GG![z^d,0,0,z^d];
        u := Sqrt(z^d);
        uz:= GG![u,0,0,u];
        if p mod 8 in [1,7] then
            S:=Append(S,ConjugateToRationalSubgroup(sub<GG|s,t,zz>));
            if u in GF(p) then
                S:=Append(S,ConjugateToRationalSubgroup(sub<GG|s,uz*t,zz>));
            end if;
        else
            if u*a in GF(p) then
                S:=Append(S,ConjugateToRationalSubgroup(sub<GG|s,uz*t,zz>));
            end if;
        end if;
    end for;
    return S;
end function;

//. Implements Thm 5.11 of Flannery-O'Brien
A5Subgroups := function(p)
    local F,GG,w,b,s,t,v,z,S;
    if p^2 mod 5 ne 1 then return []; end if;
    F:=GF(p^2);
    GG:=GL(2,F);
    w:=RootOfUnity(4,F);
    b := Sqrt(F!5);
    s:=GG![(w-1)/2,(w-1)/2,(w+1)/2,-(w+1)/2];
    t:=GG![w,0,0,-w];
    v:=GG![w/2,(1-b)/4-w*(1+b)/4,(-1+b)/4-w*(1+b)/4,-w/2];
    z:=F!PrimitiveRoot(p);
    S:=[];
    for d in Divisors(p-1) do
        if ExactQuotient(p-1,d) mod 2 eq 1 then continue; end if;
        S:=Append(S,ConjugateToRationalSubgroup(sub<GG|s,t,v,[z^d,0,0,z^d]>));
    end for;
    return S;
end function;

// returns a complete list of solvable non-conjugate subgroups of GL(2,p)
GL2Subgroups := function(p)
    assert IsPrime(p);
    if p eq 2 then G:=GL(2,2); return [sub<G|[1,0,0,1]>,sub<G|[1,1,0,1]>,sub<G|[1,1,1,0]>,G]; end if;
    return SL2Subgroups(p) cat BSubgroups(p) cat CsSubgroups(p) cat CnSubgroups(p) cat NCsSubgroups(p) cat NCnSubgroups(p) cat A4Subgroups(p) cat S4Subgroups(p) cat A5Subgroups(p);
end function;

// returns a complete list of non-conjugate subgroups of GL(2,p)
GL2SolvableSubgroups := function(p)
    assert IsPrime(p);
    if p le 3 then return GL2Subgroups(p); end if;
    return BSubgroups(p) cat CsSubgroups(p) cat CnSubgroups(p) cat NCsSubgroups(p) cat NCnSubgroups(p) cat A4Subgroups(p) cat S4Subgroups(p);
end function;

// returns label of the unique non-conjugate group locally conjugate to the specified group (or lable), if any, otherwise returns the empty string
GL2LocalConjugate := function(label)
    local id,H1,H2;
    if Type(label) ne MonStgElt then label := GL2SubgroupLabel(label); end if;
    id := GroupIdFromLabel(label);
    if id[3] ne 1 then return ""; end if;
    H1:=GL2SubgroupFromLabel(id);
    H2:=sub<GL(2,#BaseRing(H1))|[[g[2][2],g[1][2],0,g[1][1]]:g in Generators(H1)]>;
    if H1 eq H2 then return ""; end if;
    return GL2SubgroupLabel(H2);
end function;

// returns set of one or two labels of groups locally conjugate to the given group
GL2LocalConjugates := function(label)
    local id,H1,H2;
    if Type(label) ne MonStgElt then label := GL2SubgroupLabel(label); end if;
    id := GroupIdFromLabel(label);
    if id[3] ne 1 then return {label}; end if;
    H1:=GL2SubgroupFromLabel(id);
    H2:=sub<GL(2,#BaseRing(H1))|[[g[2][2],g[1][2],0,g[1][1]]:g in Generators(H1)]>;
    return {label,GL2SubgroupLabel(H2)};
end function;

// returns a list of non-conjugate locally conjugate pairs of group labels in GL(2,p)
GL2LocallyConjugatePairs := function(p)
    local S,L,H1,H2,s1,s2;
    S:=BSubgroups(p);
    L:={};
    for H1 in S do
        H2:=sub<Borel(p)|[[g[2][2],g[1][2],0,g[1][1]]:g in Generators(H1)]>;
        if H2 ne H1 then L := L join {{GL2SubgroupLabel(H1),GL2SubgroupLabel(H2)}}; end if;
    end for;
    return L;
end function;

// returns a list of labels of non-conjugate quadratic twists of a specified group (or label)
GL2Twists := function(H)
    local G,S,T,GL2,label;
    if Type(H) eq MonStgElt then H:=GL2SubgroupFromLabel(H); end if;
    GL2:=GL(2,#BaseRing(H));
    G:=sub<GL2|H,-Identity(H)>;
    if #G mod 2 ne 0 then return {GL2SubgroupLabel(H)}; end if;  // necessary check to work around magma bug in Subgroups
    S:=[K`subgroup:K in Subgroups(G:IndexEqual:=2)|not -Identity(K`subgroup) in K`subgroup];
    S:=[S[i]:i in [1..#S]| #{j:j in [1..i-1]|IsConjugate(GL2,S[i],S[j])} eq 0];
    S:=S cat [G];
    return {GL2SubgroupLabel(G):G in S};
end function;

// returns a list of twist orbits in GL(2,p), specified by labels
GL2TwistOrbits := function (p)
    return {GL2Twists(H):H in GL2Subgroups(p)};
end function;
    
// returns the index of a specifed subgroup of label in GL2
GL2SubgroupIndex := function(G)
    if Type(G) eq MonStgElt then G:=GL2SubgroupFromLabel(G); end if;
    return Index(GL(2,BaseRing(G)),G);
end function;
    
// returns true if  H1 contains a GL2-conjugate of H2
GL2SubgroupContains := function(H1,H2)
    if Type(H1) eq MonStgElt then H1:=GL2SubgroupFromLabel(H1); end if;
    if Type(H2) eq MonStgElt then H2:=GL2SubgroupFromLabel(H2); end if;
    if BaseRing(H1) ne BaseRing(H2) then return false; end if;
    return #{G:G in Conjugates(GL(2,BaseRing(H2)),H2)|G subset H1} gt 0;
end function;
    
// returns true if H is exceptional for the local/global principal for isogenies (this means every element of H has an eigenvalue but H is not in a Borel)
GL2SubgroupExceptional := function(H)
    if Type(H) eq MonStgElt then H:=GL2SubgroupFromLabel(H); end if;
    if GL2SubgroupContains(Borel(#BaseRing(H)),H) then return false; end if;
    if #{1:h in H|chi(h) eq -1} gt 0 then return false; end if;
    return true;
end function;
    
// returns a list of the labels of subgroups of a given subgroup of GL2
GL2SubgroupLabels := function(G: IndexEqual:=0)
    if Type(G) eq MonStgElt then G:=GL2SubgroupFromLabel(G); end if;
    if IndexEqual gt 0 then
        return {GL2SubgroupLabel(H`subgroup):H in Subgroups(G:IndexEqual:=IndexEqual)|Index(G,H`subgroup) eq IndexEqual};   // check index to deal with magma bug
    else
        return {GL2SubgroupLabel(H`subgroup):H in Subgroups(G)};   // check index to deal with magma bug
    end if;
end function;
        
// determine whether two lists of subgroups of GL(2,p) are the same up to ordering and conjugacy
IsConjugateList := function(S,T)
    if #S ne #T then return false; end if;
    G:=GL(2,#BaseRing(S[1]));
    Smap := Sort([<#S[i],i>:i in [1..#S]]);
    Tmap := Sort([<#T[i],i>:i in [1..#T]]);
    i:=1;
    while i le #Smap do
        j:=i+1;
        while j le #Smap and Smap[j][1] eq Smap[i][1] do j+:=1; end while;
        for k in [i..j-1] do
            if #[m:m in [i..j-1]|IsConjugate(G,S[Smap[k][2]],T[Tmap[m][2]])] ne 1 then return false; end if;
        end for;
        i := j;
    end while;
    return true;
end function;

// Test that GL2Subgroups agrees with magma
test1 :=procedure(p:verify:=true)
    local S,T,s1,s2,t1,t2;
    s1:= Cputime();
    S:=GL2Subgroups(p);
    t1 := Cputime()-s1;
    s2 := Cputime();
    T:=[H`subgroup:H in Subgroups(GL(2,p))];
    t2 := Cputime()-s2;
    if verify then
        assert IsConjugateList(S,T);
        printf "p=%o, success, Algorithm 2 time %o Magma time %o\n",p,t1,t2;
    else
        printf "p=%o, Algorithm 2 time %o Magma time %o\n",p,t1,t2;
    end if;
end procedure;

// Test that we can determine the label and the reconstruct from the label every group in GL2Subgroups
test2:=procedure(p:verify:=true)
    local G,S,Hs,id,label,s1,s2,t1,t2;
    assert IsPrime(p);
    G:=GL(2,p);
    M:=MinReps(p);
    for H in GL2Subgroups(p) do
        s1 := Cputime();
        label := GL2SubgroupLabel(H:M:=M);
        t1 := Cputime()-s1;
        s2:=Cputime();
        Hs := GL2SubgroupFromLabel(label:M:=M);
        t2:=Cputime()-s2;
        if verify and not IsConjugate(G,H,Hs) then print p, label, H,Hs; assert false; end if;
//        if t1 gt 0.1 or t2 gt 0.1 then print label,t1,t2; end if;
    end for;
    print p,"success";
end procedure;

// Test that we can determine the label and the reconstruct from the label every group in GL2Subgroups even after it has been randomly conjugated
test3:=procedure(p)
    local G,Hr,Hs,label;
    assert IsPrime(p);
    G:=GL(2,p);
    for H in GL2Subgroups(p) do
        Hr:=Random(Conjugates(G,H));
        label :=GL2SubgroupLabel(Hr);
        Hs:=GL2SubgroupFromLabel(label);
        if not IsConjugate(G,Hr,Hs) then print p, Hr,Hs; assert false; end if;
        assert GL2SubgroupLabel(Hs) eq label;
    end for;
    print p,"success";
end procedure;

// Test that GL2Subgroups agrees with GL2Counts
test4 :=procedure(p:verify:=true)
    local S,T,s1,s2,t1,t2;
    s1:= Cputime();
    S:=GL2Subgroups(p);
    t1 := Cputime()-s1;
    assert #S eq GL2Counts(p);
    printf "p=%o, time %o\n",p,t1;
end procedure;

