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

l1:=GL(27,F)![[2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,2,0,0,0,0,0,0,0,0,0,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,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,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,2,0,0,0,0,0,0,0,0,0,0,0,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,2,0,0,0,0,0,0,0],
[0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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,2,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,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,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,2],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,0,0,0,0,0,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,2,0,0,0,0,0,0,0,0,0,0],
[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([z^i:i in [0,3,2,12,9,4,8,1,5,11,15,10,7,14,13,6,0,16,9,13,2,16,1,15,8,0,4]]);

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

h1logs:=[[26904200,0,26582934,0,32081703,31921515,23897923,6673595,0,0,23302161,0,0,0,8655544,0,0,42224248,11262672,42300242,38219058,40424828,5771170,29808284,27409837,0,23661388],
[0,11853340,33699719,25538019,1078781,17894215,36986547,22737539,6239715,23034229,11259473,38815014,6562000,21032811,9079298,8787429,10876855,10154917,2638513,3236343,17290028,23659592,36153595,27053049,5529792,14739599,34347840],
[16875381,24034709,2502295,4030777,22520655,25294195,18247451,14934686,10273306,28889655,20477139,22737539,42203262,30954885,27183365,34244965,30902515,33164348,28765022,41197290,19636379,32288431,32670121,40462443,4278356,28875866,32618084],
[0,22397555,25763895,22998300,36276993,21989261,14067754,36709383,41658650,38971718,2050947,16891265,6054385,39856170,28822392,23878896,663001,10318295,25386679,17773342,13519498,28143914,26828985,18221458,15041690,10694450,4314015],
[1779944,5597143,14455051,10635925,30499015,29529000,34332499,25633955,13509535,33778419,4036979,18852085,15539175,14250106,5008285,17773342,322840,22383050,42424617,1976873,32515154,14006055,27183365,15054052,14679351,31038928,32893409],
[20742734,12092331,25574254,18427727,12414155,33758815,7593435,2502295,683649,25686627,27599772,33699719,36027567,41945373,32515154,12236461,11559652,1695859,35661238,30973051,20033017,15223893,19636379,7011335,17893263,6670032,36754869],
[42272975,13310351,33758815,2862655,35086899,36322366,22882855,25294195,3570805,40379102,33669211,17894215,15303125,29279017,32893409,15692499,31564225,5064810,19318621,17415548,36754869,36990427,32618084,7809465,1561332,25345605,20323047],
[22085144,6877645,12414155,39931053,5272014,35086899,3281000,22520655,16791429,17733879,26401850,1078781,14374227,12328849,42424617,34102687,7571108,12107371,13513579,10014491,35661238,39013354,28765022,2411665,19005747,30112987,19318621],
[0,2336979,32137245,18777370,6366314,4660597,8036193,6152841,13092700,22422576,1287821,33744145,21332385,31740860,25775416,38971718,8930552,20383685,13439887,33778419,7789784,22833287,9163154,39489173,5899119,4202711,6063037],
[0,5512080,19038158,19804585,1726575,4003063,11266285,15946043,2370265,31740860,25657429,16706490,27387420,29616304,34692780,39856170,38240713,18104335,27448721,14250106,17873871,36156532,29811786,26791881,27624137,8424209,32859836],
[315928,575525,1695859,13981791,12107371,5064810,25668259,33164348,22901383,3863463,16757175,10154917,15683261,11590389,8165521,34013161,38448749,12140126,16091400,10107795,31439606,35046424,25230457,10382639,15690214,1096957,3855960],
[0,36040141,31907775,42593942,32450971,31551883,42750318,1970431,16840304,13092700,29127087,20586620,21218250,2370265,3441240,41658650,20219499,15241817,15883745,13509535,42666926,14695992,42375996,41004048,18072281,10461977,16706066],
[0,976381,8587965,31351024,35455598,34840389,1743611,22921063,42593942,18777370,8562209,42892490,31253180,19804585,38345210,22998300,41999632,35043866,39442538,10635925,40781207,13486452,5783598,24797293,25027558,8068402,28529726],
[0,31535371,5230833,31253180,20313645,11613463,33798855,33134597,21218250,21332385,15902439,14298598,8266416,27387420,22949570,6054385,6710689,24108579,23407982,15539175,6753979,37332855,4790031,8790998,24373503,33569612,7201090],
[30806377,24474151,23129976,13841290,34653449,12249028,42413282,12744611,33632085,19210885,8236559,4792237,10644959,3818462,16210815,16813045,30783135,32505182,15471995,23685497,7480680,16757175,955299,39612654,29260570,16449910,7665451],
[0,18661686,12092331,11442745,6877645,13310351,6771215,24034709,22881225,30554090,5911293,11853340,14331770,6371702,8675424,34797820,40266886,575525,8050517,25165897,11224429,37496326,28588142,1123074,39215046,16375482,32369993],
[0,31184628,13686267,29486791,33832676,31215272,15164731,19326265,9266237,19917028,8376136,26155822,27746866,42813912,9721023,22119697,29056536,41684136,32071776,29059897,9732470,18767895,22539895,16449910,41041110,9685512,27439583],
[15595396,4791409,2536570,1743611,40198812,15780731,20598657,15262731,42750318,8036193,4876715,5179725,33798855,11266285,14679351,14067754,36517957,21674415,19005747,34332499,17893263,19727316,4278356,35166784,7289661,8895852,1561332],
[28441875,16737907,6674484,31135202,33740409,10653264,29224485,19838696,1501158,9974029,23685497,32675761,6431132,31910576,4401406,4424143,35308550,6889229,30499015,7173282,14455051,10107795,25633955,16210815,22872899,9721023,29529000],
[40482547,5911293,27599772,10712415,26401850,33669211,29260570,20477139,40528605,35492375,39612654,11259473,42157514,2678731,21280857,9629885,30851306,16757175,40483043,16210815,42465138,10382639,13013996,27155997,28984584,15486193,20488393],
[6231727,31693100,13477740,23893361,38975830,11973176,29893883,39747086,11046869,17780007,8234189,9835095,11032190,17075238,25633955,42507724,32261344,27100987,22520655,1558748,2502295,24279400,14934686,955299,15262731,22539895,25294195],
[16988911,5597595,31038938,16041267,35127363,41838847,6244222,29893883,20571754,10314785,37788766,16838891,21536647,24157993,34332499,3463771,26963687,24664559,3281000,1503103,7593435,21674415,18247451,4876715,20598657,10554976,22882855],
[39383063,10603464,13661768,35710889,1970054,20751582,35127363,38975830,42903493,952854,2932487,24463475,16046540,14776799,30499015,23297408,29777927,19510685,5272014,29726760,12414155,8597691,22520655,15471995,40198812,32071776,35086899],
[16063965,18299502,42232545,16241295,13743853,39341531,30144612,20544783,6719517,30529129,32505182,22981776,21305368,25576314,10107795,22161211,5955744,1219851,8597691,6889229,25942339,12140126,24279400,16757175,25668259,18767895,843715],
[14298370,4014497,14362682,37531842,27256460,33072582,21003750,43028257,6507306,15959590,12249028,14567839,8379027,40154324,17415548,27271965,8901344,39341531,8224419,10653264,12886419,5064810,22829130,33669211,22882855,13564295,13201033],
[0,42515942,15971616,23842613,9709112,3424149,14504031,30647497,31457474,2010995,11208315,17785659,29098050,31855267,30365664,41997240,9685512,7751453,30672175,13433611,31722647,18909030,1498790,1395767,26963687,3228504,8901344],
[3721668,33350049,20691825,18510748,13661768,42098615,31038938,13477740,38282572,18169899,2130697,13453551,14680159,16117685,14455051,22979751,35985215,5616887,12414155,29992530,25574254,25942339,2502295,7480680,2536570,9732470,33758815]];

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

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

coeffs5:=[a,1,n1,0,-a*n1-1,zz^37376627*a*n1^2+zz^15853267*n1,zz^27193453*a*n1+zz^5670093,-a*n1^2-n1];
sca:=Evaluate(scal,coeffs5)*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 overgroups H of L. We first determine the possible matrices that conjugate H into E6,
// then that the possible matrices factorize as C_GL(H) C_E6(L). Thus all overgroups are C_E6(L)-conjugate.

function CheckDetermination()
// Note that a,b are non-zero.
ProdRel([2,2,17]) eq zz^25321600*(n1+n2+zz^27193453*m2-m4)*(b+m1+zz^15853267*m3)*(b-m1+zz^37376627*m3);
// This yields three options. As we know that a=b=n1=m1=m4=1, m2=m3=n2=0 is a solution, the second might not work. The first doesn't either.

badc1a:=[a,b,n1,n2,m1,m2,m3,n1+n2+zz^27193453*m2];
badc1b:=[a,b,n1,n2,-(b+zz^15853267*m3),m2,m3,m4];
// First kill off branch b.
Evaluate(ProdRel([2,2,26]),badc1b) eq zz^3798240*b^2*n1;
Evaluate(ProdRel([2,4,16]),badc1b) eq zz^18385427*b^2*n2;
// Since not both of n1 and n2 can be zero without det(scal) being zero, we obtain a contradiction.

// Now kill branch a. This takes two steps:
Evaluate(ProdRel([2,4,13]),badc1a) eq zz^8256787*n2*(m1+zz^15853267*m3)^2;
badc2a:=ChangeCoefficient(badc1a,m3,-m1/zz^15853267);

badc2a eq [a,b,n1,n2,m1,m2,zz^5670093*m1,n1+n2+zz^27193453*m2];

Evaluate(ProdRel([2,2,26]),badc2a) eq zz^3798240*b^2*n2;
Evaluate(ProdRel([2,4,16]),badc2a) eq zz^18385427*b^2*n1;
// Again, we have the same contradiction.

// Thus we know the answer:
c1:=[a,b,n1,n2,b+zz^37376627*m3,m2,m3,m4];
Evaluate(ProdRel([2,4,16]),c1) eq zz^18385427*b^2*n2; // Since b ne 0, n2=0.
c2:=ChangeCoefficient(c1,n2,0);

Evaluate(ProdRel([2,2,26]),c2) eq zz^3798240*b^2*(n1+zz^27193453*m2-m4);
c3:=ChangeCoefficient(c2,m4,n1+zz^27193453*m2);

// At this point, like in characteristic 0, it seems easier to set b=1 and then put it back in afterwards.
Evaluate(ProdRel([1,2,26]),c3) eq zz^34184160*b*(b*m2-zz^10183174*n1*m3);
c4:=ChangeCoefficients(c3,[b,m2],[1,zz^10183174*n1*m3]);

Evaluate(ProdRel([1,2,12]),c4) eq zz^1980921*n1*(zz^27193453*a*n1+zz^5670093-m3);
c5:=ChangeCoefficient(c4,m3,zz^27193453*a*n1+zz^5670093);

// This kills off all coefficients

// Now we check that the centralizer 

// First find the centralizer in GL(27,k) of H.

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

scaents[170] eq zz^6085923*b*(n1-1);
{Evaluate(i,[a,b,1,n2,m1,m2,m3,m4]):i in scaents} eq {0};

// So we see that an element of sca centralizes H if and only if n1=1.

// Now find which elements of sca lie in E6.

ChRel([2,14,26]) eq zz^21523360*(b^3*n1-1);
ChRel([2,6,12]) eq zz^37376627*n1*b^3*(a*n1-1);

// This shows that a=b^3 and that n1=b^-3. Since the characteristic is 3, if n1=1 then b=a=1, so
// the centralizer of L in GL(27,F) (sca) modulo the centralizer of H in GL(27,F) (n1=1) is the
// centralizer of L in E6 (a=b^3, n1=b^-3). Thus all copies of H in E6 containing L are
// conjugate by an element of C_E6(L).

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