// This magma script is for enumerating subgroups GL_2(Z_2)
// we'll start with those that contain -I.

genbound := 1;

// Each subgroup G is stored in a cartesian product.
// Each entry in the catersian product has the form
// (i) The index
// (ii) The level 2^n
// * WARNING! To trick Magma, into working, we declare that the level of
// GL_2(Z_2) itself is equal to 2.
// (iii) The image of that subgroup in GL_2(Z/2^n Z)
// (iv) The maximal subgroups of G, given as pairs
// <M,n> where 2^n is the level of M. 
// (v) The genus of the corresponding modular curve
// (vi) The image of the subgroup in a permutation group
// (vii) The subgroup in the list that covers this subgroup with minimal degree.
// We arrange the information here so that the group is actually a subgroup of that one.

/* Pre-computations */

nbound := 8;

// We start by pre-computing
// GL(2,Z/2^n Z), permutations groups GL(2,Z/2^n Z)/{ +- 1 }
// and maps phi between them
// we also precompute kernels of the maps from GL(2,Z/2^m Z) ->
// to GL(2,Z/2^n Z) for various m and n

printf "Pre-computing!\n";
matgps := [ GL(2,Integers(2^n)) : n in [1..nbound]];
printf "Done!\n";
permgps := <>;
mathoms := <>;
for n in [1..nbound] do
  printf "Computing permutation representation %o.\n",n;
  G := matgps[n];
  K := sub<G | G![3,0,0,1], G![5,0,0,1], G![1,1,0,1], G![-1,0,0,-1]>;
  phi, B := CosetAction(G,K);
  Append(~permgps,B);
  Append(~mathoms,phi);
end for;
sl2list := < mathoms[i](SL(2,Integers(2^i))) : i in [1..nbound] >;
kerlist := <>;
for n in [1..nbound] do
  tempkerlist := <>;
  for n2 in [n+1..nbound] do
    phi := hom< matgps[n2] -> matgps[n] | [ matgps[n]!matgps[n2].i : i in [1..#Generators(matgps[n2])]]>;
    printf "Storing kernel of GL(2,Z/%oZ) -> GL(2,Z/%oZ).\n",2^n2,2^n;
    printf "Storing in slot %o, %o.\n",n,#tempkerlist+1;
    K := Kernel(phi);
    Append(~tempkerlist,K);
  end for;
  Append(~kerlist,tempkerlist);
end for;
printf "Pre-computations done!\n";

/* HELPER FUNCTIONS */

// This function takes a subgroup of GL_2(Z/2^n Z) and lifts it
// to a subgroup of GL_2(Z/2^m Z) for m > n.

function liftsub(G,n,m)
  H2 := GL(2,Integers(2^m));
  H3 := GL(2,Integers(2^n));
  genlist := [];
  for g in Generators(G) do
    Append(~genlist,H2!g);
  end for;
  // Add generators for the kernel of GL_2(Z/2^m Z) -> GL_2(Z/2^n Z)
  return sub<H2 | genlist, kerlist[n][m-n]>;
end function;

// Write our own version of a maximal subgroups routine
// because it's (perhaps) faster in a permutation group than in a matrix
// group. Only return those up to conjugacy in GL(2,Z/2^m Z), and
// only those with surjective determinant

function maxsub(G,m)
  permG := mathoms[m](G);
  permsubs := MaximalSubgroups(permG);
  addlist := [];
  for mm in [1..#permsubs] do
    found := false;
    done := false;
    // check determinant condition
    if Index(permsubs[mm]`subgroup,permsubs[mm]`subgroup meet sl2list[m]) ne 2^(m-1) then
      done := true;
      found := true;
    end if;
    it := 1;
    while (done eq false) do
      if (it gt #addlist) then
        done := true;
      else
        if IsConjugate(permgps[m],addlist[it],permsubs[mm]`subgroup) then
          found := true;
          done := true;
        end if;
      end if;
      it := it + 1;
    end while;
    if found eq false then
      Append(~addlist,permsubs[mm]`subgroup);
    end if;
  end for;
  retlist := [ addlist[n] @@ mathoms[m] : n in [1..#addlist]];
  // List subgroups in increasing order of index
  return Reverse(retlist);
end function;

// Given a subgroup G of GL_2(Z/2^n Z), determine the level 2^m of G
// and return the corresponding subgroup of GL_2(Z/2^m Z)

function reduce(G,n);
  lev := n;
  for m in [1..n-1] do
    if (kerlist[m][n-m] subset G) then
      lev := Min(lev,m);
    end if;
  end for;
  newG := sub<matgps[lev] | [matgps[lev]!G.i : i in [1..#Generators(G)]]>;
  return lev, newG;
end function;

// Given a subgroup G of GL_2(Z/2^n Z) determine
// the genus of the corresponding modular curve.
// This is more or less a copy of Cummins and Pauli's
// genus function and uses the Hurwitz formula

function genus(G,n)
  H := SL(2,Integers(2^n));
  S := H![0,-1,1,0];
  T := H![1,1,0,1];
  phi, perm := CosetAction(H,G meet H);
  lst := [phi(S),phi(T),phi(S*T)];
  cs := [CycleStructure(lst[i]) : i in [1..3]];
  ans := -2*Degree(perm) + 2;
  for j in [1..3] do
    for i in [1..#cs[j]] do
      ans := ans + (cs[j][i][1]-1)*cs[j][i][2];
    end for;
  end for;
  ans := ans div 2;
  return ans;
end function;

// Given a subgroup G of GL_2(Z/2^n Z) determine its genus and
// if it is conjugate to a subgroup we already have. If not,
// add it to sublist. The function returns whether or not
// it was added.

function addgroup(G,n,sublist);
  didweadd := false;
  // If we get here, G is good and the level of G is 2^n.
  // printf "Computing genus. It is ";
  gen := genus(G,n);
  //printf "%o.\n",gen;
  //printf "Computing index.\n";
  ind := Index(matgps[n],G);
  found := false;
  done := false;
  it := 1;
  //printf "Checking conjugacy.\n";
  permG := mathoms[n](G);
  while (done eq false) do
    if (it gt #sublist) then
      done := true;
      didweadd := true;
    else
      if (sublist[it][1] eq ind) and (sublist[it][2] eq 2^n) and (sublist[it][5] eq gen) then
        if IsConjugate(permgps[n],permG,sublist[it][6]) then
          //printf "Group not added! It is a conjugate of one we have!\n";
          found := true;
          done := true;
        end if;
      end if;
      if (found eq false) then
        it := it + 1;
      end if;
    end if;
  end while;
  adddata := <>;
  adddata4 := <>;
  // Check to see if there is an element in G with characteristic polynomial
  // x^2 - 1. In order for there to be Q-rational non-cuspidal points
  // such a thing must exist.
  if (didweadd eq true) then
    found := false;
    C := ConjugacyClasses(G);
    for m in [1..#C] do
      if (Trace(C[m][3]) eq 0) and (Determinant(C[m][3]) eq -1) then
        found := true;
      end if;
    end for;
    if found eq false then
      didweadd := false;
    end if;
  end if;
  if (didweadd eq true) then
    m2 := Max([3,n+1]);
    G2 := liftsub(G,n,m2);
    Mlist := maxsub(G2,m2);
    //Mlist := MaximalSubgroups(G2);
    //printf "There were %o.\n",#Mlist;
    numgood := 0;
    //printf "Found %o maximal subgroups!\n",#Mlist;
    for m in [1..#Mlist] do
      //printf "Working on maximal subgroup %o...",m;
      //M := Mlist[m]`subgroup;
      M := Mlist[m];
      newlev, newM := reduce(M,m2);
      Append(~adddata4,<newM,newlev>);
    end for;
    printf "Found %o maximal subgroups.\n",#Mlist;
    adddata := <ind,2^n,G,adddata4,gen,mathoms[n](G),0>;
  end if;
  return didweadd, adddata;
end function;

// The main loop

sublist := <>;
start := GL(2,Integers(2));
added, newdata := addgroup(start,1,sublist);
if added eq true then
  Append(~sublist,newdata);
  printf "We added a group! We have #sublist = %o.\n",#sublist;
end if;
curind := 1;
reallydone := false;
while reallydone eq false do
  printf "Adding maximal subgroups of group %o.\n",curind;
  genn := sublist[curind][5];
  if genn le genbound then
    // If the genus is less than or equal to 1, work with the subgroups
    looplen := #sublist[curind][4];
    for m in [1..looplen] do
      printf "Working on subgroup %o of %o.\n",m,looplen;
      added, newdata := addgroup(sublist[curind][4][m][1],sublist[curind][4][m][2],sublist);
      if added eq true then
        Append(~sublist,newdata);
        printf "We added a group! We have #sublist = %o.\n",#sublist;
      end if;
    end for;
  end if;
  if (curind lt #sublist) then
    curind := curind + 1;
  else
    reallydone := true;
  end if;
end while;

// Determine which groups on the list are maximal subgroups of other groups
// and determine which groups are contained in which maximal subgroups

sublistmaxsub := [ [] : i in [1..#sublist]];
sublistminsuper := [ [] : i in [1..#sublist]];

printf "Building subgroup lattice. This takes a while.\n";
for n in [1..#sublist] do
  printf "Working on n = %o.\n",n;
  for m in [1..#sublist[n][4]] do
    for j in [1..#sublist] do
      // Is the mth maximal subgroup of sublist[n] conjugate to
      // sublist[j]?
      if Index(sublistmaxsub[n],j) eq 0 then
        if sublist[j][2] eq (2^sublist[n][4][m][2]) then
          nn := sublist[n][4][m][2];
          if #sublist[j][3] eq #sublist[n][4][m][1] then
            if IsConjugate(permgps[nn],sublist[j][6],mathoms[nn](sublist[n][4][m][1])) then
              Append(~sublistmaxsub[n],j);
            end if;
          end if;
        end if;
      end if;
    end for;
  end for;
  Sort(~sublistmaxsub[n]);
end for;
printf "Finding minimal supergroup inclusions.\n";
for n in [1..#sublist] do
  for m in [1..#sublistmaxsub[n]] do
    ind := sublistmaxsub[n][m];
    Append(~sublistminsuper[ind],n);
  end for;
end for;

// Remove any subgroups that are covered by a subgroup of genus >= 2.

printf "Removing groups contained in a group of genus >= 2.\n";
shouldinclude := [];
for n in [1..#sublist] do
  genlist := [ sublist[m][5] : m in sublistminsuper[n]];
  if (#genlist gt 0) and (Max(genlist) ge 2) then
    Append(~shouldinclude,0);
  else
    Append(~shouldinclude,1);
  end if;
end for;

// Subgroup lattice tower has height hopefully less than or equal to 16
for it in [1..16] do
  for n in [2..#sublist] do
    if Min({ shouldinclude[m] : m in sublistminsuper[n]}) eq 0 then
      shouldinclude[n] := 0;
    end if;
  end for;
end for;

printf "Rewriting subgroup data.\n";
processlist := [];
newsublist := <>;
found := 0;
for n in [1..#sublist] do
  if shouldinclude[n] eq 0 then
    Append(~processlist,0);
  else
    found := found + 1;
    Append(~processlist,found);
    Append(~newsublist,sublist[n]);
  end if;
end for;

newsublistmaxsub := [];
for n in [1..#sublist] do
  if shouldinclude[n] eq 1 then
    maxsubset := {};
    done := false;
    queue := [ sublistmaxsub[n][i] : i in [1..#sublistmaxsub[n]] ];
    while #queue gt 0 do
      node := queue[1];
      Remove(~queue,1);
      if shouldinclude[node] eq 1 then
        // A maximal subgroup is still a maximal subgroup
        Include(~maxsubset,processlist[node]);
      else
        // We need to check maximal subgroups of sublist[node] now!
        for m in [1..#sublistmaxsub[node]] do
          Append(~queue,sublistmaxsub[node][m]);
        end for;
      end if;
    end while;
    Append(~newsublistmaxsub,Sort(SetToSequence(maxsubset)));
  end if;
end for;
newsublistminsuper := [ [] : i in [1..#newsublist]];
for n in [1..#newsublist] do
  for m in [1..#newsublistmaxsub[n]] do
    ind := newsublistmaxsub[n][m];
    Append(~newsublistminsuper[ind],n);
  end for;
end for;

// Sort the info by the index, then the genus, then level

printf "Sorting by index, genus, and level.\n";
newsubinddata := [ <newsublist[i][1], newsublist[i][5], newsublist[i][2]> :
i in [1..#newsublist]];
G := Sym(#newsublist);
p := G!1;
Sort(~newsubinddata,~p);
pinv := p^(-1);
rebuildnewsublist := <>;
rebuildmaxsub := [ [] : i in [1..#newsublist]];
rebuildminsuper := [ [] : i in [1..#newsublist]];
for n in [1..#newsublist] do
  Append(~rebuildnewsublist,newsublist[Image(p,n)]);
  for m in [1..#newsublistmaxsub[Image(p,n)]] do
    Append(~rebuildmaxsub[n],Image(pinv,newsublistmaxsub[Image(p,n)][m]));
  end for;
  for m in [1..#newsublistminsuper[Image(p,n)]] do
    Append(~rebuildminsuper[n],Image(pinv,newsublistminsuper[Image(p,n)][m]));
  end for;
  Sort(~rebuildmaxsub[n]);
  Sort(~rebuildminsuper[n]);
end for;
newsublist := rebuildnewsublist;
newsublistmaxsub := rebuildmaxsub;
newsublistminsuper := rebuildminsuper;

// Find the smallest degree proper genus zero or one cover and note that
// Then, arrange newsublist[n] so that it is a subgroup of that one.

for n in [2..#newsublist] do
  // find cover
  printf "Finding optimal cover of group %o.\n",n;
  queue := newsublistminsuper[n];
  allsuper := {};
  done := false;
  while done eq false do
  if #queue gt 0 then
    cur := queue[1];
    Remove(~queue,1);
    Include(~allsuper,cur);
    for m in [1..#newsublistminsuper[cur]] do
      toadd := newsublistminsuper[cur][m];
      if (not (toadd in allsuper)) and Index(queue,toadd) eq 0 then
        Append(~queue,toadd);
      end if;
    end for;
  else
    done := true;
  end if;
  end while;
  allsuper := Sort(SetToSequence(allsuper));
  minind := Max([ newsublist[m][1] : m in allsuper ]);
  allsuperminind := [ nn : nn in allsuper | newsublist[nn][1] eq minind ];
  allsupergensmall := [ nn : nn in allsuperminind | newsublist[nn][5] lt 2 ];
  cover := Min(allsupergensmall);
  newsublist[n][7] := cover;
  printf "The cover is group %o.\n",cover;
  // Find a matrix g so that g^(-1)*newsublist[n][3]*g is a subset of newsublist[cover][3];
  lev := newsublist[n][2];
  sub := newsublist[n][3];
  G := GL(2,Integers(lev));
  NN := Normalizer(G,newsublist[n][3]);
  BB := RightTransversal(G,NN);
  if newsublist[cover][2] ne newsublist[n][2] then
    covergp := liftsub(newsublist[cover][3],Valuation(newsublist[cover][2],2),Valuation(lev,2));
  else
    covergp := newsublist[cover][3];
  end if;
  found := false;
  for g in BB do
    if found eq false then
      if &and[ g^(-1)*gen*g in covergp : gen in Generators(sub)] then
        goodg := g;
        found := true;
      end if;
    end if;
  end for;
  if (found eq false) then
    printf "Error! We weren't able to find a conjugate of group %o inside group %o.\n",n,cover;
  end if;
  printf "Rewriting subgroup data so that the subgroup is contained in group %o.\n",cover;
  newsub := sub<G | [ goodg^(-1)*gen*goodg : gen in Generators(sub)]>;
  maxsubdata := <>;
  for m in [1..#newsublist[n][4]] do
    maxlev := newsublist[n][4][m][2];
    G := GL(2,Integers(2^maxlev));
    goodg2 := G!goodg;
    newmaxsub := sub<G | [goodg2^(-1)*gen*goodg2 : gen in Generators(newsublist[n][4][m][1])] >;
    Append(~maxsubdata,<newmaxsub,maxlev>);
  end for;
  newsublist[n][3] := newsub;
  newsublist[n][4] := maxsubdata;
  val := Valuation(newsublist[n][2],2);
  newsublist[n][6] := mathoms[val](newsub);
end for;

// Write data to a file

datafile := "gl2data.txt";
System("rm " cat datafile);
printf "Writing data to gl2data.txt.\n";
PrintFile(datafile,"newsublist:=");
PrintFileMagma(datafile,newsublist);
PrintFile(datafile,";");
PrintFile(datafile,Sprintf("newsublistmaxsub:=%o;\n",newsublistmaxsub));
PrintFile(datafile,Sprintf("newsublistminsuper:=%o;\n",newsublistminsuper));