/* In this case we let H denote a copy of PSL(2,11), L denote the subgroup 11.5, and let G denote E6(2^20).
   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:=2^20;
F<zz>:=GF(q);
z:=zz^19065;

l1:=GL(27,F)![[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,zz^209715,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[zz^838860,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,zz^838860,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,zz^209715,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,zz^838860,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,zz^838860,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,zz^209715,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,zz^209715,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,zz^209715,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,zz^209715,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,zz^838860,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,zz^209715,0],
[0,zz^838860,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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,zz^838860,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,zz^838860,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,zz^209715,0,0,0,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,zz^419430,0,0,0],
[0,0,0,0,0,0,zz^629145,0,0,0,0,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,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,0,0,0,0,0,1,0,0,0,0,0,0,0,0]];

l2:=GL(27,F)!DiagonalMatrix([zz^762600,zz^381300,zz^571950,1,1,zz^190650,zz^667275,zz^190650,zz^667275,zz^667275,zz^95325,zz^857925,zz^285975,zz^285975,zz^476625,
zz^476625,zz^953250,zz^953250,zz^95325,zz^953250,zz^571950,zz^571950,zz^762600,zz^762600,zz^190650,zz^381300,zz^857925]);

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

h1:=GL(27,F)![[zz^169125,0,zz^514140,0,0,0,0,0,0,zz^676500,0,0,0,0,0,0,zz^547965,0,0,0,0,0,0,0,zz^818565,0,0],
[0,zz^744150,0,0,0,zz^1046037,zz^1018977,0,0,0,zz^372075,zz^40590,0,zz^920040,zz^649440,0,0,zz^734847,0,0,zz^119232,0,zz^592782,0,0,0,0],
[zz^94710,0,zz^338250,0,0,0,0,0,0,zz^1007985,0,0,0,0,0,0,zz^608850,0,0,0,0,0,0,0,zz^676500,0,0],
[0,0,0,0,zz^947281,0,0,zz^1043448,zz^1043448,0,0,0,zz^645153,0,0,zz^16008,0,0,zz^225723,zz^1043448,0,zz^1043448,0,zz^204588,0,zz^645153,zz^435438],
[0,0,0,zz^101294,0,0,0,zz^893084,zz^473654,0,0,0,zz^704504,0,0,zz^704504,0,0,zz^494789,zz^263939,0,zz^683369,0,zz^263939,0,zz^75359,zz^285074],
[0,zz^983463,0,0,0,zz^372075,zz^40590,0,0,0,zz^476088,zz^43128,0,zz^956403,zz^110778,0,0,zz^230010,0,0,zz^290895,0,zz^324720,0,0,0,0],
[0,zz^536973,0,0,0,zz^669735,zz^744150,0,0,0,zz^672273,zz^983463,0,zz^949638,zz^476088,0,0,zz^1001220,0,0,zz^859155,0,zz^81180,0,0,0,0],
[0,0,0,zz^5127,zz^155491,0,0,zz^44075,zz^587735,0,0,0,zz^603950,0,0,zz^533840,0,0,zz^201740,zz^386015,0,zz^76055,0,zz^513935,0,zz^884390,zz^118715],
[0,0,0,zz^5127,zz^574921,0,0,zz^1007165,zz^723650,0,0,0,zz^533840,0,0,zz^411455,0,0,zz^328430,zz^463505,0,zz^805445,0,zz^705200,0,zz^1023380,zz^45530],
[zz^676500,0,zz^378840,0,0,0,0,0,0,zz^608850,0,0,0,0,0,0,zz^514140,0,0,0,0,0,0,0,zz^547965,0,0],
[0,zz^372075,0,0,0,zz^538662,zz^105702,0,0,0,zz^710325,zz^649440,0,zz^40590,zz^953865,0,0,zz^802497,0,0,zz^389832,0,zz^626607,0,0,0,0],
[0,zz^669735,0,0,0,zz^734847,zz^1046037,0,0,0,zz^230010,zz^372075,0,zz^744150,zz^710325,0,0,zz^119232,0,0,zz^383067,0,zz^180117,0,0,0,0],
[0,0,0,zz^403422,zz^344071,0,0,zz^142250,zz^701285,0,0,0,zz^176300,0,0,zz^253790,0,0,zz^1007165,zz^422690,0,zz^915305,0,zz^369185,0,zz^914915,zz^933365],
[0,zz^500610,0,0,0,zz^599547,zz^1012212,0,0,0,zz^669735,zz^744150,0,zz^439725,zz^372075,0,0,zz^626607,0,0,zz^315417,0,zz^958092,0,0,0,0],
[0,zz^230010,0,0,0,zz^802497,zz^538662,0,0,0,zz^534435,zz^710325,0,zz^372075,zz^879450,0,0,zz^389832,0,0,zz^416892,0,zz^315417,0,0,0,0],
[0,0,0,zz^1032567,zz^344071,0,0,zz^701285,zz^159470,0,0,0,zz^882935,0,0,zz^797450,0,0,zz^723650,zz^142250,0,zz^422690,0,zz^76445,0,zz^176300,zz^705200],
[zz^128535,0,zz^608850,0,0,0,0,0,0,zz^94710,0,0,0,0,0,0,zz^676500,0,0,0,0,0,0,0,zz^169125,0,0],
[0,zz^43128,0,0,0,zz^649440,zz^791505,0,0,0,zz^110778,zz^895518,0,zz^354318,zz^117543,0,0,zz^710325,0,0,zz^534435,0,zz^669735,0,0,0,0],
[0,0,0,zz^822852,zz^553786,0,0,zz^369185,zz^76445,0,0,0,zz^587735,0,0,zz^723650,0,0,zz^705200,zz^911000,0,zz^351965,0,zz^842120,0,zz^44075,zz^176300],
[0,0,0,zz^5127,zz^784636,0,0,zz^1015160,zz^673220,0,0,0,zz^464960,0,0,zz^603950,0,0,zz^324125,zz^705200,0,zz^933365,0,zz^378020,0,zz^957575,zz^1040600],
[0,zz^895518,0,0,0,zz^81180,zz^20295,0,0,0,zz^117543,zz^530208,0,zz^462558,zz^564033,0,0,zz^953865,0,0,zz^879450,0,zz^372075,0,0,0,0],
[0,0,0,zz^5127,zz^365206,0,0,zz^285770,zz^595730,0,0,0,zz^538145,0,0,zz^464960,0,0,zz^394235,zz^513935,0,zz^797450,0,zz^463505,0,zz^830885,zz^114410],
[0,zz^320493,0,0,0,zz^115005,zz^290895,0,0,0,zz^354318,zz^327258,0,zz^56658,zz^462558,0,0,zz^40590,0,0,zz^372075,0,zz^439725,0,0,0,0],
[0,0,0,zz^843987,zz^784636,0,0,zz^933365,zz^705200,0,0,0,zz^201740,0,0,zz^328430,0,0,zz^45530,zz^168305,0,zz^673220,0,zz^176300,0,zz^743555,zz^1023380],
[zz^399135,0,zz^676500,0,0,0,0,0,0,zz^128535,0,0,0,0,0,0,zz^169125,0,0,0,0,0,0,0,zz^304425,0,0],
[0,0,0,zz^403422,zz^973216,0,0,zz^3260,zz^771395,0,0,0,zz^495485,0,0,zz^176300,0,0,zz^44075,zz^495875,0,zz^788615,0,zz^491570,0,zz^723650,zz^797450],
[0,0,0,zz^613137,zz^763501,0,0,zz^286160,zz^842120,0,0,0,zz^513935,0,0,zz^705200,0,0,zz^176300,zz^578900,0,zz^72140,0,zz^771395,0,zz^797450,zz^44075]];

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

// 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[4,4]:=a; scal[5,5]:=b;
a1:=[2,12,11,15,14];
a2:=[26,27,19,16,13];
a3:=[3,1,10,25,17];
a4:=[21,23,7,6,18];
a5:=[22,24,9,8,20];
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[12,27]:=zz^838860*m2;
scal[14,13]:=zz^629145*m2;
scal[15,16]:=zz^838860*m2;

scal[16,15]:=zz^838860*m3;
scal[19,11]:=zz^629145*m3;
scal[26,2]:=zz^629145*m3;
scal[27,12]:=zz^838860*m3;

scal[3,21]:=zz^838860*n2;
scal[10,7]:=zz^629145*n2;
scal[25,6]:=zz^209715*n2;

scal[3,22]:=zz^419430*n3;
scal[17,20]:=zz^209715*n3;

scal[7,10]:=zz^629145*n4;
scal[18,17]:=zz^209715*n4;
scal[21,3]:=zz^419430*n4;
scal[23,1]:=zz^209715*n4;

scal[7,9]:= zz^629145*n6;
scal[18,20]:= zz^419430*n6;
scal[21,22]:= zz^838860*n6;
scal[23,24]:= zz^209715*n6;


scal[20,17]:=zz^838860*n7;
scal[22,3]:=zz^629145*n7;
scal[9,7]:=zz^419430*n8;
scal[20,18]:=zz^629145*n8;
scal[22,21]:=zz^209715*n8;
scal[24,23]:=zz^838860*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.

// Although we were able to do this in characteristic 0, it seems quite hard to work this out by hand for
// p=2. We will instead use Groebner bases to do the determination. This of course means we are dependent
// on the algorithms in Magma, but Groebner bases are a well-known concept, and it shouldn't be too bad.

// We are able to do the first few steps by hand, anyway, but since we rely on Groebner bases eventually,
// there seems little point.

// We first need all of the relations.

rels:=[ProdRel(i):i in seqs];

// We want to know which variables we will be able to set to be 1, using C_GL(H). We have to check that they
// are non-zero first, of course.

homs:=AHom(GModule(H),GModule(H));
scal[1,1] eq n1 and (homs.1)[1,1] eq 1;
scal[6,6] eq n5 and (homs.2)[6,6] eq 1;
scal[4,4] eq a and (homs.3)[4,4] eq 1;

// So each of n1,n5 and a may be set to be 1 if they are non-zero. Since a is always non-zero, we can immediately
// set it to be 1.

coeffs1:=[1,b,m1,m2,m3,m4,n1,n2,n3,n4,n5,n6,n7,n8,n9];
rels1:=[Evaluate(i,coeffs1):i in rels];
B:=GroebnerBasis(rels1);

// We first want to prove that m4 is non-zero. 

badcoeffs2:=[1,b,m1,m2,m3,m4,0,n2,n3,n4,n5,n6,n7,n8,n9];

badB2:=[Evaluate(i,badcoeffs2):i in B];

badB2[75] eq m3*n2*n3;
badB2[124]/zz^489335 eq n2*n3*n7;

// Note that not both of n2 and n3 can be zero. Thus one of the following holds:
// n2=0, n3 ne 0
// n2 ne 0, n3=0.
// n2,n3 ne 0, m3=n7=0.

badcoeffs3a:=[1,b,m1,m2,m3,m4,0,0,n3,n4,n5,n6,n7,n8,n9];
badcoeffs3b:=[1,b,m1,m2,m3,m4,0,n2,0,n4,n5,n6,n7,n8,n9];
badcoeffs3c:=[1,b,m1,m2,0,m4,0,n2,n3,n4,n5,n6,0,n8,n9];
badB3a:=[Evaluate(i,badcoeffs3a):i in B];
badB3b:=[Evaluate(i,badcoeffs3b):i in B];
badB3c:=[Evaluate(i,badcoeffs3c):i in B];

badB3c[103] eq m4*n3; // But n3 ne 0, and not both of m3 and m4 can be 0. Contradiction.

badB3a[7]/zz^5127 eq n3*n7*n8 and badB3a[39]/zz^590527 eq n3*n4*n5 and badB3a[146]/zz^629145 eq n3^2*n4*n8^2;
// Since n3 ne 0 always: one of n7,n8 is zero, one of n4,n5 is zero, one of n4,n8 is zero. Since n1,n2=0, we
// cannot have n4=n7=0 or n5=n8=0, so must have n4=n8=0, n5,n7 ne 0. But...

badB3a[50]/zz^714554 eq n3*(n4*n8+n5*n7)^2;
// Thus branch 3a fails. Finally, branch 3b. This is easier:

badB3b[28]/zz^62574 eq n2*n4*n6 and badB3b[145] eq n2^2*n6^2*n7 and badB3b[139] eq n2^3*n7*n9^3 and badB3b[143] eq n2^2*n4*n9^2;
// Since n2 ne 0, have n4=0 or n6=0, n7=0 or n9=0, so since n4 and n7 cannot both be 0, must be n4=n9=0 or n6=n7=0.
// The other two equations kill off these possibilities.

// Thus n1 ne 0.

coeffs2:=[1,b,m1,m2,m3,m4,1,n2,n3,n4,n5,n6,n7,n8,n9];
rels2:=[Evaluate(i,coeffs2):i in rels];
B2:=GroebnerBasis(rels2);

B2[42]/zz^489335 eq m4-zz^559240*m3*n2*n3;

coeffs3:=ChangeCoefficient(coeffs2,m4,zz^559240*m3*n2*n3);
rels3:=[Evaluate(i,coeffs3):i in rels];
B3:=GroebnerBasis(rels3);

// Now we check that n5 is non-zero.

badcoeffs3:=ChangeCoefficient(coeffs3,n5,0);
brels3:=[Evaluate(i,badcoeffs3):i in rels];
badB3:=GroebnerBasis(brels3);
badB3[20] eq m2^2;
// Since m2=0, m4 and m1 are non-zero. Thus n2,n3 are also non-zero (as m4 ne 0).
badB3[46] eq n2*n6^2; // Thus n6=0. Thus n4 ne 0 as n5=n6=0.

badcoeffs4:=ChangeCoefficients(badcoeffs3,[m2,n6],[0,0]);
brels4:=[Evaluate(i,badcoeffs4):i in rels];
badB4:=GroebnerBasis(brels4);
badB4[17] eq m3*n2*n3*n4;
// but m3,n2,n3,n4 are all non-zero. A contradiction.

coeffs3:=ChangeCoefficient(coeffs3,n5,1);
rels3:=[Evaluate(i,coeffs3):i in rels];
B3:=GroebnerBasis(rels3);

B3[53]/zz^699050 eq (n9+zz^699050*n3*n7 + zz^349525*n6*n8)^2;

coeffs4:=ChangeCoefficient(coeffs3,n9,zz^699050*n3*n7 + zz^349525*n6*n8);
rels4:=[Evaluate(i,coeffs4):i in rels];
B4:=GroebnerBasis(rels4);

B4[2] eq b*m3*n3 + zz^494462*n4^2*n6*n8^2 + zz^494462*n6*n7^2;
// Suppose that n6=0. Then from this, we see that b*m3*n3=0. Since b ne 0, and m3 ne 0 (since then m4=0),
// we have n3=0. But then n9=0, a contradiction. Thus n6 ne 0.

B4[21] eq n6*(m3 + zz^677915*n4*n8 + zz^677915*n7)^2;
// Thus m3=zz^677915*n4*n8 + zz^677915*n7.

coeffs5:=ChangeCoefficient(coeffs4,m3,zz^677915*n4*n8 + zz^677915*n7);
rels5:=[Evaluate(i,coeffs5):i in rels];
B5:=GroebnerBasis(rels5);

B5[7] eq n6*(n7+n4*n8)*(b+zz^772382*n6);
B5[37] eq n6^2*(n7+n4*n8)*(n4*n8 + zz^629145*n7);

// Thus either n7=n4*n8 or b+zz^772382*n6 and n4*n8 + zz^629145*n7. But if n7=n4*n8 then m3=m4=0. Thus:

coeffs6:=ChangeCoefficients(coeffs5,[b,n7],[zz^772382*n6,n4*n8/zz^629145]);
rels6:=[Evaluate(i,coeffs6):i in rels];
B6:=GroebnerBasis(rels6);

// Note that if n8=0 then n7=n9=0, so n8 ne 0. 

B6[8] eq n8*(m2 + zz^31287*n3);
B6[17] eq n8*(n3*n4 + zz^908765*n6);

coeffs7:=ChangeCoefficients(coeffs6,[m2,n6],[zz^31287*n3,n3*n4/zz^908765]);
rels7:=[Evaluate(i,coeffs7):i in rels];
B7:=GroebnerBasis(rels7);

B7[1] eq n3*(m1+zz^520622*n4);
B7[4] eq n3*(n3+zz^1001505*n4*n8);
// Since b ne 0, n3 ne 0. Thus m1=zz^520622*n4 and n3=zz^1001505*n4*n8.

coeffs8:=ChangeCoefficients(coeffs7,[m1,n3],[zz^520622*n4,zz^1001505*n4*n8]);
rels8:=[Evaluate(i,coeffs8):i in rels];
B8:=GroebnerBasis(rels8);

B8 eq [ n4^2*n8*(n2*n4 + zz^279620),n4*n8^2*(n2*n4 + zz^279620)];

// Thus n2,n8 are arbitrary (non-zero), and n4=zz^279620/n2.

// Give one specialization for checking inclusion into E6.

coeffs9:=ChangeCoefficients(coeffs8,[n2,n8,n4],[1,1,zz^279620]);
sca1:=GL(27,F)!Evaluate(scal,coeffs9);
h11:=h1^sca1;
{CheckRel(i,h11):i in seqs} eq {0};

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

// scasub is now the matrices sca8 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 m1*m4*n1+1;
// Thus m1,m4,n1 are all non-zero.

ChRel([2,3,11]) eq m1*m2*n1; // so m2=0.
ChRel([10,13,19]) eq m3*m4*n1; // so m3=0.
ChRel([1,5,13]) eq b*m4*n2; // so n2=0.
ChRel([1,4,13]) eq a*m4*n3; // so n3=0.
ChRel([2,6,16]) eq m1*m4*n4; // so n4=0.
ChRel([4,6,27]) eq a*m4*n6; // so n6=0.
ChRel([2,8,16]) eq m1*m4*n7; // so n7=0.
ChRel([5,8,27]) eq b*m4*n8; // so n8=0.


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

Evaluate(ChRel([2,12,27]),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([5,6,27]),co2) eq b*m4*n5-1;
Evaluate(ChRel([1,21,24]),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, m2 and m3, 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:=[Evaluate(i,[0,0,0,0,0,0,0,m1,0,m2,0,0,0,m3,0]):i in coeffs8];
ctest2:=[Evaluate(i,[0,0,n1,0,0,0,0,0,0,0,n2,0,0,0,n3]):i in co3];
// So recall that n1*n2*n3 eq 1 and m1*m2=zz^279620.

// 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=n1/n2, m3=n1^3 is the solution

{Evaluate(i,[0,0,n1/n2,zz^279620*n2/n1,n1^3,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,n1^3) 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 e be a cube root of y. Note that

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

// so as long as y has a cube root in F_q, (x,y) lies 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 55^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|[i^3:i in Mats]>;
E:=sub<NGT|MatsE>;

l1 in NGT and l2 in NGT;

repeat y1:=Random(NGT); y2:=Random(NGT); until NGT eq sub<NGT|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.";