/*This function takes a set X, consisting of a set N of integers, and a permutation g in Sym(N), and returns the set of
representatives of orbits of elements of X under the action of g.
*/

function OrbitReps(X,g)

Orbs:=[];
SpanOfOrbs:=[];
for i in X do
  if(not(Sort(i) in SpanOfOrbs)) then
    Append(~Orbs,i);
    SpanOfOrbs cat:=[Sort(i^(g^j)):j in [1..Order(g)-1]];
  end if;
end for;
return Orbs;
end function;

g:=Sym(7)!(2,3,4,5,6,7);

function FindOrb(Orbs,i)
for Orb in Orbs do if(i in Orb) then return Orb; end if; end for;
return {};
end function;

G2:=PSp(4,8);
CC:=ConjugacyClasses(G2);
load "Subgroups/functions.txt";
load "Traces/E8/TracesE8.txt";
load "Traces/E8/Tr21E8Eigs";
Tr21E8:={EigsToTracePowers(i,21):i in Tr21E8Eigs};
delete Tr21E8Eigs;
p:=2;
Tr21:={[i[1]]:i in Tr21E8};
Tr:=<{[248]},{},Tr3E8,{},Tr5E8,{},Tr7E8,{},Tr9E8,{},{},{},Tr13E8,{},{},{},{},{},{},{},Tr21E8>;
lim:=#[i:i in [j[1]:j in CC]|i le #Tr];
phi:=PowerMap(G2);
CPM:=[[phi(j,i):i in Prune(Divisors(CC[j,1]))]:j in [1..lim]];

X1:=OrderByDimension(CompositionFactors(PermutationModule(G2,GF(8))));
S41:=X1[10];
S43:=Socle(SymmetricSquare(S41)); S45:=Socle(SymmetricSquare(S43));
S42:=SocleFactors(SymmetricSquare(S41))[3];
S44:=Socle(SymmetricSquare(S42)); S46:=Socle(SymmetricSquare(S44));

S1612:=TensorProduct(S41,S42);
S1613:=TensorProduct(S41,S43);
S1614:=TensorProduct(S41,S44);
S1615:=TensorProduct(S41,S45);
S1616:=TensorProduct(S41,S46);
S1623:=TensorProduct(S42,S43);
S1624:=TensorProduct(S42,S44);
S1625:=TensorProduct(S42,S45);
S1626:=TensorProduct(S42,S46);
S1634:=TensorProduct(S43,S44);
S1635:=TensorProduct(S43,S45);
S1636:=TensorProduct(S43,S46);
S1645:=TensorProduct(S44,S45);
S1646:=TensorProduct(S44,S46);
S1656:=TensorProduct(S45,S46);

Irr:=[X1[1],S41,S42,S43,S44,S45,S46,S1612,S1613,S1614,S1615,S1616,S1623,S1624,S1625,S1626,S1634,S1635,S1636,S1645,S1646,S1656];
for i in [2..5] do for j in [i+1..6] do for k in [j+1..7] do Append(~Irr,TensorProduct(TensorProduct(Irr[i],Irr[j]),Irr[k])); end for; end for; end for;

/* This is all of the cases, without the cases excluded theoretically in the text
AnsAll:=[MultisetToSequence(i):i in
[{*1^^8,4^^12,16^^4,64^^2*},
{*1^^12,4^^15,16^^3,64^^2*},
{*1^^16,4^^18,16^^2,64^^2*},
{*1^^20,4^^21,16,64^^2*},
{*1^^8,4^^16,16^^7,64*},
{*1^^12,4^^19,16^^6,64*},
{*1^^16,4^^22,16^^5,64*},
{*1^^20,4^^25,16^^4,64*},
{*1^^24,4^^28,16^^3,64*},
{*1^^28,4^^31,16^^2,64*},
{*1^^32,4^^34,16,64*},
{*1^^36,4^^37,64*},
{*1^^8,4^^20,16^^10*},
{*1^^12,4^^23,16^^9*},
{*1^^16,4^^26,16^^8*},
{*1^^20,4^^29,16^^7*},
{*1^^24,4^^32,16^^6*},
{*1^^28,4^^35,16^^5*},
{*1^^32,4^^38,16^^4*},
{*1^^36,4^^41,16^^3*},
{*1^^40,4^^44,16^^2*},
{*1^^44,4^^47,16*},
{*1^^48,4^^50*}]];
*/

AnsAll:=[MultisetToSequence(i):i in 
[{*1^^8,4^^12,16^^4,64^^2*},
{*1^^12,4^^15,16^^3,64^^2*},
{*1^^16,4^^18,16^^2,64^^2*},
{*1^^8,4^^16,16^^7,64*},
{*1^^12,4^^19,16^^6,64*},
{*1^^16,4^^22,16^^5,64*},
{*1^^20,4^^25,16^^4,64*},
{*1^^24,4^^28,16^^3,64*},
{*1^^28,4^^31,16^^2,64*},
{*1^^8,4^^20,16^^10*},
{*1^^12,4^^23,16^^9*},
{*1^^16,4^^26,16^^8*},
{*1^^20,4^^29,16^^7*},
{*1^^24,4^^32,16^^6*},
{*1^^28,4^^35,16^^5*},
{*1^^32,4^^38,16^^4*}]];

Br:=[BrauerCharacter(i):i in Irr];
dims:=[Dimension(i):i in Irr];

function TestCase(iii)

Blocks1:=[i:i in AnsAll[iii]|i eq 1]; Blocks2:=[i:i in AnsAll[iii]|i eq 4];
Blocks3:=[i:i in AnsAll[iii]|i eq 16]; Blocks4:=[i:i in AnsAll[iii]|i eq 64];

Bits:=[];

for jj in [Blocks1,Blocks2,Blocks3,Blocks4] do
Ans:=[jj];
Ans2:=[];
for CurrentDegreeSeq in Ans do
  NewAns:=[];
  for j in SequenceToSet(CurrentDegreeSeq) do
m:=[]; for k in [1..#dims] do if(dims[k] eq j) then Append(~m,k); end if; end for;
n:=Multiplicity(SequenceToMultiset(CurrentDegreeSeq),j);
// n:=0; for k in CurrentDegreeSeq do if(k eq j) then n+:=1; end if; end for;
ThingsToAppend:=DecreasingSequences(n,m);
if(NewAns eq []) then NewAns:=ThingsToAppend; delete ThingsToAppend;
else
  NewAns2:=[];
  for alpha in NewAns do
for beta in ThingsToAppend do Append(~NewAns2,alpha cat beta);
end for;
  end for;
  delete NewAns; NewAns:=NewAns2; delete NewAns2; delete ThingsToAppend;
end if;
  end for;
  Ans2 cat:=NewAns;
end for;

Append(~Bits,Ans2);
delete Ans2;
end for;

// To reduce the number to check, do a partial automorphism collapse here. Could just do OrbitReps(Bits[2],g), but this is faster

Bits[2]:=Sort([Sort(i):i in Bits[2]]);
BB2:=[i:i in Bits[2]|Maximum([Multiplicity(i,j):j in [2..7]]) eq Multiplicity(i,2)];
BB2a:=[i:i in BB2|Multiplicity(i,2) gt Maximum([Multiplicity(i,j):j in [3..7]])];
BB2b:=[i:i in BB2|Multiplicity(i,2) eq Maximum([Multiplicity(i,j):j in [3..7]])];
BB2bb:=OrbitReps(BB2b,g);
Bits[2]:=Sort(BB2a cat BB2bb);
delete BB2; delete BB2a; delete BB2b; delete BB2bb;


al:=33;
al2:=48;
al3:=12;

Tral:=[ [ &+[Br[j,al]:j in k]:k in jj]:jj in Bits];
TrSetal:=[SequenceToSet(i):i in Tral];

// If there are no factors of a given dimension, add a zero for that entry to make the loop work
for i in [1..4] do if(TrSetal[i] eq {}) then Include(~TrSetal[i],0); end if; end for;

bigno:=&*[#i:i in TrSetal];
Ans3:=[]; nn:=0;
for i in TrSetal[1] do for j in TrSetal[2] do for k in TrSetal[3] do for l in TrSetal[4] do
  nn+:=1; if(nn mod 10000 eq 0) then "We have tested",nn,"of",bigno,"and so far found",#Ans3,"that are conspicuous"; end if;
  if([i+j+k+l] in Tr13E8) then
    a1:=[Bits[1,ii]:ii in [1..#Tral[1]]|Tral[1,ii] eq i];
    a2:=[Bits[2,ii]:ii in [1..#Tral[2]]|Tral[2,ii] eq j];
    a3:=[Bits[3,ii]:ii in [1..#Tral[3]]|Tral[3,ii] eq k];
    a4:=[Bits[4,ii]:ii in [1..#Tral[4]]|Tral[4,ii] eq l];
    // We now have to repair any of these that are empty
    if(a3 eq []) then a3:=[[]]; end if;
    if(a4 eq []) then a4:=[[]]; end if;
    Ans2:=[a cat b cat c cat d:a in a1,b in a2, c in a3, d in a4];
    for ia in [1..#Ans2] do
      if (ia mod 100000) eq 0 then nn,ia,#Ans2,#Ans3; end if;
      test:=&+[Br[kk,al2]:kk in Ans2[ia]];
      if([test] in Tr21) then delete test;
        test2:=&+[Br[kk,al3]:kk in Ans2[ia]];
        if([test2] in Tr7E8) then delete test2;
          test3:=&+[Br[kk,al3+3]:kk in Ans2[ia]];
          if([test3] in Tr7E8) then delete test3;
            test:=&+[Br[kk]:kk in Ans2[ia]];
            for jj in [2..#CPM] do
              if(CC[jj,1] le #Tr and GCD(CC[jj,1],p) eq 1) then
                ord:=CC[jj,1]; divs:=Prune(Divisors(ord));
                tt:=[test[kk]:kk in CPM[jj]];
                if(not(tt in Tr[ord])) then delete test; continue ia; end if;
              end if;
            end for;
          delete test; Append(~Ans3,Ans2[ia]);
          end if;
        end if;
      end if;
    end for;
    delete Ans2;
  end if;
end for; end for; end for; end for;
Ans3:=Sort([Sort(ii):ii in Ans3]);

delete Tral; delete TrSetal; delete Bits;

return Ans3;

end function;

/*
1 has solutions
[
    [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 4, 4, 4, 4, 5, 5, 9, 9, 10, 13, 26, 35 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 5, 5, 5, 5, 6, 6, 10, 10, 11, 14, 29, 40 ]
]
2 has solutions
[
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 7, 7, 7, 7, 8, 12, 16, 32, 32 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 4, 4, 4, 4, 5, 7, 7, 7, 7, 9, 12, 19, 38, 38 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 5, 5, 5, 5, 6, 7, 7, 7, 7, 10, 12, 21, 41, 41 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6, 6, 6, 6, 7, 9, 11, 18, 37, 37 ]
]
5 has solution
[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 6, 6, 6, 6, 7, 7, 7, 7, 11, 11, 12, 12, 15, 16, 32 ]
12 has solutions
[
[1,1,1,1,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 11, 11, 11, 11, 12, 12, 15, 15 ],
[1,1,1,1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 7, 10, 10, 10, 10, 11, 11, 14, 14 ]
]
14 has solution
[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 7, 7,
7, 7, 7, 7, 7, 7, 12, 12, 12, 12, 16, 16 ]
16 has solution
[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 8, 8, 8, 8 ]

3,4,6,7,8,9,10,11,13,15 have none
*/

// Use this to determine which 64 is which.
//j:=32; [i-1:i in [2..7]|#CompositionFactors(TensorProduct(Irr[i],Irr[j])) gt 1];

// Running TestCase above on all options yields the following.

Ans3:=[
    [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 4, 4, 4, 4, 5, 5, 9, 9, 10, 13, 26, 35 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 5, 5, 5, 5, 6, 6, 10, 10, 11, 14, 29, 40 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 7, 7, 7, 7, 8, 12, 16, 32, 32 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 4, 4, 4, 4, 5, 7, 7, 7, 7, 9, 12, 19, 38, 38 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 5, 5, 5, 5, 6, 7, 7, 7, 7, 10, 12, 21, 41, 41 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6, 6, 6, 6, 7, 9, 11, 18, 37, 37 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 6, 6, 6, 6, 7, 7, 7, 7, 11, 11, 12, 12, 15, 16, 32 ],
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,3,3,3,3,4,6,6,6,6,6,6,6,6,7,7,7,7,11,11,11,11,12,12,15,15],
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,4,5,5,5,5,5,5,5,5,6,6,6,6,7,10,10,10,10,11,11,14,14],
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,4,4,4,7,7,7,7,7,7,7,7,12,12,12,12,16,16],
    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,8,8,8,8]
];
Mods:=[DirectSum([Irr[i]:i in j]):j in Ans3];

q:=2^12;
m:=q-1;
F<w>:=GF(q);
w63:=w^65; w65:=w^63;
//n:=(q-1) div m;
G:=GroupOfLieType("E8",q);
V:=VectorSpace(GF(q),8);
rho1:=StandardRepresentation(G);
Over1:=GL(248,q);
g1:=elt<G|V![w,1,1,1,1,1,1,1]>;
g2:=elt<G|V![1,w,1,1,1,1,1,1]>;
g3:=elt<G|V![1,1,w,1,1,1,1,1]>;
g4:=elt<G|V![1,1,1,w,1,1,1,1]>;
g5:=elt<G|V![1,1,1,1,w,1,1,1]>;
g6:=elt<G|V![1,1,1,1,1,w,1,1]>;
g7:=elt<G|V![1,1,1,1,1,1,w,1]>;
g8:=elt<G|V![1,1,1,1,1,1,1,w]>;
W:=Reflections(G);
Mats1:=[rho1(i):i in [g1,g2,g3,g4,g5,g6,g7,g8]];
MatsE:=[rho1(i):i in W];
NGT:=sub<Over1|Mats1 cat MatsE>;
T:=sub<NGT|Mats1>;
E:=sub<NGT|MatsE>;

T7:=sub<T|[T.i^(9*65):i in [1..8]]>;
T21:=sub<T|[T.i^(3*65):i in [1..8]]>;
T63:=sub<T|[T.i^65:i in [1..8]]>;
T13:=sub<T|[T.i^(5*63):i in [1..8]]>;
T65:=sub<T|[T.i^63:i in [1..8]]>;
T195:=sub<T|[T.i^21:i in [1..8]]>;
T1365:=sub<T|[T.i^3:i in [1..8]]>;

R:=ResidueClassRing(1365);
X:=[[R!(Log((T1365.j)[i,i]) div 3):i in [1..248]]:j in [1..8]];


for i in [72..83] do if(EigenvaluesOfElement(ChangeRing(S41,F),CC[i,3]) eq { <w^4032, 1>, <w^3591, 1>, <w^504, 1>, <w^63, 1> }) then y:=CC[i,3]; end if; end for;
Eigs1:=[EigenvaluesOfElement(ChangeRing(i,F),y):i in Mods];
Eigs2:=[EigenvaluesOfElement(ChangeRing(i,F),y^5):i in Mods];


/* If you want to find your own representatives use this code. I used it to obtain El1 below.
El1:=[];
for i in [1..11] do repeat x1:=Random(T13); until Eigenvalues(x1) eq Eigs2[i]; "Found one"; Append(~El1,ExpressElementInGenerators(T13,x1)); El1[#El1]; end for;
*/

El1:=[
    [ 9, 8, 6, 7, 7, 7, 9, 6 ],
    [ 1, 11, 12, 8, 12, 1, 5, 4 ],
    [ 9, 0, 8, 1, 11, 5, 12, 11 ],
    [ 4, 3, 10, 0, 8, 3, 3, 2 ],
    [ 10, 4, 7, 12, 0, 8, 11, 12 ],
    [ 11, 3, 0, 4, 8, 4, 3, 8 ],
    [ 10, 6, 7, 5, 7, 3, 9, 7 ],
    [ 2, 8, 1, 4, 8, 5, 8, 7 ],
    [ 2, 1, 5, 7, 4, 4, 10, 0 ],
    [ 1, 8, 9, 0, 2, 5, 5, 12 ],
    [ 4, 6, 2, 6, 1, 4, 1, 7 ]
];

El2:=[
    [ 735, 987, 945, 966, 420, 420, 189, 126 ],
    [ 567, 504, 1071, 987, 798, 294, 105, 84 ],
    [ 1281, 819, 168, 21, 504, 1197, 252, 504 ],
    [ 357, 609, 1029, 1092, 168, 63, 63, 42 ],
    [ 483, 630, 966, 798, 0, 168, 777, 252 ],
    [ 777, 63, 0, 630, 714, 357, 882, 168 ],
    [ 756, 1218, 693, 651, 693, 609, 189, 147 ],
    [ 588, 168, 21, 1176, 168, 105, 441, 147 ],
    [ 42, 840, 105, 693, 84, 84, 210, 0 ],
    [ 567, 441, 189, 0, 588, 924, 1197, 798 ],
    [ 1176, 1218, 42, 126, 21, 1176, 21, 147 ]
];

NextEls1All:=[];
for OneToCheck in [1..#El1] do
  Els1:=ConstructPreimages([105*i:i in El1[OneToCheck]],1365,5);
  e22:={*i[1]^^i[2]:i in Eigs1[OneToCheck]*};
  e2:={*R!(Log(i) div 3):i in e22*};
  Exclude(~e2,0^^8);
  NextEls1:=[];
  for ii in [1..#Els1] do
    hh:={*&+[Els1[ii,i]*X[i,j]:i in [1..8]]:j in [1..120]*};
    if e2 eq hh join {*-i:i in hh*} then Append(~NextEls1,Els1[ii]); end if;
    if(ii mod 10000 eq 0) then "At stage",ii,"of",#Els1,"for case",OneToCheck,"there are",#NextEls1,"elements with the correct eigenvalues"; end if;
  end for;
  Append(~NextEls1All,NextEls1);
  El2[OneToCheck] eq NextEls1[1];
  delete NextEls1;
end for;

EE:=[
[E.2*E.5*E.4*E.3*E.1*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.6*E.5*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.1*E.7*E.6*E.5*E.4*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.1*E.7*E.6*E.5*E.4*E.2*E.3*E.4*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1,
E.3*E.4*E.2*E.5*E.4*E.3*E.1*E.6*E.5*E.4*E.2*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.6*E.5*E.7*E.6*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.7],
[E.1*E.4*E.2*E.3*E.1*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.7*E.6*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.7,
E.3*E.1*E.4*E.2*E.5*E.4*E.3*E.1*E.7*E.6*E.5*E.4*E.2*E.3*E.4*E.5*E.6*E.7*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.7],
[E.2*E.4*E.3*E.5*E.4*E.2*E.3*E.1*E.6*E.5*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.7*E.6*E.5*E.4*E.2*E.8,
E.2*E.4*E.2*E.3*E.1*E.4*E.5*E.4*E.2*E.3*E.4*E.5*E.6*E.5*E.4*E.3*E.7*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4],
[E.3*E.1*E.4*E.5*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.7*E.6*E.5*E.8,
E.2*E.3*E.4*E.2*E.3*E.4*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.6*E.5*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.1*E.7*E.6*E.5*E.4*E.2*E.3*E.4*E.5*E.6*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.6*E.5],
[E.1*E.3*E.1*E.4*E.5*E.4*E.2*E.3*E.4*E.5*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.6*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.6,E.1*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.3*E.1*E.4*E.6*E.5*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.1*E.7*E.6*E.5*E.4*E.2*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.1*E.7*E.6*E.5*E.4*E.2*E.3*E.4*E.5*E.6*E.7*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.6,E.3*E.1*E.4*E.2*E.3*E.1*E.4*E.5*E.4*E.2*E.3*E.4*E.6*E.5*E.4*E.2*E.3*E.4*E.5*E.6*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.6*E.5],
[E.2*E.4*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.6*E.5*E.4*E.2*E.3*E.1*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.7*E.6*E.5*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.7*E.6*E.5,
E.2*E.3*E.1*E.4*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.6*E.5*E.4*E.2*E.3*E.1*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.6*E.5*E.7*E.6*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.7*E.6*E.5],
[E.4*E.2*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.7*E.6*E.5*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4,
E.1*E.4*E.3*E.5*E.4*E.2*E.7*E.6*E.5*E.4*E.3*E.1*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.1*E.7*E.6*E.5*E.4*E.8*E.7],
[E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.2*E.3,
E.6*E.5*E.7*E.6],
[E.2*E.3*E.4*E.2*E.5*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.1*E.7*E.6*E.5*E.4*E.2*E.3*E.4*E.5*E.8*E.7*E.6*E.5*E.4*E.2*E.3,
E.1*E.4*E.2*E.3*E.5*E.4*E.2*E.3*E.4*E.5*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.1*E.7*E.6*E.5*E.4*E.2*E.3*E.4*E.5*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2],
[],
[]];

for OneToCheck in [1..#El1] do
  NextEls1:=NextEls1All[OneToCheck];
  if #NextEls1 le 1 then "Case",OneToCheck,"Nothing to check here"; continue OneToCheck; end if;
  gg1:=ProduceGroupElements(T1365,NextEls1);
  Orbs:={{i}:i in [1..#NextEls1]};
  for e in EE[OneToCheck] do
    for i in [1..#NextEls1] do
      j:=Position(gg1,gg1[i]^e);
      Orb1:=FindOrb(Orbs,i);
      if(not(j in Orb1)) then
        Orb2:=FindOrb(Orbs,j);
        Exclude(~Orbs,Orb1);
        Exclude(~Orbs,Orb2);
        Include(~Orbs,Orb1 join Orb2); //printf "Merged %o and %o.\n",Orb1,Orb2;
        delete Orb1; delete Orb2; if(#Orbs eq 1) then "Case",OneToCheck,"All merged"; break e; end if;
      end if;
    end for;
//  "So far we have merged elements into",#Orbs,"orbits";
  end for;
end for;

EigsToTest1:=[EigenvaluesOfElement(ChangeRing(Irr[i],F),y):i in [1..#Irr]];
EigsToTest2:=[{R!(Log(i[1]) div 3):i in j}:j in EigsToTest1];

for OneToCheck in [1..#El1] do
  Els2:=ConstructPreimages(El2[OneToCheck],1365,3);
  EigEls2:=[{*&+[ii[i]*X[i,j]:i in [1..8]]:j in [1..248]*}:ii in Els2];
  EigMults2:=[{*i*3:i in MultisetToSet(j)*}:j in EigEls2];
  WhatMults:=[{i:i in j|Multiplicity(j,i) eq 1}:j in EigMults2];
  OnesToTest:=Sort(SetToSequence(SequenceToSet(Ans3[OneToCheck])));
  dd:=[[i:i in [1..#WhatMults]|EigsToTest2[j] subset WhatMults[i]]:j in OnesToTest];
  "For case",OneToCheck,"the numbers of elements stabilizing the eigenspaces of each irreducible module are",[#i-1:i in dd];
  "There are",#[i:i in dd[1]|&and[i in dd[j]:j in [2..#OnesToTest]]]-1,"elements stabilizing all eigenspaces simultaneously";
  if(OneToCheck eq 2) then "There are",#[i:i in dd[2]|i in dd[4]],"elements stabilizing 4.1 and 4.4 simultaneously in Case 2"; end if;
  if(OneToCheck eq 9) then
    "There are",#[i:i in dd[2]|i in dd[5]]-1,"elements stabilizing 4.1 and 4.4 simultaneously in Case 9";
    "There are",#[i:i in dd[4]|i in dd[7]]-1,"elements stabilizing 4.3 and 4.6 simultaneously in Case 9";
  end if;
end for;

// Now test the radicals
"Now we look at the radicals from the paper.";
labs:=["1","4_1","4_2","4_3","4_4","4_5","4_6","16_{12}","16_{13}","16_{14}","16_{15}","16_{16}",
"16_{23}","16_{24}","16_{25}","16_{26}"];
S1:=Irr[1];

"Case 1:";
A1:=S1613; for i in [S41,S43,S1,S42,S1] do E,rho:=Ext(A1,i); A1:=MaximalExtension(A1,i,E,rho); DescribeLayers(Dual(A1),Irr,labs); end for;
for i in [S1,S41,S42,S43,S44] do Dimension(Ext(A1,i)) eq 0; end for;

"Case 2:";
A2:=S1614; for i in [S41,S44,S1] do E,rho:=Ext(A2,i); A2:=MaximalExtension(A2,i,E,rho); DescribeLayers(Dual(A2),Irr,labs); end for;
for i in [S1,S41,S42,S44,S45] do Dimension(Ext(A2,i)) eq 0; end for;
for i in [S1615,S1624] do E,rho:=Ext(A2,i); A2:=MaximalExtension(A2,i,E,rho); DescribeLayers(Dual(A2),Irr,labs); end for;
for i in [S41,S42,S44,S45,S1] do E,rho:=Ext(A2,i); A2:=MaximalExtension(A2,i,E,rho); DescribeLayers(Dual(A2),Irr,labs); end for;
for i in [S1,S41,S42,S44,S45] do Dimension(Ext(A2,i)) eq 0; end for;

A3:=S42; for i in [S1,S41,S1] do E,rho:=Ext(A3,i); A3:=MaximalExtension(A3,i,E,rho); DescribeLayers(Dual(A3),Irr,labs); end for;
for i in [S1,S41,S44,S1614] do Dimension(Ext(A3,i)) eq 0; end for;
for i in [S1615,S1624] do E,rho:=Ext(A3,i); A3:=MaximalExtension(A3,i,E,rho); DescribeLayers(Dual(A3),Irr,labs); end for;
for i in [S41,S44,S1614,S1,S41,S44] do E,rho:=Ext(A3,i); A3:=MaximalExtension(A3,i,E,rho); DescribeLayers(Dual(A3),Irr,labs); end for;
for i in [S1,S41,S44,S1614] do Dimension(Ext(A3,i)) eq 0; end for;

"Case 6:";
A4:=BuildLargestModule(S41,[S1,S43,S45]); DescribeLayers(A4,Irr,labs);
"The others are images under the field automorphism.";

"Case 9:";

A5:=S1615;
for i in [S41,S45,S1,S1614,S41,S44,S46,S1] do E,rho:=Ext(A5,i); A5:=MaximalExtension(A5,i,E,rho); DescribeLayers(Dual(A5),Irr,labs); end for;
for i in [S1,S41,S42,S43,S44,S45,S46,S1614,S1624] do Dimension(Ext(A5,i)) eq 0; end for;

// A5a is the image of A5 under an automorphism, so we don't need to check that it is the radical.
A5a:=S1624;
for i in [S44,S42,S1,S1614,S44,S41,S43,S1] do E,rho:=Ext(A5a,i); A5a:=MaximalExtension(A5a,i,E,rho); DescribeLayers(Dual(A5a),Irr,labs); end for;

A6:=S42; for i in [S1,S41,S1] do E,rho:=Ext(A6,i); A6:=MaximalExtension(A6,i,E,rho); DescribeLayers(Dual(A6),Irr,labs); end for;
for i in [S1,S41,S42,S44,S45,S1614] do Dimension(Ext(A6,i)); end for;
"The other radical is the image under the field automorphism of this one.";
