/* In this case we let H denote a copy of PSL(2,19), L denote the subgroup 19.9, and let G denote E6(191).
   We give matrices defining L, an element of order 19 and an element of order 9. 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:=191;
F:=GF(q);
zz:=PrimitiveElement(F);
z:=zz^10;

l1:=GL(27,F)![[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,-1,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,-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,-1,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,-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,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,1],
[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,1,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,0,-1,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,-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,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,1,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,-1,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,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,1,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,-1,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,-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,1,0,0,0,0,0,0,0,0,0,0]];

l2:=GL(27,F)!DiagonalMatrix([153,107,150,125,5,177,32,30,177,25,52,160,121,180,32,125,36,136,25,36,180,160,6,69,153,154,30]);

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

h1:=GL(27,F)![[147,0,0,0,0,0,0,0,185,154,0,142,0,0,123,136,89,0,0,0,11,0,0,0,0,0,21],
[0,141,130,50,164,171,149,8,0,0,113,0,59,19,0,0,0,111,26,118,0,116,147,13,162,144,0],
[0,80,43,98,157,33,81,96,0,0,148,0,141,44,0,0,0,65,71,80,0,137,18,43,148,82,0],
[0,110,188,21,9,45,106,71,0,0,172,0,42,46,0,0,0,91,35,109,0,68,124,116,129,100,0],
[0,148,126,157,184,67,19,7,0,0,22,0,113,108,0,0,0,28,89,114,0,10,29,152,24,114,0],
[0,120,108,103,174,46,35,131,0,0,102,0,165,171,0,0,0,176,26,71,0,76,163,108,133,79,0],
[0,33,62,44,3,116,146,175,0,0,67,0,171,49,0,0,0,46,145,22,0,28,83,94,174,33,0],
[0,137,93,73,158,28,123,187,0,0,10,0,116,96,0,0,0,101,115,174,0,97,144,167,8,85,0],
[11,0,0,0,0,0,0,0,17,63,0,84,0,0,8,13,3,0,0,0,154,0,0,0,0,0,150],
[125,0,0,0,0,0,0,0,171,150,0,19,0,0,47,151,161,0,0,0,21,0,0,0,0,0,52],
[0,126,89,42,127,145,91,53,0,0,163,0,80,103,0,0,0,27,176,145,0,90,88,35,123,171,0],
[57,0,0,0,0,0,0,0,155,188,0,150,0,0,177,144,133,0,0,0,102,0,0,0,0,0,30],
[0,109,189,75,142,158,100,77,0,0,77,0,47,74,0,0,0,171,79,187,0,106,98,118,77,86,0],
[0,98,27,85,188,44,170,16,0,0,157,0,50,23,0,0,0,149,88,77,0,73,71,152,23,116,0],
[33,0,0,0,0,0,0,0,92,174,0,121,0,0,78,5,155,0,0,0,6,0,0,0,0,0,20],
[32,0,0,0,0,0,0,0,78,183,0,85,0,0,177,52,177,0,0,0,68,0,0,0,0,0,144],
[44,0,0,0,0,0,0,0,158,11,0,103,0,0,159,153,134,0,0,0,147,0,0,0,0,0,125],
[0,173,170,120,119,108,124,3,0,0,162,0,44,123,0,0,0,88,163,173,0,47,100,181,92,98,0],
[0,95,147,175,135,16,71,138,0,0,184,0,183,137,0,0,0,53,131,137,0,4,3,24,46,77,0],
[0,148,176,23,16,174,62,145,0,0,24,0,162,3,0,0,0,68,58,130,0,8,99,76,164,114,0],
[38,0,0,0,0,0,0,0,5,178,0,27,0,0,52,26,144,0,0,0,55,0,0,0,0,0,40],
[0,80,170,77,131,22,82,54,0,0,114,0,118,106,0,0,0,46,120,40,0,174,18,18,130,4,0],
[0,43,68,27,56,62,3,44,0,0,126,0,130,162,0,0,0,102,83,170,0,93,21,127,176,2,0],
[0,34,135,3,47,188,9,135,0,0,7,0,27,27,0,0,0,127,174,60,0,33,119,98,175,142,0],
[0,44,162,23,164,49,145,54,0,0,108,0,19,72,0,0,0,88,20,106,0,96,68,172,3,117,0],
[0,43,127,152,93,94,75,167,0,0,152,0,13,172,0,0,0,156,83,18,0,167,10,143,76,73,0],
[88,0,0,0,0,0,0,0,121,107,0,160,0,0,85,27,150,0,0,0,49,0,0,0,0,0,172]];

h2:=GL(27,F)![[147,0,0,0,0,0,0,0,185,154,0,142,0,0,123,136,89,0,0,0,11,0,0,0,0,0,21],
[0,148,157,86,99,76,64,77,0,0,14,0,65,29,0,0,0,128,86,143,0,45,171,71,75,69,0],
[0,126,89,39,20,168,40,9,0,0,94,0,148,176,0,0,0,52,74,143,0,101,124,136,79,113,0],
[0,41,111,100,104,36,57,67,0,0,186,0,142,30,0,0,0,36,135,67,0,88,141,190,45,103,0],
[0,94,27,7,137,99,65,71,0,0,187,0,14,34,0,0,0,183,165,181,0,170,185,17,138,85,0],
[0,9,65,6,136,73,135,63,0,0,2,0,155,137,0,0,0,173,146,38,0,142,125,139,124,64,0],
[0,57,37,157,41,120,155,13,0,0,37,0,94,106,0,0,0,117,118,21,0,47,176,110,155,31,0],
[0,140,70,120,75,47,103,14,0,0,160,0,121,152,0,0,0,105,49,126,0,89,49,22,121,84,0],
[11,0,0,0,0,0,0,0,17,63,0,84,0,0,8,13,3,0,0,0,154,0,0,0,0,0,150],
[125,0,0,0,0,0,0,0,171,150,0,19,0,0,47,151,161,0,0,0,21,0,0,0,0,0,52],
[0,139,183,13,29,7,86,137,0,0,8,0,63,108,0,0,0,1,148,22,0,163,128,79,113,43,0],
[57,0,0,0,0,0,0,0,155,188,0,150,0,0,177,144,133,0,0,0,102,0,0,0,0,0,30],
[0,78,174,104,190,170,2,183,0,0,106,0,122,112,0,0,0,43,68,175,0,54,183,13,17,73,0],
[0,3,23,180,84,157,91,89,0,0,74,0,36,48,0,0,0,190,185,21,0,120,5,126,100,183,0],
[33,0,0,0,0,0,0,0,92,174,0,121,0,0,78,5,155,0,0,0,6,0,0,0,0,0,20],
[32,0,0,0,0,0,0,0,78,183,0,85,0,0,177,52,177,0,0,0,68,0,0,0,0,0,144],
[44,0,0,0,0,0,0,0,158,11,0,103,0,0,159,153,134,0,0,0,147,0,0,0,0,0,125],
[0,67,136,126,157,4,114,50,0,0,6,0,20,79,0,0,0,128,97,58,0,127,79,106,26,183,0],
[0,14,54,102,115,178,67,97,0,0,68,0,141,186,0,0,0,84,63,43,0,177,92,190,86,161,0],
[0,153,177,100,45,155,146,105,0,0,40,0,138,105,0,0,0,6,67,174,0,121,189,70,144,175,0],
[38,0,0,0,0,0,0,0,5,178,0,27,0,0,52,26,144,0,0,0,55,0,0,0,0,0,40],
[0,11,12,21,186,21,124,148,0,0,58,0,11,135,0,0,0,13,153,32,0,126,69,98,174,60,0],
[0,89,44,108,186,99,85,62,0,0,27,0,157,151,0,0,0,8,110,156,0,146,55,182,9,17,0],
[0,171,5,54,4,40,15,158,0,0,54,0,92,12,0,0,0,29,49,65,0,171,157,187,179,190,0],
[0,87,41,48,146,106,161,5,0,0,32,0,61,39,0,0,0,168,54,135,0,152,38,54,105,153,0],
[0,136,182,110,4,93,13,13,0,0,17,0,71,129,0,0,0,112,103,128,0,95,85,38,146,178,0],
[88,0,0,0,0,0,0,0,121,107,0,160,0,0,85,27,150,0,0,0,49,0,0,0,0,0,172]];

// Here is the hom space for L, even as a subgroup of NGT. This makes computation incredibly easy

R<a,m1,m2,m3,m4>:=PolynomialRing(F,5);
mats:=MatrixRing(R,27);
scal:=mats!DiagonalMatrix([a:i in [1..27]]);
a1:=[1,12,10,15,21,17,27,9,16];
a2:=[25,22,19,7,14,20,8,6,4];
for i in [1..9] 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[16,4]:=-m2;
scal[12,22]:=-m2;
scal[10,19]:=-m2;
scal[15,7]:=-m2;
scal[4,16]:=-m3;
scal[22,12]:=-m3;
scal[19,10]:=-m3;
scal[7,15]:=-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);
g1:=mats!h1;
g2:=mats!h2;

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,g)

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 now check that there is a unique H over L for each module action.

function CheckDetermination()

ProdRel([1,1,1],g1) eq 64*(m1+m2)*(m1+88*m2)*(m1+102*m2);
ProdRel([1,1,1],g2) eq 64*(m1+m2)*(m1+88*m2)*(m1+102*m2);
// Since not both of m1 and m2 can be zero, we see that neither m1 nor m2 is zero.


ProdRel([1,1,8],g1)/91-51*ProdRel([1,4,16],g1) eq m1*m2*(a+75*m4);
ProdRel([1,1,8],g2)/52-175*ProdRel([1,4,16],g2) eq m1*m2*(a+66*m4);

c11:=[a,m1,m2,m3,a/-75];
c21:=[a,m1,m2,m3,a/-66];

Evaluate(ProdRel([2,2,2],g1),c11) eq 62*(a+10*m3)^2*(a+52*m3);
Evaluate(ProdRel([2,2,3],g1),c11) eq 190*(a+10*m3)^2*(a+16*m3);
// Thus a+10*m3=0.
Evaluate(ProdRel([2,2,2],g2),c21) eq 182*(a+78*m3)^2*(a+75*m3);
Evaluate(ProdRel([2,2,3],g2),c21) eq 104*(a+78*m3)^2*(a+154*m3);
c12:=[Evaluate(i,[a,m1,m2,-a/10,0]):i in c11];
c22:=[Evaluate(i,[a,m1,m2,-a/78,0]):i in c21];

Evaluate(ProdRel([1,2,4],g1),c12) eq 56*a^2*(m1+102*m2);
Evaluate(ProdRel([1,2,4],g2),c22) eq 32*a^2*(m1+88*m2);
c13:=[Evaluate(i,[a,m1,-m1/102,0,0]):i in c12];
c23:=[Evaluate(i,[a,m1,-m1/88,0,0]):i in c22];

c14:=[Evaluate(i,[1,1,0,0,0]):i in c13];
c24:=[Evaluate(i,[1,1,0,0,0]):i in c23];

h11:=h1^(GL(27,F)!Evaluate(scal,c14));
h22:=h2^(GL(27,F)!Evaluate(scal,c24));

// As the two 18-dimensional factors here are not Aut(H)-conjugate, one obtains two classes of H
// containing the given L. (Note that the two dual 9s are Aut(H)-conjugate.)

H1:=sub<GL(27,F)|l1,l2,h11>;
H2:=sub<GL(27,F)|l1,l2,h22>;

// Now check that H1 and H2 lie in E6.

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

// This proves that both of the copies of PSL(2,19) we constructed are contained in E6.

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.";
