/* In this case we let H denote a copy of PSL(2,25), L denote the subgroup 5^2.12, and let G denote F4(16).
   We give matrices defining L, an element of order 5 and an element of order 12. We give an element x that
   conjugates our version of L into the normalizer of a torus of F4(16), as constructed using the GroupOfLieType
   command. We give H as a pair of matrices, one from 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.
*/

// First we set up L, H and give the matrix x that does the conjugation.
F<ww>:=GF(16);

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,ww^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,ww^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,ww^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,ww^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,ww^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,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,ww^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,ww^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,ww^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,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,ww^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,ww^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,ww^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,ww^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,ww^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,ww^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,ww^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,ww^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,ww^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,ww^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,ww^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,ww^9]];

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

l3x:=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,0,0,ww^3,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,ww^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,0,0,0,0,0,0,0,0,0,0,0,ww^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,ww^6,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,ww^12,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,ww^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,0,0,0,0,ww^9,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,ww^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,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,ww^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,ww^3,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,1,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,ww^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,ww^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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,ww^9,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,ww^6,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,ww^6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,ww^3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,ww^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,ww^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,0,0,0,0,0,0,0,0,0,0,0,0,ww^12,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,ww^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]];

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

h2:=GL(26,F)![[0,1,ww^11,ww^7,ww^12,ww^8,ww^9,ww^4,ww,ww^2,ww^6,ww^3,ww^13,ww^14,ww^10,ww^5,ww^2,ww^3,ww^7,ww^3,ww^3,ww^3,ww^4,ww^14,ww^11,ww],
[1,0,ww^6,ww^7,ww^2,ww^3,ww^14,ww^4,ww,ww^12,ww^11,ww^8,ww^13,ww^9,ww^5,ww^10,ww^7,ww^3,ww^2,ww^3,ww^3,ww^3,ww^14,ww^4,ww,ww^11],
[ww^10,1,ww^11,ww^7,ww^7,0,ww^9,0,0,0,ww^6,ww^8,ww^8,ww^4,1,0,ww^2,ww^8,ww^2,ww^3,ww^13,ww^8,ww^14,ww^14,0,ww],
[ww^5,ww^5,ww,ww^12,0,0,ww^14,ww^4,ww^11,0,ww,0,ww^3,ww^14,1,1,ww^2,ww^8,ww^2,0,0,ww^8,ww^14,ww^14,ww,ww],
[1,ww^10,ww,0,ww^7,ww^3,0,ww^14,ww,ww^12,0,ww^8,0,ww^14,ww^5,ww^5,ww^7,ww^13,0,ww^8,ww^8,ww^3,0,ww^9,ww^11,ww^11],
[ww^10,1,0,0,ww^12,ww^13,ww^4,ww^4,ww^11,ww^2,ww^11,ww^3,0,0,ww^5,ww^5,0,ww^13,ww^12,ww^8,ww^8,ww^3,ww^4,0,ww^11,ww^11],
[1,ww^10,ww^6,ww^2,0,ww^13,ww^14,0,0,ww^2,ww,0,ww^13,ww^9,0,ww^10,ww^2,ww^8,ww^2,ww^3,ww^13,ww^8,ww^14,ww^14,ww^6,0],
[ww^5,ww^5,0,ww^7,ww^2,ww^13,0,ww^9,ww^6,ww^2,0,ww^13,ww^8,0,ww^5,ww^5,ww^7,0,ww^7,ww^8,ww^8,0,ww^9,ww^9,ww^11,ww^11],
[ww^5,ww^5,0,ww^2,ww^7,ww^8,0,ww^9,ww^6,ww^7,0,ww^8,ww^13,0,ww^5,ww^5,ww^12,0,ww^12,ww^8,ww^8,0,ww^4,ww^4,ww^11,ww^11],
[ww^10,1,0,0,ww^12,ww^8,ww^14,ww^14,ww,ww^7,ww,ww^3,0,0,ww^5,ww^5,0,ww^3,ww^7,ww^8,ww^8,ww^13,ww^9,0,ww^11,ww^11],
[1,ww^10,ww^6,ww^7,0,ww^8,ww^4,0,0,ww^7,ww^11,0,ww^8,ww^9,0,1,ww^2,ww^8,ww^2,ww^13,ww^3,ww^8,ww^14,ww^14,ww,0],
[1,ww^10,ww^11,0,ww^2,ww^3,0,ww^4,ww^11,ww^12,0,ww^13,0,ww^4,ww^5,ww^5,ww^12,ww^3,0,ww^8,ww^8,ww^13,0,ww^4,ww^11,ww^11],
[ww^5,ww^5,ww^11,ww^12,0,0,ww^4,ww^14,ww,0,ww^11,0,ww^3,ww^4,ww^10,ww^10,ww^2,ww^8,ww^2,0,0,ww^8,ww^14,ww^14,ww^6,ww^6],
[ww^10,1,ww,ww^2,ww^2,0,ww^9,0,0,0,ww^6,ww^13,ww^13,ww^14,ww^10,0,ww^2,ww^8,ww^2,ww^13,ww^3,ww^8,ww^14,ww^14,0,ww^6],
[ww^2,ww^7,ww^8,ww^14,ww^4,ww^10,0,ww,ww^13,ww^4,0,ww^10,1,ww^6,0,ww^2,ww^9,ww^10,ww^4,0,0,ww^5,ww^11,ww^6,ww^13,0],
[ww^7,ww^2,0,ww^14,ww^4,ww^10,ww^6,ww,ww^13,ww^4,ww^8,ww^10,1,0,ww^2,0,ww^4,ww^5,ww^9,0,0,ww^10,ww^6,ww^11,0,ww^13],
[ww^7,ww^2,ww^13,ww^4,ww^9,0,ww,ww^6,ww^8,0,ww^13,ww^5,ww^10,ww,ww^12,ww^7,0,0,ww^4,ww^10,ww^5,0,ww^11,0,ww^8,ww^3],
[ww^12,ww^12,ww^13,ww^4,ww^9,1,ww,0,0,ww^14,ww^13,ww^5,ww^10,ww,ww^7,ww^2,0,ww^5,0,1,1,ww^10,0,0,ww^13,ww^8],
[ww^2,ww^7,ww^13,ww^4,0,ww^5,ww,ww^6,ww^8,ww^9,ww^13,0,ww^10,ww,ww^7,ww^12,ww^4,0,0,ww^5,ww^10,0,0,ww^11,ww^3,ww^8],
[ww^12,ww^12,ww^8,0,ww^4,ww^10,ww^11,ww,ww^13,ww^4,ww^3,ww^10,0,ww^6,0,0,ww^4,1,ww^14,ww^10,ww^5,1,ww,ww^11,0,0],
[ww^12,ww^12,ww^3,0,ww^4,ww^10,ww^6,ww,ww^13,ww^4,ww^8,ww^10,0,ww^11,0,0,ww^14,1,ww^4,ww^5,ww^10,1,ww^11,ww,0,0],
[ww^12,ww^12,ww^13,ww^4,ww^14,ww^5,ww,0,0,ww^9,ww^13,1,ww^10,ww,ww^2,ww^7,0,ww^10,0,1,1,ww^5,0,0,ww^8,ww^13],
[ww^2,ww^7,ww^13,ww^4,0,1,ww,ww^11,ww^3,ww^14,ww^13,0,ww^10,ww,ww^2,ww^12,ww^14,0,0,ww^10,ww^5,0,0,ww,ww^3,ww^13],
[ww^7,ww^2,ww^13,ww^4,ww^14,0,ww,ww^11,ww^3,0,ww^13,1,ww^10,ww,ww^12,ww^2,0,0,ww^14,ww^5,ww^10,0,ww,0,ww^13,ww^3],
[ww^7,ww^2,0,ww^9,ww^4,ww^10,ww^11,ww,ww^13,ww^4,ww^3,ww^10,ww^5,0,ww^7,0,ww^14,ww^10,ww^9,0,0,ww^5,ww^6,ww,0,ww^8],
[ww^2,ww^7,ww^3,ww^9,ww^4,ww^10,0,ww,ww^13,ww^4,0,ww^10,ww^5,ww^11,0,ww^7,ww^9,ww^5,ww^14,0,0,ww^10,ww,ww^6,ww^8,0]];

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

x:=GL(26,F)![[0,0,0,0,0,0,0,0,0,0,0,0,ww^4,ww^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,ww^2,ww^12,0,0,0,0,0,0,0,0,0,0,0,0],
[ww^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,ww^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,ww^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,ww^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,ww^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,ww^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,ww^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,ww^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,ww^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,ww^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,ww^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,ww^4],
[0,ww,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,ww,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,ww,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,ww,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,ww,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,ww,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,ww,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,ww,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,ww,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,ww,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,ww,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,ww,0]];

// We need the indices of the trilinear form as well.
seqs:=[[i,j,k]:i,j,k in [1..26]|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:=[63,85,106,125,146,163,285,295,304,311,320,325,388,410,431,450,471,488,610,620,629,636,645,650,823,843,869,969,1050,1113,1142,
1216,1336,1408,1418,1468,1545,1551,1638,1691,1789,1923,2008,2047,2112,2274,2330,2392,2546,2657,2751,2850,2974,3035,3098,3126];

vals:=[1,ww^10,ww^5,1,ww^5,ww^10,ww^14,ww^4,ww^4,ww^9,ww^14,ww^9,ww^3,ww^8,ww^13,ww^3,ww^13,ww^8,ww^7,ww^2,ww^2,ww^12,ww^7,ww^12,ww^7,
ww^10,ww^7,ww^4,ww^10,ww^7,ww^7,ww^4,ww^10,ww^7,ww^7,ww^4,ww^7,ww^10,ww^7,ww^4,ww^7,ww^4,ww^7,ww^7,ww^4,ww^4,ww^7,ww^4,ww^4,ww^4,
ww^4,ww^4,ww,ww,ww,ww];

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.

R<a,b,c,d>:=PolynomialRing(F,4);
scal:=DiagonalMatrix([a,b,c,c,c,c,c,c,c,c,c,c,c,c,d,d,d,d,d,d,d,d,d,d,d,d]);

// 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
mats:=MatrixRing(R,26);

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!h2;

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,3,6]) eq ww^10*c*(c+ww^5*d)^2;

// So c=ww^5*d
coeffs1:=[a,b,ww^5*d,d];

// Then we obtain the other two coefficients:
Evaluate(ProdRel([1,3,16]),coeffs1) eq d^2*(b+ww^12*d);
Evaluate(ProdRel([1,3,4]),coeffs1) eq d^2*(a+ww^5*d);

coeffs2:=[ww^5*d,ww^12*d,ww^5*d,d];

// Now we check that this element lies in F4. First construct the appropriate conjugate.
scal2:=Evaluate(scal,coeffs2);
sca:=GL(26,F)!Evaluate(scal2,[1:i in [1..4]]);

gg:=h2^sca;

// Then that gg does lie in F4.
{CheckRel(i,gg):i in seqs} eq {0};

// Now check that the normalizer of L normalizes H1.
H1:=sub<GL(26,F)|l1,l2,gg>;
NH:=sub<GL(26,F)|l1,l2,gg,l3x^(x^-1)>;
Index(NH,H1) eq 2 and IsIsomorphic(H1,PSL(2,25));
// We will check that <L,l3> is the normalizer of L in NGT, which is the normalizer of L in G, in CheckF4Form.

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:=16;
GG:=GroupOfLieType("F4",q);
z:=PrimitiveElement(F)^3;
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:

L.1^x in NGT and L.2^x in NGT;
Lx:=sub<NGT|l1^x,l2^x>;

Normalizer(NGT,Lx) eq sub<NGT|Lx,l3x>;

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