/* In this case we let H denote a copy of PSL(2,17), L denote the subgroup 17.8, and let G denote E6(137).
   We give matrices defining L, an element of order 17 and an element of order 8. We construct L to lie in
   normalizer of a torus NGT, as constructed using the GroupOfLieType command. We give H as L and a matrix
   h1, chosen at random. A posteriori, one could produce H as a subgroup of E6 and show that no other
   overgroup of L preserves the E6 trilinear form, but we prefer a less synthetic proof.
*/

q:=137;
F:=GF(q);
zz:=PrimitiveElement(F);
z:=zz;

l1:=GL(27,F)![[136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,136,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0],
[0,0,136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,136,0,0,0,0,0,0,0],
[0,136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,136,0,0,0,0,0,0,0,0],
[0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0],
[0,0,0,0,0,0,0,0,0,0,136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,136,0,0,0,0],
[0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,136],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,136,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,136,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,136,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,136,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0]];

l2:=GL(27,F)!DiagonalMatrix([1,122,60,72,50,38,74,34,59,123,16,56,115,73,119,88,1,133,50,119,60,133,34,16,74,1,38]);

L:=sub<GL(27,F)|l1,l2>;

h1:=GL(27,F)![[102,0,10,0,106,58,66,45,0,0,95,0,0,0,23,0,0,82,124,40,97,54,94,106,127,0,42],
[0,97,76,130,31,32,93,18,37,130,82,47,125,115,77,127,14,1,124,15,30,133,72,83,39,29,128],
[97,9,108,37,114,53,47,39,97,38,123,18,64,102,32,109,28,63,69,60,50,124,101,86,14,75,61],
[0,44,28,49,36,122,48,28,12,84,125,10,118,93,77,35,4,1,7,15,112,133,112,48,82,67,77],
[66,102,43,37,65,2,24,77,31,14,111,122,129,135,94,15,101,41,47,103,8,126,32,1,65,21,93],
[52,118,33,16,129,53,38,108,34,13,123,76,134,4,8,109,29,1,39,94,69,33,50,93,40,128,60],
[55,11,53,114,111,16,4,53,130,63,55,32,34,66,93,15,107,122,59,135,60,129,61,28,132,86,12],
[58,4,129,17,64,111,29,114,81,12,65,31,68,74,47,101,131,131,112,72,39,74,69,54,39,68,59],
[0,98,124,35,125,74,126,99,136,55,9,7,58,76,81,84,64,107,89,14,85,120,122,101,44,113,22],
[0,135,133,46,63,71,65,35,97,76,37,22,58,52,8,93,104,44,115,135,121,98,3,126,14,49,10],
[125,83,1,92,131,122,58,63,129,30,116,1,16,93,11,136,118,109,63,96,104,69,13,63,37,91,8],
[0,67,103,124,56,7,7,40,101,136,133,100,129,97,13,12,96,8,87,31,1,105,23,16,109,119,28],
[0,31,121,76,120,23,123,100,124,35,56,7,76,46,126,49,86,45,69,37,73,94,126,50,56,1,92],
[0,133,3,76,69,103,89,73,129,58,2,12,51,58,32,118,60,121,2,129,12,64,18,129,55,46,1],
[55,49,100,14,89,36,112,91,136,105,101,48,69,112,26,134,51,81,72,77,14,116,14,67,69,107,82],
[0,31,118,106,4,11,49,9,70,39,78,97,4,2,134,93,111,83,16,35,61,79,36,99,78,22,44],
[0,87,8,77,107,91,2,116,121,70,76,96,117,74,108,90,33,117,17,110,32,80,53,107,129,37,90],
[13,88,99,123,108,133,74,90,7,126,68,44,89,65,65,48,27,79,39,24,40,100,14,126,76,105,132],
[52,94,69,112,80,100,73,2,42,72,77,99,135,68,34,38,128,112,65,114,43,96,77,26,113,108,2],
[83,78,123,81,65,55,69,123,4,128,67,82,135,100,136,12,70,116,83,26,44,63,51,109,11,120,109],
[10,105,123,112,30,98,119,32,127,78,46,73,16,94,77,130,7,95,114,135,108,74,39,14,90,53,53],
[106,115,66,72,126,26,39,119,101,37,25,126,12,18,24,125,96,109,29,64,38,79,47,68,74,8,4],
[54,1,11,107,7,115,126,30,123,3,48,42,17,87,65,128,67,81,64,57,129,6,114,72,108,17,111],
[134,82,62,114,56,68,28,42,2,61,81,34,133,11,96,103,39,30,6,112,136,109,74,116,58,80,15],
[48,100,99,40,22,102,111,39,36,87,36,129,60,52,135,99,76,68,26,100,84,122,84,55,4,47,121],
[0,53,41,91,86,19,113,101,6,4,47,119,38,75,9,103,37,44,70,32,27,98,130,86,96,104,76],
[124,98,30,4,11,38,66,123,77,106,37,19,102,1,43,130,110,75,129,68,33,136,108,14,99,32,53]];

// Here is the centralizer in GL(27,F) for L, even as a subgroup of NGT. This makes computation incredibly easy

R<a,b,n1,n2,m1,m2,m3,m4>:=PolynomialRing(F,8);
mats:=MatrixRing(R,27);
scal:=mats!DiagonalMatrix([a:i in [1..27]]);
scal[17,17]:=n1; scal[26,26]:=n1; scal[26,17]:=n2; scal[17,26]:=-n2;
for i in [2,4,9,10,12,13,14,16] do scal[i,i]:=b; end for;
a1:=[3,6,7,18,11,20,5,8];
a2:=[21,27,25,22,24,15,19,23];
for i in [1..#a1] do scal[a1[i],a1[i]]:=m1; scal[a1[i],a2[i]]:=m2; scal[a2[i],a1[i]]:=m3; scal[a2[i],a2[i]]:=m4; end for;
scal[20,15]:=-m2;
scal[11,24]:=-m2;
scal[18,22]:=-m2;
scal[7,25]:=-m2;

scal[27,6]:=-m3;
scal[21,3]:=-m3;
scal[23,8]:=-m3;
scal[19,5]:=-m3;


scal*mats!l1 eq mats!l1*scal and scal*mats!l2 eq mats!l2*scal;

seqs:=[[i,j,k]:i,j,k in [1..NumberOfRows(l1)]|i le j and j le k];

ents:=[300,322,340,354,359,609,637,661,681,697,901,917,947,1004,1020,1164,1182,1246,1283,1311,1399,1419,1521,1558,1577,1710,1744,1787,
1821,1922,1958,1989,2029,2184,2227,2260,2357,2388,2427,2586,2609,2710,2749,2859,2991];

vals:=[1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,1,-1,1,-1,-1,1,-1,1,-1,1,-1,1,1,-1,-1,1];

f27:=[0:i in [1..#seqs]];
for i in [1..#ents] do f27[ents[i]]:=vals[i]; end for;

function Matricise(v)
return Matrix(1,NumberOfColumns(v),[R!v[i]:i in [1..NumberOfColumns(v)]]);
end function;

V:=GModule(L);
g:=mats!h1;

intseqs1:=[seqs[i]:i in ents];
intseqs:=&join{{[i[1],i[2],i[3]],[i[1],i[3],i[2]],[i[2],i[1],i[3]],[i[2],i[3],i[1]],[i[3],i[1],i[2]],[i[3],i[2],i[1]]}:i in intseqs1};

function ProdRel(seq)

u:=V.seq[1]; v:=V.seq[2]; w:=V.seq[3];
u1:=Matricise(u)*scal; v1:=Matricise(v)*scal; w1:=Matricise(w)*scal;
u2:=Matricise(u)*g*scal; v2:=Matricise(v)*g*scal; w2:=Matricise(w)*g*scal;
e1:=&+[u1[1,i[1]]*v1[1,i[2]]*w1[1,i[3]]*f27[Position(seqs,Sort([i[1],i[2],i[3]]))]:i in intseqs];
e2:=&+[u2[1,i[1]]*v2[1,i[2]]*w2[1,i[3]]*f27[Position(seqs,Sort([i[1],i[2],i[3]]))]:i in intseqs];

return e1-e2;
end function;

function CheckRel(seq,gg)
u:=V.seq[1]; v:=V.seq[2]; w:=V.seq[3];
u2:=u^gg; v2:=v^gg; w2:=w^gg;
e1:=&+[u[i[1]]*v[i[2]]*w[i[3]]*f27[Position(seqs,Sort([i[1],i[2],i[3]]))]:i in intseqs];
e2:=&+[u2[i[1]]*v2[i[2]]*w2[i[3]]*f27[Position(seqs,Sort([i[1],i[2],i[3]]))]:i in intseqs];
return e1-e2;
end function;

// We need this to understand C_E6(L) at the end of the proof.

c6:=[a,1,n1,0,69*a*n1+69,135*a*n1^2+2*n1,120*a*n1+17,69*a*n1^2+69*n1];
sca:=Evaluate(scal,c6)*DiagonalMatrix([b:i in [1..27]]);

function ChRel(seq)

u:=V.seq[1]; v:=V.seq[2]; w:=V.seq[3];
u1:=Matricise(u); v1:=Matricise(v); w1:=Matricise(w);
u2:=Matricise(u)*sca; v2:=Matricise(v)*sca; w2:=Matricise(w)*sca;
e1:=&+[u1[1,i[1]]*v1[1,i[2]]*w1[1,i[3]]*f27[Position(seqs,Sort([i[1],i[2],i[3]]))]:i in intseqs];
e2:=&+[u2[1,i[1]]*v2[1,i[2]]*w2[1,i[3]]*f27[Position(seqs,Sort([i[1],i[2],i[3]]))]:i in intseqs];

return e1-e2;
end function;

coeffs0:=[a,b,n1,n2,m1,m2,m3,m4];

function ChangeCoefficient(coeffs,coeff,target)
coeffs2:=coeffs;
nn:=Position(coeffs0,coeff);
coeffs2[nn]:=target;
return [Evaluate(i,coeffs2):i in coeffs];
end function;

function ChangeCoefficients(coeffs,coeff,target)
coeffs2:=coeffs;
for i in [1..#coeff] do
  nn:=Position(coeffs0,coeff[i]);
  coeffs2[nn]:=target[i];
end for;
return [Evaluate(i,coeffs2):i in coeffs];
end function;

// We now determine the possible copies of H lying above L.

function CheckDetermination()
// Note that a,b are non-zero. We cannot have n1+37*n2=0 or n1+100*n2=0, and cannot have m1*m4+m2*m3=0,
// from the determinant of scal being non-zero.

ProdRel([2,2,2])/76-ProdRel([2,4,17])/67 eq 6*(m2+4*m4)*(b-m1+4*m3)^2;
// Thus m2+4*m4=0 or m1=b+4*m3. The latter in fact holds, so we proceed with the former, arriving at a contradiction.

badcoeffs1:=[a,b,n1,n2,m1,-4*m4,m3,m4];

Evaluate(ProdRel([2,2,16]),badcoeffs1) eq 16*b^2*(n1-n2);
badcoeffs2:=ChangeCoefficient(badcoeffs1,n2,n1);

Evaluate(87*ProdRel([2,12,17])-18*ProdRel([3,17,21]),badcoeffs2) eq n1*m1*m3;
Evaluate(ProdRel([2,10,14])/29-ProdRel([3,17,24])/41,badcoeffs2) eq 54*n1*(m1^2+11*m1*m3+16*m3^2);
// The first shows that m1=0 or m3=0. The second then implies that m1=m3=0. But then the determinant of scal is zero.
// Thus m1=b+4*m3

coeffs1:=[a,b,n1,n2,b+4*m3,m2,m3,m4];

Evaluate(ProdRel([2,4,16]),coeffs1) eq 15*b^2*n2;
// Since b ne 0, we have n2=0.
coeffs2:=ChangeCoefficient(coeffs1,n2,0);

Evaluate(ProdRel([2,2,12]),coeffs2) eq 99*b^2*(n1+34*m2-m4);
coeffs3:=ChangeCoefficient(coeffs2,m4,n1+34*m2);

// At this point it seems as though we need to scale scal so that b=1, which of course we may do.

coeffs4:=ChangeCoefficient(coeffs3,b,1);

Evaluate(ProdRel([1,2,26]),coeffs4) eq 109*(121*n1*m3-m2);
coeffs5:=ChangeCoefficient(coeffs4,m2,121*n1*m3);

Evaluate(ProdRel([1,2,12]),coeffs5) eq 49*n1*(120*a*n1+17-m3);
// Thus 120*a*n1+17=m3
coeffs6:=ChangeCoefficient(coeffs5,m3,120*a*n1+17);
coeffs6 eq c6;
// Now all coefficients are zero.

// This is the complete subgroup of GL(26,F) that conjugates H inside E6. (Note that H<E6 to begin with.)

sca:=Evaluate(scal,coeffs6)*DiagonalMatrix([b:i in [1..27]]);

// Now find the centralizer of H in GL(27,k):

scadiff:=sca*mats!h1-mats!h1*sca;
scaents:=Sort([scadiff[i,j]:i,j in [1..27]]);
scaents[170] eq 2*b*(n1-1);
{Evaluate(i,[a,b,1,n2,m1,m2,m3,m4]):i in scaents} eq {0};

// Since a,b,n1 are all non-zero, this means the centralizer of H is all matrices sca with n1=1, a,b arbitrary.

// Now obtain C_E6(L).

Determinant(sca) eq a^9*b^27*n1^18;
// Thus a,b,n1 cannot be 0.
ChRel([2,14,26]) eq 1-b^3*n1;
ChRel([2,6,12]) eq 2*b^3*n1*(1-a*n1);
// Thus a eq b^3 and n1=b^-3. In fact, this is enough to 
{Evaluate(ChRel(i),[b^3,b,1/b^3,0,0,0,0,0]):i in seqs} eq {0};

// Notice that C_GL(H)*C_E6(L) contains all those elements sca where n1 lies inside the subset {x^3:x in F} of F.
// (This is because they intersect in Z(E6).) Of course, this means that if F is algebraically closed we obtain
// that C_GL(H)*C_E6(L) is exactly all matrices sca. Over the finite group we obtain three orbits, permuted by
// the diagonal automorphism of E6(q).

// Finally, check that h1 lies in E6.

{CheckRel(i,h1):i in seqs} eq {0};

return "";
end function;

// We now give code that can prove that the trilinear form we gave is the correct one. Choose two random elements that
// generate the subgroup 2^6.W(E6), and then build the trilinear form relations. We also check that l1 and l2 lie in
// NGT.

function CheckE6Form()

GG:=GroupOfLieType("E6",q);
W:=VectorSpace(F,6);
rho:=StandardRepresentation(GG);
Over:=GL(27,q);
g1:=elt<GG|W![z,1,1,1,1,1]>;
g2:=elt<GG|W![1,z,1,1,1,1]>;
g3:=elt<GG|W![1,1,z,1,1,1]>;
g4:=elt<GG|W![1,1,1,z,1,1]>;
g5:=elt<GG|W![1,1,1,1,z,1]>;
g6:=elt<GG|W![1,1,1,1,1,z]>;
Refs:=Reflections(GG);
Mats:=[rho(i):i in [g1,g2,g3,g4,g5,g6]];
MatsE:=[rho(i):i in Refs];
NGT:=sub<Over|Mats cat MatsE>;
T:=sub<NGT|Mats>;
E:=sub<NGT|MatsE>;

l1 in NGT and l2 in NGT;

repeat y1:=Random(E); y2:=Random(E); until E eq sub<E|y1,y2>;

mat:=[];
for h in [y1,y2] do
  for nn in [1..#seqs] do aa:=seqs[nn,1]; bb:=seqs[nn,2]; cc:=seqs[nn,3];
    val:=[F!0:i in [1..#seqs]];
    for i in [1..NumberOfRows(l1)] do if(h[aa,i] ne 0) then for j in [1..NumberOfRows(l1)] do if(h[bb,j] ne 0) then for k in [1..NumberOfRows(l1)] do if(h[cc,k] ne 0) then
      val[Position(seqs,Sort([i,j,k]))]+:=h[aa,i]*h[bb,j]*h[cc,k];
    end if; end for; end if; end for; end if; end for;
    val[nn]-:=1; Append(~mat,val); delete val;
  end for;
end for;

ftest:=Nullspace(Transpose(Matrix(F,mat))).1;
return &and[ftest[i] eq f27[i]:i in [1..#seqs]];
end function;

"The function CheckE6Form() checks that the trilinear form f27 from E6 is correct";
"The function CheckDetermination() checks that the determination of the possible conjugates into E6 is correct.";
