labs:=["1","5","5^*","10","10^*","24","40_1","40_1^*","40_2","40_2^*","74","160","160^*","280","280^*","1024"];
labsL:=["1","4","4^*","6","14","20","20^*"];
labsM:=["1","1","3","3^*","3","3^*","8","8"];


load "Subgroups/functions.txt";

G:=PSL(5,2);
Irr:=IrreducibleModules(G,GF(2));
L:=Subgroups(G:OrderEqual:=Order(PSL(4,2)))[1]`subgroup;
CC:=ConjugacyClasses(G);
z:=CC[4,3];
M:=Centralizer(G,z);
u:=CC[6,3];
v:=CC[8,3];
w:=CC[14,3];
repeat g:=Random(G); until u^g in M; u:=u^g;

Irr[4]:=ExteriorSquare(Irr[3]);
Irr[5]:=Dual(Irr[4]);
a:=IdentifyFactors(ExteriorSquare(Irr[4]),Irr)[2];
Irr:=Swap(Irr,7,a);
b:=IdentifyFactors(TensorProduct(Irr[2],Irr[5]),Irr)[2];
Irr:=Swap(Irr,9,b);
Irr[8]:=Dual(Irr[7]);
Irr[10]:=Dual(Irr[9]);

LE8:=DirectSum([Irr[i]:i in [2,2,3,3,4,5,6,6,7,8,9,10]]);
EigenvaluesOfElement(ChangeRing(LE8,GF(4)),z);

IrrL:=[Restriction(Irr[1],L)];
for i in [2,3,4,6,8,7] do 
  IrrL:=AppendWithoutDuplicates(IrrL,CompositionFactors(Restriction(Irr[i],L)));
end for;
IrrM:=[Restriction(Irr[1],M)];
for i in [2,3,5,4,6,7] do 
  IrrM:=AppendWithoutDuplicates(IrrM,CompositionFactors(Restriction(Irr[i],M)));
end for;
IrrM:=OrderByDimension(IrrM);

Bool,al,be:=IsIrreducible(ChangeRing(IrrM[2],GF(4)));
IrrMTheta:=[TensorProduct(ChangeRing(IrrM[i],GF(4)),al):i in [1,3,4,7]];


function ComputeThetaEigenspace(V)
W:=DirectSum([i:i in IndecomposableSummands(Restriction(ChangeRing(V,GF(4)),M))|SequenceToSet(IdentifyFactors(i,IrrMTheta)) ne {0}]);
return ActionOfUnipotent([W],u);
end function;

function ComputeOneEigenspace(V)
W:=DirectSum([i:i in IndecomposableSummands(Restriction(V,M))|SequenceToSet(IdentifyFactors(i,IrrM)) subset {1,3,4,7}]);
return ActionOfUnipotent([W],u);
end function;



"Step 3 calculations.";

M1:=Irr[10];
for i in [6,2,7,6,8,3,2,5] do
  E,rho:=Ext(M1,Irr[i]); M1:=MaximalExtension(M1,Irr[i],E,rho); DescribeLayers(Dual(M1),Irr,labs);
end for;

for i in [2,3,4,5,6,7,8] do Dimension(Ext(M1,Irr[i])) eq 0; end for;

// Now add the 40.2d:

E,rho:=Ext(M1,Irr[9]); M1:=MaximalExtension(M1,Irr[9],E,rho); DescribeLayers(Dual(M1),Irr,labs);

"This is the radical as stated in the paper:";

M1:=RemoveFromTop(Dual(M1),Irr,Irr,[1..8]); DescribeLayers(M1,Irr,labs);
"w acts with blocks";
ActionOfUnipotent([M1],w);
// M1a will be the module M1'
M1a:=M1;
for i in [3,4,5] do Dimension(Ext(M1a,Irr[i])) eq 0; end for;
"This module is M_1'";
j:=2; E,rho:=Ext(M1a,Irr[j]); M1a:=MaximalExtension(M1a,Irr[j],E,rho); DescribeLayers(Dual(M1a),Irr,labs);
// M1b will be the module M1''

M1b:=Dual(M1a);
for i in [4,5] do
  E,rho:=Ext(M1b,Irr[i]); M1b:=MaximalExtension(M1b,Irr[i],E,rho); DescribeLayers(M1b,Irr,labs);
end for;
"That is the module M_1''";
"w acts with blocks";
ActionOfUnipotent([M1b],w);

// Now make M2

M2:=Dual(RemoveFromTop(Dual(M1b),Irr,Irr,[4,5]));
j:=2; E,rho:=Ext(M2,Irr[j]); M2:=MaximalExtension(M2,Irr[j],E,rho); 
"M2 is the module";
DescribeLayers(M2,Irr,labs);
"and w acts with blocks";
ActionOfUnipotent([M2],w);

// M2a and M2b are the result of adding 10 and 10d

j:=5; E,rho:=Ext(M2,Irr[j]); M2a:=MaximalExtension(M2,Irr[j],E,rho); 
j:=4; E,rho:=Ext(M2,Irr[j]); M2b:=MaximalExtension(M2,Irr[j],E,rho); 
"The two extensions are";
DescribeLayers(Dual(M2a),Irr,labs);
DescribeLayers(Dual(M2b),Irr,labs);

"Actions of v on former module and 10 are";
ActionOfUnipotent([M2a,Irr[5]],v);

j:=4; E,rho:=Ext(Dual(M2b),Irr[j]); M2b2:=MaximalExtension(Dual(M2b),Irr[j],E,rho);
"Action of w on extension by 10s";
DescribeLayers(Dual(M2b2),Irr,labs);

"Actions of w on latter module and extension by 10s";
ActionOfUnipotent([M2b,M2b2],w);



"Step 4 calculations.";

N1:=Irr[8];
j:=6; E,rho:=Ext(N1,Irr[j]); N1:=MaximalExtension(N1,Irr[j],E,rho); DescribeLayers(Dual(N1),Irr,labs);
j:=7; E,rho:=Ext(N1,Irr[j]); N1:=MaximalExtension(N1,Irr[j],E,rho); DescribeLayers(Dual(N1),Irr,labs); 

for i in [9,10] do Dimension(Ext(N1,Irr[i])) eq 0; end for;

N2:=DirectSum([N1,Irr[9],Irr[10],Irr[6]]);
"theta-eigenspace of N2 is";
ComputeThetaEigenspace(N2);



"Step 5 calculations.";

A1:=Irr[8];

for j in [6,3,5,9,4,2,3] do
  E,rho:=Ext(A1,Irr[j]); A1:=MaximalExtension(A1,Irr[j],E,rho); DescribeLayers(Dual(A1),Irr,labs);
end for;

for i in [2,3,4,5,6,9,10] do Dimension(Ext(A1,Irr[i])) eq 0; end for;

j:=7; E,rho:=Ext(A1,Irr[j]); A1:=MaximalExtension(A1,Irr[j],E,rho);

"The radical has structure:";
DescribeLayers(Dual(A1),Irr,labs);
"and the quotient by the socle is the sum of modules";
for i in IndecomposableSummands(JacobsonRadical(A1)) do DescribeLayers(i,Irr,labs); end for;

SubsA1:=Submodules(Dual(A1):CodimensionLimit:=40);
for i in SubsA1 do if(6 in IdentifyFactors(Socle(Dual(i)),Irr)) then A1a:=Dual(i); break i; end if; end for;

"The next radical in the paper is as follows:";
j:=6; E,rho:=Ext(A1a,Irr[j]); A1a:=MaximalExtension(A1a,Irr[j],E,rho); DescribeLayers(Dual(A1a),Irr,labs);

for i in [2,3,6,9,10] do Dimension(Ext(A1a,Irr[i])) eq 0; end for;

A2:=RemoveFromTop(Dual(A1),Irr,Irr,[2,8,10]);

for j in [9,10,3] do
  E,rho:=Ext(A2,Irr[j]); A2:=MaximalExtension(A2,Irr[j],E,rho); DescribeLayers(Dual(A2),Irr,labs);
end for;

for i in [2,3,9,10] do Dimension(Ext(A2,Irr[i])) eq 0; end for;



"We next set up the big radicals we need in Step 7.";
"We will use a submodule of them for Step 6.";

B1:=Irr[9];

for j in [6,3] do
  E,rho:=Ext(B1,Irr[j]); B1:=MaximalExtension(B1,Irr[j],E,rho); DescribeLayers(Dual(B1),Irr,labs);
end for;

for i in [2,3,4,5,6,7,10] do Dimension(Ext(B1,Irr[i])) eq 0; end for;

B2:=Dual(A1);
B2:=RemoveFromTop(B2,Irr,Irr,[2,3,4,5,6,7,9]);
DescribeLayers(B2,Irr,labs);

B3:=Irr[10];
for j in [6,2] do
  E,rho:=Ext(B3,Irr[j]); B3:=MaximalExtension(B3,Irr[j],E,rho); DescribeLayers(Dual(B3),Irr,labs);
end for;

for i in [2,3,4,5,6] do Dimension(Ext(B3,Irr[i])) eq 0; end for;

j:=7; E,rho:=Ext(B3,Irr[j]); B3:=MaximalExtension(B3,Irr[j],E,rho);
DescribeLayers(Dual(B3),Irr,labs); DescribeLayers(B3,Irr,labs);

B3a:=RemoveFromTop(Dual(B3),Irr,Irr,[6]);
B3a:=DirectSum(B3a,Dual(B3a)); // This is the final possibility for M1 in this step

"Step 6 calculations.";

"The unique extension of M1 by a 5^\\pm top and bottom";
j:=2; E,rho:=Ext(B3a,Irr[j]); B4:=Dual(MaximalExtension(B3a,Irr[j],E,rho));
j:=2; E,rho:=Ext(B4,Irr[j]); B4:=MaximalExtension(B4,Irr[j],E,rho);
for i in IndecomposableSummands(B4) do DescribeLayers(i,Irr,labs); end for;
"Proof that the 10^* remains a submodule";
j:=5; E,rho:=Ext(B4,Irr[j]); B4:=MaximalExtension(B4,Irr[j],E,rho);
DescribeLayers(B4,Irr,labs);

"The case there M2 is M1+10+10^*.";

XB4:=IndecomposableSummands(B4);
ToAdd1:=XB4[1];
ToAdd2:=BuildLargestModule(Irr[2],[Irr[5]]);

B4a:=DirectSum([ToAdd1,Irr[5]]);
B4c:=Dual(DirectSum([ToAdd1,ToAdd2]));
SubsB4c:=Submodules(B4c:CodimensionLimit:=5);
B4c:=Dual([i:i in SubsB4c|not(IsDecomposable(i))][1]);
"First case is (plus dual):";
for i in IndecomposableSummands(B4a) do DescribeLayers(i,Irr,labs); end for;
"Second case is (plus dual):";
DescribeLayers(B4c,Irr,labs);

"Action on 1-eigenspace of the two modules";
ComputeOneEigenspace(DirectSum(B4a,B4a)),ComputeOneEigenspace(DirectSum(B4c,B4c));

"Add the 24s on to the bottom";
j:=6; E,rho:=Ext(B4a,Irr[j]); B4a1:=MaximalExtension(B4a,Irr[j],E,rho);
j:=6; E,rho:=Ext(B4c,Irr[j]); B4c1:=MaximalExtension(B4c,Irr[j],E,rho);
j:=6; E,rho:=Ext(Dual(B4a),Irr[j]); B4a1d:=MaximalExtension(Dual(B4a),Irr[j],E,rho);
j:=6; E,rho:=Ext(Dual(B4c),Irr[j]); B4c1d:=MaximalExtension(Dual(B4c),Irr[j],E,rho);
B4a2:=DirectSum(B4a1,B4a1d);
B4c2:=DirectSum(B4c1,B4c1d);
"The theta-eigenspaces:";
ComputeThetaEigenspace(B4a2),ComputeThetaEigenspace(B4c2);

"Now move on to 10/10^* and 10^*/10.";

// 10/10* and 10*/10 use exactly the same commands as above, with a slightly different ToAdd.
// We do the case where there are two summands first. First build 10/10^* and 10^*/10.

TT1:=BuildLargestModule(Irr[4],[Irr[5]]);
TT2:=BuildLargestModule(Irr[5],[Irr[4]]);

"The direct sum case for M_2'. We have two options depending on whether it's 10/10^* or 10^*/10";

B4aq:=DirectSum([ToAdd1,Dual(ToAdd1),TT1]);
B4ar:=DirectSum([ToAdd1,Dual(ToAdd1),TT2]);
"Action on 1-eigenspace of the two modules";
ComputeOneEigenspace(B4aq),ComputeOneEigenspace(B4ar);
B4a1q1:=IndecomposableSummands(B4a1)[2];
B4a1q2:=IndecomposableSummands(B4a1d)[2];
B4a1q:=DirectSum([B4a1q1,B4a1q2,TT1]);
B4a1r:=DirectSum([B4a1q1,B4a1q2,TT2]);
"The theta-eigenspaces:";
ComputeThetaEigenspace(B4a1q),ComputeThetaEigenspace(B4a1r);

"The non-direct sum case for M_2'. Again, two cases.";

j:=2; E,rho:=Ext(TT1,Irr[j]); TT1a:=MaximalExtension(TT1,Irr[j],E,rho);
j:=2; E,rho:=Ext(TT2,Irr[j]); TT2a:=MaximalExtension(TT2,Irr[j],E,rho);

B4as:=DirectSum(ToAdd1,TT1a);
B4at:=DirectSum(ToAdd1,TT2a);
SubsB4as:=Submodules(Dual(B4as):CodimensionLimit:=5);
SubsB4at:=Submodules(Dual(B4at):CodimensionLimit:=5);
B4as:=Dual([i:i in SubsB4as|not(IsDecomposable(i))][1]);
B4at:=Dual([i:i in SubsB4at|not(IsDecomposable(i))][1]);
B4as:=DirectSum(B4as,Dual(ToAdd1));
B4at:=DirectSum(B4at,Dual(ToAdd1));
j:=2; E,rho:=Ext(Dual(B4as),Irr[j]); B4as:=MaximalExtension(Dual(B4as),Irr[j],E,rho);
j:=2; E,rho:=Ext(Dual(B4at),Irr[j]); B4at:=MaximalExtension(Dual(B4at),Irr[j],E,rho);
SubsB4as:=Submodules(Dual(B4as):CodimensionLimit:=5);
SubsB4at:=Submodules(Dual(B4at):CodimensionLimit:=5);
B4as:=Dual([i:i in SubsB4as|not(IsDecomposable(i))][1]);
B4at:=Dual([i:i in SubsB4at|not(IsDecomposable(i))][1]);
"Action on 1-eigenspace of the two modules";
ComputeOneEigenspace(B4as),ComputeOneEigenspace(B4at);

j:=6; E,rho:=Ext(B4as,Irr[j]); B4a1s:=MaximalExtension(B4as,Irr[j],E,rho);
j:=6; E,rho:=Ext(B4at,Irr[j]); B4a1t:=MaximalExtension(B4at,Irr[j],E,rho);
"The theta-eigenspaces:";
ComputeThetaEigenspace(B4a1s),ComputeThetaEigenspace(B4a1t);

"Now the case where M_1 is not a summand of M_2.";

Dimension(Ext(B3a,Irr[4])) eq 0;

j:=5; E,rho:=Ext(B3a,Irr[j]); B4:=MaximalExtension(B3a,Irr[j],E,rho);
j:=2; E,rho:=Ext(B4,Irr[j]); B41:=MaximalExtension(B4,Irr[j],E,rho);
"The extensions with 5.";
for i in IndecomposableSummands(B41) do DescribeLayers(i,Irr,labs); end for;

SubsB41:=Prune(Submodules(Dual(B41):CodimensionLimit:=5));
SubsB41:=[Dual(i):i in SubsB41|not(8 in IdentifyFactors(i/JacobsonRadical(i),Irr))];

B41a:=[i:i in SubsB41|IsDecomposable(i)][1];
B41b:=[i:i in SubsB41|not(IsIsomorphic(i,B41a))][1];

j:=2; E,rho:=Ext(Dual(B41a),Irr[j]); B41a1:=MaximalExtension(Dual(B41a),Irr[j],E,rho);
"Adding on the 5s.";
for i in IndecomposableSummands(Dual(B41a1)) do DescribeLayers(i,Irr,labs); end for;
j:=5; E,rho:=Ext(B41a1,Irr[j]); B41a2:=MaximalExtension(B41a1,Irr[j],E,rho);
"And the 10s.";
for i in IndecomposableSummands(Dual(B41a2)) do DescribeLayers(i,Irr,labs); end for;

SubsB41a2:=Prune(Prune(Submodules(Dual(B41a2):CodimensionLimit:=10))); // Remove the module itself and the submodule of codim 5.
SubsB41a2:=[i:i in SubsB41a2|IsSelfDual(i)];
"The 1-eigenspace for M_2''";
[ComputeOneEigenspace(i):i in SubsB41a2];

SubsB41a3:=[];
for i in SubsB41a2 do
  j:=6; E,rho:=Ext(i,Irr[j]); Append(~SubsB41a3,MaximalExtension(i,Irr[j],E,rho));
end for;

"The 1-eigenspace for M_2'', extended by the 24s";
[ComputeOneEigenspace(i):i in SubsB41a3];
Subs2B41a3:=[Prune(Prune(Prune(Prune(Submodules(Dual(i):CodimensionLimit:=24))))):i in SubsB41a3];

// Now obtain the precise one we need in both cases.
SubsB41a4:=[];
SomeEigs:=[ComputeOneEigenspace(i):i in Subs2B41a3[1]];
for i in [1..7] do if("$4^{18}" in SomeEigs[i]) then Append(~SubsB41a4,Subs2B41a3[1,i]); end if; end for;
SomeEigs:=[ComputeOneEigenspace(i):i in Subs2B41a3[2]];
for i in [1..7] do if("$4^{18}" in SomeEigs[i]) then Append(~SubsB41a4,Subs2B41a3[2,i]); end if; end for;
"The unique extension by two 24s with no extra blocks has 1-eigenspace:";
[ComputeOneEigenspace(i):i in SubsB41a4];
"and theta-eigenspace";
[ComputeThetaEigenspace(i):i in SubsB41a4];

"So now we have to use the diagonal 5 for M_1'."; //This was B41b above.

"The 1-eigenspace of M_2' is";
ComputeOneEigenspace(B41b);
"The 1-eigenspaces of 24 and 5^*/10 are";
ComputeOneEigenspace(Irr[6]); ComputeOneEigenspace(BuildLargestModule(Irr[4],[Irr[3]]));

"Compute the extensions of M_1 and M_2' by 24.";
j:=6; E,rho:=Ext(B3a,Irr[j]); E; Test1:=MaximalExtension(B3a,Irr[j],E,rho); DescribeLayers(Dual(Test1),Irr,labs);
j:=6; E,rho:=Ext(Dual(B41b),Irr[j]); E; Test2:=MaximalExtension(Dual(B41b),Irr[j],E,rho); DescribeLayers(Dual(Test2),Irr,labs);
Test1a:=SocleSeries(Dual(Test1))[3]; // Remove the errant 24.
XTest1a:=IndecomposableSummands(Test1a);
// I want to fix which summand is which in this pair.
if(IdentifyModule(Socle(XTest1a[1]),Irr) eq 9) then Reverse(~XTest1a); end if;

// Make the three modules with the 24s

B5a:=DirectSum([XTest1a[1],Dual(XTest1a[1])]);
B5b:=DirectSum([XTest1a[2],Dual(XTest1a[2])]);

"First possibility for M_3. Adding on 10^* gives:";

j:=5; E,rho:=Ext(B5a,Irr[j]); E; B5a1:=MaximalExtension(B5a,Irr[j],E,rho);
for i in IndecomposableSummands(B5a1) do DescribeLayers(i,Irr,labs); end for;
"with theta-eigenspace";
ComputeThetaEigenspace(B5a1);

"Second possibility for M_3. 1-eigenspace is";
ComputeOneEigenspace(B5b);

// Third case:

SubsTest1a:=Submodules(Test1a:CodimensionLimit:=24);
B5c:=[i:i in SubsTest1a|not(IsDecomposable(i))][1];

"Third possibility for M_3. 1-eigenspace is";
ComputeOneEigenspace(B5c);



"Calculations for Step 7.";

// On to Step 7. We have already checked the first few claims:
"First radical:";
DescribeLayers(Dual(B1),Irr,labs);

"Second and third modules:";
DescribeLayers(B2,Irr,labs);
DescribeLayers(RemoveFromTop(Dual(B3),Irr,Irr,[6]),Irr,labs);


// On to the later ones. Start with where 40.1d lies in fourth layer.

B6a:=RemoveFromTop(B2,Irr,Irr,[8,3,6]);

"The 24 cannot stop 14 falling into the socle on restriction to L:";
E,rho:=Ext(B6a,Irr[6]); B6atest:=MaximalExtension(B6a,Irr[6],E,rho);
DescribeLayers(B6atest,Irr,labs); DescribeLayers(Restriction(B6atest,L),IrrL,labsL);

"And clearly neither can 5:";
E,rho:=Ext(B6a,Irr[2]); B6a1:=MaximalExtension(B6a,Irr[2],E,rho); DescribeLayers(B6a1,Irr,labs);

"Now for the diagonal submodule possibility for 40_1^*.";

B7:=DirectSum(B2,RemoveFromTop(Dual(B3),Irr,Irr,[6]));
// This is the sum of the second and third radicals.
// This has three copies of 40_1^* on the top. We can take the one from the fourth layer off.
B7:=SocleSeries(B7)[3];
//and then remove any summands other than 40_1^*.
B7:=RemoveFromTop(B7,Irr,Irr,[2,3,4,5]);

SubsB7:=Prune(Submodules(B7:CodimensionLimit:=40));
for i in SubsB7 do if(not(IsDecomposable(i))) then B7:=i; end if; end for;
"It has structure";
DescribeLayers(B7,Irr,labs);
"This is its theta-eigenspace.";
ComputeThetaEigenspace(B7);



"Step 8 calculations.";

E,rho:=Ext(Irr[6],Irr[7]); C1a:=Extension(Irr[6],Irr[7],E.1,rho);
E,rho:=Ext(Irr[6],Irr[8]); C1b:=Extension(Irr[6],Irr[8],E.1,rho);
E,rho:=Ext(Irr[6],Irr[9]); C1c:=Extension(Irr[6],Irr[9],E.1,rho);
E,rho:=Ext(Irr[6],Irr[10]); C1d:=Extension(Irr[6],Irr[10],E.1,rho);
SubsC1ab:=Submodules(DirectSum(C1a,C1b):CodimensionLimit:=24);
SubsC1cd:=Submodules(DirectSum(C1c,C1d):CodimensionLimit:=24);
C1ab:=[i:i in SubsC1ab|not(IsDecomposable(i))][1];
C1cd:=[i:i in SubsC1cd|not(IsDecomposable(i))][1];
SubsC1:=Submodules(DirectSum(C1ab,C1cd):CodimensionLimit:=24);
C1:=[i:i in SubsC1|not(IsDecomposable(i))][1];
E,rho:=Ext(C1,Irr[6]); C1:=MaximalExtension(C1,Irr[6],E,rho);

"The restriction of M1 to L has structure:";
DescribeLayers(Restriction(C1,L),IrrL,labsL);

C2:=C1;
for j in [2,3] do
  E,rho:=Ext(C2,Irr[j]); C2:=MaximalExtension(C2,Irr[j],E,rho); DescribeLayers(Dual(C2),Irr,labs);
end for;
C2:=Dual(C2);
for j in [2,3] do
  E,rho:=Ext(C2,Irr[j]); C2:=MaximalExtension(C2,Irr[j],E,rho); DescribeLayers(Dual(C2),Irr,labs);
end for;
"The module above is the 'unique such module' from the paper. u acts on it with blocks";
ActionOfUnipotent([C2],u);
"1-eigenspace of u on M2 and 10";
ComputeOneEigenspace(C2); ComputeOneEigenspace(Irr[4]);

"Case where 10+10^* is a summand of M3";
C3a:=DirectSum([C1,Irr[4],Irr[5]]);
ComputeOneEigenspace(C3a); ComputeThetaEigenspace(C3a);

"Now the extension with 5, and its theta-eigenspace";
j:=3; E,rho:=Ext(C3a,Irr[j]); C3a:=MaximalExtension(C3a,Irr[j],E,rho); DescribeLayers(Dual(C3a),Irr,labs);
ComputeThetaEigenspace(C3a);

"Case where 10^*/10 is a summand";

C3b:=DirectSum([C1,BuildLargestModule(Irr[4],[Irr[5]])]);
j:=3; E,rho:=Ext(C3b,Irr[j]); C3b1:=MaximalExtension(C3b,Irr[j],E,rho); DescribeLayers(Dual(C3b1),Irr,labs);
SubsC3b:=Submodules(Dual(C3b1):CodimensionLimit:=5);

SomeEigs:=[ComputeOneEigenspace(i):i in SubsC3b];
"We choose the extension with sixteen blocks of size 4.";
for i in [1..3] do if("$4^{16}" in SomeEigs[i]) then C3b2:=SubsC3b[i]; end if; end for;
"And it has 1-eigenspace";
ComputeOneEigenspace(C3b2);
"and summands";
for i in IndecomposableSummands(C3b2) do DescribeLayers(i,Irr,labs); end for;

"Now the other option for M3.";
j:=4; E,rho:=Ext(C1,Irr[j]); C3c:=MaximalExtension(C1,Irr[j],E,rho); DescribeLayers(C3c,Irr,labs);

"The restrictions of M3 to M. First the 1-eigenspace and second the theta-eigenspace.";
for i in IndecomposableSummands(Restriction(C3c,M)) do if(SequenceToSet(IdentifyFactors(i,IrrM)) subset {1,3,4,7}) then DescribeLayers(i,IrrM,labsM); end if; end for;

for i in IndecomposableSummands(Restriction(C3c,M)) do if(SequenceToSet(IdentifyFactors(i,IrrM)) subset {2,5,6,8}) then DescribeLayers(i,IrrM,labsM); end if; end for;

j:=4; E,rho:=Ext(Dual(C3c),Irr[j]); C3d:=MaximalExtension(Dual(C3c),Irr[j],E,rho);

"Compare with the action of M on the extended module as given below.";

for i in IndecomposableSummands(Restriction(Dual(C3d),M)) do if(SequenceToSet(IdentifyFactors(i,IrrM)) subset {1,3,4,7}) then DescribeLayers(i,IrrM,labsM); end if; end for;

for i in IndecomposableSummands(Restriction(Dual(C3d),M)) do if(SequenceToSet(IdentifyFactors(i,IrrM)) subset {2,5,6,8}) then DescribeLayers(i,IrrM,labsM); end if; end for;

"Now place both 5^*s underneath M_1', and compute its theta-eigenspace.";
j:=3; E,rho:=Ext(C3c,Irr[j]); C3e:=MaximalExtension(C3c,Irr[j],E,rho); DescribeLayers(C3e,Irr,labs);
ComputeThetaEigenspace(C3e);

"Placing the two 5^*s on top of M_1', along with the 10^*s, yields";
j:=2; E,rho:=Ext(C3d,Irr[j]); C3f:=MaximalExtension(C3d,Irr[j],E,rho);

C3f:=Dual(C3f);
DescribeLayers(C3f,Irr,labs);
"and it has 1-eigenspace";
ComputeOneEigenspace(C3d);



"Step 9 calculations.";

j:=7; E,rho:=Ext(Irr[j],Irr[6]); D401:=MaximalExtension(Irr[j],Irr[6],E,rho);
j:=8; E,rho:=Ext(Irr[j],Irr[6]); D401d:=MaximalExtension(Irr[j],Irr[6],E,rho);
j:=9; E,rho:=Ext(Irr[j],Irr[6]); D402:=MaximalExtension(Irr[j],Irr[6],E,rho);

D1a:=DirectSum([D401,Dual(D401),Irr[9],Irr[10]]);
D1b:=DirectSum([D402,Dual(D402),Irr[7],Irr[8]]);

SubsD4012:=Submodules(Dual(DirectSum(D401,D402)):CodimensionLimit:=24);
SubsD401d2:=Submodules(Dual(DirectSum(D401d,D402)):CodimensionLimit:=24);
D4012:=Dual([i:i in SubsD4012|not(IsDecomposable(i))][1]);
D401d2:=Dual([i:i in SubsD401d2|not(IsDecomposable(i))][1]);
D1c:=DirectSum(D4012,Dual(D4012));
D1d:=DirectSum(D401d2,Dual(D401d2));

"First theta-eigenspace is";
ComputeThetaEigenspace(D1a);

"Now the 1- and theta-eigenspaces of the three modules";
ComputeOneEigenspace(D1b),ComputeOneEigenspace(D1c),ComputeOneEigenspace(D1d);
ComputeThetaEigenspace(D1b),ComputeThetaEigenspace(D1c),ComputeThetaEigenspace(D1d);

"The action of w on the three modules";
ActionOfUnipotent([D1b,D1c,D1d],w);

"Now the cases where M_1 is a summand of M_2";


Add1:=DirectSum([Irr[4],Irr[5]]);
Add2:=BuildLargestModule(Irr[4],[Irr[5]]);
Add3:=BuildLargestModule(Irr[5],[Irr[4]]);

D1bA1:=DirectSum(D1b,Add1);
D1cA1:=DirectSum(D1c,Add1);
D1dA1:=DirectSum(D1d,Add1);
"The theta-eigenspaces for 10+10^* case";
ComputeThetaEigenspace(D1b),ComputeThetaEigenspace(D1c),ComputeThetaEigenspace(D1d);
j:=2; E,rho:=Ext(D1bA1,Irr[j]); D1bA1a:=MaximalExtension(D1bA1,Irr[j],E,rho);
j:=2; E,rho:=Ext(D1cA1,Irr[j]); D1cA1a:=MaximalExtension(D1cA1,Irr[j],E,rho);
j:=2; E,rho:=Ext(D1dA1,Irr[j]); D1dA1a:=MaximalExtension(D1dA1,Irr[j],E,rho);
"And with the 5s on the bottom";
ComputeThetaEigenspace(D1bA1a),ComputeThetaEigenspace(D1cA1a),ComputeThetaEigenspace(D1dA1a);

"Next the 10/10^* and 10^*/10 cases";
D1bA2:=DirectSum(D1b,Add2);
D1cA2:=DirectSum(D1c,Add2);
D1dA2:=DirectSum(D1d,Add2);
D1bA3:=DirectSum(D1b,Add3);
D1cA3:=DirectSum(D1c,Add3);
D1dA3:=DirectSum(D1d,Add3);
"The 1-eigenspace of all six modules M_2 is";
{ComputeOneEigenspace(i):i in [D1bA2,D1cA2,D1dA2,D1bA3,D1cA3,D1dA3]};
"and the theta-eigenspace is";
{ComputeThetaEigenspace(i):i in [D1bA2,D1cA2,D1dA2,D1bA3,D1cA3,D1dA3]};



j:=2; E,rho:=Ext(D1bA2,Irr[j]); D1bA2a:=MaximalExtension(D1bA2,Irr[j],E,rho);
j:=2; E,rho:=Ext(D1cA2,Irr[j]); D1cA2a:=MaximalExtension(D1cA2,Irr[j],E,rho);
j:=2; E,rho:=Ext(D1dA2,Irr[j]); D1dA2a:=MaximalExtension(D1dA2,Irr[j],E,rho);
j:=2; E,rho:=Ext(D1bA3,Irr[j]); D1bA3a:=MaximalExtension(D1bA3,Irr[j],E,rho);
j:=2; E,rho:=Ext(D1cA3,Irr[j]); D1cA3a:=MaximalExtension(D1cA3,Irr[j],E,rho);
j:=2; E,rho:=Ext(D1dA3,Irr[j]); D1dA3a:=MaximalExtension(D1dA3,Irr[j],E,rho);
"The theta-eigenspace of the extensions by 5 are:";
{ComputeThetaEigenspace(i):i in [D1bA2a,D1cA2a,D1dA2a,D1bA3a,D1cA3a,D1dA3a]};


D1bA2b:=[Dual(i):i in Submodules(Dual(D1bA2a):CodimensionLimit:=5)|"4^{16}" in ComputeOneEigenspace(i)][1];
D1cA2b:=[Dual(i):i in Submodules(Dual(D1cA2a):CodimensionLimit:=5)|"4^{16}" in ComputeOneEigenspace(i)][1];
D1dA2b:=[Dual(i):i in Submodules(Dual(D1dA2a):CodimensionLimit:=5)|"4^{16}" in ComputeOneEigenspace(i)][1];
D1bA3b:=[Dual(i):i in Submodules(Dual(D1bA3a):CodimensionLimit:=5)|"4^{16}" in ComputeOneEigenspace(i)][1];
D1cA3b:=[Dual(i):i in Submodules(Dual(D1cA3a):CodimensionLimit:=5)|"4^{16}" in ComputeOneEigenspace(i)][1];
D1dA3b:=[Dual(i):i in Submodules(Dual(D1dA3a):CodimensionLimit:=5)|"4^{16}" in ComputeOneEigenspace(i)][1];

j:=2; E,rho:=Ext(Dual(D1bA2b),Irr[j]); D1bA2c:=MaximalExtension(Dual(D1bA2b),Irr[j],E,rho);
j:=2; E,rho:=Ext(Dual(D1cA2b),Irr[j]); D1cA2c:=MaximalExtension(Dual(D1cA2b),Irr[j],E,rho);
j:=2; E,rho:=Ext(Dual(D1dA2b),Irr[j]); D1dA2c:=MaximalExtension(Dual(D1dA2b),Irr[j],E,rho);
j:=2; E,rho:=Ext(Dual(D1bA3b),Irr[j]); D1bA3c:=MaximalExtension(Dual(D1bA3b),Irr[j],E,rho);
j:=2; E,rho:=Ext(Dual(D1cA3b),Irr[j]); D1cA3c:=MaximalExtension(Dual(D1cA3b),Irr[j],E,rho);
j:=2; E,rho:=Ext(Dual(D1dA3b),Irr[j]); D1dA3c:=MaximalExtension(Dual(D1dA3b),Irr[j],E,rho);

D1bA2d:=[Dual(i):i in Submodules(Dual(D1bA2c):CodimensionLimit:=5)|"4^{16}" in ComputeOneEigenspace(i)][1];
D1cA2d:=[Dual(i):i in Submodules(Dual(D1cA2c):CodimensionLimit:=5)|"4^{16}" in ComputeOneEigenspace(i)][1];
D1dA2d:=[Dual(i):i in Submodules(Dual(D1dA2c):CodimensionLimit:=5)|"4^{16}" in ComputeOneEigenspace(i)][1];
D1bA3d:=[Dual(i):i in Submodules(Dual(D1bA3c):CodimensionLimit:=5)|"4^{16}" in ComputeOneEigenspace(i)][1];
D1cA3d:=[Dual(i):i in Submodules(Dual(D1cA3c):CodimensionLimit:=5)|"4^{16}" in ComputeOneEigenspace(i)][1];
D1dA3d:=[Dual(i):i in Submodules(Dual(D1dA3c):CodimensionLimit:=5)|"4^{16}" in ComputeOneEigenspace(i)][1];
"The 1-eigenspaces of the extensions by 5 and 5^* are";
{ComputeOneEigenspace(i):i in [D1bA2d,D1cA2d,D1dA2d,D1bA3d,D1cA3d,D1dA3d]};

"The descriptions of the unique extension by the 5s and 5^*s in the six cases";
for j in [D1bA2d,D1cA2d,D1dA2d,D1bA3d,D1cA3d,D1dA3d] do
  [DescribeLayers(i,Irr,labs):i in IndecomposableSummands(j)];
end for;

"These have theta-eigenspace";
{ComputeThetaEigenspace(i):i in [D1bA2d,D1cA2d,D1dA2d,D1bA3d,D1cA3d,D1dA3d]};


j:=3; E,rho:=Ext(D1bA2d,Irr[j]); D1bA2e:=MaximalExtension(D1bA2d,Irr[j],E,rho);
j:=3; E,rho:=Ext(D1cA2d,Irr[j]); D1cA2e:=MaximalExtension(D1cA2d,Irr[j],E,rho);
j:=3; E,rho:=Ext(D1dA2d,Irr[j]); D1dA2e:=MaximalExtension(D1dA2d,Irr[j],E,rho);
j:=3; E,rho:=Ext(D1bA3d,Irr[j]); D1bA3e:=MaximalExtension(D1bA3d,Irr[j],E,rho);
j:=3; E,rho:=Ext(D1cA3d,Irr[j]); D1cA3e:=MaximalExtension(D1cA3d,Irr[j],E,rho);
j:=3; E,rho:=Ext(D1dA3d,Irr[j]); D1dA3e:=MaximalExtension(D1dA3d,Irr[j],E,rho);

//Check that there are three copies
{Dimension(i):i in [D1bA2e,D1cA2e,D1dA2e,D1bA3e,D1cA3e,D1dA3e]} eq {253};

"Putting all of the 5^*s on the bottom yields theta-eigenspaces";
{ComputeThetaEigenspace(i):i in [D1bA2e,D1cA2e,D1dA2e,D1bA3e,D1cA3e,D1dA3e]};


"Now we assume that M_1 is not a summand of M_2.";

D1s:=[D1b,D1c,D1d];
D1t:=[];

for i in D1s do
  ii:=i;
  for j in [2,3] do 
    E,rho:=Ext(ii,Irr[j]); ii:=MaximalExtension(ii,Irr[j],E,rho);
  end for;
  Append(~D1t,ii);
end for;

D1u:=[];
for i in D1t do
  for j in [4,5] do 
    E,rho:=Ext(i,Irr[j]); Append(~D1u,MaximalExtension(i,Irr[j],E,rho));
  end for;
end for;

"1- and theta-eigenspaces of M_3':";

{ComputeOneEigenspace(i):i in D1u};
{ComputeThetaEigenspace(i):i in D1u};

"Remove the 5 and 5^*s, and obtain 1-eigenspace";
{ComputeOneEigenspace(RemoveFromTop(Dual(i),Irr,Irr,[2,3])):i in D1u};
"Remove the 5^*s:";
{ComputeOneEigenspace(RemoveFromTop(Dual(i),Irr,Irr,[2])):i in D1u};
"Remove the 5:";
{ComputeOneEigenspace(RemoveFromTop(Dual(i),Irr,Irr,[3])):i in D1u};


"Thus the {5,10}-radical is not semisimple.";

"The 1-eigenspace of 10^*/5 is";
ComputeOneEigenspace(BuildLargestModule(Irr[2],[Irr[5]]));

"First case, two 5s in socle";

E,rho:=Ext(D1b,Irr[5]); D3bq:=MaximalExtension(D1b,Irr[5],E,rho);
E,rho:=Ext(D1c,Irr[5]); D3cq:=MaximalExtension(D1c,Irr[5],E,rho);
E,rho:=Ext(D1d,Irr[5]); D3dq:=MaximalExtension(D1d,Irr[5],E,rho);

E,rho:=Ext(D3bq,Irr[2]); D3b1q:=MaximalExtension(D3bq,Irr[2],E,rho);
E,rho:=Ext(D3cq,Irr[2]); D3c1q:=MaximalExtension(D3cq,Irr[2],E,rho);
E,rho:=Ext(D3dq,Irr[2]); D3d1q:=MaximalExtension(D3dq,Irr[2],E,rho);

E,rho:=Ext(Dual(D3b1q),Irr[5]); D3b2q:=MaximalExtension(Dual(D3b1q),Irr[5],E,rho);
E,rho:=Ext(Dual(D3c1q),Irr[5]); D3c2q:=MaximalExtension(Dual(D3c1q),Irr[5],E,rho);
E,rho:=Ext(Dual(D3d1q),Irr[5]); D3d2q:=MaximalExtension(Dual(D3d1q),Irr[5],E,rho);

E,rho:=Ext(D3b2q,Irr[2]); D3b3q:=Dual(MaximalExtension(D3b2q,Irr[2],E,rho));
E,rho:=Ext(D3c2q,Irr[2]); D3c3q:=Dual(MaximalExtension(D3c2q,Irr[2],E,rho));
E,rho:=Ext(D3d2q,Irr[2]); D3d3q:=Dual(MaximalExtension(D3d2q,Irr[2],E,rho));
"The 1- and theta-eigenspaces of the modules built up as in the paper are";
{ComputeOneEigenspace(i):i in [D3b3q,D3c3q,D3d3q]};
{ComputeThetaEigenspace(i):i in [D3b3q,D3c3q,D3d3q]};

"Second case, the radical is 5^*+(10/5^*)";

E,rho:=Ext(D1b,Irr[4]); D3br:=MaximalExtension(D1b,Irr[4],E,rho);
E,rho:=Ext(D1c,Irr[4]); D3cr:=MaximalExtension(D1c,Irr[4],E,rho);
E,rho:=Ext(D1d,Irr[4]); D3dr:=MaximalExtension(D1d,Irr[4],E,rho);

E,rho:=Ext(D3br,Irr[3]); D3b1r:=MaximalExtension(D3br,Irr[3],E,rho);
E,rho:=Ext(D3cr,Irr[3]); D3c1r:=MaximalExtension(D3cr,Irr[3],E,rho);
E,rho:=Ext(D3dr,Irr[3]); D3d1r:=MaximalExtension(D3dr,Irr[3],E,rho);

DescribeLayers(D3d1r,Irr,labs);

E,rho:=Ext(Dual(D3b1r),Irr[4]); D3b2r:=MaximalExtension(Dual(D3b1r),Irr[4],E,rho);
E,rho:=Ext(Dual(D3c1r),Irr[4]); D3c2r:=MaximalExtension(Dual(D3c1r),Irr[4],E,rho);
//E,rho:=Ext(Dual(D3d1r),Irr[4]); D3d2r:=MaximalExtension(Dual(D3d1r),Irr[4],E,rho);

E,rho:=Ext(D3b2r,Irr[3]); D3b3r:=Dual(MaximalExtension(D3b2r,Irr[3],E,rho));
E,rho:=Ext(D3c2r,Irr[3]); D3c3r:=Dual(MaximalExtension(D3c2r,Irr[3],E,rho));
//E,rho:=Ext(D3d2r,Irr[3]); D3d3r:=Dual(MaximalExtension(D3d2r,Irr[3],E,rho));


"The 1- and theta-eigenspaces of the modules built up as in the paper are";
{ComputeOneEigenspace(D3b3r),ComputeOneEigenspace(D3c3r)};
{ComputeThetaEigenspace(D3b3r),ComputeThetaEigenspace(D3c3r)};
"Removing the 5^*s from the bottom yields 1-eigenspace";
ComputeOneEigenspace(RemoveFromTop(Dual(D3b3r),Irr,Irr,[2]));
ComputeOneEigenspace(RemoveFromTop(Dual(D3c3r),Irr,Irr,[2]));

"Third case, the radical is 5+(10/5^*)";

E,rho:=Ext(D1b,Irr[2]); D3bs:=MaximalExtension(D1b,Irr[2],E,rho);
E,rho:=Ext(D1c,Irr[2]); D3cs:=MaximalExtension(D1c,Irr[2],E,rho);
E,rho:=Ext(D1d,Irr[2]); D3ds:=MaximalExtension(D1d,Irr[2],E,rho);

E,rho:=Ext(Dual(D3bs),Irr[2]); D3b1s:=MaximalExtension(Dual(D3bs),Irr[2],E,rho);
E,rho:=Ext(Dual(D3cs),Irr[2]); D3c1s:=MaximalExtension(Dual(D3cs),Irr[2],E,rho);
E,rho:=Ext(Dual(D3ds),Irr[2]); D3d1s:=MaximalExtension(Dual(D3ds),Irr[2],E,rho);

// Don't need to dualize again as these modules are self-dual.

E,rho:=Ext(D3b1s,Irr[4]); D3b2s:=(MaximalExtension(D3b1s,Irr[4],E,rho));
E,rho:=Ext(D3c1s,Irr[4]); D3c2s:=(MaximalExtension(D3c1s,Irr[4],E,rho));
E,rho:=Ext(D3d1s,Irr[4]); D3d2s:=(MaximalExtension(D3d1s,Irr[4],E,rho));

E,rho:=Ext(D3b2s,Irr[3]); D3b3s:=(MaximalExtension(D3b2s,Irr[3],E,rho));
E,rho:=Ext(D3c2s,Irr[3]); D3c3s:=(MaximalExtension(D3c2s,Irr[3],E,rho));
E,rho:=Ext(D3d2s,Irr[3]); D3d3s:=(MaximalExtension(D3d2s,Irr[3],E,rho));

"Note that both 10s are submodules of the third case:";
DescribeLayers(D3d3s,Irr,labs);

D3b4s:=RemoveFromTop(Dual(D3b3s),Irr,Irr,[5]);
D3c4s:=RemoveFromTop(Dual(D3c3s),Irr,Irr,[5]);
//D3d4s:=RemoveFromTop(Dual(D3d3s),Irr,Irr,[5]);

"The theta-eigenspaces of the two modules are";
{ComputeThetaEigenspace(D3b4s),ComputeThetaEigenspace(D3c4s)};

"Fourth case, the radical is 5^*+(10^*/5)";

E,rho:=Ext(D1b,Irr[3]); D3bt:=MaximalExtension(D1b,Irr[3],E,rho);
E,rho:=Ext(D1c,Irr[3]); D3ct:=MaximalExtension(D1c,Irr[3],E,rho);
E,rho:=Ext(D1d,Irr[3]); D3dt:=MaximalExtension(D1d,Irr[3],E,rho);

E,rho:=Ext(Dual(D3bt),Irr[3]); D3b1t:=MaximalExtension(Dual(D3bt),Irr[3],E,rho);
E,rho:=Ext(Dual(D3ct),Irr[3]); D3c1t:=MaximalExtension(Dual(D3ct),Irr[3],E,rho);
E,rho:=Ext(Dual(D3dt),Irr[3]); D3d1t:=MaximalExtension(Dual(D3dt),Irr[3],E,rho);

E,rho:=Ext(D3b1t,Irr[5]); D3b2t:=(MaximalExtension(D3b1t,Irr[5],E,rho));
E,rho:=Ext(D3c1t,Irr[5]); D3c2t:=(MaximalExtension(D3c1t,Irr[5],E,rho));
E,rho:=Ext(D3d1t,Irr[5]); D3d2t:=(MaximalExtension(D3d1t,Irr[5],E,rho));

E,rho:=Ext(D3b2t,Irr[2]); D3b3t:=(MaximalExtension(D3b2t,Irr[2],E,rho));
E,rho:=Ext(D3c2t,Irr[2]); D3c3t:=(MaximalExtension(D3c2t,Irr[2],E,rho));
E,rho:=Ext(D3d2t,Irr[2]); D3d3t:=(MaximalExtension(D3d2t,Irr[2],E,rho));

D3b4t:=RemoveFromTop(Dual(D3b3t),Irr,Irr,[4]); 
D3c4t:=RemoveFromTop(Dual(D3c3t),Irr,Irr,[4]); 
D3d4t:=RemoveFromTop(Dual(D3d3t),Irr,Irr,[4]); 

"The theta-eigenspaces of the modules are";
{ComputeThetaEigenspace(i):i in [D3b4t,D3c4t,D3d4t]};


"Step 10 calculations.";

// In Step 9 we already proved this assertion. We added 5^\pm first, then the 10s.

"Note that all 5 and 5^* appear in the socle of these modules.";
[DescribeLayers(i,Irr,labs):i in D1u];

// And that's it, done!
"Everything is proved.";