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

// Make a list of the subgroups that correspond to 9B0-9a, 9H0-9b,
// 9I0-9a, 9I0-9b, 9I0-9c, 9J0-9a, and 27A0-27a. These generators
// come from the Table 1 on page 23 of the arXiv version of
// Andrew Sutherland and David 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.

imlist := <
sub<GL9 | GL9![1,1,0,1], GL9![2,0,0,5], GL9![1,0,0,2]>, // 9B0-9a
sub<GL9 | GL9![1,3,0,1], GL9![5,0,3,2], GL9![2,1,0,1]>, // 9H0-9b
sub<GL9 | GL9![2,1,0,5], GL9![1,2,3,2]>, // 9I0-9a
sub<GL9 | GL9![2,1,0,5], GL9![4,0,3,5]>, // 9I0-9b
sub<GL9 | GL9![2,2,0,5], GL9![2,2,3,1]>, // 9I0-9c
sub<GL9 | GL9![1,3,0,1], GL9![2,2,3,8], GL9![1,2,0,2]>, // 9J0-9a
sub<GL27 | GL27![1,1,0,1], GL27![2,1,9,5], GL27![1,2,3,2]> // 27A0-27a
>;

// 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;

// What are the torsion orbits?

function torsorbits(G,r);
  a := Valuation(r,2);
  b := Valuation(r,3);
      ord := 2^a*3^b;
      GL2 := GL(2,Integers(ord));
      redgp := sub<GL2 | [ Transpose(GL2!Eltseq(g)) : g in Generators(G)]>;
      R := RSpace(redgp);
      // R is the module that Magma thinks that redgp naturally acts on.
      // That is, row vectors. That's why we take the transpose in the definition
      // of redgp.
      orbs := [];
      for a in R do
        if (ord*a eq 0) and ((ord div 2)*a ne 0) and ((ord div 3)*a ne 0) then
          if not a in (&join orbs) then
	    Append(~orbs,Orbit(redgp,a));
	  end if;
        end if;
      end for;
      retlist := Sort([ #o : o in orbs]);
      //printf "Orbits of action on points of order %o is %o.\n",ord,retlist;
  return retlist;
end function;

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

for curim in [1..#imlist] do
printf "Working on 3-adic image %o of %o.\n",curim,#imlist;
G := imlist[curim];
if Characteristic(BaseRing(G)) eq 9 then
  G27 := liftlevel(G,9,27);
else
  G27 := G;
end if;
G2 := liftlevel(G,Characteristic(BaseRing(G)),54);
L := LowIndexSubgroups(G2,6);

for i in [1..#L] do
  curG := L[i];
  hom1 := hom<curG -> GL2 | x :-> GL2!Eltseq(x)>;
  hom2 := hom<curG -> GL27 | x :-> GL27!Eltseq(x)>;  
  if Index(G2,curG) eq 6 then
    if Image(hom2) eq G27 then
      K1 := Kernel(hom1);
      K2 := Kernel(hom2);
      if Index(curG,K1) eq 6 then
        if (K2 subset K1) then
          tors54 := torsorbits(curG,54);
          if (6 in tors54) or (18 in tors54) or (54 in tors54) then
            printf "Subgroup %o of %o is a suitable mod 54 image.\n",i,#L;
            orbs := torsorbits(curG,18);
	    // Because our subgroups contain -I, the degree on X_1(18) is
	    // half the numbers in the list orbs.
	    if 6 in orbs then
	      printf "But subgroup %o gives a degree 3 point on X_1(18).\n",i;
	    else
	      Append(~goodsubs,curG);
	    end if;
          end if;
        end if;
      end if;  
    end if;
  end if;
end for;

end for;

// We find that there are four good subgroups.
IsConjugateSubgroup(GL54,goodsubs[1],goodsubs[2]);
IsConjugateSubgroup(GL54,goodsubs[1],goodsubs[3]);
IsConjugateSubgroup(GL54,goodsubs[1],goodsubs[4]);
// The latter three are conjugate in GL(2,Z/54Z) to subgroups of the first.

K := goodsubs[1];
// So there's only one modular curve we need to handle: X_K.
genus(K);

K2 := MinimalOvergroups(GL54,K)[1];
genus(K2);

v := Matrix(Integers(54),[[27],[0]]);
stabgenlist := [];
// The orbit of v has size 3.
done := false;
while done eq false do
  x := Random(K);
  if x*v eq v then
    if #sub<K | stabgenlist cat [x]> gt #sub<K | stabgenlist> then
      Append(~stabgenlist,x);
      if #sub<K | stabgenlist> eq (#K/3) then
        done := true;
      end if;
    end if;
  end if;
end while;
stabgenlist;

// We see that the stabilizer of v is contained in Gamma_0(54).