load "Subgroups/E8/ComputeFactorsE8.txt";

function OrbitReps(X,frob)

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

G2:=PSU(3,4);
CC:=ConjugacyClasses(G2);
Irr:=IrreducibleModules(G2,GF(16));
S31:=Irr[2];
S31d:=Dual(S31);
S32:=GModule(G2,FrobeniusImage(S31,1));
S32d:=Dual(S32);
S912:=TensorProduct(S31,S32);
S912d:=Dual(S912);

// We call it S921 rather than S912 bar here.

S921:=TensorProduct(S31,S32d);
S921d:=Dual(S921);

S81:=TensorProduct(S31,S31d); S81:=S81/Fix(S81);
S82:=GModule(G2,FrobeniusImage(S81,1));

S2412:=TensorProduct(S31,S82);
S2421:=TensorProduct(S32,S81);
S2412d:=Dual(S2412);
S2421d:=Dual(S2421);
S64:=TensorProduct(S81,S82);

Irr:=[Irr[1],S31,S31d,S32,S32d,S81,S82,S912,S912d,S921,S921d,S2412,S2412d,S2421,S2421d,S64];
labs:=["1","3_1","3_1^*","3_2","3_2^*","8_1","8_2","9_{1,2}","9_{1,2}^*","\\bar 9_{1,2}","\\bar 9_{1,2}^*","24_{1,2}","24_{1,2}^*","24_{2,1}","24_{2,1}^*","64"];

frob:=Sym(#Irr)![IdentifyModule(GModule(G2,FrobeniusImage(i,1)),Irr):i in Irr];

Ans3,Mods,I:=ComputeConspicuousFactorsE8(Irr,CC,Tr,[13,15],[]);
Ans3a:=[IdentifyFactors(i,Irr):i in Mods];
Ans4:=OrbitReps(Ans3a,frob);

cohom:=[i:i in [1..14]|Dimension(Ext(Irr[i],Irr[1])) eq 1];
pres:=[&+[Multiplicity(i,j):j in cohom]-Multiplicity(i,1):i in Ans4];

Ans5:=[Ans4[i]:i in [1..#Ans4]|not(1 in Ans4[i]) or pres[i] gt 0];

"The sets of factors with positive pressure or no trivial factors are";
[{*labs[j]:j in i*}:i in Ans5];

"Radical for Case 1";

P912d:=ProjectiveIndecomposableModule(S912d);
P921d:=ProjectiveIndecomposableModule(S921d);

JA1:=RemoveFromTop(P912d,Irr,Irr,[1,2,3,4,5,6,7,8,9]);
A1:=Dual(P912d/JA1);
#[i:i in CompositionFactors(A1)|Dimension(i) eq 1] eq 2;

"Radical for Case 3";

// Note that there are no 8_2s in A1, so it is also this radical.
A1a:=RemoveFromTop(A1,Irr,Irr,[1..7]);
DescribeLayers(A1a,Irr,labs);

"and the smaller module is";
A1b:=RemoveFromTop(A1a,Irr,Irr,[2..9]);
DescribeLayers(A1b,Irr,labs);

"Radicals for Case 4";

"Check the number of trivial factors.";
JA2:=RemoveFromTop(P912d,Irr,Irr,[1..5] cat [7..11] cat [14,15]);
A2:=Dual(P912d/JA2);
#[i:i in CompositionFactors(A2)|Dimension(i) eq 1] eq 5;

JA3:=RemoveFromTop(P921d,Irr,Irr,[1..5] cat [7..11] cat [14,15]);
A3:=Dual(P921d/JA3);
#[i:i in CompositionFactors(A3)|Dimension(i) eq 1] eq 3;


A2a:=RemoveFromTop(A2,Irr,Irr,[1..8] cat[10..15]);
#[i:i in CompositionFactors(A2a)|Dimension(i) eq 1] eq 3;

"and the smaller module is";
A2b:=RemoveFromTop(A2a,Irr,Irr,[2..15]);
DescribeLayers(A2b,Irr,labs);

facts:=Reverse(IdentifyLayers(A2,Irr));
"Check the locations of the trivials.";
[i:i in [1..14]|1 in facts[i]] eq [2,4,5,6,8];


"Check that there is a 1+1+1 subquotient.";

J1:=RemoveFromTop(JacobsonRadical(P912d),Irr,Irr,[1..5] cat [7,10,11,14,15]);
J2:=RemoveFromTop(J1,Irr,Irr,[8]);
J3:=RemoveFromTop(J2,Irr,Irr,[1..5] cat [7,10,11,14,15]);
J4:=RemoveFromTop(J3,Irr,Irr,[9]);
J5:=RemoveFromTop(J4,Irr,Irr,[1..5] cat [7,10,11,14,15]);
JJ:=Dual(P912d/J5);
JJ2:=RemoveFromTop(JJ,Irr,Irr,[2..15]);
Dimension(AHom(JJ2,Irr[1])) eq 3;


K2:=RemoveFromTop(J1,Irr,Irr,[9]);
K3:=RemoveFromTop(K2,Irr,Irr,[1..5] cat [7,10,11,14,15]);
K4:=RemoveFromTop(K3,Irr,Irr,[8]);
K5:=RemoveFromTop(K4,Irr,Irr,[1..5] cat [7,10,11,14,15]);
KK:=Dual(P912d/K5);
KK2:=RemoveFromTop(KK,Irr,Irr,[2..15]);
Dimension(AHom(JJ2,Irr[1])) eq 3 and Dimension(AHom(KK2,Irr[1])) eq 3;




"From now we need to introduce the subgroup L, the centralizer of an element of order 5.";

z:=CC[5,3];
L:=Centralizer(G2,z);
IrrL:=[Restriction(Irr[1],L),OrderByDimension(CompositionFactors(Restriction(Irr[2],L)))[1]];
for i in [2,3,4] do Insert(~IrrL,i+1,TensorPower(IrrL[2],i)); end for;
Append(~IrrL,SocleFactors(Restriction(S82,L))[2]);
for i in [2,3,4,5] do Append(~IrrL,TensorProduct(IrrL[i],IrrL[6])); end for;
Append(~IrrL,SocleFactors(Restriction(S81,L))[2]);
for i in [2,3,4,5] do Append(~IrrL,TensorProduct(IrrL[i],IrrL[11])); end for;
Append(~IrrL,TensorProduct(IrrL[6],IrrL[11]));
for i in [2,3,4,5] do Append(~IrrL,TensorProduct(IrrL[i],IrrL[16])); end for;


// I never got on with this numbering, so now I use the new one:
IrrL:=[IrrL[i]:i in [1,6,11,16,2,7,12,17,3,8,13,18,4,9,14,19,5,10,15,20]];

labsL:=["1","2_1","2_2","4","1a","2_1a","2_2a","4a","1b","2_1b","2_2b","4b","1c","2_1c","2_2c","4c","1d","2_1d","2_2d","4d"];

"Labels for L-modules are 1,2_1,2_2,4 (labels for SL(2,4)) followed by nothing, a, b, c or d, depending on whether it is the 1-, z-, z^2-, z^3- or z^4-eigenspace.";

Ld:=DerivedSubgroup(L);
IrrLd:=[Restriction(IrrL[i],Ld):i in [1..4]];


"Here are the summands of the restriction of modules to L:";
for i in Irr do [DescribeLayers(j,IrrL,labsL):j in IndecomposableSummands(Restriction(i,L))]; end for;

"This verifies the table.";


"Now verify the statements about the permutation module P_L.";

PL:=PermutationModule(G2,L,GF(16));
"Its structure is";
DescribeLayers(PL,Irr,labs);

JPL:=RemoveFromTop(PL,Irr,Irr,[1..5] cat [8..15]);
Dimension(PL)-Dimension(JPL) eq 1;
J2PL:=RemoveFromTop(PL,Irr,Irr,[1..6] cat [8..15]);
PLquot:=PL/J2PL;
"The quotient is";
DescribeLayers(PLquot,Irr,labs);

"Verify the statement about 1-dimensional fixed points.";

// Construct all submodules of P_L with only 9s as quotients. We then dualize them.

AllSubs:=[RemoveFromTop(PL,Irr,Irr,[1..7] cat [12..16])];

ToCheckSubs:=AllSubs;

repeat 
  MM:=ToCheckSubs[1];
  Remove(~ToCheckSubs,1);
  for j in [8,9,10,11] do
    homs:=AHom(MM,Irr[j]);
    if(Dimension(homs) gt 1) then "The claim from the paper is false"; end if;
    if(Dimension(homs) eq 0) then continue j; end if;
    MM2:=Kernel(homs.1); MM2:=RemoveFromTop(MM2,Irr,Irr,[1..7] cat [12..16]);
    if(IsIsomorphicToList(MM2,AllSubs)) then continue j;
      else Append(~AllSubs,MM2); Append(~ToCheckSubs,MM2);
    end if;
  end for;
until ToCheckSubs eq [];

Dimension(AllSubs[#AllSubs]) eq 0; Prune(~AllSubs); #AllSubs eq 34;
AllQuots:=[Dual(i):i in AllSubs];

&and[Dimension(Fix(Restriction(i,L))) eq 1:i in AllQuots];
&and[#SocleFactors(i) in {2,6}:i in AllQuots];
TestSocle:=DirectSum([Irr[i]:i in [1,6,7]]);
&and[IsIsomorphic(i/JacobsonRadical(i),TestSocle):i in AllQuots|#SocleFactors(i) eq 6];
&and[{6,7} subset IdentifyFactors(i/SocleSeries(i)[4],Irr):i in AllQuots|#SocleFactors(i) eq 6];

"Case 2";

"First check the composition factors on M(A_8).";

Poss9s:=[[4,4,1],[4,2,2,1],[4,2,3,1],[4,3,3,1],[4,2,1,1,1],[4,3,1,1,1],[4,1,1,1,1,1],
[2,2,2,2,1],[2,2,2,3,1],[2,2,3,3,1],[2,3,3,3,1],[3,3,3,3,1],
[2,2,2,1,1,1],[2,2,3,1,1,1],[2,3,3,1,1,1],[3,3,3,1,1,1],
[2,2,1,1,1,1,1],[2,3,1,1,1,1,1],[3,3,1,1,1,1,1],
[2,1,1,1,1,1,1,1],[3,1,1,1,1,1,1,1]];

MCase2:=DirectSum([Irr[i]:i in [1] cat Ans5[2]]);
ToObtain:=IdentifyFactors(Restriction(MCase2,Ld),IrrLd);

"The correct factors are:";
for i in Poss9s do
  MM:=DirectSum([IrrLd[j]:j in i]);
  MM2:=DirectSum([ExteriorPower(MM,3),ExteriorPower(Dual(MM),3),TensorProduct(MM,Dual(MM))]);
  if(IdentifyFactors(MM2,IrrLd) eq ToObtain) then {*labsL[j]:j in i*}; end if;
end for;

"Now to split it between A1 and A6.";

// There are only three options.
M11:=DirectSum([IrrLd[1],IrrLd[1]]); M21:=DirectSum([IrrLd[i]:i in [1,2,2,3]]);
M12:=IrrLd[2]; M22:=DirectSum([IrrLd[i]:i in [1,1,1,2,3]]);
M13:=IrrLd[3]; M23:=DirectSum([IrrLd[i]:i in [1,1,1,2,2]]);
M1:=DirectSum([TensorProduct(M21,Dual(M21)),TensorProduct(M11,Dual(M11))]);
M2:=DirectSum([TensorProduct(M22,Dual(M22)),TensorProduct(M12,Dual(M12))]);
M3:=DirectSum([TensorProduct(M23,Dual(M23)),TensorProduct(M13,Dual(M13))]);

SequenceToMultiset(IdentifyFactors(Restriction(MCase2,L),IrrL)) eq
{* 1^^15, 2^^7, 3^^8, 4^^2, 5^^11, 6^^7, 7^^6, 8^^3, 9^^9, 10^^7, 11^^5, 12^^4,
13^^9, 14^^7, 15^^5, 16^^4, 17^^11, 18^^7, 19^^6, 20^^3 *};

// So we are looking for {* 1^^15, 2^^7, 3^^8, 4^^2 *} for the 1-eigenspace

not SequenceToMultiset(IdentifyFactors(M1,IrrLd)) eq {* 1^^15, 2^^7, 3^^8, 4^^2 *};
SequenceToMultiset(IdentifyFactors(M2,IrrLd)) eq {* 1^^15, 2^^7, 3^^8, 4^^2 *};
not SequenceToMultiset(IdentifyFactors(M3,IrrLd)) eq {* 1^^15, 2^^7, 3^^8, 4^^2 *};

// Thus it is as in the paper.

"Now the radicals from that case.";

JP9121:=RemoveFromTop(P912d,Irr,Irr,[1,2,3,4,5,8,9,10,11,14,15]);
JP9122:=RemoveFromTop(P912d,Irr,Irr,[1,2,3,4,5,8,9,10,11]);
JP9122:=RemoveFromTop(JP9122,Irr,Irr,[6,7]);

JP912:=JP9121 meet JP9122;
JP912:=RemoveFromTop(JP912,Irr,Irr,[1,2,3,4,5,8,9,10,11]);
Mod1:=RemoveFromTop(Dual(P912d/JP912),Irr,Irr,[1..7] cat [14,15]);

JP9211:=RemoveFromTop(P921d,Irr,Irr,[1,2,3,4,5,8,9,10,11,14,15]);
JP9212:=RemoveFromTop(P921d,Irr,Irr,[1,2,3,4,5,8,9,10,11]);
JP9212:=RemoveFromTop(JP9212,Irr,Irr,[6,7]);

JP921:=JP9211 meet JP9212;
JP921:=RemoveFromTop(JP921,Irr,Irr,[1,2,3,4,5,8,9,10,11]);
Mod2:=RemoveFromTop(Dual(P921d/JP921),Irr,Irr,[1..7] cat [14,15]);

Dimension(Mod1) eq 211 and Dimension(Mod2) eq 211;


aa1:=IdentifyFactors(JacobsonRadical(JacobsonRadical(Mod1))/SocleSeries(Mod1)[2],Irr);
aa2:=IdentifyFactors(JacobsonRadical(JacobsonRadical(Mod2))/SocleSeries(Mod2)[2],Irr);
&and[not(i in j):i in [1,6,7],j in [aa1,aa2]];

"And the improved radicals:";

JP9122imp:=RemoveFromTop(JacobsonRadical(P912d),Irr,Irr,[6,7]);

JP912imp:=JP9121 meet JP9122imp;
JP912imp:=RemoveFromTop(JP912imp,Irr,Irr,[8,9,10,11]);
Mod1imp:=RemoveFromTop(Dual(P912d/JP912imp),Irr,Irr,[1..7] cat [14,15]);

JP9212imp:=RemoveFromTop(JacobsonRadical(P921d),Irr,Irr,[6,7]);

JP921imp:=JP9211 meet JP9212imp;
JP921imp:=RemoveFromTop(JP921imp,Irr,Irr,[8,9,10,11]);
Mod2imp:=RemoveFromTop(Dual(P921d/JP921imp),Irr,Irr,[1..7] cat [14,15]);

DescribeLayers(Mod1imp,Irr,labs);
DescribeLayers(Mod2imp,Irr,labs);



"Case 5";

MCase5:=DirectSum([Irr[i]:i in [1] cat Ans5[5]]);
ToObtain:=IdentifyFactors(Restriction(MCase5,Ld),IrrLd);

"The correct factors are:";
for i in Poss9s do
  MM:=DirectSum([IrrLd[j]:j in i]);
  MM2:=DirectSum([ExteriorPower(MM,3),ExteriorPower(Dual(MM),3),TensorProduct(MM,Dual(MM))]);
  if(IdentifyFactors(MM2,IrrLd) eq ToObtain) then {*labsL[j]:j in i*}; end if;
end for;

"Only one way to split it between A1 and A6.";

"The semisimple case:";

M12:=IrrLd[2]; M22:=DirectSum([IrrLd[i]:i in [1,2,2,2]]);
M2:=DirectSum([TensorProduct(M22,Dual(M22)),TensorProduct(M12,Dual(M12))]);

{*DescribeLayers(i,IrrLd,labsL):i in IndecomposableSummands(M2)*} diff {*"1"*};

M22:=DirectSum([IrrLd[2],IrrLd[2],BuildLargestModule(IrrLd[1],[IrrLd[2]])]);
M2:=DirectSum([TensorProduct(M22,Dual(M22)),TensorProduct(M12,Dual(M12))]);

"The 1-eigenspace of the non-semisimple case:";
{*DescribeLayers(i,IrrLd,labsL):i in IndecomposableSummands(M2)*} diff {*"1"*};

"and the other modules, in the order they are given in the paper.";

{*DescribeLayers(i,IrrLd,labsL):i in IndecomposableSummands(ExteriorPower(M22,3))*};
{*DescribeLayers(i,IrrLd,labsL):i in IndecomposableSummands(TensorProduct(M12,ExteriorSquare(M22)))*};
{*DescribeLayers(i,IrrLd,labsL):i in IndecomposableSummands(TensorProduct(M12,M22))*};
{*DescribeLayers(i,IrrLd,labsL):i in IndecomposableSummands(TensorProduct(M12,ExteriorSquare(M22)))*};
{*DescribeLayers(i,IrrLd,labsL):i in IndecomposableSummands(M22)*};

"Finally for this case, the claims about the quotient of the permutation module.";

PLq:=PLquot/Fix(PLquot);

"The socle of the whole quotient contains 2_2 in both z and z^4-eigenspaces:";
{7,19} subset IdentifyFactors(Socle(Restriction(PLq,L)),IrrL);

homs:=AHom(Irr[2],PLq);
PLq1:=PLq/Image(homs.1);
homs:=AHom(Irr[3],PLq);
PLq2:=PLq/Image(homs.1);
PLq1:=Dual(RemoveFromTop(Dual(PLq1),Irr,Irr,[1,4,5,8,9,10,11]));
PLq2:=Dual(RemoveFromTop(Dual(PLq2),Irr,Irr,[1,4,5,8,9,10,11]));
"The smaller modules";
DescribeLayers(PLq1,Irr,labs);
DescribeLayers(PLq2,Irr,labs);
"2_2 lies in the eigenspaces stated.";
7 in IdentifyFactors(Socle(Restriction(PLq1,L)),IrrL);
19 in IdentifyFactors(Socle(Restriction(PLq2,L)),IrrL);



"Case 6 now.";

"Check table of eigenspace actions:";

MCase6:=DirectSum([Irr[i]:i in Ans5[6]]);
SequenceToMultiset(IdentifyFactors(Restriction(MCase6,L),IrrL)) eq 
{* 1^^20, 2^^8, 3^^8, 4^^4, 5^^7, 6^^9, 7^^4, 8^^3, 9^^7, 10^^4, 11^^9, 12^^3,
13^^7, 14^^4, 15^^9, 16^^3, 17^^7, 18^^9, 19^^4, 20^^3 *};

"Now find the action on M(D_6).";

// Note that the action on the 1-eigenspace is field-invariant, so we may work up to field automorphism.
// As 4^4 lies in eigenspace, #4s x #1s le 4.
// as 1^16 lies in eigenspace, #1s le 5

Poss12s:=[[4,4,4],[4,4,2,2],[4,4,2,3],[4,4,2,1,1],
[4,2,2,2,2],[4,2,2,2,3],[4,2,2,3,3],[4,2,2,2,1,1],[4,2,2,3,1,1],
[4,2,2,1,1,1,1],[4,2,3,1,1,1,1],
[2,2,2,2,2,2],[2,2,2,2,2,3],[2,2,2,2,3,3],[2,2,2,3,3,3],
[2,2,2,2,2,1,1],[2,2,2,2,3,1,1],[2,2,2,3,3,1,1],
[2,2,2,2,1,1,1,1],[2,2,2,3,1,1,1,1],[2,2,3,3,1,1,1,1]];

for i in Poss12s do
  MM:=DirectSum([IrrLd[j]:j in i]);
  if(SequenceToMultiset(IdentifyFactors(ExteriorSquare(MM),IrrLd)) eq {* 1^^18, 2^^8, 3^^8, 4^^4 *}) then
    {*labsL[j]:j in i*};
  end if;
end for;

load "Traces/TracesD6.txt";

MM1:=DirectSum([IrrLd[i]:i in [1,1,2,2,2,2,3]]);
BrMM1:=BrauerCharacter(MM1);

// There is only one of these.
for i in Tr5D6 do if(i[1] eq [BrMM1[4]]) then TestD6:=i; end if; end for;

4*TestD6[1,1]+TestD6[4,1]+2*TestD6[2,1]+2*TestD6[3,1]+6 ne BrauerCharacter(Restriction(MCase6,Ld))[4];
//Check the other case, to be sure.

MM2:=DirectSum([IrrLd[i]:i in [2,2,3,3,4]]);
BrMM2:=BrauerCharacter(MM2);

// There is only one of these up to swapping L(\lambda_2) and L(\lambda_3).
for i in Tr5D6 do if(i[1] eq [BrMM2[4]]) then TestD6:=i; end if; end for;

4*TestD6[1,1]+TestD6[4,1]+2*TestD6[2,1]+2*TestD6[3,1]+6 eq BrauerCharacter(Restriction(MCase6,Ld))[4];

"Action of L' on exterior square is";
{*DescribeLayers(i,IrrLd,labsL):i in IndecomposableSummands(ExteriorSquare(MM2))*};

"Action of L on 1,8_1/\\bar 9_{1,2}^*. Check we have the right module:";
E1:=BuildLargestModule(S921d,[S81,Irr[1]]);
DescribeLayers(E1,Irr,labs);
"and its restriction is";
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(Restriction(E1,L))*};

"The other module is";
E2:=BuildLargestModule(S921d,[S32d]);
DescribeLayers(E2,Irr,labs);
"and its restriction is";
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(Restriction(E2,L))*};

"The radical is";
E3:=BuildLargestModule(S921,[Irr[i]:i in [1,6,7,8,9,10,11,12,13,14,15]]);
DescribeLayers(E3,Irr,labs);

"The 1 is a quotient:";
not(1 in IdentifyFactors(RemoveFromTop(E3,Irr,Irr,[1,10]),Irr));



"Case 8, from the Difficult Cases section.";

LE8:=DirectSum([Irr[i]:i in Ans5[8]]);
"The eigenvalues of z are";
EigenvaluesOfElement(LE8,z);
"so the trace is 2. The restriction to L has factors";
{*labsL[i]:i in IdentifyFactors(Restriction(LE8,L),IrrL)*};

SequenceToMultiset(IdentifyFactors(Restriction(LE8,L),IrrL)) eq {* 1^^8, 2^^6, 3^^6, 4^^4, 5^^10, 6^^7, 7^^7, 8^^3, 9^^10, 10^^7, 11^^7, 12^^3, 13^^10, 14^^7, 15^^7, 16^^3, 17^^10, 18^^7, 19^^7, 20^^3 *};

"Test the claim about the factors on M(A_4). The only compatible set of factors is";
Poss5s:=[[4,1],[2,2,1],[2,3,1],[3,3,1],[2,1,1,1],[3,1,1,1]];

for i in Poss5s do for j in Poss5s do
  MM:=DirectSum(TensorPower(DirectSum([IrrLd[k]:k in i]),2),TensorPower(DirectSum([IrrLd[k]:k in j]),2));
  if(SequenceToMultiset(IdentifyFactors(MM,IrrLd))
    eq {* 1^^10, 2^^6, 3^^6, 4^^4*}) then {*labs[k]:k in i*},"and",{*labs[k]:k in j*}; // Add 1^2 as L+0=MxM
  end if;
end for; end for;

// Make the nine modules:

M1:=DirectSum([IrrL[1],IrrL[2],IrrL[3]]);
M2:=DirectSum([IrrL[3],BuildLargestModule(IrrL[2],[IrrL[1]])]);
M3:=Dual(M2);
M4:=GModule(L,FrobeniusImage(M2,1));
M5:=Dual(M4);
M7:=BuildLargestModule(IrrL[1],[IrrL[2],IrrL[3]]);
M6:=Dual(M7);
M8:=SocleSeries(BuildLargestModule(IrrL[3],[IrrL[1],IrrL[2]]))[3];
M9:=Dual(M8);

"Summands of L(A_4) restricted to L' in the nine cases.";

SequenceToMultiset(&cat[[DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(j,Dual(j)))]:j in [M1,M2,M3,M4,M5,M6,M7,M8,M9]]);

// Make the radical:


JB1:=RemoveFromTop(P912d,Irr,Irr,[2..15]);
B1:=Dual(P912d/JB1);
B1a:=RemoveFromTop(B1,Irr,Irr,[2..7] cat[12..15]);
"The first radical:";
DescribeLayers(B1a,Irr,labs);

IsIrreducible(Socle(Dual(RemoveFromTop(B1a,Irr,Irr,[2..7] cat[9..15]))));
B1b:=RemoveFromTop(B1a,Irr,Irr,[8]);
B1b:=RemoveFromTop(B1b,Irr,Irr,[2..7] cat[12..15]);
"The later modules:";
DescribeLayers(B1b,Irr,labs);
DescribeLayers(FrobeniusImage(B1b,3),Irr,labs);
B1bd:=Dual(B1b);
JB1c:=RemoveFromTop(JacobsonRadical(B1bd),Irr,Irr,[2,3,4,5,7]);
DescribeLayers(Dual(B1bd/JB1c),Irr,labs);
JB1d:=RemoveFromTop(JB1c,Irr,Irr,[8,9]);
DescribeLayers(Dual(B1bd/JB1d),Irr,labs);
"To check that the 9_{1,2} does not lie in the 8_2, we take the {9_{1,2}}-residual, which is";
DescribeLayers(RemoveFromTop(Dual(B1bd/JB1d),Irr,Irr,[2..7] cat [9..15]),Irr,labs);
"and for the other claim,";
DescribeLayers(RemoveFromTop(Dual(B1bd/JB1d),Irr,Irr,[2..8] cat [10..15]),Irr,labs);



"Step 1.";

"This should be doubled for V_1:";
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(M1,Dual(M1)))*};

"The other V_i:";
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(M1,ExteriorSquare(M1)))*};



"Step 2.";

"V_1 in case (4,4) is";
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(M4,Dual(M4)))*};

"The other V_i:";
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(M4,ExteriorSquare(M4)))*};
"and (2,2)";
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(Dual(M2),ExteriorSquare(M2)))*};



"Step 3.";

"Exterior square of M_1 is";
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(ExteriorSquare(M8))*};

"The next command shows that the tensor product of it and all possibilities for M_2 is the same";

StripDuplicates([TensorProduct(i,ExteriorSquare(M8)):i in [M1,M2,M3,M4,M5,M6,M7,M8,M9]]);
"and this module is";
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(M1,ExteriorSquare(M8)))*};

"Now prove that there is no 2_1 or 2_2 in the zeta-eigenspace.";

not &or[&or[IsIsomorphic(i,j):i in [IrrL[2],IrrL[3]],j in IndecomposableSummands(TensorProduct(M8,ExteriorSquare(k)))]:k in [M1,M2,M3,M4,M5]];

"Now check that L(A_4)^L is non-zero. This gives the contradiction.";

&and[Dimension(Fix(TensorProduct(i,Dual(i)))) gt 1:i in [M1,M2,M3,M4,M5]]; //Subtract 1 as L+0=MxM^*.



"Step 4.";

"The tensor product in the article:";
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(Dual(M6),ExteriorSquare(M1)))*};



"Step 5.";

"We check the lack of the 2_i stated.";

not &or[&or[IsIsomorphic(i,j):i in [IrrL[2],IrrL[3]],j in IndecomposableSummands(TensorProduct(Dual(k),ExteriorSquare(M2)))]:k in [M6,M7]];

not &or[&or[IsIsomorphic(i,j):i in [IrrL[2]],j in IndecomposableSummands(TensorProduct(M2,ExteriorSquare(k)))]:k in [M4,M5]];

not &or[&or[IsIsomorphic(i,j):i in [IrrL[3]],j in IndecomposableSummands(TensorProduct(Dual(k),ExteriorSquare(M1)))]:k in [M4]];



"Step 6.";

"The module V_{z^2}:";
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(ExteriorSquare(M1),Dual(M2)))*};

"Check that 24_{2,1} has no extensions with other factors.";
for i in [3..15] do Dimension(Ext(S2421,Irr[i])) eq 0; end for;

"L-action of 3_1/24_{2,1}";
E,rho:=Ext(S31,S2421); BB:=Extension(S31,S2421,E.1,rho);
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(Restriction(BB,L))*};



"Step 7.";

"V_z has structure";
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(ExteriorSquare(M7),M4))*};

"For the other three cases the V_{z^2} is";
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(ExteriorSquare(M4),Dual(M6)))*};
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(ExteriorSquare(M6),Dual(M6)))*};
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(ExteriorSquare(M6),Dual(M9)))*};



"Step 8.";

"The next two modules are the same, and are V_z and V_{z^2}:";
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(ExteriorSquare(M8),M8))*};
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(ExteriorSquare(M8),Dual(M8)))*};



"Step 9.";

{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(ExteriorSquare(M8),M6))*};
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(TensorProduct(ExteriorSquare(M6),Dual(M8)))*};


"Case 24_{2,1}^*. Radical is";

C1:=BuildLargestModule(S2421d,[S31,S31d,S32,S32d,S2412,S2412d,S2421,S2421d]);
DescribeLayers(C1,Irr,labs);

"Case 24_{2,1}^*+3_1^*. Radical is";

C2:=BuildLargestModule(S31d,[S31,S31d,S32,S32d,S2412,S2412d]);
DescribeLayers(C2,Irr,labs);

"Guaranteed submodule is";
C3:=RemoveFromTop(C2,Irr,Irr,[2,3,4,5]);
DescribeLayers(C3,Irr,labs);

Wsub:=DirectSum(C3,S2421d);
Wpyx:=[RemoveFromTop(C2,Irr,Irr,[5,3,4]),C1];
"The action of L on W' is";
{*DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(Restriction(Wsub,L))*};

"No 9_{1,2}^* on top of the pyx for W:";
{Dimension(Ext(S912d,i)):i in Wpyx} eq {0};

"Now eliminate the 9_{1,2} only in socle of W_0 possibility.";
Wpyxext:=[BuildLargestModule(i,[S912]):i in Wpyx];
"{9_{1,2}}'-radicals in the paper:";
[DescribeLayers(RemoveFromTop(i,Irr,Irr,[2,3,4,5]),Irr,labs):i in Wpyxext];

"No 3_1^*/3_2/3_1^* when we add the 3s on the extended pyx.";

D1:=DirectSum(Wpyxext);
D1:=BuildLargestModule(D1,[S31d,S32]);
{*DescribeLayers(i,Irr,labs):i in IndecomposableSummands(D1)*};


"The socle being 9_{1,2}^*, try to add 3_2^*/3_1/9_{1,2}^* on top of the pyx.";

Wtest:=DirectSum(Wpyx cat [S912d]);
Wtest2:=BuildLargestModule(Wtest,[S31]);
Wtest3:=BuildLargestModule(Wtest2,[S32d]);
{*DescribeLayers(i,Irr,labs):i in IndecomposableSummands(Wtest3)*};
Wtest4:=RemoveFromTop(Wtest3,Irr,Irr,[2,3,4]);
"The restrictions to L of the modules";
[DescribeLayers(i,Irr,labs):i in IndecomposableSummands(Wtest4)];
"have socles";
[DescribeLayers(Socle(Restriction(i,L)),IrrL,labsL):i in IndecomposableSummands(Wtest4)];


"Now we have the final case, of 9_{1,2}+9_{1,2}^* as the socle of W_0.";

W:=[RemoveFromTop(Wpyx[1],Irr,Irr,[2,4]),Wpyx[2]];

"Add each of 9_{1,2},9_{1,2}^*,\\bar 9_{1,2},\\bar 9_{1,2}^*";
Wext1:=[BuildLargestModule(i,[S912]):i in W];
Wext2:=[BuildLargestModule(i,[S912d]):i in W];
Wext3:=[BuildLargestModule(i,[S921]):i in W];
Wext4:=[BuildLargestModule(i,[S921d]):i in W];

"The socles of the summands upon restriction to L for \\bar 9_{1,2}^\\pm:";
[DescribeLayers(Socle(Restriction(i,L)),IrrL,labsL):i in Wext3];
[DescribeLayers(Socle(Restriction(i,L)),IrrL,labsL):i in Wext4];


"The summand of W_0 we are trying to add:";
W0a:=GModule(G2,FrobeniusImage(B1b,1));
DescribeLayers(W0a,Irr,labs);

"Contradiction for \\bar 9_{1,2}^* on top of W";
Wext4:=[BuildLargestModule(i,[S921d]):i in W];
Wext4a:=[BuildLargestModule(i,[Irr[5]]):i in Wext4];
Wext4b:=[BuildLargestModule(i,[Irr[3]]):i in Wext4a];
"The built up module is";
[DescribeLayers(i,Irr,labs):i in Wext4b];
"and the z^4-eigenspace is";
for i in IndecomposableSummands(Restriction(DirectSum([Wext4b[1]]),L)) do
  if(SequenceToSet(IdentifyFactors(i,IrrL)) subset {17,18,19,20}) then DescribeLayers(i,IrrL,labsL); end if;
end for;

"Contradiction for \\bar 9_{1,2}^* a submodule of W";

WW:=W cat [Irr[11]];
WWa:=[BuildLargestModule(i,[Irr[5]]):i in WW];
WWb:=[BuildLargestModule(i,[Irr[3]]):i in WWa];
"The built up module is the sum of";
[DescribeLayers(i,Irr,labs):i in WWb];
"and their restrictions to L have socles";
[DescribeLayers(Socle(Restriction(i,L)),IrrL,labsL):i in WWb];



"Now the case where 3.1 is the socle of W. The radical first."; 

F1:=BuildLargestModule(S31,[S31,S31d,S32,S32d,S2412,S2412d,S2421,S2421d]);
DescribeLayers(F1,Irr,labs);
"and its dual is";
DescribeLayers(Dual(F1),Irr,labs);


F1a:=SocleSeries(F1)[4];
F1arem1:=RemoveFromTop(F1a,Irr,Irr,[2,3,5,13,14]);
F1arem2:=RemoveFromTop(F1a,Irr,Irr,[2,3,4,5,13]);
"The {3_2}'-residual of the Jacobson radical of the original W' still contains the 3_1/3_2^*/3_1:";
DescribeLayers(F1arem1,Irr,labs);
"and has simple top:";
IsIrreducible(F1arem1/JacobsonRadical(F1arem1));

"The {24_{2,1}}'-residual of the Jacobson radical of the original W' still contains the 3_1/3_2^*/3_1:";
DescribeLayers(F1arem2,Irr,labs);

// Remove the 3_2 and the 24
F1b:=RemoveFromTop(F1a,Irr,Irr,[4,14]);
// Remove the 3_1 that forms the 3_1/3_2^*/3_1.
F1c:=RemoveFromTop(F1b,Irr,Irr,[2]);
"The substitute for W' is therefore";
DescribeLayers(F1c,Irr,labs);

"Now (try to) attach 9s to the top of this. Cannot attach 9_{1,2}:";
Dimension(Ext(S912,F1c)) eq 0;
"Can attach a single 9_{1,2}^*";
E,rho:=Ext(S912d,F1c); Dimension(E) eq 1; F1d:=Extension(S912d,F1c,E.1,rho);

"Remove any 3-dimensional factors so we can understand what is going on.";
F1e:=RemoveFromTop(F1d,Irr,Irr,[2,3,4,5]);
DescribeLayers(F1e,Irr,labs);
"The socle of the restriction to L is";
DescribeLayers(Socle(Restriction(F1e,L)),IrrL,labsL);

"Now try to add the modules suggested to the pyx. This doesn't work:";

for j in [S31d,S32,S912] do Dimension(Ext(j,F1c)) eq 0; end for;

F1g1:=RemoveFromTop(F1e,Irr,Irr,[9]);
F1g2:=BuildLargestModule(S912,[S31d,S32]);
"The module";
DescribeLayers(F1g1,Irr,labs);
"restricts to L as";
DescribeLayers(Restriction(F1g1,L),IrrL,labsL);
"The module";
DescribeLayers(F1g2,Irr,labs);
"restricts to L as";
DescribeLayers(Restriction(F1g2,L),IrrL,labsL);



"Now the case where 3.1+3_1^* is the socle of W."; 

J1a:=F1c;
J1b:=GModule(G2,FrobeniusImage(J1a,2));
J1:=[J1a,J1b];
"The module W' is";
[DescribeLayers(i,Irr,labs):i in J1];

// J2,J2a, etc., will be the case with a 24_{2,1}, and J3,J3a, etc., for the case with 24_{2,1}^*.

"The two options for a submodule of W are";
J2:=[RemoveFromTop(i,Irr,Irr,[2,3,4,5,12,15]):i in J1]; J2[2]:=S31d;
[DescribeLayers(i,Irr,labs):i in J2];
"and";
J3:=[RemoveFromTop(i,Irr,Irr,[2,3,4,5,12,14]):i in J1];
[DescribeLayers(i,Irr,labs):i in J3];

"The z^4-eigenspaces of these modules are:";
[[DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(Restriction(j,L))|SequenceToSet(IdentifyFactors(i,IrrL)) subset {17..20}]:j in J2];
"and";
[[DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(Restriction(j,L))|SequenceToSet(IdentifyFactors(i,IrrL)) subset {17..20}]:j in J3];

"Adding the 3_2 yields";
J2a:=[J2[1],BuildLargestModule(J2[2],[S32])];
[[DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(Restriction(j,L))|SequenceToSet(IdentifyFactors(i,IrrL)) subset {17..20}]:j in J2a];

J3a:=[J3[1],BuildLargestModule(J3[2],[S32])];
[[DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(Restriction(j,L))|SequenceToSet(IdentifyFactors(i,IrrL)) subset {17..20}]:j in J3a];
"Add on the \\bar 9_{1,2} to both summands of these modules:";
J2b:=[BuildLargestModule(i,[S921]):i in J2a];
J3b:=[BuildLargestModule(i,[S921]):i in J3a];
[DescribeLayers(i,Irr,labs):i in J2b];
[DescribeLayers(i,Irr,labs):i in J3b];

"And their restrictions to L have z^4-eigenspace";
[[DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(Restriction(j,L))|SequenceToSet(IdentifyFactors(i,IrrL)) subset {17..20}]:j in J2b];
[[DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(Restriction(j,L))|SequenceToSet(IdentifyFactors(i,IrrL)) subset {17..20}]:j in J3b];
"The duals in the last case are";
[[DescribeLayers(Dual(i),IrrL,labsL):i in IndecomposableSummands(Restriction(j,L))|SequenceToSet(IdentifyFactors(i,IrrL)) subset {17..20}]:j in J3b];


"Now suppose that W_0 has \\bar 9_{1,2}^* in the socle. The module we use is";
J4:=[BuildLargestModule(RemoveFromTop(i,Irr,Irr,[12]),[S31,S31d,S32,S32d]):i in J1];
[DescribeLayers(i,Irr,labs):i in J4];
"Add on the \\bar 9_{1,2}^*";
J4a:=[BuildLargestModule(i,[S921d]):i in J4];
[DescribeLayers(i,Irr,labs):i in J4a];
"Add on the 3_2^*s";
J4b:=[BuildLargestModule(i,[S32d]):i in J4a];
[DescribeLayers(i,Irr,labs):i in J4b];
"However, there are now no 3_1^* to add.";
[Dimension(Ext(S31d,i)) eq 0:i in J4b];

"Take the {3_i^\\pm}-residual to obtain";
J4ares:=[RemoveFromTop(i,Irr,Irr,[2,3,4,5]):i in J4a];
[DescribeLayers(i,Irr,labs):i in J4ares];

"Removing the 24_{2,1} and 24_{2,1}^* and examining the z^3-eigenspace:";
J4ares1:=[RemoveFromTop(i,Irr,Irr,[14]):i in J4ares];
&cat[[DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(Restriction(j,L))|SequenceToSet(IdentifyFactors(i,IrrL)) subset {13..16}]:j in J4ares1];

J4ares2:=[RemoveFromTop(i,Irr,Irr,[15]):i in J4ares];
&cat[[DescribeLayers(i,IrrL,labsL):i in IndecomposableSummands(Restriction(j,L))|SequenceToSet(IdentifyFactors(i,IrrL)) subset {13..16}]:j in J4ares2];

"And that is it.";