// Converted from subsystems.tex by MagmaCode 1.2 on Wed Jan 11 10:10:08 2012
//
// Reflection subgroups of complex reflection groups
// Author: Don Taylor
// Created: 30 April 2008
// Revised: 22 December 2011
// Revised: 10 January 2012

vector_reps := function(S)
  T := {}; // representatives
  U := {}; // subspaces
  for v in S do
    X := sub<Parent(v)|v>;
    if X notin U then Include(~T,v); Include(~U,X); end if;
  end for;
  return T;
end function;

orbit_reps := function(G,T)
  reps := [];
  while #T gt 0 do
    t := Rep(T);
    Append(~reps,t);
    S := t^G;
    T := { x : x in T | x notin S };
  end while;
  return reps;
end function;

components := function(R)
  seq := [];
  X := R;
  while not IsEmpty(X) do
    r := X[1];
    Exclude(~X,r);
    C := [r];
    ndx := 0;
    while ndx lt #C do
      ndx +:= 1;
      a := C[ndx];
      extn := [ x : x in X | a*x ne x*a ];
      C cat:= extn;
      X := [ x : x in X | x notin extn ];
    end while;
    Append(~seq,C);
  end while;
  return seq;
end function;

extendGrp := function(W,refreps, H)
  N := Normaliser(W,H);
  X := [ r : r in refreps | r notin H ];
  orbreps := orbit_reps(N,X);

  return [ G : r in orbreps |
        not exists{ E : E in Self() | IsConjugate(W,E,G) }
                where G is sub< W | H, r > ];
end function;

reduceList := func< W,mainlst,extn |
  [ G : G in extn | not exists{ M : M in mainlst | IsConjugate(W,M,G) } ] >;

fix := func< G | &meet[Eigenspace(r,1) : r in Generators(G)] >;

rank := func< G | Dimension(G) - Dimension(Fix(GModule(G))) >;

parabolicClosure := function(G,K)
  F := fix(K);
  R := &join[ Class(G,r) : r in Generators(G)];
  T := {G| r : r in R | r notin K and F subset Eigenspace(r,1)};
  L := K;
  while #T gt 0 do
    x := Rep(T);
    L := sub<G | L, x >;
    T := { t : t in T | t notin L };
  end while;
  return L;
end function;

isParabolic := func< G,P | P eq parabolicClosure(G,P) >;

 groupName := AssociativeArray(CartesianProduct(Integers(),Integers()));
 groupName[<2,1>] := "A1";           // ShephardTodd(1,1,2)
 groupName[<3,1>] := "L1";           // ShephardTodd(3,1,1)
 groupName[<4,1>] := "Z4";           // ShephardTodd(4,1,1)
 groupName[<5,1>] := "Z5";           // ShephardTodd(5,1,1)
 groupName[<6,3>] := "A2";           // ShephardTodd(1,1,3)
 groupName[<8,4>] := "B2";           // ShephardTodd(2,1,2)
 groupName[<10,5>] := "D2(5)";       // ShephardTodd(5,5,2)
 groupName[<12,6>] := "D2(6)";       // ShephardTodd(6,6,2)
 groupName[<16,6>] := "B2(4)";       // ShephardTodd(4,2,2)
 groupName[<16,8>] := "G(8,8,2)";    // ShephardTodd(8,8,2)
 groupName[<18,5>] := "B2(3)";       // ShephardTodd(3,1,2)
 groupName[<20,10>] := "G(10,10,2)"; // ShephardTodd(10,10,2)
 groupName[<24,4>] := "L2";          // ShephardTodd(4)
 groupName[<24,6>] := "A3";          // ShephardTodd(1,1,4)
 groupName[<32,6>] := "G(4,1,2)";    // ShephardTodd(4,1,2)
 groupName[<32,10>] := "G(8,4,2)";   // ShephardTodd(8,4,2)
 groupName[<36,8>] := "G(6,2,2)";    // ShephardTodd(6,2,2)
 groupName[<48,9>] := "B3";          // ShephardTodd(2,1,3)
 groupName[<48,10>] := "G6";         // ShephardTodd(6)
 groupName[<48,12>] := "G12";        // ShephardTodd(12)
 groupName[<50,7>] := "G(5,1,2)";    // ShephardTodd(5,1,2)
 groupName[<54,9>] := "D3(3)";       // ShephardTodd(3,3,3)
 groupName[<64,10>] := "G(8,2,2)";   // ShephardTodd(8,2,2)
 groupName[<72,8>] := "G5";          // ShephardTodd(5)
 groupName[<96,6>] := "G8";          // ShephardTodd(8)
 groupName[<96,12>] := "D3(4)";      // ShephardTodd(4,4,3)
 groupName[<96,18>] := "G13";        // ShephardTodd(13)
 groupName[<100,12>] := "G(10,2,2)"; // ShephardTodd(10,2,2)
 groupName[<120,10>] := "A4";        // ShephardTodd(1,1,5)
 groupName[<120,15>] := "H3";        // ShephardTodd(23)
 groupName[<144,14>] := "G7";        // ShephardTodd(7)
 groupName[<144,20>] := "G14";       // ShephardTodd(14)
 groupName[<162,12>] := "B3(3)";     // ShephardTodd(3,1,3)
 groupName[<192,12>] := "D4";        // ShephardTodd(2,2,4)
 groupName[<192,15>] := "B3(4)";     // ShephardTodd(4,2,3)
 groupName[<192,18>] := "G9";        // ShephardTodd(9)
 groupName[<240,30>] := "G22";       // ShephardTodd(22)
 groupName[<288,14>] := "G10";       // ShephardTodd(10)
 groupName[<288,26>] := "G15";       // ShephardTodd(15)
 groupName[<336,21>] := "J3(4)";     // ShephardTodd(24)
 groupName[<360,20>] := "G20";       // ShephardTodd(20)
 groupName[<384,16>] := "B4";        // ShephardTodd(2,1,4)
 groupName[<576,26>] := "G11";       // ShephardTodd(11)
 groupName[<600,12>] := "G16";       // ShephardTodd(16)
 groupName[<648,12>] := "L3";        // ShephardTodd(25)
 groupName[<648,18>] := "D4(3)";     // ShephardTodd(3,3,4)
 groupName[<720,15>] := "A5";        // ShephardTodd(1,1,6)
 groupName[<720,50>] := "G21";       // ShephardTodd(21)
 groupName[<1152,24>] := "F4";       // ShephardTodd(28)
 groupName[<1200,42>] := "G17";      // ShephardTodd(17)
 groupName[<1296,21>] := "M3";       // ShephardTodd(26)
 groupName[<1536,24>] := "D4(4)";    // ShephardTodd(4,4,4)
 groupName[<1800,32>] := "G18";      // ShephardTodd(18)
 groupName[<1920,20>] := "D5";       // ShephardTodd(2,2,5)
 groupName[<2160,45>] := "J3(5)";    // ShephardTodd(27)
 groupName[<3072,28>] := "B4(4)";    // ShephardTodd(4,2,4)
 groupName[<3600,62>] := "G19";      // ShephardTodd(19)
 groupName[<5040,21>] := "A6";       // ShephardTodd(1,1,7)
 groupName[<7680,40>] := "N4";       // ShephardTodd(29)
 groupName[<9720,30>] := "D5(3)";    // ShephardTodd(3,3,5)
 groupName[<14400,60>] := "H4";      // ShephardTodd(30)
 groupName[<23040,30>] := "D6";      // ShephardTodd(2,2,6)
 groupName[<40320,28>] := "A7";      // ShephardTodd(1,1,8)
 groupName[<46080,60>] := "O4";      // ShephardTodd(31)
 groupName[<51840,36>] := "E6";      // ShephardTodd(35)
 groupName[<51840,45>] := "K5";      // ShephardTodd(33)
 groupName[<155520,40>] := "L4";     // ShephardTodd(32)
 groupName[<174960,45>] := "D6(3)";  // ShephardTodd(3,3,6)
 groupName[<322560,42>] := "D7";     // ShephardTodd(2,2,7)
 groupName[<362880,36>] := "A8";     // ShephardTodd(1,1,9)
 groupName[<2903040,63>] := "E7";    // ShephardTodd(36)
 groupName[<5160960,56>] := "D8";    // ShephardTodd(2,2,8)
 groupName[<39191040,126>] := "K6";  // ShephardTodd(34)
 groupName[<696729600,120>] := "E8"; // ShephardTodd(37)

name := func< n, r |
  IsDefined(groupName,<n,r>) select groupName[<n,r>]
    else "<"*IntegerToString(n)*"|"*IntegerToString(r)*">" >;

setup := function(n)
  roots, coroots, rho, W, J := ComplexRootDatum(n);
  K := BaseRing(J);

  V := VectorSpace(K,Nrows(J),J);
  roots := ChangeUniverse(roots,V);
  coroots := ChangeUniverse(coroots,V);
  rho := map< roots->coroots | a:->rho(a) >;
  Phi := vector_reps(roots);
  R := {@ W!PseudoReflection(a,rho(a)) : a in Phi @};
  R join:= {@ r^2 : r in R | Order(r) eq 4 @};
  return W, R;
end function;

getTag := function(W,refreps,H)
  sform := [];
  R := { r : r in refreps | r in H };
  R := Setseq( R diff { r^2 : r in R | r^2 in R } );
  for C in components(R) do
    if Generators(W) subset C then
      sform := [name(#W,#C)]; break;
    else
      Append(~sform, name(Order(sub<W|C>),#C));
    end if;
  end for;
  return &cat Sort(sform);
end function;

sName := function(W,refreps,refgroup,H)
  tag := getTag(W,refreps,H);
  error if tag notin Keys(refgroup), "sName: group not in list";
  base := tag;
  ndx := 1;
  while not IsConjugate(W,refgroup[tag],H) do
    ndx +:= 1;
    tag := base*"-"*IntegerToString(ndx);
  end while;
  return tag;
end function;

newName := function(W,refreps,refgroup,H)
  tag := getTag(W,refreps,H);
  base := tag;
  ndx := 1;
  while tag in Keys(refgroup) do
    ndx +:= 1;
    tag := base*"-"*IntegerToString(ndx);
  end while;
  return tag;
end function;

rankOne := function(W,refreps)
  refgroup := AssociativeArray(Parent(""));

  rank1 := [ sub<W|r> : r in orbit_reps(W,refreps) ];

  names := [];
  for H in rank1 do
    tag := newName(W,refreps,refgroup,H);
    refgroup[tag] := H;
    Append(~names,tag);
  end for;
  return names, refgroup;
end function;

extension := AssociativeArray(Parent(""));

simpleExtensions := procedure(~subtypes,~extension,~refgroup,W,refreps)
  ndx := 0;
  subgroups := [refgroup[tag] : tag in subtypes];
  while ndx lt #subtypes do
    ndx +:= 1;
    tag := subtypes[ndx];
    print "Extending:",tag;
    H := refgroup[tag];
    extn := extendGrp(W,refreps,H);
    extn2 := reduceList(W,subgroups,extn);
    if #extn2 gt 0 then
      for H in extn2 do
        Append(~subgroups,H);
        etag := newName(W,refreps,refgroup,H);
        refgroup[etag] := H;
        Append(~subtypes,etag);
      end for;
    end if;
    extension[tag] := [ sName(W,refreps,refgroup,H) : H in extn ];
    print " ",extension[tag];
  end while;
end procedure;

testext := procedure(names,refgroup,extension,W)
  for tag in names do
    H := refgroup[tag];
    if not isParabolic(W,H) then
      m := rank(H);
      for extn in extension[tag] do
        K := refgroup[extn];
        if rank(K) gt m then
          flag := isParabolic(W,K);
          print tag, extn, flag;
          if flag then print "~~~COUNTER EXAMPLE~~~~~~~~~"; end if;
        end if;
      end for;
    end if;
  end for;
end procedure;

extends := function(A, B, refgroup, extension, W)
  flag, G := IsDefined(refgroup,B);
  if flag then
    flag, E := IsDefined(extension,A);
    if flag then
       return exists(X){ X : X in E | IsConjugate(W,X,G)};
    else
      return "Extensions not defined";
    end if;
  else
    return "Identity not known";
  end if;
end function;

getList := function(n)
  error if n notin [4..37], "n must be in the range 4..37";
  W, refreps := setup(n);
  names, refgroup := rankOne(W,refreps);
  extension := AssociativeArray(Parent(""));
  simpleExtensions(~names,~extension,~refgroup,W,refreps);
  return names, refgroup, extension, W,refreps;
end function;

printTable := procedure(n)
  names, refgroup, extension, W, _ := getList(n);
  for X in names do
    H := refgroup[X];
    para := isParabolic(W,H) select "P" else "N";
    print para, "|", X, "|", extension[X];
  end for;
end procedure;