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

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

l2:=GL(27,F)!DiagonalMatrix([z^i:i in [17,7,2,10,9,5,9,4,8,12,16,7,11,11,16,6,1,14,15,13,4,3,18,17,6,1,5]]);

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

h1vals:=[[1831746,61284,794925,1103245,936154,1857423,1616519,823787,477639,1713113,938273,226435,1166231,53596,1442696,1013964,1453291,1915017,444681,415975,1951156,325027,1316747,1278839,122133,1045720,658334],
[1911612,871798,633755,222231,1339900,86873,835477,778065,1000449,1318447,1350495,1697553,537161,833358,1919486,1728950,19337,374843,835659,692129,555538,237389,1289741,1100201,1176043,1887730,1903924],
[1617293,1558919,1640658,722526,1925613,1272542,460770,1755330,278140,699182,1635676,1070638,811726,1300007,1147395,1238343,1650834,1751762,1851174,1243016,290487,991332,1358800,152450,1726624,185991,784261],
[1617225,839007,414138,390794,273969,441674,762250,1342446,355708,1659506,470380,350726,1797186,332343,1935223,991619,964154,1911190,73850,1922876,1830727,566224,466812,152382,1479900,1452435,1906517],
[1244542,1751084,1411633,68377,1346234,1407083,534823,1761679,126751,1609977,155613,270099,1587227,1163572,306420,361984,212687,1700925,648573,786027,1354108,1729631,1246843,1924907,498057,1943284,345790],
[1343443,1628813,1889318,1366838,584715,421698,1496402,125458,1291048,1428502,440484,66210,1902898,1915873,1741421,1286369,1020266,1275890,1077546,1395544,125977,1745788,710852,835610,12928,689017,1400855],
[1924907,1246661,1899914,556658,534823,365646,1842634,772142,615032,145134,835296,1244066,286290,610665,155613,498057,1657978,236082,1136854,1274308,785117,264788,1735124,947826,1865140,1189249,1407083],
[104215,161289,213390,108894,733719,1872990,1697306,1813742,93736,1683726,314248,630018,216106,218701,1355685,483425,1211454,871954,1222910,505260,1423221,184684,1540196,783898,1791236,1535829,896947],
[580435,1206041,1511692,1897648,1874283,1907824,409440,916104,1500006,979516,351066,717760,1007052,1495333,1815909,64689,59196,639228,271760,1032190,1404385,1932734,1131602,1068716,552970,547477,1419543],
[685153,393283,801978,117566,273629,914522,761910,1375338,1801884,846,777068,1858126,563566,1051847,288787,743527,852858,1357498,1854558,766892,1863619,1461596,1345812,1173434,1231808,1341139,426241],
[732681,1247699,607716,1703932,1594757,748872,321316,828228,42678,1599436,1143998,1234724,352092,20843,1651831,1410927,1351160,101052,698664,1385934,893103,244582,1736890,1706648,433846,960639,96685],
[123639,1697553,145474,1687074,1812039,477394,832882,1246794,512168,830166,1337520,1400674,536642,537161,373933,199481,1907586,1839686,347378,203848,778065,1702232,801460,1431450,851668,19337,1063435],
[1885803,1359529,708930,49654,45287,1183326,697474,1655250,1623828,357974,1277256,1359010,776830,841705,1922569,1600433,269918,1548038,556318,193184,1657845,1944404,675972,1299762,1092600,566797,219739],
[773168,1655726,1197211,537935,1574756,1196301,1021849,1657845,158985,846255,946007,1359529,841705,756974,1757418,1377906,566797,83195,1044599,681465,1185706,479561,1164253,909241,623871,1694166,1765292],
[1237104,1816690,119435,1215651,1745564,96685,1594757,1869665,1507521,1111155,1651831,271137,997405,832254,452742,730562,1937201,1565895,210383,897653,1784934,1709425,1248609,732681,1410927,824566,649592],
[705576,1523358,107587,169251,1698332,1491961,1834405,894609,1606629,1463099,1308131,1947013,572473,349946,627766,1713894,629885,1112855,428163,31797,729458,794857,993201,1766869,1144903,157746,1642768],
[116943,738909,1445242,1066950,521075,197898,13242,594678,573176,544470,220404,674034,167122,464001,806445,1555049,1256214,865046,1899982,1430084,919053,1582696,1231740,1080530,575892,1256733,843211],
[1298241,1813987,312618,780434,775761,1173094,1264042,974750,1872780,1768682,1642992,1325706,211690,699971,1154711,804467,1584618,528750,1236506,1388592,1463031,1178076,236600,1786522,1292748,119775,684813],
[33497,527271,617622,1101810,1882125,1180342,417282,1531298,1710904,518210,493072,38990,1378686,1866967,4791,325367,872022,1442098,21150,442192,66455,9464,1383668,521778,813648,1360303,692061],
[1546731,1925681,1551404,586528,1608395,1087156,143552,402464,107026,972484,769158,1437400,604368,1092649,280877,1470941,1944064,1183000,31008,690626,890745,1030914,323158,81888,6098,479221,598875],
[1231584,1891886,701671,597175,326148,1873509,1710281,1423221,582017,218883,379123,161289,218701,1699686,1270954,318274,1535829,1360235,1711191,993541,310586,672965,75353,1080777,1459987,1112174,1193144],
[1352987,1368145,1196924,1080204,496079,1334604,984360,1932216,904774,1564392,1478134,879864,299668,787949,989853,178081,40756,869688,1448608,928118,467373,4230,1460294,1841268,666362,529037,846323],
[1522339,1598129,742024,158424,1144047,1430424,1632328,512236,1234398,626240,194950,1109848,161992,650273,1659793,1507181,820556,1058968,47320,1351118,1000517,637926,105750,57496,42338,1308837,942143],
[1278839,1202997,1283206,1591526,1616519,1349590,639438,1503470,965920,248270,1912240,1534246,580190,189669,938273,122133,463754,450174,932962,904256,1800349,813308,1805028,1931026,1033820,476729,1857423],
[1766869,970451,595868,657532,1834405,218520,1248364,249296,141786,1951380,331050,646076,64640,1549035,1308131,1144903,1603852,1601136,916444,520078,1871171,1283138,1481482,725432,155366,1606447,1491961],
[1662496,654178,1933523,1555231,298548,1819773,1497637,919053,1061457,1032751,1783007,738909,464001,1591370,1646934,1082910,1256733,1353327,435139,1918365,495398,117853,1720021,1093505,578487,1552930,678060],
[144354,1492740,1401037,878557,1476546,1400855,584715,1102539,802767,940221,1741421,652251,939311,531740,341204,1437176,1665579,787609,589265,907263,1398736,1257507,222571,1343443,1286369,1500428,1905278]];

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

// 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:=[6,8,12,24,25,17,11,7,13];
a2:=[27,21,2,1,16,26,15,5,14];
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[27,6]:=-m3;
scal[17,26]:=-m2;
scal[25,16]:=m2;
scal[24,1]:=m2;
scal[12,2]:=-m2;
scal[8,21]:=-m2;
scal[6,27]:=m2;
scal[13,14]:=-m2;
scal[5,7]:=-m3;
scal[15,11]:=-m3;
scal[26,17]:=m3;
scal[16,25]:=-m3;
scal[1,24]:=-m3;
scal[2,12]:=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,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,4,1,1,4,1,4,4,1,4,1,4,1,4,1,1,4,4,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;

coeffs0:=[a,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 check that there is a unique H over L.

function CheckDetermination()

ProdRel([3,4,22]) eq zz^1619037*a*(a+m1+2*m3)*(a-m1-2*m3);

// Kill off bad branch first.
badc1:=ChangeCoefficient(coeffs0,m3,2*a+2*m1);
Evaluate(ProdRel([3,3,4]),badc1) eq zz^436883*a^2*(m2+3*m4);
badc2:=ChangeCoefficient(badc1,m4,3*m2);
Evaluate(ProdRel([1,9,20]),badc2) eq -a^3;

c1:=ChangeCoefficient(coeffs0,m3,3*a+2*m1);
Evaluate(ProdRel([3,18,22]),c1) eq 2*a*(a+3*m2+4*m4)^2;
c2:=ChangeCoefficient(c1,m4,a+3*m2);
Evaluate(ProdRel([1,3,10]),c2) eq zz^1002261*a^2*(a-m1+2*m2);
c3:=ChangeCoefficient(c2,m2,2*a+3*m1);

// This is the centralizer in GL(27,F) of h1:
sca:=Evaluate(scal,c3);
sca*mats!h1 eq mats!h1*sca;

// Thus H=<l1,l2,h1> is unique containing L.

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

// Now check that H lies in E6.

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

// This proves that the copy of PSL(2,19) we constructed is 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.";
