"Verbosity 0 just checks that there are elements that power to a given element of order 65 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 65 that power to it have the correct eigenvalues on L(G), and that those elements are all conjugate under the Weyl group.";
printf("\n");
"Verbosity 2 in addition checks that the 36 elements of order 65 are the only ones with the correct eigenvalues on L(G)";
printf("\n");
"Verbosity 3 in addition constructs PSU(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, and that the eigenvalues we are testing for our elements of orders 65 and 13 are correct.";
printf("\n");
"Verbosity 4 in addition computes the radicals as given in the text.";
printf("\n");
"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:=2^12;
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>;

T1365:=sub<T|[T.i^3:i in [1..8]]>;

R:=ResidueClassRing(1365);
X:=[[R!((Log((T1365.j)[i,i])*1365) div 4095):i in [1..248]]:j in [1..8]];


EigsToTest:=[
    { <w^3843, 1>, <w^3276, 1>, <w^1008, 1>, <w^63, 1> },
    { <w^4032, 1>, <w^252, 1>, <w^819, 1>, <w^3087, 1> },
    { <w^3591, 1>, <w^126, 1>, <w^2016, 1>, <w^2457, 1> },
    { <w^3969, 1>, <w^1638, 1>, <w^2079, 1>, <w^504, 1> },
    { <w^3906, 1>, <w^756, 1>, <w^3339, 1>, <w^1071, 1>, <w^3024, 1>, <w^189, 1> },
    { <w^3717, 1>, <w^2583, 1>, <w^1512, 1>, <w^2142, 1>, <w^378, 1>, <w^1953, 1> }
];

Eigs1:={ <w^2835, 5>, <w^1638, 3>, <w^1071, 3>, <w^2016, 4>, <w^1197, 3>, <w^819, 3>, <w^2772, 3>, <w^3276, 3>, <w^2457,
3>, <w^2205, 5>, <w^3780, 5>, <w^3591, 4>, <w^4032, 4>, <w^2709, 3>, <w^1890, 5>, <w^2646, 3>, <w^3969, 4>, <w^693,
3>, <w^882, 4>, <w^3150, 5>, <w^378, 3>, <w^2898, 3>, <w^3402, 3>, <w^1953, 3>, <w^3087, 4>, <w^2268, 4>, <w^2520,
5>, <w^1575, 5>, <w^756, 3>, <w^441, 4>, <w^1827, 4>, <w^2961, 4>, <w^315, 5>, <w^1323, 3>, <w^3465, 5>, <w^3906,
3>, <w^3528, 4>, <w^1386, 3>, <w^3843, 4>, <w^1449, 3>, <w^567, 4>, <w^3717, 3>, <w^189, 3>, <w^1008, 4>, <w^3654,
4>, <w^2394, 3>, <w^1134, 4>, <w^504, 4>, <w^252, 4>, <w^126, 4>, <w^3213, 4>, <w^630, 5>, <w^945, 5>, <w^63, 4>,
<w^1764, 4>, <w^3339, 3>, <w^3024, 3>, <w^2142, 3>, <w^1701, 3>, <w^1512, 3>, <w^2079, 4>, <w^2331, 4>, <w^2583,
3>, <1, 8>, <w^1260, 5> };
Eigs2:={ <w^1575, 19>, <w^2205, 19>, <w^3150, 19>, <w^3780, 19>, <w^2835, 19>, <w^1260, 19>, <w^945, 19>, <w^3465, 19>,
<1, 20>, <w^630, 19>, <w^315, 19>, <w^1890, 19>, <w^2520, 19> };


if(Verbosity ge 3) then
  "Checking that eigenvalues of simple modules on PSU(4,4) are correct";
  G2:=PSU(4,4);
  CC:=ConjugacyClasses(G2);
  for i in CompositionFactors(PermutationModule(G2,GF(16))) do if(Dimension(i) eq 4) then S41:=i; break i; end if; end for;
  S1:=TrivialModule(G2,GF(16));
  S41d:=Dual(S41);
  S42:=Socle(SymmetricSquare(S41));
  S42d:=Dual(S42);
  S61:=ExteriorSquare(S41);
  S62:=ExteriorSquare(S42);
  Irr:=[S1,S41,S41d,S42,S42d,S61,S62];
  Irr:=AppendWithoutDuplicates(Irr,CompositionFactors(ExteriorSquare(S61)));
  Irr:=AppendWithoutDuplicates(Irr,CompositionFactors(ExteriorSquare(S62)));
  S141:=Irr[8]; S142:=Irr[9];
  S1612:=TensorProduct(S41,S42);
  S1621:=TensorProduct(S41,S42d);
  S2412:=TensorProduct(S41,S62);
  S2421:=TensorProduct(S61,S42);
  S1612d:=Dual(S1612);
  S1621d:=Dual(S1621);
  S2412d:=Dual(S2412);
  S2421d:=Dual(S2421);

  Irr cat:=[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}^*"];

  for i in [79..94] do
    if(EigenvaluesOfElement(ChangeRing(S41,F),CC[i,3]) eq { <w^3843, 1>, <w^3276, 1>, <w^1008, 1>, <w^63, 1> }) then x:=CC[i,3]; end if;
  end for;
  [EigenvaluesOfElement(ChangeRing(i,F),x):i in [S41,S41d,S42,S42d,S61,S62]] eq EigsToTest;

  MM:=DirectSum([S1,S1,S41,S41,S42,S42,S61,S62,S1612,S1621,S2412,S2421]);
  MM:=DirectSum([MM,Dual(MM),S141,S142]);
  EigenvaluesOfElement(ChangeRing(MM,F),x) eq Eigs1;
  EigenvaluesOfElement(ChangeRing(MM,F),x^5) eq Eigs2;
end if;


El1:=[1155,1155,630,840,630,525,840,630];

if(Verbosity eq 1) then
  "Check our element of order 13 has correct eigenvalues";
  Eigenvalues(ProduceGroupElement(T1365,El1)) eq Eigs2;
end if;

Els1:=ConstructPreimages(El1,1365,5);
e11:={*i[1]^^i[2]:i in Eigs1*};
e1:={*R!((Log(i)*1365) div 4095):i in e11*};

NextEls1:=[[1323,231,399,168,672,105,168,126],[1323,1050,945,168,672,105,168,126],[231,231,399,1260,672,378,168,126],
[231,1050,945,1260,672,378,168,126],[777,504,126,987,945,378,1260,126],[777,1050,945,987,945,378,1260,126],
[231,504,126,168,945,1197,1260,126],[231,1050,945,168,945,1197,1260,126],[777,504,126,987,1218,105,987,399],
[777,1050,945,987,1218,105,987,399],[504,504,126,1260,1218,1197,987,399],[504,1050,945,1260,1218,1197,987,399],
[504,231,399,987,945,378,1260,399],[504,1050,945,987,945,378,1260,399],[1323,231,399,168,945,1197,1260,399],
[1323,1050,945,168,945,1197,1260,399],[777,504,126,168,672,105,168,672],[777,231,399,168,672,105,168,672],
[1050,504,126,1260,672,378,168,672],[1050,231,399,1260,672,378,168,672],[504,231,399,987,1218,105,987,672],
[504,1050,945,987,1218,105,987,672],[231,231,399,1260,1218,1197,987,672],[231,1050,945,1260,1218,1197,987,672],
[1323,504,126,987,945,378,1260,945],[1323,231,399,987,945,378,1260,945],[777,504,126,168,945,1197,1260,945],
[777,231,399,168,945,1197,1260,945],[231,504,126,168,672,105,168,1218],[231,1050,945,168,672,105,168,1218],
[504,504,126,1260,672,378,168,1218],[504,1050,945,1260,672,378,168,1218],[1323,504,126,987,1218,105,987,1218],
[1323,231,399,987,1218,105,987,1218],[1050,504,126,1260,1218,1197,987,1218],[1050,231,399,1260,1218,1197,987,1218]];

if(Verbosity ge 2) then
  NextElsa:=[];
  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(~NextElsa,Els1[ii]); end if;
    if(ii mod 10000 eq 0) then "At stage",ii,"of",#Els1,"there are",#NextElsa,"elements with the correct eigenvalues"; end if;
  end for;
  NextEls1 eq NextElsa;
end if;

if(Verbosity ge 1) then 
  gg1:=ProduceGroupElements(T1365,NextEls1);
  g:=gg1[1]^5;
  Orbs:={{i}:i in [1..#NextEls1]};
  EltsToMerge:=[E.1*E.3*E.1*E.4*E.2*E.3*E.1*E.4*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.5*E.6*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.1*E.7*E.8*E.7*E.6*E.5*E.4*E.3,E.3*E.4*E.2*E.3*E.5*E.4*E.2*E.3*E.6*E.5*E.7*E.6*E.5*E.8*E.7*E.6*E.5*E.4*E.3];
  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 1) then "All merged"; break e; end if;
      end if;
    end for;
    "So far we have merged elements into",#Orbs,"orbits";
  end for;

  El2:=NextEls1[1];
  Els2:=ConstructPreimages(El2,1365,3);
  Remove(~Els2,2176); //Removes the element of order 65.
  E195:=[{*&+[ii[i]*X[i,j]:i in [1..8]]:j in [1..248]*}:ii in Els2];
  E195Powers:=[{*3*i:i in MultisetToSet(j)*}:j in E195];
  Eigs195Powers:=[{*w^(3*Integers()!i):i in j*}:j in E195Powers];

  "For each of S41,S41d,S42,S42d,S61,S62, we find the number of elements of order 195 stabilizing their constituent eigenspaces";
  for i in EigsToTest do
    nn:=[j:j in [1..6560]|&*[Multiplicity(Eigs195Powers[j],k[1]):k in i] eq 1]; nn ne [];
    if(Position(EigsToTest,i) eq 1) then 2754 in nn; end if;
    if(Position(EigsToTest,i) eq 5) then 57 in nn; end if;
  end for;

  El31:=[843,596,1189,593,877,525,593,71];
  El32:=[1363,11,19,1308,32,5,8,6];
  [R!7*i:i in El31] eq Els2[2754];
  [R!7*i:i in El32] eq Els2[57];
  #Eigenvalues(ProduceGroupElement(T1365,El31)) eq #Eigenvalues(ProduceGroupElement(T1365,Els2[2754]));
  #Eigenvalues(ProduceGroupElement(T1365,El32)) eq #Eigenvalues(ProduceGroupElement(T1365,Els2[57]));
end if;

if(Verbosity ge 4) then
  // Start with Case 5
  A1:=S141;
  for i in [S62,S1,S2421,S2421d,S141,S62,S1,S141] do
    E,rho:=Ext(A1,i); A1:=MaximalExtension(A1,i,E,rho); DescribeLayers(Dual(A1),Irr,labs);
  end for;
  for i in [1,7,8,16,17] do Dimension(Ext(A1,Irr[i])) eq 0; end for;  

  A2:=RemoveFromTop(A1,Irr,Irr,[1..16]); DescribeLayers(Dual(A2),Irr,labs);
  for i in [1,7,8] do Dimension(Ext(A2,Irr[i])) eq 0; end for;  

  A3:=JacobsonRadical(A1); A3:=RemoveFromTop(A3,Irr,Irr,[1,7,16,17]);
  A3:=Dual(A1/A3); DescribeLayers(A3,Irr,labs);
  for i in [1,16,17] do Dimension(Ext(Irr[i],A3)) eq 0; end for;  

  // Now Case 6.
  B1:=S2412d;
  for i in [S1621,S61,S1,S141,S142,S42,S42d,S61,S1] do
    E,rho:=Ext(B1,i); B1:=MaximalExtension(B1,i,E,rho); DescribeLayers(Dual(B1),Irr,labs);
  end for;
  for i in [1..13] cat [16,17] do Dimension(Ext(B1,Irr[i])) eq 0; end for;

  B2:=Dual(B1);
  E,rho:=Ext(S2412d,B2); B2:=Extension(S2412d,B2,E.1,rho); DescribeLayers(B2,Irr,labs);
  B2a:=RemoveFromTop(B2,Irr,Irr,[1..13]); DescribeLayers(B2a,Irr,labs);

  E,rho:=Ext(S41d,S1612); B3:=Extension(S41d,S1612,E.1,rho);
  E,rho:=Ext(S42,S1621); B4:=Extension(S42,S1621,E.1,rho);
  for i in [2,3,4,5,6,7,10,11,12,13] do
    Dimension(Ext(Irr[i],B3)) eq 0;
    Dimension(Ext(Irr[i],B4)) eq 0;
  end for;

end if;