GL9 := GL(2,Integers(9));
GL27 := GL(2,Integers(27));

// Much of this code comes from the file "S3entanglement" associated
// to the paper "Odd degree isolated points on X_1(N) with rational
// j-invariant" by Bourdon, Gill, Rouse, and Watson. The group
// generators come from Table 1 on page 23 of the arXiv version of
// Sutherland and Zywina's "Modular curves of prime-power level
// with infinitely many rational points."

// These matrices represent the action of Galois on E[3^k] = (Z/3^k Z)^2.
// The elements of (Z/3^k Z)^2 are represented as column vectors.

H9:=sub<GL9 | GL9![2,0,0,5], GL9![4,2,3,4], GL9![1,0,0,2]>; // 9C0-9a

// A slight modification of Cummins and Pauli's code for using Riemann-Hurwitz
// to compute the genus of X_G.

function genus2(G)
  GG := sub<GL(2,BaseRing(G)) | SetToSequence(Generators(G)) cat [GL(2,BaseRing(G))![-1,0,0,-1]]>;
  md := Modulus(BaseRing(GG));
  H := SL(2,Integers(md));
  S := H![0,-1,1,0];
  T := H![1,1,0,1];
  phi, perm := CosetAction(H,GG meet H);
  lst := [phi(S),phi(T),phi(S*T)];
  //printf "Permutation for S = %o.\n",phi(S);
  //printf "Permutation for T = %o.\n",phi(T);
  //printf "Permutation for S*T = %o.\n",phi(S*T);
  cs := [CycleStructure(lst[i]) : i in [1..3]];
  gen := -2*Degree(perm) + 2;
  einfty := #Orbits(sub<perm | lst[2]>);
  e2 := #Fix(lst[1]);
  e3 := #Fix(lst[3]);
  ind := Degree(perm);
  for j in [1..3] do
    for i in [1..#cs[j]] do
      gen := gen + (cs[j][i][1]-1)*cs[j][i][2];
    end for;
  end for;
  gen := gen div 2;
  //printf "The genus = %o.\n",gen;
  genhur := 1 + (ind/12) - (e2/4) - (e3/3) - (einfty/2);
  //printf "The Hurwitz formula is %o = 1 + %o/12 - %o/4 - %o/3 - %o/2.\n",genhur, ind,e2,e3,einfty;
  return gen, ind, einfty, e2, e3;
end function;

// Take a subgroup G of GL(2,Z/rZ) and for a multiple s of r
// compute the full preimage of G in GL(2,Z/sZ).

function liftlevel(G,r,s)
  bigG := GL(2,Integers(s));
  redhom := hom<bigG -> GL(2,Integers(r)) | x :-> GL(2,Integers(r))!Eltseq(x)>;
  newgens := [];
  for k in Generators(Kernel(redhom)) do
    Append(~newgens,bigG!k);
  end for;
  F := Factorization(s);
  coprime := 1;
  if not (&and [ GCD(F[i][1],r) ne 1 : i in [1..#F]]) then
    coprime := &*[ F[i][1]^F[i][2] : i in [1..#F] | GCD(F[i][1],r) eq 1];
  end if;
  for k in Generators(G) do
    // Take k and make it a matrix mod s with determinant that's a unit
    // mod coprime
    newelt := bigG![CRT([Integers()!k[1][1],1],[r,coprime]),
    CRT([Integers()!k[1][2],0],[r,coprime]),
    CRT([Integers()!k[2][1],0],[r,coprime]),
    CRT([Integers()!k[2][2],1],[r,coprime])];
    Append(~newgens,newelt);
  end for;
  return sub<bigG | newgens>;
end function;

goodsubs := <>;
GL2 := GL(2,Integers(2));
GL54 := GL(2,Integers(54));

H27 := liftlevel(H9,9,27);
H54 := liftlevel(H9,9,54);
L := LowIndexSubgroups(H54,6);

for i in [1..#L] do
  curH := L[i];
  hom1 := hom<curH -> GL2 | x :-> GL2!Eltseq(x)>;
  hom2 := hom<curH -> GL27 | x :-> GL27!Eltseq(x)>;  
  if Index(H54,curH) eq 6 then
    if Image(hom2) eq H27 then
      K1 := Kernel(hom1);
      K2 := Kernel(hom2);
      if Index(curH,K1) eq 3 then
        if (K2 subset K1) then
       	      Append(~goodsubs,curH);
        end if;
      end if;  
    end if;
  end if;
end for;

// We find that there are four good subgroups (up to conjugacy).

function reduce(G,n)
  H:=GL(2,Integers(n));
  R:=sub<H|Generators(G)>;
  return R;
end function;

R1:=reduce(goodsubs[1],18);
R2:=reduce(goodsubs[2],18);
R3:=reduce(goodsubs[3],18);
R4:=reduce(goodsubs[4],18);

genus2(R1);
genus2(R2);
genus2(R3);
genus2(R4);

GL18 := GL(2,Integers(18));

H7:=sub<GL18|[11,16,6,17],[13,3,9,14]>;
H8:=sub<GL18|[5,2,12,11],[17,5,15,16]>;
H9:=sub<GL18|[11,1,3,4],[7,15,9,14]>;
H10:=sub<GL18|[16,1,3,11],[5,15,9,16]>;  

IsConjugateSubgroup(GL18,H7,R1);
IsConjugateSubgroup(GL18,H8,R2);
IsConjugateSubgroup(GL18,H9,R3);
IsConjugateSubgroup(GL18,H10,R4);