/* In this case we let H denote a copy of PSL(2,11), L denote the subgroup 11.5, and let G denote E6(331).
   We give matrices defining L, an element of order 11 and an element of order 5. 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 work form there, but we
   prefer a less synthetic proof.
*/

q:=331;
F:=GF(q);
zz:=PrimitiveElement(F);
z:=zz^6;

l1:=GL(27,F)![[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,330,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,267,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,267,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,181,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,8],
[0,0,0,0,0,0,0,0,0,0,0,181,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[207,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,8,0,0,0],
[0,0,0,0,0,0,0,0,64,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,267,0,0],
[0,0,0,0,0,0,0,0,0,0,150,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,267,0],
[0,0,0,0,181,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,150,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,207,0,0,0,0],
[0,0,0,0,0,0,0,181,0,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,207,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,0,0,0,0,0,0,0,267,0,0,0,0,0],
[0,0,0,323,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,207,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,323,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,207,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,330,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,330,0,0,0,0,0,0,0,0,0,0,0,0],
[0,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]];

l2:=GL(27,F)![[180,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,80,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,74,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,74,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,293,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,111,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,120,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,274,0,0,0,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,274,0,0,0,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,270,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,167,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,274,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,120,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,80,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,85,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,80,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,293,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,111,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,293,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,270,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,180,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,167,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,167,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,85,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,74]];

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

h1:=GL(27,F)![[249,246,0,0,141,0,288,0,0,0,0,0,315,125,0,0,176,0,0,278,0,194,0,0,0,0,190],
[120,22,0,0,6,0,6,0,0,0,0,0,82,13,0,0,283,0,0,282,0,112,0,0,0,0,162],
[0,0,198,0,0,41,0,55,144,0,27,128,0,0,8,138,0,0,180,0,0,0,236,277,0,233,0],
[0,0,0,122,0,0,0,0,0,17,0,0,0,0,0,0,0,273,0,0,25,0,0,0,96,0,0],
[278,82,0,0,36,0,234,0,0,0,0,0,247,162,0,0,263,0,0,283,0,316,0,0,0,0,6],
[0,0,204,0,0,153,0,21,211,0,215,25,0,0,151,98,0,0,310,0,0,0,58,61,0,76,0],
[170,190,0,0,300,0,22,0,0,0,0,0,275,18,0,0,303,0,0,93,0,93,0,0,0,0,141],
[0,0,55,0,0,233,0,180,172,0,38,207,0,0,290,101,0,0,153,0,0,0,31,260,0,98,0],
[0,0,317,0,0,83,0,250,0,0,95,248,0,0,2,81,0,0,217,0,0,0,16,234,0,31,0],
[0,0,0,195,0,0,0,0,0,131,0,0,0,0,0,0,0,186,0,0,241,0,0,0,136,0,0],
[0,0,146,0,0,200,0,156,261,0,0,276,0,0,25,185,0,0,76,0,0,0,55,185,0,55,0],
[0,0,98,0,0,121,0,174,120,0,151,76,0,0,94,287,0,0,184,0,0,0,231,160,0,65,0],
[81,82,0,0,10,0,270,0,0,0,0,0,303,82,0,0,33,0,0,325,0,283,0,0,0,0,161],
[263,227,0,0,107,0,316,0,0,0,0,0,6,249,0,0,93,0,0,97,0,6,0,0,0,0,84],
[0,0,282,0,0,151,0,127,316,0,65,241,0,0,55,311,0,0,310,0,0,0,175,95,0,178,0],
[0,0,220,0,0,128,0,185,159,0,304,122,0,0,136,55,0,0,90,0,0,0,24,143,0,207,0],
[251,135,0,0,125,0,303,0,0,0,0,0,52,87,0,0,53,0,0,161,0,6,0,0,0,0,83],
[0,0,0,90,0,0,0,0,0,96,0,0,0,0,0,0,0,233,0,0,314,0,0,0,276,0,0],
[0,0,180,0,0,98,0,153,279,0,78,139,0,0,98,237,0,0,76,0,0,0,1,109,0,83,0],
[278,214,0,0,59,0,325,0,0,0,0,0,87,31,0,0,43,0,0,303,0,247,0,0,0,0,52],
[0,0,0,276,0,0,0,0,0,58,0,0,0,0,0,0,0,98,0,0,106,0,0,0,98,0,0],
[103,83,0,0,18,0,48,0,0,0,0,0,59,244,0,0,238,0,0,10,0,36,0,0,0,0,125],
[0,0,168,0,0,71,0,95,205,0,180,226,0,0,277,183,0,0,270,0,0,0,180,311,0,222,0],
[0,0,175,0,0,330,0,273,246,0,304,314,0,0,31,215,0,0,131,0,0,0,136,153,0,212,0],
[0,0,0,186,0,0,0,0,0,236,0,0,0,0,0,0,0,25,0,0,98,0,0,0,71,0,0],
[0,0,21,0,0,76,0,310,298,0,180,142,0,0,178,174,0,0,195,0,0,0,200,308,0,198,0],
[6,107,0,0,238,0,278,0,0,0,0,0,43,321,0,0,112,0,0,33,0,263,0,0,0,0,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,m1,m2,m3,m4,n1,n2,n3,n4,n5,n6,n7,n8,n9>:=PolynomialRing(F,15);
mats:=MatrixRing(R,27);
scal:=mats!DiagonalMatrix([R!0:i in [1..27]]);
scal[9,9]:=a; scal[11,11]:=b;
a1:=[1,20,22,17,7];
a2:=[23,6,12,26,15];
a3:=[27,2,14,13,5];
a4:=[3,16,8,24,19];
a5:=[4,18,10,25,21];
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;
for i in [1..#a1] do
  scal[a3[i],a3[i]]:=n1; scal[a3[i],a4[i]]:=n2; scal[a3[i],a5[i]]:=n3;
  scal[a4[i],a3[i]]:=n4; scal[a4[i],a4[i]]:=n5; scal[a4[i],a5[i]]:=n6;
  scal[a5[i],a3[i]]:=n7; scal[a5[i],a4[i]]:=n8; scal[a5[i],a5[i]]:=n9;
end for;
scal[1,23]:=207*m2;
scal[7,15]:=207*m2;
scal[17,26]:=124*m2;
scal[22,12]:=323*m2;
scal[6,20]:=207*m3;
scal[12,22]:=181*m3;
scal[26,17]:=-m3;

scal[27,3]:=207*n2;
scal[2,16]:=-n2;
scal[14,8]:=8*n2;
scal[13,24]:=124*n2;

scal[27,4]:=-n3;
scal[2,18]:=-64*n3;
scal[14,10]:=207*n3;
scal[13,25]:=150*n3;

scal[3,27]:=8*n4;
scal[16,2]:=-n4;
scal[8,14]:=207*n4;
scal[24,13]:=323*n4;

scal[3,4]:=-8*n6;
scal[16,18]:=64*n6;
scal[8,10]:=-181*n6;
scal[24,25]:=124*n6;

scal[4,27]:=-n7;
scal[18,2]:=-150*n7;
scal[10,14]:=8*n7;
scal[25,13]:=64*n7;

scal[4,3]:=124*n8;
scal[18,16]:=150*n8;
scal[10,8]:=64*n8;
scal[25,24]:=323*n8;

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(l2)]|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;

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)*scal; v2:=Matricise(v)*scal; w2:=Matricise(w)*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;

coeffs0:=[a,b,m1,m2,m3,m4,n1,n2,n3,n4,n5,n6,n7,n8,n9];

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;

// Now we prove that all possible images h1 under scal that lie in E6 are conjugate under C_E6(L).

function CheckDetermination()

// Our first goal is to determine the restriction on scal so that h1^scal lies in E6.
// To do this, we will work modulo C_GL(H), which of course we are allowed to do.

// We then choose a particular matrix that works, called sca1, and set h11=h1^sca1 in E6. We then
// let sca10 denote all matrices (mod C_GL(H)) that map h1 into E6, and write scasub for sca1^-1*sca10.

// Then scasub is (mod C_GL(H)) the set of all elements that map h11 into E6.

// We then compute C_E6(L). Our final goal is to prove that scasub = C_E6(L). This is only true modulo
// C_GL(H), so we aim to prove that there is a bijection f between scasub and C_E6(L) such that x * f(x)
// lies in C_GL(H). because C_E6(H) has order 3, this is not true on the level of finite fields, as
// scasub and C_E6(L) both have order (q-1)^3 but C_E6(L) contains C_E6(H) of order 3. We show that, at
// the level of algebraic groups though, it is true.

// Start computing the restrictions on scal.

// Suppose that n7 eq 0. This is not true.

badcoeffs1:=[a,b,m1,m2,m3,m4,n1,n2,n3,n4,n5,n6,0,n8,n9];

Evaluate(ProdRel([2,21,21]),badcoeffs1) eq 163*n8*n9*(m1+286*n1);
Evaluate(ProdRel([1,4,10]),badcoeffs1) eq 174*n8*n9*(m1+196*n1);
// Thus either n8=0 or n9=0.

badcoeffs2a:=[a,b,m1,m2,m3,m4,n1,n2,n3,n4,n5,n6,0,0,n9];
badcoeffs2b:=[a,b,m1,m2,m3,m4,n1,n2,n3,n4,n5,n6,0,n8,0];

Evaluate(ProdRel([2,5,18]),badcoeffs2a) eq 177*n2*n9*(m1+66*n1);// n9 cannot be 0
Evaluate(ProdRel([1,10,13]),badcoeffs2a) eq 328*n2*n9*(m1+99*n1);// so n2=0.

badcoeffs3a:=[a,b,m1,m2,m3,m4,n1,0,n3,n4,n5,n6,0,0,n9];

Evaluate(ProdRel([2,9,10]),badcoeffs3a) eq 122*n5*n9*(m1+37*n1);// n5,n9 cannot be 0
Evaluate(ProdRel([2,9,25]),badcoeffs3a) eq 236*n5*n9*(m1+45*n1);// so contradiction.

// Thus we cannot have both n7 and n8 eq to 0, so bad branch 2a closes. We now close bad branch 2b,
// and then n7 is always non-zero.

Evaluate(ProdRel([1,7,25]),badcoeffs2b) eq 33*n3*n8*(m1+11*n1);// n8 cannot be 0 by branch a
Evaluate(ProdRel([2,2,21]),badcoeffs2b) eq 306*n3*n8*(m1+111*n1);// so n3=0.

badcoeffs3b:=[a,b,m1,m2,m3,m4,n1,n2,0,n4,n5,n6,0,n8,0];

Evaluate(ProdRel([2,11,18]),badcoeffs3b) eq 220*n6*n8*(m1+101*n1);// n8,n6 cannot be 0
Evaluate(ProdRel([1,11,18]),badcoeffs3b) eq 170*n6*n8*(m1+103*n1);// so contradiction.

// Thus n7 ne 0 in any solution. Can therefore scale by a scalar element and assume that n7=1.

// We can actually scale by any element of C_GL(H). This is 3-dimensional.
H:=sub<GL(27,F)|l1,l2,h1>;
homs:=AHom(GModule(H),GModule(H));
// Check what the bijective maps here are:
scc:=m1*mats!homs.1+m2*mats!homs.2+m3*mats!homs.3;
Determinant(scc) eq m1^10*m2^12*m3^5;
// Note that
scal[9,9] eq a and scc[9,9] eq m2;
scal[25,13] eq 64*n7 and scc[25,25] eq m3;

// Thus we may assume that a=n7=1, and we still have one variable to play with.

coeffs1:=ChangeCoefficients(coeffs0,[a,n7],[1,1]);

194*Evaluate(ProdRel([4,6,21]),coeffs1) eq 305*m3*n8*n9-m4;

coeffs2:=ChangeCoefficient(coeffs1,m4,305*m3*n8*n9);
coeffs2 eq [1,b,m1,m2,m3,305*m3*n8*n9,n1,n2,n3,n4,n5,n6,1,n8,n9];

// We now need to confirm that neither n8 nor n9 is equal to zero. We start by assuming that both of them are zero.
bcoeffs3a:=ChangeCoefficients(coeffs2,[n8,n9],[0,0]);
Evaluate(ProdRel([1,4,10]),bcoeffs3a) eq 265*m2;
// If m2=m4=0 then det(scal)=0, so this isn't allowed. Thus exactly one of n8 and n9 is zero (or neither).

f1:=n6*n8+181*n5*n9;
Evaluate(ProdRel([3,25,25]),coeffs2)/241 eq f1;
// Thus if n8=0 then n5=0, and if n9=0 then n6=0.
Evaluate(ProdRel([1,4,4]),coeffs2) eq 247*(n1*n8*n9+n2*n9+n3*n8);
// Thus if n8=0 then n2=n5=0, and if n9=0 then n3=n6=0. Thus det(scal)=0 in both cases, and again we have a contradiction.

// Thus n8,n9 are both non-zero.

(Evaluate(ProdRel([3,4,4]),coeffs2)/78-f1)/86 eq n9*(n5-274*n4*n8);
// Thus n5=274*n4*n8.
coeffs3:=ChangeCoefficient(coeffs2,n5,274*n4*n8);

Evaluate(ProdRel([3,4,4]),coeffs3) eq 78*n8*(n6-56*n4*n9);
171*Evaluate(ProdRel([1,10,10]),coeffs3) eq 123*m1*n8*n9 + 169*n2*n9 + 107*n3*n8 - m2;
// This gives us n6 and m2.
coeffs4:=ChangeCoefficients(coeffs3,[m2,n6],[123*m1*n8*n9 + 169*n2*n9 + 107*n3*n8,56*n4*n9]);

// Note that if m3=0 then m4=0, so m3 cannot equal 0 either as then det(scal)=0.
// So, n8, n9, m3 cannot be 0.

Evaluate(ProdRel([4,15,23])*176-ProdRel([3,6,21])*15,coeffs4) eq 27*m3*n8*n9*(m3-49*n4);
// Thus m3=49*n4. In particular, n4 ne 0 as well.
coeffs5:=ChangeCoefficient(coeffs4,m3,49*n4);

Evaluate(ProdRel([3,11,11]),coeffs5) eq 98*n4^2*n8*n9*(n9-264*n4);
coeffs6:=ChangeCoefficient(coeffs5,n9,264*n4);

Evaluate(ProdRel([3,3,3]),coeffs6) eq 320*n4^3*n8*(n4-261*b*n8);
coeffs7:=ChangeCoefficient(coeffs6,n4,261*b*n8);

Evaluate(173*ProdRel([2,3,24])-72*ProdRel([2,6,23]),coeffs7) eq b^2*n8^3*(n3-133*b*n2);
coeffs8:=ChangeCoefficient(coeffs7,n3,133*b*n2);

Evaluate(ProdRel([2,6,12]),coeffs8)/304 eq b^3*n8^4*(m1-299*n1);
Evaluate(ProdRel([1,4,4]),coeffs8)/12 eq b*n8*(n2-61*n1*n8);
coeffs9:=ChangeCoefficients(coeffs8,[m1,n2],[299*n1,61*n1*n8]);
// In particular, we see that if n1=0 then m1=m2=0, so n1 ne 0. We now set it equal to 1 because of

scal[2,2] eq n1 and scc[2,2] eq m1;

// Thus by killing off our last variable from C_GL(H), we obtain:
coeffs10:=ChangeCoefficient(coeffs9,n1,1);

// In fact, coeffs9 already kills off all coefficients, but coeffs10 is what we shall use.
// We need a particular case from this to check whether H lies in E6.
coeffs11:=ChangeCoefficients(coeffs10,[b,n8],[1,1]);
sca1:=GL(27,F)!Evaluate(scal,coeffs11);
h11:=h1^sca1;
{CheckRel(i,h11):i in seqs} eq {0};

sca10:=Evaluate(scal,coeffs10);
scasub:=mats!sca1^-1*sca10;

// scasub is now the matrices sca10 transported so they map h11 to the other possible options for h11.


// Need that scasub*C_GL(H)=C_E6(L)*C_GL(H).
// Alternatively that (any element of scasub)*(some element of C_E6(L)) lies in C_GL(H) and vice versa.

// So now we need to compute C_E6(L).

ChRel([1,15,27]) eq 330*m1*m4*n1+1;
// Thus m1,m4,n1 are all non-zero.

ChRel([1,1,22]) eq 248*m1^2*m2; //so m2=0.
ChRel([1,3,15]) eq 323*m1*n4*m4; //so n4=0.
ChRel([1,4,15]) eq m1*n7*m4; //so n7=0.
ChRel([1,26,26]) eq 329*m1*m3*m4; //so m3=0.
ChRel([2,9,15]) eq 64*a*m4*n3; //so n3=0.
ChRel([2,11,15]) eq -b*m4*n2; //so n2=0.
ChRel([3,9,26]) eq 8*a*m4*n6; //so n6=0.
ChRel([4,11,26]) eq 124*b*m4*n8; //so n8=0.

co1:=[a,b,m1,0,0,m4,n1,0,0,0,n5,0,0,0,n9];

Evaluate(ChRel([1,17,26]),co1) eq m1^2*m4-1;
-Evaluate(ChRel([1,15,27]),co1) eq m1*m4*n1-1;

// Thus m1=n1.
co2:=[a,b,m1,0,0,m4,m1,0,0,0,n5,0,0,0,n9];

// Now we need to obtain the final relations.

-Evaluate(ChRel([1,15,27]),co2) eq m1^2*m4-1;
Evaluate(ChRel([3,11,26]),co2) eq b*m4*n5-1;
-Evaluate(ChRel([1,19,25]),co2) eq m1*n5*n9-1;
-Evaluate(ChRel([4,9,26]),co2) eq a*m4*n9-1;
//Then:
//a=m4^-1*n9^-1=m1^2*m1*n5=m1^3*n5
//b=m4^-1*n5^-1=m1^2*m1*n9=m1^3*n9
//m4=m1^-2=(n5*n9)^2.

co3:=[m1^3*n5,m1^3*n9,m1,0,0,n5^2*n9^2,m1,0,0,0,n5,0,0,0,n9];

// The remaining coefficients satisfy m1*n5*n9=1.

// Now we show that any element of scasub can be producted with some element of C_E6(L) to lie in C_GL(H) and
// vice versa, to make the bijection referred to at the start of this function.
// First, we label our elements of scasub with m1 and m2, and C_E6(L) with n1, n2 and n3, to distinguish them.
// Then we force the product of them to commute with h11 and find the relations.
ctest1:=ChangeCoefficients(coeffs10,[b,n8],[m1,m2]);
ctest2:=ChangeCoefficients(co3,[m1,n5,n9],[n1,n2,n3]);

// This is all elements X of C_GL(L) that map h11 into E6, modulo C_GL(H)
scat1:=mats!sca1^-1*Evaluate(scal,ctest1);

// This is all elements of C_E6(L), modulo C_GL(H).
scat2:=Evaluate(scal,ctest2);

scat12:=scat1*scat2;

scadiff:=scat12*mats!h11-mats!h11*scat12;
scaents:=Sort([scadiff[i,j]:i,j in [1..27]]);

// Here is a solution for m1 and m2 in terms of n1,n2,n3.
// Claim that m1=n2/n3=n1*n2^2, m2=n1/n2 is the solution

{Evaluate(i,[0,0,n1*n2^2,n1/n2,0,0,n1,n2,1/n1/n2,0,0,0,0,0,0]):i in scaents} eq {0};

// Now note that the map (n1,n2) -> (n1*n2^2,n1/n2) is a homomorphism on F_q^x x F_q^x, with kernel (z,z) for z a cube root of 1.
// (This assumes that z lies in F_q. If it does not, then the map is bijective.)

// Thus the image is exactly one third of F_q^x x F_q^x, of course a subgroup.

// If x,y in F_q^x, let d and e be cube roots of them. Note that

// (d*e^2,d*e^-1) -> (d*e^2*d^2*e^-2=d^3=x,d*e^2/d/e^-1=e^3=y),

// so as long as x,y have cube roots in F_q, they lie in the image. In particular, this map is surjective on the algebraic closure.
// (One also sees that because the kernel is finite.)

// Therefore X=C_E6(L) mod C_GL(H), as needed.

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(l2)] do if(h[aa,i] ne 0) then for j in [1..NumberOfRows(l2)] do if(h[bb,j] ne 0) then for k in [1..NumberOfRows(l2)] 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.";