/*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(15)!(2,3,4,5)(6,9,11,8)(7,10)(12,15,14,13);


G2:=PSp(4,4);
CC:=ConjugacyClasses(G2);
load "Subgroups/E8/ComputeFactorsE8.txt";
load "Subgroups/functions.txt";

load "Traces/E8/Tr17E8Eigs";
Tr17E8:={[&+i]:i in Tr17E8Eigs};
delete Tr17E8Eigs;

Tr cat:=<{},Tr17E8>;

X:=OrderByDimension(CompositionFactors(PermutationModule(G2,GF(4))));
S1:=X[1];
S41:=X[6];
S43:=Socle(SymmetricSquare(S41));
S42:=SocleFactors(ExteriorSquare(S41))[2];
S44:=Socle(SymmetricSquare(S42));

S1612:=TensorProduct(S41,S42);
S1613:=TensorProduct(S41,S43);
S1614:=TensorProduct(S41,S44);
S1623:=TensorProduct(S42,S43);
S1624:=TensorProduct(S42,S44);
S1634:=TensorProduct(S43,S44);
Irr:=[S1,S41,S42,S43,S44,S1612,S1613,S1614,S1623,S1624,S1634];

for i in [2..3] do for j in [i+1..4] do for k in [j+1..5] do Append(~Irr,TensorProduct(TensorProduct(Irr[i],Irr[j]),Irr[k])); end for; end for; end for;

AnsA:=RestrictedPartitions(248,{1,4,16,64});
AnsB:=[i:i in AnsA|Multiplicity(i,1)+Multiplicity(i,16)-Multiplicity(i,4)-Multiplicity(i,64) in {-2,23}];

Ans3,Mods,Irr:=ComputeConspicuousFactorsE8LowMemory(Irr,CC,Tr,[17,15,5],AnsB);
Ans3a:=OrbitReps(Ans3,g);
[SequenceToMultiset(i):i in Ans3a] eq [
    {* 1^^68, 2^^32, 3^^12, 4 *},
    {* 1^^32, 2^^16, 3^^16, 4^^6, 6^^4 *},
    {* 1^^32, 2^^9, 3, 5^^8, 8, 11^^8 *},
    {* 1^^32, 2^^8, 3, 4^^8, 5, 7, 10^^8 *},
    {* 1^^24, 2^^12, 3^^9, 4^^3, 5^^8, 8^^4, 10^^2 *},
    {* 1^^16, 2^^9, 3^^4, 4^^9, 5^^4, 7^^4, 8^^2, 9^^2 *},
    {* 1^^12, 2^^7, 3^^5, 4^^5, 5^^2, 7^^2, 8^^3, 10, 15 *},
    {* 1^^12, 2^^6, 3^^4, 4^^5, 5^^4, 7^^2, 8^^2, 9, 10, 14 *},
    {* 1^^12, 2^^5, 3^^5, 4, 5^^4, 6, 8, 10, 14^^2 *},
    {* 1^^8, 2^^4, 3^^2, 4^^4, 5^^2, 7^^2, 8, 9, 13, 15 *},
    {* 1^^8, 2^^4, 3^^2, 4^^4, 5^^2, 8^^2, 9^^2, 12, 14 *}
];

// Now choose those of positive pressure

Ans4:=[Ans3a[i]:i in [2,5,6,7,8,9,10,11]];
Mods:=[DirectSum([Irr[i]:i in j]):j in Ans4];

Z:=Subgroups(G2:OrderEqual:=3600);
if(IsIrreducible(Restriction(S42,Z[1]`subgroup))) then LL1:=Z[1]`subgroup; LL2:=Z[2]`subgroup; else LL1:=Z[2]`subgroup; LL2:=Z[1]`subgroup; end if;
CCLL1:=ConjugacyClasses(LL1);
CCLL2:=ConjugacyClasses(LL2);
y1:=CCLL1[5,3];
y2:=CCLL2[5,3];
C1:=Centralizer(LL1,y1);
C2:=Centralizer(LL2,y2);
L1:=DerivedSubgroup(C1);
L2:=DerivedSubgroup(C2);
for i in [8,9,10,11] do if(DerivedSubgroup(Centralizer(LL1,CCLL1[i,3])) eq L1) then z1:=CCLL1[i,3]; break i; end if; end for;
for i in [8,9,10,11] do if(DerivedSubgroup(Centralizer(LL2,CCLL2[i,3])) eq L2) then z2:=CCLL2[i,3]; break i; end if; end for;
D1:=Centralizer(LL1,z1);
D2:=Centralizer(LL2,z2);

IrrD1:=[TrivialModule(D1,GF(4))];
XRes:=IndecomposableSummands(Restriction(S41,D1));
if(IsAbsolutelyIrreducible(XRes[1])) then
  Append(~IrrD1,XRes[1]);
  Append(~IrrD1,GModule(D1,FrobeniusImage(IrrD1[2],1)));
  Append(~IrrD1,TensorProduct(IrrD1[2],IrrD1[3]));
  Append(~IrrD1,XRes[2]);
else
  Append(~IrrD1,XRes[2]);
  Append(~IrrD1,GModule(D1,FrobeniusImage(IrrD1[2],1)));
  Append(~IrrD1,TensorProduct(IrrD1[2],IrrD1[3]));
  Append(~IrrD1,XRes[1]);
end if;
for i in [2..4] do Append(~IrrD1,TensorProduct(IrrD1[i],IrrD1[5])); end for;
Append(~IrrD1,OrderByDimension(CompositionFactors(TensorPower(IrrD1[5],2)))[3]);
for i in [2..4] do Append(~IrrD1,TensorProduct(IrrD1[i],IrrD1[9])); end for;

IrrD2:=[TrivialModule(D2,GF(4))];
XRes:=IndecomposableSummands(Restriction(S42,D2));
if(IsAbsolutelyIrreducible(XRes[1])) then
  Append(~IrrD2,XRes[1]);
  Insert(~IrrD2,2,GModule(D2,FrobeniusImage(IrrD2[2],1)));
  Append(~IrrD2,TensorProduct(IrrD2[2],IrrD2[3]));
  Append(~IrrD2,XRes[2]);
else
  Append(~IrrD2,XRes[2]);
  Insert(~IrrD2,2,GModule(D2,FrobeniusImage(IrrD2[2],1)));
  Append(~IrrD2,TensorProduct(IrrD2[2],IrrD2[3]));
  Append(~IrrD2,XRes[1]);
end if;
for i in [2..4] do Append(~IrrD2,TensorProduct(IrrD2[i],IrrD2[5])); end for;
Append(~IrrD2,OrderByDimension(CompositionFactors(TensorPower(IrrD2[5],2)))[3]);
for i in [2..4] do Append(~IrrD2,TensorProduct(IrrD2[i],IrrD2[9])); end for;
for i in [1..4] do IrrD2:=Swap(IrrD2,4+i,8+i); end for;

IrrL1:=[Restriction(IrrD1[i],L1):i in [1..4]];
IrrL2:=[Restriction(IrrD2[i],L2):i in [1..4]];

labs:=["1","4_1","4_2","4_3","4_4","16_{1,2}","16_{1,3}","16_{1,4}","16_{2,3}","16_{2,4}","16_{3,4}","64","64","64","64"];
labsL2:=["1","2_1","2_2","4"];

// Case 2

// Action of L1 on M(D5).

Tens:=[DirectSum([IrrL1[i]:i in j]):j in [[4,4,2],[4,4,1,1],[4,2,2,2],[4,2,2,3],
[4,2,2,1,1],[4,2,3,1,1],[4,2,1,1,1,1],[2,2,2,2,2],[2,2,2,2,3],[2,2,2,3,3],
[2,2,2,2,1,1],[2,2,2,3,1,1],[2,2,3,3,1,1],[2,2,2,1,1,1,1],[2,2,3,1,1,1,1],
[2,2,1,1,1,1,1,1],[2,3,1,1,1,1,1,1],[2,1,1,1,1,1,1,1,1]]];
aa:=[SequenceToMultiset(IdentifyFactors(ExteriorSquare(i),IrrL1)):i in Tens];
[i:i in [1..#aa]|Multiplicity(aa[i],1) lt 24 and Minimum(Multiplicity(aa[i],2),Multiplicity(aa[i],3)) le 3 and Multiplicity(aa[i],4) eq 0];


// This yields 14 and 16, but unless L(A2) is trivial, we have a copy of both 2.1 and 2.2 in it. So it must be trivial.



// Case 2 radicals
A1:=S1624; for i in [S42,S44,S1,S41,S43,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,S1614,S1624] do Dimension(Ext(A1,i)); end for;

A2:=S44; for i in [S1624,S1, S41,S42,S43,S44, S1, S41,S42,S43,S44, S1, S41,S43, S1614,S1, S41,S44, S1, S42, S1, S41, S1, S44] 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,S43,S44,S1614,S1624] do Dimension(Ext(A2,i)); end for;
A2a:=RemoveFromTop(Dual(A2),Irr,Irr,[1,2,3,4]);
DescribeLayers(A2a,Irr,labs);

A3:=S1614; for i in [S41,S1,S42,S1,S41,S43,S1,S42,S44,S1,S41,S1614,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,S42,S43,S44,S1614,S1624] do Dimension(Ext(A3,i)); end for;
A3a:=RemoveFromTop(Dual(A3),Irr,Irr,[1,2,3,4]);
DescribeLayers(A3a,Irr,labs);

// Case 3
"First find action of L2. As there is no 4, cannot be 4 or 2_1,2_2 in factors. As 1-eigenspace is balanced, must be different actions.";
M1:=DirectSum([IrrL2[i]:i in [1,2,2]]); M2:=DirectSum([IrrL2[i]:i in [1,3,3]]);
{*labsL2[j]:j in IdentifyFactors(TensorProduct(M1,ExteriorSquare(M2)),IrrL2)*};
"This is the right way round";
M2:=DirectSum([IrrL2[i]:i in [1,2,2]]); M1:=DirectSum([IrrL2[i]:i in [1,3,3]]);
{*labsL2[j]:j in IdentifyFactors(TensorProduct(M1,ExteriorSquare(M2)),IrrL2)*};
"This is the wrong way round";

"Restrictions of modules mentioned are semisimple:";
//for i in [S41,S42,S43,S44,S1613] do DescribeLayers(Restriction(i,L2),IrrL2,labsL2); end for;
&and[IsSemisimple(Restriction(i,L2)):i in [S41,S42,S43,S44,S1613]];

"Set up a copy of SL(2,16), which functions as our version of an A1 subgroup.";
TestG:=SL(2,16);
TestH:=Subgroups(TestG:OrderEqual:=60)[1]`subgroup;
S21:=GModule(TestG);
S22:=GModule(TestG,FrobeniusImage(S21,1));
S23:=GModule(TestG,FrobeniusImage(S21,2));
S24:=GModule(TestG,FrobeniusImage(S21,3));
S412:=TensorProduct(S21,S22);
S413:=TensorProduct(S21,S23);
S414:=TensorProduct(S21,S24);
S423:=TensorProduct(S22,S23);
S424:=TensorProduct(S22,S24);
S434:=TensorProduct(S23,S24);
IrrTest:=[TrivialModule(TestG,GF(16)),S21,S22,S23,S24,S412,S413,S414,S423,S424,S434];
S221:=SymmetricSquare(S21);
S231:=GModule(TestG,FrobeniusImage(S221,1));
S211:=GModule(TestG,FrobeniusImage(S221,3));
IrrHTest:=[Restriction(IrrTest[i],TestH):i in [1,2,3,6]];

"For 1+1+0 and 0/2+2, we check that the factors in the socle of each eigenspace for A1 and that for SL(2,4) are the same";

M1:=DirectSum([S21,S21,TrivialModule(TestG,GF(16))]); M2:=DirectSum([S221,S22]);
Z0:=DirectSum([TensorProduct(M1,Dual(M1)),TensorProduct(M2,Dual(M2))]);
Z1:=TensorProduct(M1,ExteriorSquare(M2));
Z2:=TensorProduct(ExteriorSquare(M1),Dual(M2));
Z3:=Dual(Z2); Z4:=Dual(Z1);

for ZZ in [Z0,Z1,Z2,Z3,Z4] do
Sort(Multiplicities(SequenceToMultiset(IdentifyFactors(Socle(ZZ),IrrTest)))) eq Sort(Multiplicities(SequenceToMultiset(IdentifyFactors(Socle(Restriction(ZZ,TestH)),IrrHTest))));
end for;

"And now the complete (dual of) z1-, (dual of) z2-, z3- and z4-eigenspaces for SL(2,4)";

for ZZ in [Dual(Z1),Dual(Z2),Z3,Z4] do
DescribeLayers(Restriction(ZZ,TestH),IrrHTest,labsL2);
end for;

"For 0/4+4 and 2+2+0, we check that the factors in the socle of each eigenspace for A1 and that for SL(2,4) are the same";

M1:=DirectSum([S231,S23]); M2:=DirectSum([S22,S22,TrivialModule(TestG,GF(16))]);
Z0:=DirectSum([TensorProduct(M1,Dual(M1)),TensorProduct(M2,Dual(M2))]);
Z1:=TensorProduct(M1,ExteriorSquare(M2));
Z2:=TensorProduct(ExteriorSquare(M1),Dual(M2));
Z3:=Dual(Z2); Z4:=Dual(Z1);

for ZZ in [Z0,Z1,Z2,Z3,Z4] do
Sort(Multiplicities(SequenceToMultiset(IdentifyFactors(Socle(ZZ),IrrTest)))) eq Sort(Multiplicities(SequenceToMultiset(IdentifyFactors(Socle(Restriction(ZZ,TestH)),IrrHTest))));
end for;

"And now the complete z1-, (dual of) z2-, z3- and (dual of) z4-eigenspaces for SL(2,4)";

for ZZ in [Z1,Dual(Z2),Z3,Dual(Z4)] do
DescribeLayers(Restriction(ZZ,TestH),IrrHTest,labsL2);
end for;

"For 1+1+0 and 2+2+0, we check that the factors in the socle of each eigenspace for A1 and that for SL(2,4) are the same";

M1:=DirectSum([S21,S21,TrivialModule(TestG,GF(16))]); M2:=DirectSum([S22,S22,TrivialModule(TestG,GF(16))]);
Z0:=DirectSum([TensorProduct(M1,Dual(M1)),TensorProduct(M2,Dual(M2))]);
Z1:=TensorProduct(M1,ExteriorSquare(M2));
Z2:=TensorProduct(ExteriorSquare(M1),Dual(M2));
Z3:=Dual(Z2); Z4:=Dual(Z1);

for ZZ in [Z0,Z1,Z2,Z3,Z4] do
Sort(Multiplicities(SequenceToMultiset(IdentifyFactors(Socle(ZZ),IrrTest)))) eq Sort(Multiplicities(SequenceToMultiset(IdentifyFactors(Socle(Restriction(ZZ,TestH)),IrrHTest))));
end for;

"And now the complete z1- and z2-eigenspaces for SL(2,4)";

for ZZ in [Z1,Z2] do
DescribeLayers(Restriction(ZZ,TestH),IrrHTest,labsL2);
end for;

"Now switch to y2. Check factors on M(A8).";

Nines:=[DirectSum([IrrL2[i]:i in j]):j in [[4,4,1],[4,2,2,1],[4,2,3,1],[4,2,1,1,1],
[4,1,1,1,1,1],[2,2,2,1],[2,2,2,3,1],[2,2,3,3,1],[2,2,2,1,1,1],[2,2,3,1,1,1],
[2,2,1,1,1,1,1],[2,3,1,1,1,1,1],[2,1,1,1,1,1,1]]];
for i in Nines do Exclude(SequenceToMultiset(IdentifyFactors(TensorProduct(i,Dual(i)),IrrL2)),1); end for;
"Only cases 1 and 8 work for the 1-eigenspace. These are factors for all of L(E8):";
for i in [Nines[j]:j in [1,8]] do Exclude(SequenceToMultiset(IdentifyFactors(DirectSum([TensorProduct(i,Dual(i)),ExteriorPower(i,3),ExteriorPower(Dual(i),3)]),IrrL2)),1); end for;

"Construct the five modules' actions on L(E8) and compare with those that are left";

N1:=DirectSum([IrrL2[i]:i in [1,2,2,3,3]]);
N2:=DirectSum([IrrL2[i]:i in [2,3,3]] cat [SymmetricSquare(IrrL2[3])]);
N3:=GModule(L2,FrobeniusImage(N2,1));
N4:=DirectSum([IrrL2[i]:i in [2,3]] cat [BuildLargestModule(IrrL2[1],[IrrL2[2],IrrL2[3]])]);
N5:=DirectSum([IrrL2[i]:i in [2,3]] cat [JacobsonRadical(BuildLargestModule(IrrL2[3],[IrrL2[1],IrrL2[2]]))]);

// Must add an extra trivial as we pick up two from the other action and only one from this.

NMods:=[DirectSum([ExteriorPower(i,3),Dual(ExteriorPower(i,3)),TensorProduct(i,Dual(i)),IrrL2[1]]):i in [N1,N2,N3,N4,N5]];

"Do any of them coincide?";

M12:=DirectSum(SymmetricSquare(IrrL2[3]),IrrL2[2]);
M22:=DirectSum(SymmetricSquare(IrrL2[2]),IrrL2[3]);
for i in [M12,Dual(M12)] do for j in [M22,Dual(M22)] do
  NN:=DirectSum([
  TensorProduct(i,Dual(i)),TensorProduct(j,Dual(j)),
  TensorProduct(i,ExteriorSquare(j)),
  TensorProduct(Dual(j),ExteriorSquare(i)),
  TensorProduct(j,ExteriorSquare(Dual(i))),
  TensorProduct(Dual(i),Dual(ExteriorSquare(j)))]);
  &or[IdentifyLayers(ii,IrrL2) eq IdentifyLayers(NN,IrrL2):ii in NMods];
end for; end for;


// Case 4 possibilities for A8
"Case 4: proof that there is no copy of L1 inside A8.";
t1:=SequenceToMultiset(IdentifyFactors(Restriction(Mods[4],L1),IrrL1));
for i in Nines do aa:=Exclude(SequenceToMultiset(IdentifyFactors(DirectSum([TensorProduct(i,Dual(i)),ExteriorPower(i,3),ExteriorPower(Dual(i),3)]),IrrL2)),1); aa,aa eq t1; end for;

// Case 5

"First check that Table 8.3 is correct";
S64134:=TensorProduct(S1634,S41);

for i in [S41,S42,S43,S44,S1613,S1614,S1623,S1624,S64134] do IdentifyLayers(Restriction(i,D2),IrrD2); end for;

"Action of L2 is semisimple";
M1:=DirectSum([IrrL2[i]:i in [1,3,4]]); M2:=IrrL2[3];

for i in IndecomposableSummands(ExteriorPower(M1,3)) do DescribeLayers(i,IrrL2,labsL2); end for;
for i in IndecomposableSummands(TensorProduct(M1,M2)) do DescribeLayers(i,IrrL2,labsL2); end for;

"Action of L2 is non-semisimple";
M1:=DirectSum(SymmetricSquare(IrrL2[2]),IrrL2[4]); M2:=IrrL2[3];

for i in IndecomposableSummands(ExteriorPower(M1,3)) do DescribeLayers(i,IrrL2,labsL2); end for;
for i in IndecomposableSummands(TensorProduct(M1,M2)) do DescribeLayers(i,IrrL2,labsL2); end for;

M1:=DirectSum(SymmetricSquare(S21),S412); M2:=S22;

for i in IndecomposableSummands(ExteriorPower(M1,3)) do IdentifyLayers(i,IrrTest); end for;
for i in IndecomposableSummands(TensorProduct(M1,M2)) do IdentifyLayers(i,IrrTest); end for;

LL2bar:=Normalizer(G2,LL2);
IrrLL2old:=IrreducibleModules(LL2,GF(4));
IrrLL2:=[IrrLL2old[1]];
for i in [2..5] do if(IdentifyModule(Restriction(IrrLL2old[i],L2),IrrL2) eq 2) then Append(~IrrLL2,IrrLL2old[i]); end if; end for;
Append(~IrrLL2,GModule(LL2,FrobeniusImage(IrrLL2[2],1)));
IrrLL2:=AppendWithoutDuplicates(IrrLL2,CompositionFactors(Restriction(Induction(IrrLL2[2],LL2bar),LL2)));
Append(~IrrLL2,GModule(LL2,FrobeniusImage(IrrLL2[4],1)));
Append(~IrrLL2,TensorProduct(IrrLL2[2],IrrLL2[3]));
Append(~IrrLL2,TensorProduct(IrrLL2[4],IrrLL2[5]));
Append(~IrrLL2,TensorProduct(IrrLL2[2],IrrLL2[4]));
Append(~IrrLL2,TensorProduct(IrrLL2[2],IrrLL2[5]));
Append(~IrrLL2,TensorProduct(IrrLL2[3],IrrLL2[4]));
Append(~IrrLL2,TensorProduct(IrrLL2[3],IrrLL2[5]));
Append(~IrrLL2,TensorProduct(IrrLL2[2],IrrLL2[7]));
Append(~IrrLL2,TensorProduct(IrrLL2[3],IrrLL2[7]));
Append(~IrrLL2,TensorProduct(IrrLL2[4],IrrLL2[6]));
Append(~IrrLL2,TensorProduct(IrrLL2[5],IrrLL2[6]));
Append(~IrrLL2,TensorProduct(IrrLL2[6],IrrLL2[7]));
// Order is 00,10,20,01,02,30,03,11,12,21,22,13,23,31,32,33

TestEigzeta:=Restriction(TensorProduct(IrrLL2[13],IrrLL2[15]),D2);
for i in IndecomposableSummands(TestEigzeta) do IdentifyLayers(i,IrrD2); end for;

TestEigzeta:=Restriction(TensorProduct(IrrLL2[2],IrrLL2[4]),D2);
for i in IndecomposableSummands(TestEigzeta) do IdentifyLayers(i,IrrD2); end for;

// Case 6 extensions

"Check the extensions claims for Case 6";

for i in [S1612,S1624] do Dimension(Ext(S64134,i)) eq 0; end for;
E,rho:=Ext(S64134,S1614); A4:=MaximalExtension(S64134,S1614,E,rho); #CompositionFactors(A4) eq 2;
Dimension(Ext(A4,S64134));

// Case 7 extensions

// S64124:=TensorProduct(S1624,S41);
// S64234:=TensorProduct(S1634,S42);
// Dimension(Ext(S64124,S1613)) eq 0; //This follows since untwisting it by the Frobenius yields Dimension(Ext(S64134,S1624)) eq 0;
// Dimension(Ext(S64234,S1613)) eq 0; //This follows since twisting it by the Frobenius yields Dimension(Ext(S64134,S1624)) eq 0;
// This was checked in Case 6.

// Case 8
"Now consider Case 8";

Nines:=[DirectSum([IrrL1[i]:i in j]):j in [[4,4,1],[4,2,2,1],[4,2,3,1],[4,2,1,1,1],
[4,1,1,1,1,1],[2,2,2,1],[2,2,2,3,1],[2,2,3,3,1],[2,2,2,1,1,1],[2,2,3,1,1,1],
[2,2,1,1,1,1,1],[2,3,1,1,1,1,1],[2,1,1,1,1,1,1]]];
"These are the possible composition factors for L1";
for i in Nines do Exclude(SequenceToMultiset(IdentifyFactors(TensorProduct(i,Dual(i)),IrrL1)),1); end for;

"But these the composition factors of L1 on L(E8)";
SequenceToMultiset(IdentifyFactors(Restriction(Mods[8],L1),IrrL1));
