// Extracted from refSub.tex by MagmaCode 2.2 on Wed Feb  5 12:25:09 2025
//
// Reflection subgroups of complex reflection groups
// Author: Don Taylor
// Created: 15 January 2025
// Replaces: subsystems.m (from 30 April 2008)

 toStr := IntegerToString;
 Z := Integers();
 groupName := AssociativeArray(CartesianProduct(<Z,Z,Z>));
 groupName[<2,1,1>] := "A1";             // ShephardTodd(2,1,1)
 groupName[<3,0,2>] := "Z3";             // ShephardTodd(3,1,1)
 groupName[<4,1,3>] := "Z4";             // ShephardTodd(4,1,1)
 groupName[<5,0,4>] := "Z5";             // ShephardTodd(5,1,1)
 groupName[<6,3,1>] := "A2";             // ShephardTodd(1,1,3)
 groupName[<8,4,2>] := "B2";             // ShephardTodd(2,1,2)
 groupName[<10,5,1>] := "D2(5)";         // ShephardTodd(5,5,2)
 groupName[<12,6,2>] := "D2(6)";         // ShephardTodd(6,6,2)
 groupName[<16,6,3>] := "B2(4)";         // ShephardTodd(4,2,2)
 groupName[<16,8,2>] := "D2(8)";         // ShephardTodd(8,8,2)
 groupName[<18,3,3>] := "B2(3)";         // ShephardTodd(3,1,2)
 groupName[<20,10,2>] := "D2(10)";       // ShephardTodd(10,10,2)
 groupName[<24,0,2>] := "G4";            // ShephardTodd(4) = SL(2,3)
 groupName[<24,6,1>] := "A3";            // ShephardTodd(1,1,4)
 groupName[<24,8,2>] := "B2(6)";         // ShephardTodd(6,3,2)
 groupName[<32,6,4>] := "G(4,1,2)";      // ShephardTodd(4,1,2)
 groupName[<32,10,3>] := "B2(8)";        // ShephardTodd(8,4,2)
 groupName[<36,6,4>] := "G(6,2,2)";      // ShephardTodd(6,2,2)
 groupName[<48,9,2>] := "B3";            // ShephardTodd(2,1,3)
 groupName[<48,6,3>] := "G6";            // ShephardTodd(6)
 groupName[<48,12,1>] := "G12";          // ShephardTodd(12)
 groupName[<50,5,5>] := "G(5,1,2)";      // ShephardTodd(5,1,2)
 groupName[<54,9,1>] := "D3(3)";         // ShephardTodd(3,3,3)
 groupName[<64,10,5>] := "G(8,2,2)";     // ShephardTodd(8,2,2)
 groupName[<72,0,4>] := "G5";            // ShephardTodd(5)
 groupName[<96,6,3>] := "G8";            // ShephardTodd(8)
 groupName[<96,12,1>] := "D3(4)";        // ShephardTodd(4,4,3)
 groupName[<96,18,2>] := "G13";          // ShephardTodd(13)
 groupName[<100,10,6>] := "G(10,2,2)";   // ShephardTodd(10,2,2)
 groupName[<120,10,1>] := "A4";          // ShephardTodd(1,1,5)
 groupName[<120,15,1>] := "G23";         // ShephardTodd(23)
 groupName[<144,6,5>] := "G7";           // ShephardTodd(7)
 groupName[<144,12,3>] := "G14";         // ShephardTodd(14)
 groupName[<150,15,1>] := "D3(5)";       // ShephardTodd(5,5,3)
 groupName[<162,9,3>] := "B3(3)";        // ShephardTodd(3,1,3)
 groupName[<192,12,1>] := "D4";          // ShephardTodd(2,2,4)
 groupName[<192,15,2>] := "B3(4)";       // ShephardTodd(4,2,3)
 groupName[<192,18,4>] := "G9";          // ShephardTodd(9)
 groupName[<240,30,1>] := "G22";         // ShephardTodd(22)
 groupName[<288,6,5>] := "G10";          // ShephardTodd(10)
 groupName[<288,18,4>] := "G15";         // ShephardTodd(15)
 groupName[<336,21,1>] := "G24";         // ShephardTodd(24)
 groupName[<360,0,2>] := "G20";          // ShephardTodd(20)
 groupName[<384,16,2>] := "B4";          // ShephardTodd(2,1,4)
 groupName[<576,18,6>] := "G11";         // ShephardTodd(11)
 groupName[<600,0,4>] := "G16";          // ShephardTodd(16)
 groupName[<648,0,2>] := "G25";          // ShephardTodd(25)
 groupName[<648,18,1>] := "D4(3)";       // ShephardTodd(3,3,4)
 groupName[<720,15,1>] := "A5";          // ShephardTodd(1,1,6)
 groupName[<720,30,3>] := "G21";         // ShephardTodd(21)
 groupName[<1152,24,2>] := "G28";        // ShephardTodd(28)
 groupName[<1200,30,5>] := "G17";        // ShephardTodd(17)
 groupName[<1296,9,3>] := "G26";         // ShephardTodd(26)
 groupName[<1536,24,1>] := "D4(4)";      // ShephardTodd(4,4,4)
 groupName[<1800,0,6>] := "G18";         // ShephardTodd(18)
 groupName[<1920,20,1>] := "D5";         // ShephardTodd(2,2,5)
 groupName[<2160,45,1>] := "G27";        // ShephardTodd(27)
 groupName[<3072,28,2>] := "B4(4)";      // ShephardTodd(4,2,4)
 groupName[<3600,30,7>] := "G19";        // ShephardTodd(19)
 groupName[<5040,21,1>] := "A6";         // ShephardTodd(1,1,7)
 groupName[<7680,40,1>] := "G29";        // ShephardTodd(29)
 groupName[<9720,30,1>] := "D5(3)";      // ShephardTodd(3,3,5)
 groupName[<14400,60,1>] := "G30";       // ShephardTodd(30)
 groupName[<23040,30,1>] := "D6";        // ShephardTodd(2,2,6)
 groupName[<40320,28,1>] := "A7";        // ShephardTodd(1,1,8)
 groupName[<46080,60,1>] := "G31";       // ShephardTodd(31)
 groupName[<51840,36,1>] := "G35";       // ShephardTodd(35)
 groupName[<51840,45,1>] := "G33";       // ShephardTodd(33)
 groupName[<155520,0,2>] := "G32";       // ShephardTodd(32)
 groupName[<174960,45,1>] := "D6(3)";    // ShephardTodd(3,3,6)
 groupName[<322560,42,1>] := "D7";       // ShephardTodd(2,2,7)
 groupName[<362880,36,1>] := "A8";       // ShephardTodd(1,1,9)
 groupName[<2903040,63,1>] := "G36";     // ShephardTodd(36)
 groupName[<5160960,56,1>] := "D8";      // ShephardTodd(2,2,8)
 groupName[<39191040,126,1>] := "G34";   // ShephardTodd(34)
 groupName[<696729600,120,1>] := "G37";  // ShephardTodd(37)

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

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

components := function(ctx, H)
  refs := [ r : r in &join ctx`refs | r in H ];
  rtGrps := {@ Stabiliser(H,Basis(Eigenspace(r,1))) : r in refs @};
  rtRefs := [ [ r : r in refs | r in H ] : H in rtGrps ];
  seq := [];
  X := [ 1..#rtRefs ];
  while not IsEmpty(X) do
    n := X[1];
    Exclude(~X,n);
    C := rtRefs[n];

    while exists(ndx){ x : c in C, r in rtRefs[x], x in X | c*r ne r*c } do
      C cat:= rtRefs[ndx];
      Exclude(~X,ndx);
    end while;
    Append(~seq,C);
  end while;
  return seq;
end function;

getTag := function(ctx,H)
  W := ctx`group;

  if forall{ r : r in Generators(W) | r in H } then
    refno2 := &+[ Z | c[2] : c in ctx`refrep | c[1] eq 2 ];
    sform := [name(#W,refno2,#ctx`refrep)];
  else
    sform := [];
    for C in components(ctx,H) do
      refno2 := #{ r : r in C | Order(r) eq 2 };
      H := sub<W|C>;
      refno := #class_reps(H,C);
      Append(~sform, name(#H,refno2,refno));
    end for;
  end if;
  return &cat Sort(sform);
end function;

cTag := function(H)
  refreps := [ c : c in Classes(H) | Rank(Matrix(c[3]) - H.0) eq 1 ];
  refno2 := &+[ Z | c[2] : c in refreps | c[1] eq 2 ];
  return "<"*toStr(#H)*"|"*toStr(refno2)*"|"*toStr(#refreps)*">";
end function;

baseName := function(ctx,refgroup,H)
  base := getTag(ctx,H);
  tag := base;
  ndx := 1;
  while tag in Keys(refgroup) do
    ndx +:= 1;
    tag := base*"-"*toStr(ndx);
  end while;
  return base, ndx;
end function;

context := recformat< group : GrpMat, refrep : SeqEnum, refs : SeqEnum >;

prepare := function(W)
  cc := Classes(W);
  refreps := [ c : c in cc | Rank(Matrix(c[3]) - W.0) eq 1 ];
  R := [ Class(W,c[3]) : c in refreps ];
  return rec< context | group := W, refrep := refreps, refs := R >;
end function;

setup := function(ndx)
  if Type(ndx) eq RngIntElt then
    W := ShephardTodd( ndx : NumFld);
  else
    m, p, n := Explode(ndx);
    W := ShephardTodd(m,p,n);
  end if;
  return prepare(W);
end function;

ccCheck := function(W, baseTag, limit, refgroup, H)
  tag := baseTag;
  if tag in Keys(refgroup) then
    ndx := 1;
    while tag in Keys(refgroup) and not IsConjugate(W,refgroup[tag],H) do
      ndx +:= 1;
      tag := baseTag * "-" * toStr(ndx);
    end while;
    if ndx lt limit then return true, tag; end if;
  end if;
  return false, _;
end function;

rankOne := function(ctx)
  refgroup := AssociativeArray(Parent(""));
  names := [];
  for H in [ sub< ctx`group | c[3] > : c in ctx`refrep ] do
    baseTag, n := baseName(ctx, refgroup, H);
    if n eq 1 then
      tag := baseTag;
    else
      found, theTag := ccCheck(ctx`group,baseTag,n,refgroup,H);
      if found then
        continue;
      else
        tag := baseTag * "-" * toStr(n);
      end if;
    end if;
    Append(~names,tag);
    refgroup[tag] := H;
  end for;
  return names, refgroup;
end function;

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

extendGrp := function(ctx, H)
  W := ctx`group;
  N := Normaliser(W,H);
  X := [ r : r in &join ctx`refs | r notin H ];
  ccreps := class_reps(N,X);

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

simpleExtensions := procedure(~names,~extension,~refgroup,ctx : 
  Verbose := false)
  ndx := 0;
  while ndx lt #names do
    enames := {@ @};
    ndx +:= 1;
    tag := names[ndx];
    if Verbose then print "Extending:",tag; end if;
    H0 := refgroup[tag];
    extn := extendGrp(ctx,H0);
    for H in extn do
      baseTag, n := baseName(ctx, refgroup, H);
      if baseTag in names then
        found, theTag := ccCheck(ctx`group,baseTag,n,refgroup,H);
        if not found then
          etag := n eq 1 select baseTag else baseTag * "-" * toStr(n);
          refgroup[etag] := H;
          Include(~enames,etag);
          Append(~names,etag);
          if Verbose then print "  Adding", etag; end if;
        else
          Include(~enames,theTag);
          if Verbose then print "  Skipping", theTag; end if;
        end if;
      else
        refgroup[baseTag] := H;
        Include(~enames,baseTag);
        Append(~names,baseTag);
        if Verbose then print "  Adding", baseTag; end if;         
      end if;
    end for;
    extension[tag] := Setseq(enames);
    if Verbose then print extension[tag]; end if;
  end while;
end procedure;

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

parabolicClosure := function(ctx,K)
  W := ctx`group;
  F := fix(K);
  R := &join ctx`refs;
  T := {W| 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<W | 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) >;

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

testExt := procedure(n)
  ctx := setup(n);
  names, refgroup := rankOne(ctx);
  extension := AssociativeArray(Parent(""));
  simpleExtensions(~names,~extension,~refgroup,ctx);
  for tag in names do
    H := refgroup[tag];
    if not isParabolic(ctx,H) then
      m := rank(H);
      for extn in extension[tag] do
        K := refgroup[extn];
        if rank(K) gt m then
          flag := isParabolic(ctx,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 : Verbose := false, Sorted := false)
  ctx := setup(n);
  names, refgroup := rankOne(ctx);
  extension := AssociativeArray(Parent(""));
  simpleExtensions(~names,~extension,~refgroup,ctx : Verbose := Verbose);
  if Sorted then
    fn := function(A,B)
      H := refgroup[A]; h := Ngens(H); rH := rank(H);
      K := refgroup[B]; k := Ngens(K); rK := rank(K);
      return (rH eq rK) select ((h eq k) select #H - #K else h - k) 
             else rH - rK;
    end function;
    Sort(~names, fn);
  end if;
  return names, refgroup, extension, ctx;
end function;

printTable := procedure(n : Verbose := false, Sorted := false)
  names, refgroup, extension, ctx := getList(n : Verbose := Verbose, 
         Sorted := Sorted);
  maxlen := Max([ #nm : nm in names ]);
  W := ctx`group;
  for X in names do
    H := refgroup[X];
    para := isParabolic(ctx,H) select "P" else " ";
    if Sorted then
      rk := rank(H);
      ngens := Ngens(H);
      nconj := Index(W,Normaliser(W,H));
      h := maxlen - #X;
      pad := " "^h * "|";
      print para, rk, ngens, nconj, "|", X, pad, extension[X];
    else
      print para, "|", X, "|", extension[X];
    end if;
  end for;
end procedure;

