/*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(6)!(2,3,4,5,6);
numofsimps:=6;

G2:=SuzukiGroup(32);
CC:=ConjugacyClasses(G2);
load "Subgroups/functions.txt";

load "Traces/E8/Tr25E8Eigs";
Tr25E8:={EigsToTracePowers(i,25):i in Tr25E8Eigs};
Tr25:={i[1]:i in Tr25E8|i[2] eq -2};
delete Tr25E8Eigs; delete Tr25E8;
load "Traces/E8/Tr31E8Eigs";
Tr31:={&+j:j in Tr31E8Eigs};
delete Tr31E8Eigs;

S1:=TrivialModule(G2,GF(32));
S41:=GModule(G2);
S42:=GModule(G2,FrobeniusImage(S41,1));
S43:=GModule(G2,FrobeniusImage(S41,2));
S44:=GModule(G2,FrobeniusImage(S41,3));
S45:=GModule(G2,FrobeniusImage(S41,4));

S1612:=TensorProduct(S41,S42);
S1613:=TensorProduct(S41,S43);
S1614:=TensorProduct(S41,S44);
S1615:=TensorProduct(S41,S45);
S1623:=TensorProduct(S42,S43);
S1624:=TensorProduct(S42,S44);
S1625:=TensorProduct(S42,S45);
S1634:=TensorProduct(S43,S44);
S1635:=TensorProduct(S43,S45);
S1645:=TensorProduct(S44,S45);

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

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 *}
]];

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


al:=11;
al2:=6;

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 100000 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 Tr31) 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 10000) eq 0 then nn,ia,#Ans2,#Ans3; end if;
      test:=&+[Br[kk,al2]:kk in Ans2[ia]];
      if(test in Tr25) then
        delete test; Append(~Ans3,Ans2[ia]);
      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 solution
[
    [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 5, 5, 6, 6, 7, 7, 10, 12,
    20, 23 ]
]
2 has solutions
[
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4,
    5, 6, 7, 8, 11, 26, 26 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 4, 4, 4, 4, 5, 5, 5,
    5, 5, 8, 9, 14, 21, 21 ]
]
6 has solution
[
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 4, 4, 4, 4, 5, 5,
    5, 5, 6, 6, 6, 6, 8, 8, 10, 10, 14, 16, 21 ]
]
15 has solution
[
    [ 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, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 8, 8, 10, 10, 10, 10, 16,
    16 ]
]
17 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, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5,
    5, 5, 5, 5, 5, 5, 8, 8, 8, 8, 14, 14 ]
]
19 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, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 9, 9, 9, 9 ]
]

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

These cases have solutions:
[ 1, 2, 6, 13, 15, 17, 19 ]
*/
AnsTot:=[
    [ 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 5, 5, 6, 6, 7, 7, 10, 12,
    20, 23 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4,
    5, 6, 7, 8, 11, 26, 26 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 4, 4, 4, 4, 5, 5, 5,
    5, 5, 8, 9, 14, 21, 21 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 4, 4, 4, 4, 5, 5,
    5, 5, 6, 6, 6, 6, 8, 8, 10, 10, 14, 16, 21 ],
    [ 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, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 8, 8, 10, 10, 10, 10, 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, 2,
    2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5,
    5, 5, 5, 5, 5, 5, 8, 8, 8, 8, 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, 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, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 9, 9, 9, 9 ]
];
Mods:=[DirectSum([Irr[i]:i in j]):j in AnsTot];

"Computing the number of elements stabilizing the submodules.";

q:=2^10;
F<w>:=GF(q);
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>;

T31:=sub<T|[T.i^33:i in [1..8]]>;
T93:=sub<T|[T.i^11:i in [1..8]]>;

R:=ResidueClassRing(93);
X:=[[R!(Log((T93.j)[i,i]) div 11):i in [1..248]]:j in [1..8]];


for i in [11..25] do if(EigenvaluesOfElement(ChangeRing(S41,F),CC[i,3]) eq { <w^726, 1>, <w^297, 1>, <w^990, 1>, <w^33, 1> }) then x:=CC[i,3]; end if; end for;
Eigs1:=[EigenvaluesOfElement(ChangeRing(i,F),x):i in Mods];

// repeat x1:=Random(T31); until Eigenvalues(x1) in Eigs1;

//Elements of T31 that have eigenvalues in Eigs1
El1:=[
[ 28, 2, 14, 10, 0, 4, 29, 14 ],
[ 14, 28, 20, 5, 27, 11, 29, 17 ],
[ 10, 6, 19, 30, 27, 6, 16, 25 ],
[ 13, 20, 25, 26, 17, 8, 21, 20 ],
[ 6, 17, 9, 6, 20, 11, 25, 14 ],
[ 22, 3, 8, 19, 4, 10, 0, 21 ],
[ 30, 1, 1, 21, 17, 2, 21, 22 ]
];

&and[Eigenvalues(ProduceGroupElement(T31,El1[i])) eq Eigs1[i]:i in [1..#Eigs1]];

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

for OneToCheck in [1..#El1] do
  Els1:=ConstructPreimages([3*i:i in El1[OneToCheck]],93,3);
  EigEls1:=[{*&+[ii[i]*X[i,j]:i in [1..8]]:j in [1..248]*}:ii in Els1];
  EigMults1:=[{*i*3:i in MultisetToSet(j)*}:j in EigEls1];
  WhatMults:=[{i:i in j|Multiplicity(j,i) eq 1}:j in EigMults1];
  OnesToTest:=Sort(SetToSequence(SequenceToSet(AnsTot[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";
end for;

labs:=["1","4_1","4_2","4_3","4_4","4_5","16_{12}","16_{13}","16_{14}","16_{15}",
"16_{23}","16_{24}","16_{25}","16_{34}","16_{35}","16_{45}"];

"Now check the statements in the proof:";

"Case 2 radical";
A1:=BuildLargestModule(S41,[S1,S43,S42,S41]); DescribeLayers(A1,Irr,labs);

"Case 4 radicals:";

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

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

A4:=S1615; for i in [S41,S45,S1613,S1,S41,S43,S1,S44,S1,S41,S1,S43] do E,rho:=Ext(A4,i); A4:=MaximalExtension(A4,i,E,rho); DescribeLayers(Dual(A4),Irr,labs); end for;
for i in [S1,S41,S43,S44,S45,S1613] do Dimension(Ext(A4,i)) eq 0; end for;
"The Ext^1 with 16_{15} is",Dimension(Ext(A4,S1615));
"That last 1 is unimportant, as mentioned in the article";

"Case 5 radicals";

// Use A3 as before, but now need to add copies of S1645 on it as well.

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

A6:=S1645; for i in [S44,S45,S1615,S1,S41,S45,S1613,S1,S41,S43,S1,S44,S1,S41,S1,S43] do E,rho:=Ext(A6,i); A6:=MaximalExtension(A6,i,E,rho); DescribeLayers(Dual(A6),Irr,labs); end for;
for i in [S1,S41,S43,S44,S45,S1613,S1615] do Dimension(Ext(A6,i)) eq 0; end for;

A5a:=RemoveFromTop(JacobsonRadical(A5),Irr,Irr,[1,2,4,5,6,10]);
A5b:=Dual(A5/A5a); DescribeLayers(A5b,Irr,labs);

A6a:=RemoveFromTop(JacobsonRadical(A6),Irr,Irr,[1,2,4,5,6,10]);
A6b:=Dual(A6/A6a); DescribeLayers(A6b,Irr,labs);


"Case 1 radicals and unipotent actions.";

Dimension(Ext(Irr[7],Irr[20])) eq 0; Dimension(Ext(Irr[7],Irr[23])) eq 0;

// We find the submodules of P44 and P45 first. To make all of the other modules in the proof we will take submodules of these.

A7:=S44; for i in [S1,S41,S42,S1612,S1,S41,S42,S45,S1,S42,S1] do E,rho:=Ext(A7,i); A7:=MaximalExtension(A7,i,E,rho); DescribeLayers(Dual(A7),Irr,labs); end for;
for i in [S1,S41,S42,S45,S1612] do Dimension(Ext(A7,i)) eq 0; end for;
A7stage1:=A7;
for i in [S1615,S1624,S41,S45,S1612,S1,S41,S42,S1] do E,rho:=Ext(A7,i); A7:=MaximalExtension(A7,i,E,rho); DescribeLayers(Dual(A7),Irr,labs); end for;
for i in [S1,S41,S42,S45,S1612] do Dimension(Ext(A7,i)) eq 0; end for;
A7stage2:=A7;
E,rho:=Ext(A7,S44); A7:=MaximalExtension(A7,S44,E,rho); DescribeLayers(Dual(A7),Irr,labs);
A7stage3:=A7;
A7:=RemoveFromTop(Dual(A7),Irr,Irr,[1,3,10,12]); DescribeLayers(A7,Irr,labs);

A8:=S45; for i in [S1,S42,S1,S44] do E,rho:=Ext(A8,i); A8:=MaximalExtension(A8,i,E,rho); DescribeLayers(Dual(A8),Irr,labs); end for;
for i in [S1,S41,S42,S44,S1612] do Dimension(Ext(A8,i)) eq 0; end for;
A8stage1:=A8;
for i in [S1615,S1624,S41,S42,S44,S1,S1612,S41] do E,rho:=Ext(A8,i); A8:=MaximalExtension(A8,i,E,rho); DescribeLayers(Dual(A8),Irr,labs); end for;
for i in [S1,S41,S42,S44,S1612] do Dimension(Ext(A8,i)) eq 0; end for;
A8stage2:=A8;
E,rho:=Ext(A8,S45); A8:=MaximalExtension(A8,S45,E,rho); DescribeLayers(Dual(A8),Irr,labs);
A8stage3:=A8;
A8:=RemoveFromTop(Dual(A8),Irr,Irr,[1,3,10,12]); DescribeLayers(A8,Irr,labs);

v:=CC[3,3];
"Action of v on radical for 4_4:";
ActionOfUnipotent([A7],v);

// Now the one for P41. As with A7stage3 and A8stage3, we keep A9stage2 for later use.

"Not the radical for 4_1:";

A9:=S41; for i in [S1612,S1,S41,S42,S44,S1,S41,S42,S44,S1,S44,S1,S41] do E,rho:=Ext(A9,i); A9:=MaximalExtension(A9,i,E,rho); DescribeLayers(Dual(A9),Irr,labs); end for;
for i in [S1,S41,S42,S44,S45,S1612] do Dimension(Ext(A9,i)) eq 0; end for;
A9stage1:=A9;
for i in [S1615,S1624,S41,S42,S45,S1,S45,S1,S42,S1,S44,S1612,S41,S42,S1,S44] do E,rho:=Ext(A9,i); A9:=MaximalExtension(A9,i,E,rho); DescribeLayers(Dual(A9),Irr,labs); end for;
for i in [S1,S41,S42,S44,S45,S1612] do Dimension(Ext(A9,i)) eq 0; end for;
A9stage2:=A9;
A9a:=RemoveFromTop(Dual(A9),Irr,Irr,[1,3,5,6,7,10,12]); DescribeLayers(A9a,Irr,labs);

"Now the smaller modules with no 4_1s in the middle, for when there is 4_1^2 in the socle";

A9b:=RemoveFromTop(JacobsonRadical(A9stage2),Irr,Irr,[1,3,5,6,7,10,12]);
A9b:=RemoveFromTop(A9b,Irr,Irr,[2]);
A9b:=Dual(A9/A9b);
A9b:=RemoveFromTop(A9b,Irr,Irr,[1,3,7,10,12]); 

A9ba:=RemoveFromTop(A9b,Irr,Irr,[1,3,5,6,7,10,12]); 

DescribeLayers(A9ba,Irr,labs);

DescribeLayers(A9b,Irr,labs);

A7b:=RemoveFromTop(A7stage3,Irr,Irr,[1,3,5,6,7,10,12]);
A7b:=RemoveFromTop(A7b,Irr,Irr,[2]);
A7b:=Dual(A7stage3/A7b);
A7b:=RemoveFromTop(A7b,Irr,Irr,[1,3,7,10,12]); DescribeLayers(A7b,Irr,labs);

A8b:=RemoveFromTop(A8stage3,Irr,Irr,[1,3,5,6,7,10,12]);
A8b:=RemoveFromTop(A8b,Irr,Irr,[2]);
A8b:=Dual(A8stage3/A8b);
A8b:=RemoveFromTop(A8b,Irr,Irr,[1,3,7,10,12]); DescribeLayers(A8b,Irr,labs);

"The action of v on these modules";
ActionOfUnipotent([A9b,A7b,A8b],v);

"Now all 4_1, 4_4 and 4_5 in top and socle:";
A7c:=RemoveFromTop(JacobsonRadical(A7stage3),Irr,Irr,[1,3,7,10,12]);
A7c:=RemoveFromTop(A7c,Irr,Irr,[2,5,6]);
A7c:=Dual(A7stage3/A7c);
A7c:=RemoveFromTop(A7c,Irr,Irr,[1,3,7,10,12]); DescribeLayers(A7c,Irr,labs);

A8c:=RemoveFromTop(JacobsonRadical(A8stage3),Irr,Irr,[1,3,7,10,12]);
A8c:=RemoveFromTop(A8c,Irr,Irr,[2,5,6]);
A8c:=Dual(A8stage3/A8c);
A8c:=RemoveFromTop(A8c,Irr,Irr,[1,3,7,10,12]); DescribeLayers(A8c,Irr,labs);

A9c:=RemoveFromTop(JacobsonRadical(A9stage2),Irr,Irr,[1,3,7,10,12]);
A9c:=RemoveFromTop(A9c,Irr,Irr,[2,5,6]);
A9c:=Dual(A9stage2/A9c);
A9c:=RemoveFromTop(A9c,Irr,Irr,[1,3,7,10,12]); DescribeLayers(A9c,Irr,labs);

"and the unipotent element:";
ActionOfUnipotent([A9c,A7c,A8c],v);

"Now 4_4 and 4_5 in the socle and top only:";
A7d:=RemoveFromTop(JacobsonRadical(A7stage3),Irr,Irr,[1,2,3,7,10,12]);
A7d:=RemoveFromTop(A7d,Irr,Irr,[5,6]);
A7d:=Dual(A7stage3/A7d);
A7d:=RemoveFromTop(A7d,Irr,Irr,[1,3,7,10,12]); DescribeLayers(A7d,Irr,labs);

A8d:=RemoveFromTop(JacobsonRadical(A8stage3),Irr,Irr,[1,2,3,7,10,12]);
A8d:=RemoveFromTop(A8d,Irr,Irr,[5,6]);
A8d:=Dual(A8stage3/A8d);
A8d:=RemoveFromTop(A8d,Irr,Irr,[1,3,7,10,12]); DescribeLayers(A8d,Irr,labs);

A9d:=RemoveFromTop(JacobsonRadical(A9stage2),Irr,Irr,[1,2,3,7,10,12]);
A9d:=RemoveFromTop(A9d,Irr,Irr,[5,6]);
A9d:=Dual(A9stage2/A9d);
A9d:=RemoveFromTop(A9d,Irr,Irr,[1,3,7,10,12]); DescribeLayers(A9d,Irr,labs);

"Take the {1}'-heart";
A7d1:=RemoveFromTop(A7d,Irr,Irr,[2..12]);
A7d1:=Dual(RemoveFromTop(Dual(A7d1),Irr,Irr,[2..12]));
A8d1:=RemoveFromTop(A8d,Irr,Irr,[2..12]);
A8d1:=Dual(RemoveFromTop(Dual(A8d1),Irr,Irr,[2..12]));
A9d1:=RemoveFromTop(A9d,Irr,Irr,[2..12]);
A9d1:=Dual(RemoveFromTop(Dual(A9d1),Irr,Irr,[2..12]));

IsIsomorphic(A8d1,A9d1);
DescribeLayers(A7d1,Irr,labs);
DescribeLayers(A8d1,Irr,labs);
"and the dual:";
DescribeLayers(Dual(A8d1),Irr,labs);

"The action of v on these two modules is";

ActionOfUnipotent([A9d1,A7d1],v);


// Here are the two 4.1 and 4.5 we need for socle their sum
"Now 4_5 in top and socle";

A8e:=RemoveFromTop(JacobsonRadical(A8stage3),Irr,Irr,[1,2,3,5,7,10,12]);
A8e:=RemoveFromTop(A8e,Irr,Irr,[6]);
A8e:=Dual(A8stage3/A8e);
A8e:=RemoveFromTop(A8e,Irr,Irr,[1,3,5,7,10,12]); DescribeLayers(A8e,Irr,labs);

A9e1:=RemoveFromTop(JacobsonRadical(A9stage2),Irr,Irr,[1,2,3,5,7,10,12]);
A9e1:=RemoveFromTop(A9e1,Irr,Irr,[6]);
A9e1:=Dual(A9stage2/A9e1);
A9e1:=RemoveFromTop(A9e1,Irr,Irr,[1,3,5,7,10,12]); DescribeLayers(A9e1,Irr,labs);

"The action of v on these two modules is";
ActionOfUnipotent([A9e1,A8e],v);

// The next two are 4.1 ad 4.4 we need for socle their sum

"Now 4_4 in top and socle";

Ords:=[[1,3,6,7],[2],[1,3,6,7],[10,12],[1,3,6,7],[2],[1,3,6,7],[2,5]];

A7e:=JacobsonRadical(A7stage3);
for i in Ords do A7e:=RemoveFromTop(A7e,Irr,Irr,i); end for;
A7e:=Dual(A7stage3/A7e);
A7e:=RemoveFromTop(A7e,Irr,Irr,[1,3,6,7,10,12]); DescribeLayers(A7e,Irr,labs);

A9e2:=JacobsonRadical(A9stage2);
for i in Ords do A9e2:=RemoveFromTop(A9e2,Irr,Irr,i); end for;
A9e2:=Dual(A9stage2/A9e2);
A9e2:=RemoveFromTop(A9e2,Irr,Irr,[1,3,6,7,10,12]); DescribeLayers(A9e2,Irr,labs);

"The action of v on these two modules is";
ActionOfUnipotent([A9e2,A7e],v);
