// This magma script is for enumerating subgroups GL_2(Z_2)
// Now we're working on those that don't contain -I.
// Each group G is stored as a list with five pieces of information:
// i) The number of the group <G,-I>, ii) The index of G, iii) the level of G
// iv) The group G itself, and v) The name of X_G. 

/* Pre-computations */

nbound := 7;

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

load "gl2data.txt";
// The following is the list of curves that have non-CM rational points.
testlist := [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,22,
23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,
48,49,50,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,73,74,75,76,77,78,
79,80,81,82,83,84,85,86,87,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,
104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,
123,124,150,153,155,156,165,166,167,177,178,181,183,185,187,188,189,190,191,
192,193,194,195,196,197,199,200,202,203,204,205,207,208,209,210,211,212,213,
214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,
233,234,235,236,237,238,239,240,241,242,243,281,284,288,289,291,295,297,302,
304,308,309,312,318,320,323,324,326,328,349,350,441,556,558,563,566,619,649];

finecount := [];
finesublist := [];
for n in testlist do
  printf "Looking for good subgroups of G_%o.\n",n;
  lev := Valuation(newsublist[n][2],2);
  G := newsublist[n][3];
  maxsubgplev := Max([3,lev+1]);
  G2 := liftsub(G,lev,maxsubgplev);
  Mlist := MaximalSubgroups(G2);
  sublist := [];
  for m in [1..#Mlist] do
    newlev, sub := reduce(Mlist[m]`subgroup,maxsubgplev);
    if (newlev ne lev) then
      G3 := liftsub(G,lev,newlev);
    else
      G3 := G;
    end if;
    minusone := G3![-1,0,0,-1];
    SLL := SL(2,Integers(2^newlev));
    if Index(G3,sub) eq 2 then
      if not (minusone in sub) then
        if Index(sub,sub meet SLL) eq 2^(newlev-1) then
          found := false;
          C := ConjugacyClasses(sub);
          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 true) then
            isnew := true;
            if #sublist gt 0 then
              for j in [1..#sublist] do
                largerlev := Max(sublist[j][2],newlev);
                newgp := sub;
                oldgp := sublist[j][1];
                if sublist[j][2] lt largerlev then
                  oldgp := liftsub(oldgp,sublist[j][2],largerlev);
                end if;
                if newlev lt largerlev then
                  newgp := liftsub(newgp,newlev,largerlev);
                end if;
                bol := IsConjugate(GL(2,Integers(2^largerlev)),oldgp,newgp);
                if (bol eq true) then
                  isnew := false;
                end if;
              end for;
            end if;
            if (isnew eq true) then
              printf "Maximal subgroup %o works and is new.\n",m;
              Append(~sublist,<sub,newlev>);
              Append(~finesublist,<n,2*newsublist[n][1],2^newlev,sub>);
            end if;
          end if;
        end if;
      end if;
    end if;
  end for;
  Append(~finecount,<n,#sublist>);
end for;

// Sort the data.

subinddata := [ <finesublist[i][1], finesublist[i][3], #ConjugacyClasses(finesublist[i][4])> : i in [1..#finesublist]];
symgp := Sym(#subinddata);
p := symgp!1;
Sort(~subinddata,~p);
rebuildfinesublist := [];
for n in [1..#finesublist] do
  covergp := finesublist[n][1];
  index := finesublist[n][2];
  level := finesublist[n][3];
  gp := finesublist[n][4];
  name := "X" cat IntegerToString(covergp);
  if (n eq 1) then
    cnt := 1;
  end if;
  if (n gt 1) then
    if (rebuildfinesublist[n-1][1] ne covergp) then
      cnt := 1;
    else
      cnt := Min([ j : j in [1..#rebuildfinesublist] | rebuildfinesublist[j][1] eq covergp]);
      cnt := (n - cnt) + 1;
    end if;
  end if;
  name := name cat CodeToString(cnt + 96);
  Append(~rebuildfinesublist,<covergp,index,level,gp,name>);
end for;
finesublist := rebuildfinesublist;

// Build lattice data

superlist := [ [] : i in [1..#finesublist]];
for n in [1..#finesublist] do
  printf "Finding supergroups for %o.",finesublist[n][5];
  ind := finesublist[n][2];
  lev := finesublist[n][3];
  gp := finesublist[n][4];
  conjlist := SetToSequence(Conjugates(GL(2,Integers(lev)),gp));
  for j in [1..n-1] do
    if (ind mod finesublist[j][2] eq 0) and (ind gt finesublist[j][2]) then
      if (lev ge finesublist[j][3]) then
        chkgp := finesublist[j][4];
        if (finesublist[j][3] lt lev) then
          chkgp := liftsub(finesublist[j][4],Valuation(finesublist[j][3],2),Valuation(lev,2));
        end if;
        found := false;
        for k in [1..#conjlist] do
          if conjlist[k] subset chkgp then
            found := true;
          end if;
        end for;
        if (found eq true) then
          Append(~superlist[n],j);
          printf "The group %o is one.\n",finesublist[j][5];
        end if;
      end if;
    end if;
  end for;
  printf "There were %o.\n",#superlist[n];
end for;
finesublistminsuper := [ [] : i in [1..#finesublist]];
for n in [1..#finesublist] do
  supers := superlist[n];
  nonminsupers := {};
  for j in supers do
    for k in superlist[j] do
      Include(~nonminsupers,k);
    end for;
  end for;
  finesublistminsuper[n] := Sort([ m : m in supers | not (m in nonminsupers)]);
end for;
finesublistmaxsub := [ [] : i in [1..#finesublist]];
for n in [1..#finesublist] do
  for j in finesublistminsuper[n] do
    Append(~finesublistmaxsub[j],n);
  end for;
end for;
for n in [1..#finesublist] do
  finesublistmaxsub[n] := Sort(finesublistmaxsub[n]);
end for;

datafile := "gl2finedata.txt";
System("rm " cat datafile);
printf "Writing data to gl2finedata.txt\n";
PrintFile(datafile,"finesublist:=");
PrintFileMagma(datafile,finesublist);
PrintFile(datafile,";");
PrintFile(datafile,Sprintf("finesublistmaxsub:=%o;\n",finesublistmaxsub));
PrintFile(datafile,Sprintf("finesublistminsuper:=%o;\n",finesublistminsuper));