"Verbosity 0 just checks that there are elements that power to a given element of order 85 with the correct eigenvalues on L(G) and stabilizing a given set of eigenspaces";
printf "\n";
"Verbosity 1 in addition checks that the element of order 13 has the correct eigenvalues on L(G), that the elements of order 85 that power to it have the correct eigenvalues on L(G), and that those elements form two orbits under the Weyl group, one given by raising the other to the power 69.";
printf "\n";
"Verbosity 2 in addition checks that the 16 elements of order 85 are the only ones with the correct eigenvalues on L(G)";
printf "\n";
"Verbosity 3 in addition constructs PSL(4,4) and checks that the eigenvalues of an element of order 85 on the various simple modules, and on L(G), are what are claimed.";
printf "\n";
"Verbosity 4 in addition computes the radicals from the paper.";
printf "\n";
"The default is verbosity 4.";
Verbosity:=4;

load "Subgroups/functions.txt";

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

q:=3571;
F:=GF(q);
w:=PrimitiveElement(F);
FF<ww>:=GF(2^8);
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>;

T17:=sub<T|[T.i^210:i in [1..8]]>;
T1785:=sub<T|[T.i^2:i in [1..8]]>;

R:=ResidueClassRing(1785);
X:=[[R!((Log((T1785.j)[i,i])*1785) div (q-1)):i in [1..248]]:j in [1..8]];


EigsToTest:=  [{ <ww^66, 1>, <ww^36, 1>, <ww^144, 1>, <ww^9, 1> },
   { <ww^111, 1>, <ww^219, 1>, <ww^246, 1>, <ww^189, 1> },
   { <ww^18, 1>, <ww^132, 1>, <ww^33, 1>, <ww^72, 1> },
   { <ww^237, 1>, <ww^183, 1>, <ww^123, 1>, <ww^222, 1> },
   { <ww^180, 1>, <ww^45, 1>, <ww^210, 1>, <ww^153, 1>, <ww^75, 1>, <ww^102, 1> },
   { <ww^51, 1>, <ww^105, 1>, <ww^204, 1>, <ww^90, 1>, <ww^150, 1>, <ww^165, 1> }];

E1:={<ww^108,4>,<ww^63,2>,<ww^120,3>,<ww^15,3>,<ww^156,4>,<ww^93,3>,<ww^189,3>,<1,8>,<ww^42,3>,<ww^57,4>,<ww^225,3>,<ww^90,2>,<ww^138,3>,<ww^162,3>,
<ww^105,2>,<ww^96,2>,<ww^219,3>,<ww^78,4>,<ww^21,3>,<ww^153,2>,<ww^3,2>,<ww^27,4>,<ww^234,3>,<ww^117,3>,<ww^210,2>,<ww^99,4>,<ww^192,2>,<ww^84,3>,<ww^228,4>,
<ww^135,3>,<ww^18,3>,<ww^102,2>,<ww^243,2>,<ww^45,2>,<ww^231,2>,<ww^216,4>,<ww^171,3>,<ww^195,3>,<ww^180,2>,<ww^201,4>,<ww^36,3>,<ww^75,2>,<ww^6,2>,
<ww^132,3>,<ww^207,2>,<ww^111,3>,<ww^66,3>,<ww^249,2>,<ww^114,4>,<ww^150,2>,<ww^69,3>,<ww^147,4>,<ww^144,3>,<ww^24,2>,<ww^159,2>,<ww^51,2>,<ww^186,3>,
<ww^141,4>,<ww^129,2>,<ww^165,2>,<ww^33,3>,<ww^252,2>,<ww^168,3>,<ww^174,3>,<ww^30,3>,<ww^60,3>,<ww^87,3>,<ww^126,2>,<ww^48,2>,<ww^54,4>,<ww^177,4>,
<ww^81,3>,<ww^237,3>,<ww^246,3>,<ww^183,3>,<ww^12,2>,<ww^204,2>,<ww^213,3>,<ww^9,3>,<ww^123,3>,<ww^240,3>,<ww^198,4>,<ww^222,3>,<ww^39,4>,<ww^72,3>};
E2:={<ww^150,15>,<ww^225,14>,<ww^135,14>,<ww^240,14>,<1,16>,<ww^195,14>,<ww^90,15>,<ww^165,15>,<ww^45,15>,<ww^60,14>,<ww^30,14>,<ww^75,15>,<ww^120,14>,<ww^15,14>,<ww^180,15>,
<ww^210,15>,<ww^105,15>};
Eigs1:={<w^(14*Log(i[1])),i[2]>:i in E1};
Eigs2:={<w^(14*Log(i[1])),i[2]>:i in E2};

if(Verbosity ge 3) then
  "Checking that eigenvalues of simple modules on PSL(4,4) are correct";
  G2:=PSL(4,4);
  for i in CompositionFactors(PermutationModule(G2,GF(4))) do if(Dimension(i) eq 4) then S41:=i; break i; end if; end for;
  S41d:=Dual(S41);
  S42:=Socle(SymmetricSquare(S41));
  S42d:=Dual(S42);
  S61:=ExteriorSquare(S41);
  S62:=ExteriorSquare(S42);
  S2412:=TensorProduct(S41,S62);
  S2412d:=TensorProduct(S41d,S62);
  S2421:=TensorProduct(S42,S61);
  S2421d:=TensorProduct(S42d,S61);
  S1612:=TensorProduct(S41,S42);
  S1612d:=TensorProduct(S41d,S42d);
  S1621:=TensorProduct(S41,S42d);
  S1621d:=TensorProduct(S41d,S42);
  for i in CompositionFactors(ExteriorSquare(S61)) do if(Dimension(i) eq 14) then S141:=i; else S1:=i; end if; end for;
  for i in CompositionFactors(ExteriorSquare(S62)) do if(Dimension(i) eq 14) then S142:=i; end if; end for;
  MM:=DirectSum([S1,S1,S1,S1,S41,S41,S41d,S41d,S42,S42,S42d,S42d,S61,S61,S62,S62,S141,S142,S1612,S1612d,S1621,S1621d,S2412,S2412d,S2421,S2421d]);
  CC:=ConjugacyClasses(G2);
  for i in [69..84] do if(EigenvaluesOfElement(ChangeRing(S41,FF),CC[i,3]) eq {<ww^66,1>,<ww^36,1>,<ww^144,1>,<ww^9,1>}) then x:=CC[i,3]; end if; end for;
  Irr:=[S1,S41,S41d,S42,S42d,S61,S62,S141,S142,S1612,S1612d,S1621,S1621d,S2412,S2412d,S2421,S2421d];
  labs:=["1","4_1","4_1^*","4_2","4_2^*","6_1","6_2","14_1","14_2","16_{1,2}","16_{1,2}^*","\\bar{16}_{1,2}","\\bar{16}_{1,2}^*","24_{1,2}","24_{1,2}^*","24_{2,1}","24_{2,1}^*"];

  [EigenvaluesOfElement(ChangeRing(i,FF),x):i in [S41,S41d,S42,S42d,S61,S62]] eq EigsToTest;
  EigenvaluesOfElement(ChangeRing(MM,FF),x) eq E1;
  EigenvaluesOfElement(ChangeRing(MM,FF),x^5) eq E2;

end if;

El1:=[ 1155, 945, 1470, 735, 1050, 1260, 1575, 1260 ];

if(Verbosity gt 0) then
  "We now test that our given element of order 13 has the correct eigenvalues on L(G).";
  TestElt1:=ProduceGroupElement(T1785,El1);
  Eigenvalues(TestElt1) eq Eigs2;
end if;

NextEls1:=[
    [ 1659, 189, 1365, 1218, 210, 252, 1743, 252 ],
    [ 588, 546, 1008, 861, 1281, 966, 1743, 252 ],
    [ 588, 903, 1365, 504, 924, 1323, 1743, 252 ],
    [ 1659, 1617, 1008, 1575, 567, 1680, 1743, 252 ],
    [ 231, 1617, 651, 1575, 567, 252, 1386, 609 ],
    [ 231, 189, 1008, 1218, 210, 609, 1386, 609 ],
    [ 945, 546, 651, 861, 1281, 1323, 1386, 609 ],
    [ 945, 903, 1008, 504, 924, 1680, 1386, 609 ],
    [ 945, 189, 294, 504, 924, 252, 1029, 966 ],
    [ 1659, 546, 1722, 147, 210, 966, 1029, 966 ],
    [ 1659, 903, 294, 1575, 1638, 1323, 1029, 966 ],
    [ 945, 1617, 1722, 861, 1281, 1680, 1029, 966 ],
    [ 1302, 1617, 1365, 861, 1281, 252, 672, 1323 ],
    [ 1302, 189, 1722, 504, 924, 609, 672, 1323 ],
    [ 231, 546, 1365, 147, 210, 1323, 672, 1323 ],
    [ 231, 903, 1722, 1575, 1638, 1680, 672, 1323 ]
];

if(Verbosity gt 1) then

  "Finding all elements that are of order 85, preimages of our given element of order 13, and have the right eigenvalues on L(G). This takes some time.";
  Els1:=ConstructPreimages(El1,1785,5);
  e11:={*i[1]^^i[2]:i in Eigs1*};
  e1:={*R!((Log(i)*1785) div (q-1)):i in e11*};

  NextEls:=[];
  for ii in [1..#Els1] do
    if e1 eq {*&+[Els1[ii,i]*X[i,j]:i in [1..8]]:j in [1..248]*} then Append(~NextEls,Els1[ii]); end if;
    if(ii mod 10000 eq 0) then "At stage",ii,"of",#Els1,"there are",#NextEls,"elements with the correct eigenvalues"; end if;
  end for;

  NextEls1 eq NextEls;
end if;

gg1:=ProduceGroupElements(T1785,NextEls1);

if(Verbosity eq 1) then
  "Testing that elements of order 85 do at least have the correct eigenvalues on L(G), but not that these are the only ones.";

  {Eigenvalues(i):i in gg1} eq {Eigs1};
end if;

if(Verbosity ge 1) then

// It takes ages to search through E to find the right elements, so here they are.

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

  g:=gg1[1]^5;
  Orbs:={{i}:i in [1..#NextEls1]};
  for e in EltsToMerge do
    for i in [1..#NextEls1] do
      j:=Position(gg1,gg1[i]^e);
      Orb1:=FindOrb(Orbs,i);
      if(not(j in Orb1)) then
        Orb2:=FindOrb(Orbs,j);
        Exclude(~Orbs,Orb1);
        Exclude(~Orbs,Orb2);
        Include(~Orbs,Orb1 join Orb2); //printf "Merged %o and %o.\n",Orb1,Orb2;
        delete Orb1; delete Orb2; if(#Orbs eq 2) then "All merged up to powers"; break e; end if;
      end if;
    end for;
    "So far we have merged elements into",#Orbs,"orbits";
  end for;
  "Orbits are",Orbs;
  Orbs:=SetToSequence(Orbs);
  Position(gg1,gg1[Random(Orbs[1])]^69) in Orbs[2];

  El2:=NextEls1[1];
  Els2:=ConstructPreimages(El2,1785,3);
  Remove(~Els2,957); //Removes the element of order 85.
  E255:=[{*&+[ii[i]*X[i,j]:i in [1..8]]:j in [1..248]*}:ii in Els2];
  E255Powers:=[{* 3*i:i in MultisetToSet(j) *}:j in E255];
  Eigs255Powers:=[{*ww^(Integers()!i div 7):i in j*}:j in E255Powers];


  "For each of S41,S41d,S42,S42d,S61,S62, we find the number of elements of order 255 stabilizing their constituent eigenspaces";
  for i in EigsToTest do
    #[j:j in [1..6560]|&*[Multiplicity(Eigs255Powers[j],k[1]):k in i] eq 1];
  end for;

  "For each of S41,S41d,S42,S42d,S61,S62, we find the number of elements of order 255 stabilizing their constituent eigenspaces";
  for i in EigsToTest do
    [j:j in [1..6560]|&*[Multiplicity(Eigs255Powers[j],k[1]):k in i] eq 1][1];
  end for;

El3:=Els2[109]; //This stabs both 6.1 and 4.1
g:=ProduceGroupElement(T1785,El3);
n:=#Eigenvalues(g);
Els3:=ConstructPreimages(El3,1785,7);
for ii in Els3 do
 nn:=#{&+[ii[i]*X[i,j]:i in [1..8]]:j in [1..248]};
 if(nn eq n) then "Here is an element of order 1785:"; ii; break ii; end if; 
end for;

// This is an element that seventh powers to Els2[109] and has the same number of eigenvlaues as Els2[109].
// [ 334, 1284, 1595, 143, 95, 12, 83, 12 ]

end if;

// Now we want to compute the radicals from the article.

if(Verbosity ge 4) then
  // Case 5 first.

  A1:=S141;
  for i in [S1,S62,S2421d,S2421,S141,S1,S62,S141] do
    E,rho:=Ext(A1,i); A1:=MaximalExtension(A1,i,E,rho); DescribeLayers(Dual(A1),Irr,labs);
  end for;
  for i in [S1,S62,S141,S2421] do Dimension(Ext(A1,i)) eq 0; end for; // Don't need other 24 as module is graph-stable.

  // Now we can check the second module:
  A2:=RemoveFromTop(Dual(A1),Irr,Irr,[7,8]); DescribeLayers(A2,Irr,labs);

  A3:=Dual(RemoveFromTop(A2,Irr,Irr,[1,16])); DescribeLayers(A3,Irr,labs);
  for i in [S1,S62,S141] do Dimension(Ext(i,A3)) eq 0; end for;

  A4:=RemoveFromTop(Dual(A3),Irr,Irr,[17]); E,rho:=Ext(S141,A4); A4:=Extension(S141,A4,E.1,rho); DescribeLayers(A4,Irr,labs);
  for i in [S1,S62,S141] do Dimension(Ext(i,A4)) eq 0; end for;

  // Now case 6  

  B1:=S2412d;
  for i in [S1621,S61,S1,S2421,S141,S142,S42d,S61,S62,S1,S41,S62,S2421] do
    E,rho:=Ext(B1,i); B1:=MaximalExtension(B1,i,E,rho); DescribeLayers(Dual(B1),Irr,labs);
  end for;
  B1:=Dual(B1);
  for i in [S1,S41,S41d,S42,S42d,S61,S62,S141,S142,S1612,S1612d,S1621,S1621d,S2421,S2421d] do Dimension(Ext(i,B1)) eq 0; end for;

  B2a:=RemoveFromTop(B1,Irr,Irr,[1..15]); B2:=Dual(B1/B2a);
  DescribeLayers(B2,Irr,labs);

  B2c:=RemoveFromTop(B1,Irr,Irr,[1] cat [6..15]); B2b:=Dual(B1/B2c);
  DescribeLayers(B2b,Irr,labs);


// Now the radicals with socle S1612 and S1621.

  B3:=S1612d;
  for i in [S41,S42] do
    E,rho:=Ext(B3,i); B3:=MaximalExtension(B3,i,E,rho); DescribeLayers(Dual(B3),Irr,labs);
  end for;
  &and[Dimension(Ext(i,Dual(B3))) eq 0:i in [S41,S41d,S42,S42d,S61,S62,S1612,S1612d,S1621,S1621d]];

  &and[Dimension(Ext(S1621,i)) eq 0:i in [S41,S41d,S42,S42d,S61,S62,S1612,S1612d,S1621,S1621d]];


end if;