/* 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:=3^20;
F<zz>:=GF(q);
z:=zz^63396080;

l1a:=DiagonalMatrix([zz^i:i in [1743392200,3138105960,1046035320,1046035320,1394713760,2789427520,348678440,1394713760,697356880,1394713760,2440749080,1046035320,2440749080,
2092070640,2440749080,2092070640,0,2092070640,0,697356880,0,2092070640,1394713760,0,0,1743392200,2789427520]]);

l1b:=PermutationMatrix(F,Sym(27)!(1,8,13,20,9)(2,10,22,11,5)(3,6,26,27,14)(4,15,21,25,17)(7,12,19,18,24));

l1:=l1a*l1b;
l2:=GL(27,F)!DiagonalMatrix([z^i:i in [35,45,5,35,15,15,50,50,30,25,5,40,40,20,50,0,30,30,10,10,40,20,0,35,10,45,25]]);

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

h1logs:=[[2607884200,3405014575,0,33169241,1753541606,0,0,1339962600,1834572610,1204715237,1308507490,0,589681593,0,376424486,0,1609152154,0,0,2408238070,70644669,3470390995,0,0,3324499687,0,0],
[0,677185400,0,1246739446,718478870,0,0,0,0,589681593,2331802953,0,0,0,2072501048,0,2299257295,0,0,0,1586020458,1974265411,0,0,49040640,0,0],
[0,0,1673536050,0,0,3027229329,1277214151,0,0,0,0,2331929001,0,1370837217,0,2139892113,0,829996781,1833460864,0,0,0,809969567,3384571708,0,872063815,1877181753],
[0,0,0,2607884200,0,0,0,0,0,0,0,0,0,0,3010200291,0,2550251400,0,0,0,850083800,0,0,0,2031556200,0,0],
[0,1752999830,0,3470854366,2550251400,0,0,0,0,3010200291,1974265411,0,0,0,2375803685,0,986433021,0,0,0,1554937615,334466699,0,0,2437363266,0,0],
[0,0,1375597721,0,0,3120024850,434958444,0,0,0,0,1029631151,0,901481015,0,2508492854,0,208300314,656758559,0,0,0,272344113,1325018332,0,506580132,1130667649],
[0,0,1312815924,0,0,355602791,1673536050,0,0,0,0,1030597015,0,1675630259,0,3282472409,0,829256897,2947498177,0,0,0,3233616233,1114686850,0,2425328491,2572110855],
[1339962600,1588025795,0,53324613,1472668842,0,0,850083800,196560531,491706327,2704937898,0,2521610230,0,2015722469,0,1708465104,0,0,1003400097,2794351608,3222065471,0,0,844600448,0,0],
[2260308690,3083909205,0,2365013936,3222065471,0,0,2684311819,2031556200,1472668842,935393988,0,1820349937,0,1708465104,0,2388743975,0,0,2155436610,740956069,1588025795,0,0,2492191952,0,0],
[0,1079366657,0,773161834,21840059,0,0,0,0,2031556200,2607884200,0,0,0,2437363266,0,2851894061,0,0,0,1653474357,2979525430,0,0,1020043745,0,0],
[0,1003400097,0,1653474357,1248624539,0,0,0,0,2607884200,850083800,0,0,0,3470854366,0,1383553051,0,0,0,2851894061,2016933430,0,0,773161834,0,0],
[0,0,1986770523,0,0,1544150402,875005535,0,0,0,0,3120024850,0,49060913,0,3329729955,0,1130667649,2103615809,0,0,0,436234478,2301804335,0,156252377,298082046],
[1079366657,935393988,0,2666701417,1685159079,0,0,2408238070,2740731113,1715478659,3405014575,0,2550251400,0,284582427,0,2912952700,0,0,2684311819,1619067695,2704937898,0,0,789788022,0,0],
[0,0,866980633,0,0,2301804335,2559004971,0,0,0,0,1325018332,0,3058612850,0,2859685436,0,1832298221,2149085151,0,0,0,809969567,571698776,0,1375597721,1072436788],
[0,0,0,21840059,0,0,0,0,0,0,0,0,0,0,850083800,0,2684311819,0,0,0,1248624539,0,0,0,1769044779,0,0],
[0,0,991919083,0,0,2214001997,2419745467,0,0,0,0,3241487080,0,2370965084,0,0,0,1839080064,2460368727,0,0,0,1121374487,947771560,0,1792126702,1789612560],
[0,0,0,2550251400,0,0,0,0,0,0,0,0,0,0,196560531,0,2031556200,0,0,0,677185400,0,0,0,850083800,0,0],
[0,0,2572110855,0,0,2069727576,1877181753,0,0,0,0,3333099001,0,180433316,0,2988744771,0,185948450,1673536050,0,0,0,678171846,489743273,0,2519590533,679253378],
[0,0,2797488934,0,0,1299855854,489743273,0,0,0,0,3202440041,0,2929555264,0,1406939980,0,1673536050,1114686850,0,0,0,296203598,2119992793,0,2626076366,1251946134],
[2521610230,497338139,0,284582427,2066885527,0,0,2331802953,1772215090,923012259,1954905625,0,196560531,0,2005773702,0,2732671912,0,0,677185400,982524169,2948275975,0,0,2572495338,0,0],
[0,0,0,850083800,0,0,0,0,0,0,0,0,0,0,1974265411,0,677185400,0,0,0,2550251400,0,0,0,2607884200,0,0],
[0,1248624539,0,1544667526,777267651,0,0,0,0,1829860870,3294141670,0,0,0,49040640,0,3162350256,0,0,0,1700591769,2607884200,0,0,2375803685,0,0],
[0,0,1540376218,0,0,2575297868,2985130009,0,0,0,0,552335164,0,1540376218,0,403463979,0,186671893,961483066,0,0,0,0,186671893,0,413245726,533009654],
[0,0,1003318334,0,0,49060913,1114686850,0,0,0,0,901481015,0,1794068166,0,1877574019,0,2947498177,920905057,0,0,0,678171846,3058612850,0,2075027319,2797488934],
[0,0,0,2031556200,0,0,0,0,0,0,0,0,0,0,3238099971,0,850083800,0,0,0,2607884200,0,0,0,677185400,0,0],
[0,0,1664192735,0,0,3413631268,3330970080,0,0,0,0,1078126084,0,3027229329,0,1695131618,0,950886397,796798260,0,0,0,2938708187,959646678,0,1114686850,2397607150],
[0,0,829256897,0,0,3333099001,829996781,0,0,0,0,790629072,0,2828406212,0,1189490091,0,1777756310,2011361935,0,0,0,2095652397,1833460864,0,2397607150,185948450]];

h1:=GL(27,F)![[zz^h1logs[i,j]-1:j in [1..27]]:i in [1..27]];

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

scal[6,5]:=zz^1394713760*m3;
scal[27,10]:=zz^1394713760*m3;

scal[15,8]:=zz^697356880*n2;
scal[17,9]:=zz^2789427520*n2;
scal[21,13]:=zz^3138105960*n2;
scal[25,20]:=zz^2092070640*n2;

scal[15,7]:=zz^2440749080*n3;
scal[17,18]:=zz^1394713760*n3;
scal[21,12]:=zz^348678440*n3;
scal[25,19]:=zz^1394713760*n3;

scal[8,15]:=zz^2789427520*n4;
scal[9,17]:=zz^697356880*n4;
scal[13,21]:=zz^348678440*n4;
scal[20,25]:=zz^1394713760*n4;

scal[8,7]:=zz^1743392200*n6;
scal[9,18]:=zz^2092070640*n6;
scal[13,12]:=zz^697356880*n6;
scal[20,19]:=zz^2789427520*n6;

scal[12,21]:=zz^2092070640*n7;
scal[18,17]:=zz^1046035320*n7;
scal[19,25]:=zz^1046035320*n7;
scal[24,4]:=zz^2440749080*n7;

scal[12,13]:=zz^1046035320*n8;
scal[18,9]:=zz^3138105960*n8;
scal[19,20]:=zz^2440749080*n8;
scal[24,1]:=zz^1743392200*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));
//homs.1+homs.3 is the identity.
scal[4,4] eq n1 and (homs.1)[4,4] eq 1;
scal[16,16] eq a and (homs.3)[16,16] eq 1;
// So if n1 ne 0 then we may set n1=1. homs.2 is a 'unipotent' element of the hom-space, so we don't yet know
// what to set that to to keep our equations nice.

badcoeffs1:=[1,b,m1,m2,m3,m4,0,n2,n3,n4,n5,n6,n7,n8,n9];
brels1:=[Evaluate(i,badcoeffs1):i in rels];
bB1:=GroebnerBasis(brels1);

bB1[1] eq b*m1*n2 and bB1[8] eq b*m3*n2;
// Since b ne 0 and not both of m1 and m3 can be zero, need n2=0. Hence n3 ne 0.

bB1[30] eq m1*n3^2 and bB1[65] eq m3*n3^2;
// This yields a contradiction. Thus n1 ne 0.

// So each of n1 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,1,n2,n3,n4,n5,n6,n7,n8,n9];
rels1:=[Evaluate(i,coeffs1):i in rels];
B:=GroebnerBasis(rels1);

B[15] eq n8*(m4-zz^187096113*n9);
B[10] eq n5*(m2-zz^1803074312*n3);
B[12] eq m3-zz^1102327458*n8*n9;
B[34] eq n8*n9*(n7-zz^435798490*n9^2);

// So to conclude that the last terms are zero, we need that n5, n8 and n9 are all non-zero.
// Start with n8. If it is zero then from the third equation, m3=0.
badcoeffs2a:=ChangeCoefficients(coeffs1,[m3,n8],[0,0]);
brels2a:=[Evaluate(i,badcoeffs2a):i in rels];
bB2a:=GroebnerBasis(brels2a);
bB2a[1] eq m1; // contradiction as m3=0.
// Thus if n9=0 then m4=0 and m3=0, a contradiction. Hence n8,n9 are both non-zero. B[34] now shows n7 ne 0.

// Now we need to consider n5, which we will later see can actually be zero. Assume it is, and we can include
// our other conditions that we have verified.
badcoeffs2b:=ChangeCoefficients(coeffs1,[m3,m4,n7,n5],[zz^1102327458*n8*n9,zz^187096113*n9,zz^435798490*n9^2,0]);
brels2b:=[Evaluate(i,badcoeffs2b):i in rels];
bB2b:=GroebnerBasis(brels2b);

bB2b[5]+zz^59682112*bB2b[9] eq n8*(m2-zz^1803074312*n3);
// So the equation holds even if n5=0.

// We also check another one while we are at it. This one will not be needed yet, but we will see it later.
// Note that it shows that (up to constant) n3=n9^-1.
B[19] eq n3*n5*(n3*n9-zz^871745660);
// We also need to check it if n5=0. (We will deal with n3=0 next.)
bB2b[9]*n9^2-zz^1187247800*bB2b[13] eq n8*n9*(n3*n9-zz^871745660);

//If n3=0 then m2=0:
badcoeffs2c:=ChangeCoefficients(coeffs1,[m2,n3],[0,0]);
brels2c:=[Evaluate(i,badcoeffs2c):i in rels];
bB2c:=GroebnerBasis(brels2c);
bB2c[1] eq m1; //a contradiction.

// Thus n3,n7,n8,n9 all non-zero. We also have n3*n9-zz^871745660.

coeffs2:=ChangeCoefficients(coeffs1,[m3,m4,n7,m2],[zz^1102327458*n8*n9,zz^187096113*n9,zz^435798490*n9^2,zz^1803074312*n3]);
rels2:=[Evaluate(i,coeffs2):i in rels];
B2:=GroebnerBasis(rels2);

B2[8] eq m1-(zz^3208196577*n5*n9+zz^854617107*n6*n8);

coeffs3:=ChangeCoefficient(coeffs2,m1,zz^3208196577*n5*n9+zz^854617107*n6*n8);
rels3:=[Evaluate(i,coeffs3):i in rels];
B3:=GroebnerBasis(rels3);

B3[19] eq n8*n9*(n4-(zz^2615038740*n6*n9+zz^1950858160));

// Instead of using this, multiply through by n3 and then use the relation n3*n9=zz^871745660. This removes the constant term.
// One obtains n3*n4+zz^207465960*n3-n6
eq1:=n3*n9-zz^871745660;
eq2:=n4-(zz^2615038740*n6*n9+zz^1950858160);
eq2*n3+zz^2615038740*n6*eq1 eq n3*n4+zz^207465960*n3-n6;

coeffs4:=ChangeCoefficient(coeffs3,n6,n3*n4+zz^207465960*n3);
rels4:=[Evaluate(i,coeffs4):i in rels];
B4:=GroebnerBasis(rels4);
B4[7] eq n3*(n2-zz^2004851470*n3^2*n8);

coeffs5:=ChangeCoefficient(coeffs4,n2,zz^2004851470*n3^2*n8);
coeffs5 eq [Evaluate(i,[0,b,0,0,0,0,0,zz^2004851470*n3^2*n8,n3,n4,n5,0,0,n8,n9]):i in coeffs4];
rels5:=[Evaluate(i,coeffs5):i in rels];
B5:=GroebnerBasis(rels5);
B5[4] eq n5*n9*(b*n5+zz^2365888910*n4+zz^568453840);

// Again, instead of the constant term, this time we replace 1 by n3*n9/zz^871745660, which is the same thing.
// We obtain b*n5+zz^2365888910*n4+zz^568453840*n3*n9/zz^871745660, and rearrange for n4.

// n4=(b*n5+zz^568453840*n3*n9/zz^871745660)/-zz^2365888910;

// Note that need to take care of the case n5=0.
bB2b[10] eq n6*n8*(n4-zz^568453840/-zz^2365888910);
// We haven't yet shown that n6 is non-zero. In general, it can be 0, as we have that element homs.2 to play with.
// But one cannot set n5 and n6 to be zero simultaneously, as n3 is non-zero, and 
bB2b[9] eq n8*(n3+zz^1187247800*n6);
// Thus the equation holds even if n5=0.

coeffs6:=ChangeCoefficient(coeffs5,n4,(b*n5+zz^568453840*n3*n9/zz^871745660)/-zz^2365888910);
rels6:=[Evaluate(i,coeffs6):i in rels];
B6:=GroebnerBasis(rels6);
B6[2] eq n8*(b*n3*n8-zz^1232683980*n9);

coeffs7:=ChangeCoefficient(coeffs6,n9,b*n3*n8/zz^1232683980);
rels7:=[Evaluate(i,coeffs7):i in rels];
B7:=GroebnerBasis(rels7);
B7 eq [n3*n5*(b*n3^2*n8-zz^2104429640),n3*n8*(b*n3^2*n8-zz^2104429640)];

// We still have our parameter from homs.2, which we can use to set n5=0.

coeffs8:=ChangeCoefficient(coeffs7,n5,0);

// remove any copies of b*n3^2*n8 using the equation b*n3^2*n8-zz^2104429640 from B7.

coeffs9:=[1,b,zz^1410761507*n3*n8,zz^1803074312*n3,zz^3356427878*b*n3*n8^2,zz^2441196533*b*n3*n8,1,
zz^2004851470*n3^2*n8,n3,zz^3432741530,0,zz^556144400*n3,zz^74860170*b*n8,n8,zz^2254100420*b*n3*n8];

// We modify the entries 3,10,12,13.
[i:i in [1..#coeffs8]|coeffs8[i] ne coeffs9[i]] eq [ 3, 10, 12, 13 ];

eq3:=b*n3^2*n8-zz^2104429640;
coeffs8[3]-LeadingCoefficient(coeffs8[3])*n3*n8*eq3 eq coeffs9[3];
coeffs8[10]-LeadingCoefficient(coeffs8[10])*eq3 eq coeffs9[10];
coeffs8[12]-LeadingCoefficient(coeffs8[12])*n3*eq3 eq coeffs9[12];
coeffs8[13]-LeadingCoefficient(coeffs8[13])*b*n8*eq3 eq coeffs9[13];
// This shows that it is just a simplification.

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

scal9:=Evaluate(scal,coeffs9);
scasub:=mats!scal10^-1*scal9;

// 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([2,12,27]) eq m1*m4*n9-1;
// Thus m1,m4,n9 are all non-zero.

ChRel([2,2,14]) eq m1*m2*m4; //so m2=0.
ChRel([3,14,18]) eq m3*m4*n9; //so m3=0.
ChRel([4,22,23]) eq -b*m1*n2; //so n2=0.
ChRel([4,6,11]) eq -m1*m4*n3; //so n3=0.
ChRel([1,16,22]) eq -a*m1*n4; //so n4=0.
ChRel([1,6,11]) eq -m1*m4*n6; //so n6=0.
ChRel([7,11,16]) eq a*m1*n7; //so n7=0.
ChRel([7,11,23]) eq b*m1*n8; //so n8=0.

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

Evaluate(ChRel([7,10,24]),co1) eq m1*n9^2-1;
Evaluate(ChRel([2,12,27]),co1) eq m1*m4*n9-1;

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

// Now we need to obtain the final relations.

Evaluate(ChRel([2,12,27]),co2) eq m1*m4^2-1;
Evaluate(ChRel([2,20,23]),co2) eq b*m1*n5-1;
Evaluate(ChRel([1,17,26]),co2) eq m4*n1*n5-1;
Evaluate(ChRel([2,16,25]),co2) eq a*m1*n1-1;
//Then:
//a=m1^-1*n1^-1=m4^2*m4*n5=m4^3*n5
//b=m1^-1*n5^-1=m4^2*m4*n1=m4^3*n1
//m1=m4^-2=(n1*n5)^2.

co3:=[m4^3*n5,m4^3*n1,n1^2*n5^2,0,0,m4,n1,0,0,0,n5,0,0,0,m4];

// The remaining coefficients satisfy m4*n1*n5=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:=[Evaluate(i,[0,m1,0,0,0,0,0,0,m2,0,0,0,0,m3,0]):i in coeffs9];
ctest2:=[Evaluate(i,[0,0,0,0,0,n1,n2,0,0,0,n3,0,0,0,0]):i in co3];

// This is all elements X of C_GL(L) that map h11 into E6, modulo C_GL(H)
scat1:=mats!scal10^-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=1/n1/n2^2, m2=n2/n1, m3=zz^2104429640*n1^3 is the solution

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

// The map (x,y) -> (x^-1*y^-2,y*x^-1) is a homomorphism from F_q^x x F_q^x to itself. The kernel is all maps such that
// x*y^2=1 and x=y, so x^3=1, so x=1 as F has characteristic 3. Thus the map is an injection, so a bijection.

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