/* In this case we let H denote a copy of PSL(3,3), L denote the parabolic subgroup, and let G denote F4(4).
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(7^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.

There are three conjugacy classes of subgroups in N_G(T) isomorphic to L. The one has the wrong Brauer character
on the minimal module, and the other is the image of it under the graph automorphism, hence has the wrong
Brauer character on the adjoint representation. Thus L is unique up to isomorphism in N_G(T), hence in G.

*/
F<w>:=GF(4);

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

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

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

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

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

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

H:=sub<GL(26,F)|[l1,l2,l3,h1]>;

j1:=GL(26,F)![[1,0,0,w^2,w^2,0,w^2,w^2,w^2,w^2,w^2,0,w^2,w^2,w^2,1,w^2,w,w^2,w,1,1,w,1,1,w^2],
[0,1,0,1,w,0,w,1,1,w,1,w^2,w,1,w,w^2,1,w^2,w,0,w,w^2,w,w^2,w,w^2],
[1,w,w^2,0,0,0,1,0,w,0,0,w^2,w,w^2,0,0,w,w,w^2,0,1,w,0,w^2,0,w],
[0,0,w^2,w^2,w^2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0],
[w^2,w,w^2,0,0,0,0,1,0,w,w,w^2,0,0,w^2,1,w^2,0,w,0,w,1,0,0,w^2,1],
[w^2,w,0,1,0,w,0,0,0,w^2,w^2,0,1,0,w,w^2,w^2,0,0,w^2,1,0,w,w^2,w,0],
[0,0,0,0,0,w,w,w,0,0,w^2,0,w^2,0,0,0,0,0,0,0,0,0,0,w^2,w^2,w^2],
[1,w,0,0,1,w,0,0,w^2,0,1,0,w^2,w,0,0,0,1,w^2,w,0,1,w^2,w,w^2,0],
[w^2,w,0,0,w,0,w^2,0,0,1,0,w^2,w^2,0,1,w,0,0,w,w,0,w^2,w^2,1,0,w],
[1,w,0,w,0,0,0,w^2,1,0,w^2,w^2,0,1,0,0,w,w^2,0,w^2,w^2,0,w,0,1,1],
[w^2,0,0,1,1,0,0,0,1,1,1,0,1,w^2,w^2,1,w^2,w,w^2,1,1,1,1,w^2,w^2,0],
[1,w,0,0,w,0,w,w^2,w^2,0,0,w,w,w,0,0,0,1,1,w^2,0,w,1,0,w,0],
[w,1,0,1,0,1,1,w^2,0,w,1,w^2,w^2,0,1,w,w^2,0,0,w^2,1,0,w,1,1,w^2],
[w^2,w,0,0,w^2,0,w,0,0,1,0,w,w,1,1,w,0,w,w^2,w^2,0,1,1,w^2,0,1],
[w^2,1,0,w,0,0,0,1,w^2,0,1,1,0,w^2,w,1,w,w,0,w^2,w^2,0,w,0,w,w],
[w,w^2,0,1,0,0,0,w^2,w,0,w^2,w^2,0,w,w^2,w,1,1,0,w,w,0,1,0,1,1],
[w^2,1,w,w,w^2,0,w,0,1,0,0,1,w^2,w,0,0,w^2,1,0,0,1,1,w,1,0,w^2],
[1,w^2,0,0,1,0,w^2,0,0,w,0,w^2,w^2,w^2,w,w^2,0,1,1,1,0,w,w,1,0,w],
[w,0,0,0,0,0,w,w,w,w,w^2,0,w^2,w^2,w^2,1,w^2,w,w^2,0,w,w,0,1,1,1],
[1,w^2,0,1,w^2,0,0,1,0,w^2,w,w^2,0,0,1,w,w^2,0,w^2,w^2,0,w^2,0,0,w^2,1],
[1,w,w,1,w,0,w^2,0,w,0,0,w,1,w^2,0,0,0,w,w^2,0,w,1,1,w,0,1],
[0,w,w,0,0,0,w^2,1,w^2,1,w,1,1,1,w,w^2,w^2,w^2,0,0,w,0,0,w,w^2,0],
[w^2,w,w,w,1,0,0,w^2,0,w,1,w,0,0,w^2,1,w^2,0,0,1,1,w,0,0,w,w^2],
[w^2,w,0,w,0,0,w^2,w,0,w^2,w,w,0,0,w,w^2,1,0,0,1,w,0,w^2,w,0,w],
[w^2,w^2,0,w,w^2,1,0,0,w^2,w,1,0,0,w,1,w,1,1,w,0,w,w^2,1,0,w^2,0],
[w^2,1,0,0,1,1,w^2,1,w,0,w^2,w^2,1,1,0,0,0,w^2,w^2,w,0,1,w^2,1,1,0]];

J:=sub<GL(26,F)|[l1,l2,l3,h1,j1]>;

/* 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.
*/

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

ents:=[61,97,118,135,138,157,175,186,190,205,227,240,255,262,263,277,304,322,330,345,386,400,408,422,432,460,479,496,500,515,529,
549,552,565,568,576,587,613,629,659,665,670,702,769,789,845,846,877,896,909,924,927,943,952,1026,1064,1135,1139,1150,1151,1195,
1201,1216,1220,1226,1234,1278,1299,1377,1389,1426,1439,1450,1455,1468,1474,1480,1492,1528,1603,1631,1651,1654,1661,1676,1696,1732,
1825,1828,1855,1859,1877,1885,1888,1905,1943,2016,2030,2043,2055,2063,2066,2077,2094,2097,2135,2181,2197,2209,2231,2262,2277,2285,
2292,2295,2300,2331,2348,2366,2374,2414,2432,2437,2444,2450,2452,2478,2509,2524,2527,2538,2543,2555,2593,2624,2634,2641,2645,2656,
2673,2676,2712,2740,2751,2759,2762,2779,2820,2884,2889,2896,2902,2960,2968,2975,2978,3011,3026,3034,3049,3074,3077,3080,3098,3128,
3133,3148,3166,3201,3228];

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


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

R<a,b,c>:=PolynomialRing(F,3);
mats:=MatrixRing(R,NumberOfRows(l1));
scal:=mats!DiagonalMatrix([a,a,b,b,b,b,b,b,b,b,c,c,c,c,c,c,c,c,c,c,c,c,c,c,c,c]);
scal*(mats!l1) eq (mats!l1)*scal and scal*(mats!l2) eq (mats!l2)*scal and scal*(mats!l3) eq (mats!l3)*scal;

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;

function CheckDetermination()

ProdRel([1,4,18]) eq w*c^2*(a+c);
ProdRel([1,4,21]) eq c^2*(b+w*c);

// Thus a=c and b=w*c

// We now double check by placing J=PSL_4(3) inside F4 to contain this copy of H.

coeffs:=[1,w,1];
sca:=GL(26,F)!Evaluate(scal,coeffs);
j2:=j1^sca;
h2:=h1^sca;

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

return "";
end function;

// We now give code that can prove that the trilinear form we gave is the correct one.

function CheckF4Form()

q:=4;
F<w>:=GF(q);
GG:=GroupOfLieType("F4",q);
W:=VectorSpace(GF(q),4);
rho:=StandardRepresentation(GG);
Over:=GL(26,q);
g1:=elt<GG|W![w,1,1,1]>;
g2:=elt<GG|W![1,w,1,1]>;
g3:=elt<GG|W![1,1,w,1]>;
g4:=elt<GG|W![1,1,1,w]>;
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:

x*l1*x^-1 in NGT and x*l2*x^-1 in NGT and x*l3*x^-1 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^-1*y1*x,x^-1*y2*x] 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^-1*rho(Random(GG))*x;
  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.";
