/* In this case we let H denote a copy of PSL(2,13), L denote the subgroup 13.6, and let G denote F4(2^12).
We give matrices defining L, an element of order 13 and an element of order 6. We give an element x that
conjugates our version of L into the normalizer of a torus of F4(2^12), as constructed using the GroupOfLieType
command. We give H as a pair of matrices, one of which lies in L and the other chosen at random. A posteriori,
one could produce H as a subgroup of F4 and show that no other overgroup of L preserves the F4 trilinear form,
but we prefer a less synthetic proof.
*/
F<w>:=GF(2^12);
y:=w^(Order(w) div 3);
z:=w^(Order(w) div 13);

l1:=GL(26,F)![[z^10,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,z,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,z^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,z^6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,z^3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,z^7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,z^11,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,z^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,z^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,0,z^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,0,z^12,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,z^9,0,0,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,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,z,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,z^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,z^5,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,z^5,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,z^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,z^9,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,z^6,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,z^10,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,z^7,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,z^11,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,z^12,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,z^3]];

l2:=GL(26,F)![[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],
[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,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,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,1,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,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,1,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,1,0],
[0,0,0,0,0,0,0,0,0,0,0,0,y,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,y^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,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,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,1,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,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,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,1],
[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]];

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

x:=IdentityMatrix(F,26); x[13,14]:=y; x[14,13]:=y;

// We need the indices of the trilinear form as well.
seqs:=[[i,j,k]:i,j,k in [1..NumberOfRows(l1)]|i le j and j le k];

/* We next give the symmetric trilinear form that comes from F4. This was found by choosing elements of F4, mapping
them to lie over our copy of L, and then using the G-invariance of the form to obtain constraints on the structure
constants of it. Code is included at the end to reconstruct this, but it can take a while, so we give it here for
ease of calculations.
*/

ents:=[260,273,295,313,327,332,571,584,597,608,628,644,840,855,883,896,926,942,1081,1098,1158,1171,1181,1209,1295,1314,
1410,1423,1433,1452,1581,1613,1640,1653,1674,1773,1807,1836,1848,1861,2012,2039,2052,2072,2167,2196,2207,2220,2361,2374,
2384,2483,2494,2507,2615,2628];

vals:=[1,y^2,y^2,y^2,y^2,y^2,y^2,y,y,y^2,y^2,y^2,y^2,y^2,y^2,1,y^2,y^2,y^2,y^2,y^2,1,y^2,y^2,y^2,y^2,y^2,1,y^2,y^2,y^2,y^2,y,y,y^2,y^2,y^2,y^2,1,y^2,y^2,y,y,y^2,y^2,y^2,1,y^2,y,y,y^2,y^2,1,y^2,y^2,1];

f26:=[F!0:i in [1..#seqs]];
for i in [1..#ents] do f26[ents[i]]:=vals[i]; end for;


// Now set up the centralizer of L in GL(26,F), which is given by all possible matrices scal below.

aa1:=[15,22,20,11,5,8];
aa2:=[2,1,12,25,26,16];
bb1:=[3,6,17,24,21,9];
bb2:=[19,23,18,7,4,10];

R<a,b,m1,m2,m3,m4,m5,m6,m7,m8>:=PolynomialRing(F,10);
mats:=MatrixRing(R,NumberOfRows(l1));
scal:=mats!1;
for i in [1..6] do
  scal[aa1[i],aa1[i]]:=m1; scal[aa1[i],aa2[i]]:=m2;
  scal[aa2[i],aa1[i]]:=m3; scal[aa2[i],aa2[i]]:=m4;
end for;
for i in [1..6] do
  scal[bb1[i],bb1[i]]:=m5; scal[bb1[i],bb2[i]]:=m6;
  scal[bb2[i],bb1[i]]:=m7; scal[bb2[i],bb2[i]]:=m8;
end for;
scal[13,13]:=a; scal[14,14]:=b;

// 12-dimensional for H:
yy:=w^65;
tt:=GL(12,F)![[yy^16, yy^23, yy^13, yy^2 , yy^58, yy^41, yy^51, yy^10, yy^56, yy^23, yy^10,     1],
[yy^23, yy^13, yy^2 , yy^58, yy^41, yy^16, yy^10, yy^56, yy^23, yy^10,     1, yy^51],
[yy^13, yy^2 , yy^58, yy^41, yy^16, yy^23, yy^56, yy^23, yy^10,     1, yy^51, yy^10],
[yy^2 , yy^58, yy^41, yy^16, yy^23, yy^13, yy^23, yy^10,     1, yy^51, yy^10, yy^56],
[yy^58, yy^41, yy^16, yy^23, yy^13, yy^2 , yy^10,     1, yy^51, yy^10, yy^56, yy^23],
[yy^41, yy^16, yy^23, yy^13, yy^2 , yy^58,     1, yy^51, yy^10, yy^56, yy^23, yy^10],
[yy^13, yy^35, yy^18, yy^48, yy^35, yy^25, yy^41, yy^16, yy^23, yy^13, yy^2 , yy^58],
[yy^35, yy^18, yy^48, yy^35, yy^25, yy^13, yy^16, yy^23, yy^13, yy^2 , yy^58, yy^41],
[yy^18, yy^48, yy^35, yy^25, yy^13, yy^35, yy^23, yy^13, yy^2 , yy^58, yy^41, yy^16],
[yy^48, yy^35, yy^25, yy^13, yy^35, yy^18, yy^13, yy^2 , yy^58, yy^41, yy^16, yy^23],
[yy^35, yy^25, yy^13, yy^35, yy^18, yy^48, yy^2 , yy^58, yy^41, yy^16, yy^23, yy^13],
[yy^25, yy^13, yy^35, yy^18, yy^48, yy^35, yy^58, yy^41, yy^16, yy^23, yy^13, yy^2 ]];

uu:=GL(14,F)![[     0 ,w^3325 ,w^1825 ,w^3190  ,w^460 ,w^1825 ,w^3190  ,w^460 ,w^3125  ,w^395 ,w^1760 ,w^3125  ,w^395 ,w^1760],
[ w^770  ,    0 ,w^2140  ,w^775 ,w^3505 ,w^2140  ,w^775 ,w^3505  ,w^710 ,w^3440 ,w^2075  ,w^710 ,w^3440 ,w^2075],
[w^2270 ,w^1955 ,w^2015 ,w^3055 ,w^3575 ,w^3835 ,w^3965 ,w^4030 ,w^2015  ,w^975  ,w^455  ,w^195   ,w^65    ,  1],
[ w^905 ,w^3320 ,w^3055 ,w^3575 ,w^3835 ,w^3965 ,w^4030 ,w^2015  ,w^975  ,w^455  ,w^195   ,w^65    ,  1 ,w^2015],
[w^3635  ,w^590 ,w^3575 ,w^3835 ,w^3965 ,w^4030 ,w^2015 ,w^3055  ,w^455  ,w^195   ,w^65    ,  1 ,w^2015  ,w^975],
[w^2270 ,w^1955 ,w^3835 ,w^3965 ,w^4030 ,w^2015 ,w^3055 ,w^3575  ,w^195   ,w^65    ,  1 ,w^2015  ,w^975  ,w^455],
[ w^905 ,w^3320 ,w^3965 ,w^4030 ,w^2015 ,w^3055 ,w^3575 ,w^3835   ,w^65     , 1 ,w^2015  ,w^975  ,w^455  ,w^195],
[w^3635  ,w^590 ,w^4030 ,w^2015 ,w^3055 ,w^3575 ,w^3835 ,w^3965     , 1 ,w^2015  ,w^975  ,w^455  ,w^195   ,w^65],
[ w^970 ,w^3385 ,w^2145 ,w^1105  ,w^585  ,w^325  ,w^195  ,w^130 ,w^4030 ,w^2015 ,w^3055 ,w^3575 ,w^3835 ,w^3965],
[w^3700  ,w^655 ,w^1105  ,w^585  ,w^325  ,w^195  ,w^130 ,w^2145 ,w^2015 ,w^3055 ,w^3575 ,w^3835 ,w^3965 ,w^4030],
[w^2335 ,w^2020  ,w^585  ,w^325  ,w^195  ,w^130 ,w^2145 ,w^1105 ,w^3055 ,w^3575 ,w^3835 ,w^3965 ,w^4030 ,w^2015],
[ w^970 ,w^3385  ,w^325  ,w^195  ,w^130 ,w^2145 ,w^1105  ,w^585 ,w^3575 ,w^3835 ,w^3965 ,w^4030 ,w^2015 ,w^3055],
[w^3700  ,w^655  ,w^195  ,w^130 ,w^2145 ,w^1105  ,w^585  ,w^325 ,w^3835 ,w^3965 ,w^4030 ,w^2015 ,w^3055 ,w^3575],
[w^2335 ,w^2020  ,w^130 ,w^2145 ,w^1105  ,w^585  ,w^325  ,w^195 ,w^3965 ,w^4030 ,w^2015 ,w^3055 ,w^3575 ,w^3835]];

h1:=Matrix(F,26,26,[0:i in [1..26^2]]);
ab1:=aa1 cat bb1;
ab2:=[13,14] cat aa2 cat bb2;
for i in [1..12] do for j in [1..12] do h1[ab1[i],ab1[j]]:=tt[i,j]; end for; end for;
for i in [1..14] do for j in [1..14] do h1[ab2[i],ab2[j]]:=uu[i,j]; end for; end for;
h1:=GL(26,F)!h1;

// We also make the .2 on top of the PSL(2,13).

j1:=GL(26,F)![[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,1,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,1,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,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,1,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,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,1,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,w^1365,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,w^2730,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,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,1,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,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,1,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,1,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,0,0,0,0,0,0,1,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]];

// We also need this matrix ring so we can multiply elements of F4 by scal, as we have to coerce them into this ring first

function Matricise(v)
return Matrix(1,NumberOfColumns(v),[R!v[i]:i in [1..NumberOfColumns(v)]]);
end function;

// We now show that g, which is the second generator of H (the first lies in L, so is fine) satisfies the equations mentioned
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]]*f26[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]]*f26[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]]*f26[Position(seqs,Sort([i[1],i[2],i[3]]))]:i in intseqs];
e2:=&+[u2[i[1]]*v2[i[2]]*w2[i[3]]*f26[Position(seqs,Sort([i[1],i[2],i[3]]))]:i in intseqs];
return e1-e2;
end function;

coeffs0:=[a,b,m1,m2,m3,m4,m5,m6,m7,m8];

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;

function CheckDetermination()

aa:=ProdRel([3,11,15]);
bb:=ProdRel([5,8,20]);
aa*w^2535+bb*w^715 eq m2*m6*(m1+w^1235*m5);
//if m2 is 0 then m1 is 0 as we see from aa:
Evaluate(aa,[a,b,m1,0,m3,m4,m5,m6,m7,m8]) eq w^1495*m1^3;
//if m6 is 0 then m5 is 0 as we see from another relation
Evaluate(ProdRel([3,6,21]),[a,b,m1,m2,m3,m4,m5,0,m7,m8]) eq w^3315*m5^3;

//Thus m1=w^1235*m5.
coeffs1:=ChangeCoefficient(coeffs0,m1,w^1235*m5);


(Evaluate(ProdRel([5,13,22]),coeffs1)/w^1105+Evaluate(ProdRel([3,13,24]),coeffs1)/w^2730)/w^2990 eq (a+w^2415*b)*(m2+w^1235*m6)^2;
// Thus either a=w^2415*b (not true in the end) or m2=w^1235*m6 (is true in the end)
badcoeffs2:=ChangeCoefficient(coeffs1,a,w^2415*b);

Evaluate(ProdRel([3,5,9]),badcoeffs2) eq w^1430*(m2+w^1235*m6)*(m2+w^2600*m6)*(m2+w^3965*m6);
// The first of these is our previous (correct) assertion, so we may assume this is not the case.
badcoeffs2a:=ChangeCoefficient(badcoeffs2,m2,w^2600*m6);
badcoeffs2b:=ChangeCoefficient(badcoeffs2,m2,w^3965*m6);

Evaluate(ProdRel([1,3,15]),badcoeffs2a) eq w^2995*b*(m5+w^390*m6)^2;
Evaluate(ProdRel([3,5,6]),badcoeffs2a) eq w^455*(m5+w^65*m6)*(m5+w^1885*m6)*(m5+w^2145*m6);
// These are clearly contradictory, so branch 2a closes.

Evaluate(ProdRel([1,3,15]),badcoeffs2b) eq w^2995*b*(m5+w^3120*m6)^2;
Evaluate(ProdRel([3,5,6]),badcoeffs2b) eq w^455*(m5+w^520*m6)*(m5+w^780*m6)*(m5+w^2795*m6);
// These are again clearly contradictory, so branch 2b, hence the branch, closes.

coeffs2:=ChangeCoefficient(coeffs1,m2,w^1235*m6);

Evaluate(ProdRel([3,5,6]),coeffs2) eq w^455*(m5+w^1430*m6)*(m5+w^3250*m6)*(m5+w^3510*m6);
// The last of these is the correct one, so we eliminate the other two first.

badcoeffs3a:=ChangeCoefficient(coeffs2,m5,w^1430*m6);
badcoeffs3b:=ChangeCoefficient(coeffs2,m5,w^3250*m6);

Evaluate(ProdRel([1,3,15]),badcoeffs3a) eq w^1295*m6^2*(a+w^790*b);
Evaluate(ProdRel([1,3,15]),badcoeffs3b) eq w^1490*m6^2*(a+w^3130*b);

// This allows us to relate a to b.
badcoeffs3a1:=ChangeCoefficient(badcoeffs3a,a,w^790*b);
badcoeffs3b1:=ChangeCoefficient(badcoeffs3b,a,w^3130*b);


Factorization(Evaluate(ProdRel([3,5,18]),badcoeffs3a1)/w^2665-Evaluate(ProdRel([1,5,9]),badcoeffs3a1)/w^845) eq [<m7+w^3120*m8,1>,<m6,2>];
Factorization(Evaluate(ProdRel([3,5,18]),badcoeffs3b1)/w^2665-Evaluate(ProdRel([1,5,9]),badcoeffs3b1)/w^845) eq [<m7+w^390*m8,1>,<m6,2>];

Factorization(Evaluate(ProdRel([1,3,9]),badcoeffs3a1)/w^3965-Evaluate(ProdRel([1,6,9]),badcoeffs3a1)/w^2795) eq [<m7+w^3055*m8,1>,<m6,2>];
Factorization(Evaluate(ProdRel([1,3,9]),badcoeffs3b1)/w^3965-Evaluate(ProdRel([1,6,9]),badcoeffs3b1)/w^2795) eq [<m7 + w^3965*m8, 1>,<m6, 2>];

// This shows that there is a contradiction with both of these branches.

// We thus pursue the case m5=w^3510*m6, which will lead to the right answer.
coeffs3:=ChangeCoefficient(coeffs2,m5,w^3510*m6);
Evaluate(ProdRel([1,3,15]),coeffs3) eq w^3245*m6^2*(a+w^595*b);
// As with the bad branches, this relates a and b
coeffs4:=ChangeCoefficient(coeffs3,a,w^595*b);
coeffs4 eq [w^595*b,b,w^1235*w^3510*m6,w^1235*m6,m3,m4,w^3510*m6,m6,m7,m8];

// Now we need to take one more difference between relations, and then everything is easy from here on.
Evaluate(ProdRel([1,5,11]),coeffs4)*w^1300+Evaluate(ProdRel([1,3,9]),coeffs4)*w^3185 eq m6^2*(m7+w^585*m8);
coeffs5:=ChangeCoefficient(coeffs4,m7,w^585*m8);

Evaluate(ProdRel([1,6,24]),coeffs5)/w^3315 eq m6^2*(w^3970*b+m3);
Evaluate(ProdRel([1,8,22]),coeffs5)/w^2275 eq m6^2*(w^3385*b+m4);
Evaluate(ProdRel([1,3,24]),coeffs5)/w^325 eq m6^2*(w^3450*b+m8);

coeffs6:=ChangeCoefficients(coeffs5,[m3,m4,m8],[w^3970*b,w^3385*b,w^3450*b]);

// This satisfies all of the relations.
// Check that this is a solution. First, that all matrices with the specialization above yield the same conjugate of h1.

scal6:=Evaluate(scal,coeffs6);
sca:=GL(26,F)!Evaluate(scal6,[1:i in [1..10]]);

gg:=h1^sca;
mats!sca^-1*scal6*mats!gg eq mats!gg*mats!sca^-1*scal6;

// Then that gg does lie in F4.

{*CheckRel(i,gg):i in seqs*} eq {*0^^3276*};

// Then that we obtain PGL(2,13) containing H. (We check that j1 lies in F4 in the next function.)

J:=sub<GL(26,F)|l1,l2,gg,j1>;
Bool:=IsIsomorphic(J,PGL(2,13)); Bool;

return "";
end function;

// We now give code that can prove that the trilinear form we gave is the correct one. The normalizer of a torus
// has a small space of trilinear forms for F4. We use two random generators to quickly reduce to that space, and
// and then choose a random element of F4 to reduce the space of trilinear forms down to 1 dimension.

function CheckF4Form()

q:=2^12;
F<w>:=GF(q);
GG:=GroupOfLieType("F4",q);
z:=PrimitiveElement(F)^(Order(w) div 13);
W:=VectorSpace(GF(q),4);
rho:=StandardRepresentation(GG);
Over:=GL(26,q);
g1:=elt<GG|W![z,1,1,1]>;
g2:=elt<GG|W![1,z,1,1]>;
g3:=elt<GG|W![1,1,z,1]>;
g4:=elt<GG|W![1,1,1,z]>;
Refs:=Reflections(GG);
Mats:=[rho(i):i in [g1,g2,g3,g4]];
MatsE:=[rho(i):i in Refs];
NGT:=sub<Over|Mats cat MatsE>;

// Confirm that L is a subgroup of F4 via conjugation by x, and so is j1 in NGT:

x^-1*l1*x in NGT and x^-1*l2*x in NGT and x^-1*j1*x in NGT;
// Now build the matrix of linear relations that the trilinear form must satisfy.

mat:=[];
repeat y1:=Random(NGT); y2:=Random(NGT); until NGT eq sub<NGT|y1,y2>;
for h in [x*y1*x^-1,x*y2*x^-1] do
  for nn in [1..#seqs] do a:=seqs[nn,1]; b:=seqs[nn,2]; c:=seqs[nn,3];
    val:=[F!0:i in [1..#seqs]];
    for i in [1..NumberOfRows(h)] do if(h[a,i] ne 0) then for j in [1..NumberOfRows(h)] do if(h[b,j] ne 0) then for k in [1..NumberOfRows(h)] do
      val[Position(seqs,Sort([i,j,k]))]+:=h[a,i]*h[b,j]*h[c,k];
    end for; end if; end for; end if; end for;
    val[nn]-:=1; Append(~mat,val); delete val;
  end for;
  delete h;
end for;

repeat
  h:=x*rho(Random(GG))*x^-1;
  for iii in [1..5] do nn:=Random([1..#seqs]); a:=seqs[nn,1]; b:=seqs[nn,2]; c:=seqs[nn,3];
    val:=[F!0:i in [1..#seqs]];
    for i in [1..NumberOfRows(h)] do if(h[a,i] ne 0) then for j in [1..NumberOfRows(h)] do if(h[b,j] ne 0) then for k in [1..NumberOfRows(h)] do
      val[Position(seqs,Sort([i,j,k]))]+:=h[a,i]*h[b,j]*h[c,k];
    end for; end if; end for; end if; end for;
    val[nn]-:=1; Append(~mat,val); delete val;
  end for;
  delete h;
until Dimension(Nullspace(Transpose(Matrix(F,mat)))) eq 1;

ftest:=Nullspace(Transpose(Matrix(F,mat))).1;
"Is the trilinear form given in this document correct?";
return &and[ftest[i] eq f26[i]:i in [1..#f26]];

end function;

"The function CheckF4Form() checks that the trilinear form f26 from F4 is correct";
"The function CheckDetermination() checks that the determination of the possible conjugates into F4 is correct.";
