/* In this case we let H denote a copy of PSL(2,17), L denote the subgroup 17.8, and let G denote F4(13^4).
We give matrices defining L, an element of order 17 and an element of order 8. We give L directly as a subgroup
of the normalizer of a torus of F4(13^4), as constructed using the GroupOfLieType command. We give H as L and a
single new matrix, which were chosen more or less at random (preserving the direct sum decomposition to enable
smaller polynomials to be found). A posteriori, one could produce H as a subgroup of F4 and show that no other
overgroup of L preserves the F4 trilinear form (other than its other conjugate), but we prefer a less synthetic
proof.*/
F<ww>:=GF(13^4);
zz:=ww^1680;

l1:=GL(26,F)!DiagonalMatrix([zz^i:i in [12,15,1,5,4,8,11,7,10,11,14,14,0,0,3,3,7,6,6,10,9,13,12,16,2,5]]);

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

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

h1:=GL(26,F)![[ww^10540,0,0,0,0,0,ww^1530,0,ww^15130,0,ww^4080,0,ww^26152,0,ww^10200,0,ww^25330,0,ww^5610,0,0,0,0,0,0,ww^8500],
[0,ww^1020,ww^15470,ww^24724,ww^4080,ww^8330,0,ww^15374,0,ww^26934,0,ww^23534,ww^16111,ww^4211,0,ww^7214,0,ww^8574,0,ww^15374,ww^8330,ww^10200,ww^22684,ww^15470,ww^13260,0],
[0,ww^15470,ww^8330,ww^8574,ww^15300,ww^10200,0,ww^24724,0,ww^7214,0,ww^1094,ww^16111,ww^4211,0,ww^1094,0,ww^9254,0,ww^22684,ww^4080,ww^13260,ww^12654,ww^8330,ww^1190,0],
[0,ww^27296,ww^11146,ww^10200,ww^3666,ww^11826,0,ww^8330,0,ww^15300,0,ww^15470,ww^17397,ww^5497,0,ww^1190,0,ww^13260,0,ww^22610,ww^9786,ww^17946,ww^18360,ww^946,ww^25256,0],
[0,ww^4080,ww^15300,ww^1094,ww^22610,ww^1190,0,ww^7214,0,ww^22684,0,ww^8574,ww^1831,ww^18491,0,ww^12654,0,ww^10444,0,ww^23534,ww^1190,ww^22610,ww^1094,ww^13260,ww^24480,0],
[0,ww^8330,ww^10200,ww^9254,ww^1190,ww^13260,0,ww^8574,0,ww^1094,0,ww^8404,ww^16111,ww^4211,0,ww^10444,0,ww^15374,0,ww^12654,ww^15300,ww^1190,ww^21494,ww^4080,ww^22610,0],
[ww^1530,0,0,0,0,0,ww^10200,0,ww^10540,0,ww^850,0,ww^11872,0,ww^25330,0,ww^22780,0,ww^4080,0,0,0,0,0,0,ww^5610],
[0,ww^17946,ww^27296,ww^8330,ww^9786,ww^11146,0,ww^15470,0,ww^4080,0,ww^27540,ww^17397,ww^5497,0,ww^15300,0,ww^10200,0,ww^1190,ww^946,ww^11826,ww^22610,ww^10976,ww^17946,0],
[ww^15130,0,0,0,0,0,ww^10540,0,ww^18360,0,ww^5610,0,ww^11872,0,ww^1530,0,ww^10200,0,ww^8500,0,0,0,0,0,0,ww^11050],
[0,ww^946,ww^9786,ww^15300,ww^25256,ww^3666,0,ww^4080,0,ww^1190,0,ww^8330,ww^3117,ww^19777,0,ww^22610,0,ww^1190,0,ww^24480,ww^17946,ww^13016,ww^27540,ww^11826,ww^25426,0],
[ww^4080,0,0,0,0,0,ww^850,0,ww^5610,0,ww^22780,0,ww^11872,0,ww^24820,0,ww^15810,0,ww^25330,0,0,0,0,0,0,ww^10200],
[0,ww^26106,ww^3666,ww^15470,ww^11146,ww^10976,0,ww^27540,0,ww^8330,0,ww^18360,ww^3117,ww^19777,0,ww^10200,0,ww^8330,0,ww^15300,ww^27296,ww^946,ww^1190,ww^17946,ww^9786,0],
[ww^28,0,0,0,0,0,ww^14308,0,ww^14308,0,ww^14308,0,6,0,ww^28,0,ww^14308,0,ww^28,0,0,0,0,0,0,ww^14308],
[ww^26208,ww^7689,ww^7689,ww^6403,ww^21969,ww^7689,ww^11928,ww^6403,ww^11928,ww^20683,ww^11928,ww^20683,8,3,ww^26208,ww^20683,ww^11928,ww^6403,ww^26208,ww^6403,ww^21969,ww^7689,ww^6403,ww^21969,ww^7689,ww^11928],
[ww^10200,0,0,0,0,0,ww^25330,0,ww^1530,0,ww^24820,0,ww^26152,0,ww^22780,0,ww^19890,0,ww^850,0,0,0,0,0,0,ww^4080],
[0,ww^9786,ww^3666,ww^1190,ww^15226,ww^13016,0,ww^15300,0,ww^22610,0,ww^10200,ww^3117,ww^19777,0,ww^18360,0,ww^22610,0,ww^27540,ww^25256,ww^25426,ww^15470,ww^17946,ww^26106,0],
[ww^25330,0,0,0,0,0,ww^22780,0,ww^10200,0,ww^15810,0,ww^11872,0,ww^19890,0,ww^18360,0,ww^24820,0,0,0,0,0,0,ww^850],
[0,ww^11146,ww^11826,ww^13260,ww^13016,ww^17946,0,ww^10200,0,ww^1190,0,ww^8330,ww^17397,ww^5497,0,ww^22610,0,ww^1190,0,ww^18360,ww^3666,ww^25256,ww^1020,ww^9786,ww^15226,0],
[ww^5610,0,0,0,0,0,ww^4080,0,ww^8500,0,ww^25330,0,ww^26152,0,ww^850,0,ww^24820,0,ww^10200,0,0,0,0,0,0,ww^1530],
[0,ww^17946,ww^25256,ww^22610,ww^26106,ww^15226,0,ww^1190,0,ww^24480,0,ww^15300,ww^17397,ww^5497,0,ww^27540,0,ww^18360,0,ww^15470,ww^25426,ww^24066,ww^8330,ww^13016,ww^17946,0],
[0,ww^8330,ww^4080,ww^7214,ww^1190,ww^15300,0,ww^26934,0,ww^15374,0,ww^24724,ww^1831,ww^18491,0,ww^22684,0,ww^1094,0,ww^22854,ww^13260,ww^1190,ww^23534,ww^10200,ww^22610,0],
[0,ww^10200,ww^13260,ww^15374,ww^22610,ww^1190,0,ww^9254,0,ww^10444,0,ww^26934,ww^16111,ww^4211,0,ww^22854,0,ww^22684,0,ww^21494,ww^1190,ww^22610,ww^15374,ww^15300,ww^18360,0],
[0,ww^25256,ww^15226,ww^18360,ww^3666,ww^24066,0,ww^22610,0,ww^27540,0,ww^1190,ww^17397,ww^5497,0,ww^15470,0,ww^1020,0,ww^8330,ww^26106,ww^17946,ww^10200,ww^25426,ww^27296,0],
[0,ww^15470,ww^8330,ww^26934,ww^13260,ww^4080,0,ww^8404,0,ww^9254,0,ww^15374,ww^1831,ww^18491,0,ww^15374,0,ww^7214,0,ww^10444,ww^10200,ww^15300,ww^22854,ww^8330,ww^1190,0],
[0,ww^13260,ww^1190,ww^22684,ww^24480,ww^22610,0,ww^15374,0,ww^22854,0,ww^7214,ww^16111,ww^4211,0,ww^23534,0,ww^12654,0,ww^15374,ww^22610,ww^18360,ww^24724,ww^1190,ww^1020,0],
[ww^8500,0,0,0,0,0,ww^5610,0,ww^11050,0,ww^10200,0,ww^11872,0,ww^4080,0,ww^850,0,ww^1530,0,0,0,0,0,0,ww^10540]];

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

// 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:=[273,295,313,327,332,571,584,597,608,628,644,840,855,883,926,942,1081,1098,1158,1181,1209,1295,1314,1410,1433,1452,1581,1613,
1640,1653,1674,1773,1807,1836,1861,2012,2039,2052,2072,2167,2196,2220,2361,2374,2384,2483,2507,2615,2718,2731];
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,-1,1,1,-2,-2];

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,m1,m2,m3,m4>:=PolynomialRing(F,7);
mats:=MatrixRing(R,26);
scal:=mats!DiagonalMatrix([c:i in [1..26]]);
scal[13,13]:=a; scal[14,14]:=a-2*b; scal[14,13]:=b;

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

for i in [1..#a2] do
  scal[a2[i],a2[i]]:=m1; scal[a2[i],a3[i]]:=m2;
  scal[a3[i],a2[i]]:=m3; scal[a3[i],a3[i]]:=m4;
end for;
scal[1,23]:=-m2;
scal[17,8]:=-m2;
scal[11,12]:=-m2;
scal[9,20]:=-m2;

scal[8,17]:=-m3;
scal[12,11]:=-m3;
scal[20,9]:=-m3;
scal[23,1]:=-m3;

mats!l1*scal eq scal*mats!l1 and mats!l2*scal eq scal*mats!l2;

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 CentRel(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)*scal; v2:=Matricise(v)*scal; w2:=Matricise(w)*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,c,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;

function CheckDetermination()

ProdRel([1,11,13]) eq m2^2*(a + ww^16688*m1);
// The branch m2=0 is the bad one.

badcoeffs1:=ChangeCoefficient(coeffs0,m2,0);

// Note that a-2b is a factor of the determinant, and so a cannot equal 2b.
Evaluate(ProdRel([1,1,2]),badcoeffs1) eq ww^1775*(a+ww^28*m1)*(a+ww^14308*m1)*(a-2*b);
Evaluate(ProdRel([2,13,13]),badcoeffs1) eq ww^1831*(a + ww^11928*m1)*(a + ww^26208*m1)*(a-2*b);
// This yields a contradiction.

coeffs1:=ChangeCoefficient(coeffs0,m1,-a/ww^16688);
Evaluate(ProdRel([1,6,8])*ww^23828-ProdRel([1,23,25])*ww^21448,coeffs1) eq a*c*(m4-ww^13186*c);
// Since a and c cannot be zero, we obtain the following.
coeffs2:=ChangeCoefficient(coeffs1,m4,ww^13186*c);

a1:=(a*c+ww^13759*a*m2+ww^1859*b*m2+ww^10642*m2*m3);
a2:=a*c + ww^8262*m2*m3;

Evaluate(ProdRel([1,14,26]),coeffs2) eq ww^24321*m2*a1;
Evaluate(ProdRel([1,2,12]),coeffs2) eq ww^15538*c*a2;
// Since c ne 0, a2 must equal zero. But a ne 0 as well, so m2 and m3 cannot be 0.
// Thus a1=0. In particular,
(a1-a2)/ww^22542 eq m2*(ww^19777*a + ww^7877*b + 12*m3);

coeffs3:=ChangeCoefficient(coeffs2,m3,ww^19777*a+ww^7877*b);
Evaluate(ProdRel([1,1,4]),coeffs3) eq ww^7821*(a-2*b)*(a + ww^4788*m2)*(a + ww^19068*m2);
// As a cannot equal 2b, we obtain two possibilites.

coeffs4a:=ChangeCoefficient(coeffs3,m2,-a/ww^4788);
coeffs4b:=ChangeCoefficient(coeffs3,m2,-a/ww^19068);
// Each of these works, and yields two different solutions.

Evaluate(ProdRel([1,1,2]),coeffs4a) eq ww^8915*a^2*(a-2*b+ww^19589*c);
Evaluate(ProdRel([1,1,2]),coeffs4b) eq ww^8915*a^2*(a-2*b+ww^5309*c);

coeffs5a:=ChangeCoefficient(coeffs4a,c,(2*b-a)/ww^19589);
coeffs5b:=ChangeCoefficient(coeffs4b,c,(2*b-a)/ww^5309);


// These satisfy all relations. Remove the centralizer of H in GL(26,k), which involves specializing a and b.

coeffs6a:=ChangeCoefficients(coeffs5a,[a,b],[1,0]);
coeffs6b:=ChangeCoefficients(coeffs5b,[a,b],[1,0]);

sca1:=GL(26,F)!Evaluate(scal,coeffs6a);
sca2:=GL(26,F)!Evaluate(scal,coeffs6b);
h11:=h1^sca1;
h12:=h1^sca2;

H1:=sub<GL(26,F)|l1,l2,h11>;
H2:=sub<GL(26,F)|l1,l2,h12>;

// We check that H1 < G.

{CheckRel(i,h11):i in seqs} eq {0};
// H2 will be shown later to be conjugate to H1, so we don't need to check that.

// Now we find the centralizer of L in G, which has order 2, and note that this element conjugates H1 to H2.

CentRel([1,13,26]) eq a*m2^2;
CentRel([3,8,21]) eq c^2*m3;
CentRel([5,14,22]) eq b*c^2;
// Thus m2=m3=b=0, since a and c cannot be 0.

CentRel([2,13,25]) eq a*c^2-1;
CentRel([5,6,26]) eq m1*c^2-1;
// Thus a=m1.

CentRel([5,17,18]) eq c*m1*m4-1;
// Since a*c^2=1 and a=m1, we see that c=m4. To go further, its quickest to evaluate another relation with this information:

c1:=[a,0,c,a,0,0,c];

// We need to remove the centre of GL(26,k), which gives us the three scalar matrices. Easiest way to do that is to
// force the determinant to be 1.

Evaluate(Determinant(scal),c1) eq a^10*c^16;
// Since a*c^2=1, this means that (with that condition) the determinant is a^2. Thus a=\pm 1.
Evaluate(CentRel([9,14,17]),c1) eq a^3-1;
// Hence a=1. Thus c=\pm 1 as a*c^2=1.
c2:=[1,0,-1,1,0,0,-1];
l3:=GL(26,F)!Evaluate(scal,c2);
h11^l3 eq h12;

// We will check that l3 lies in NGT in the other function.

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 3-space of trilinear forms for F4. We use two random generators to quickly reduce to that space, and
// then choose a random element of F4 to reduce the space of trilinear forms down to 1 dimension.

function CheckF4Form()

q:=13^4;
z:=ww^420;
GG:=GroupOfLieType("F4",q);
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]>;
W:=Reflections(GG);
Mats1:=[rho(i):i in [g1,g2,g3,g4]];
MatsE:=[rho(i):i in W];
NGT:=sub<Over|Mats1 cat MatsE>;
T:=sub<NGT|Mats1>;
T17:=sub<NGT|[i^4:i in Mats1]>;
E:=sub<NGT|MatsE>;

// Confirm that L is a subgroup of F4, and that l3 above is in NGT:

l1 in NGT and l2 in NGT;

l3:=Evaluate(scal,[1,0,-1,1,0,0,-1]); l3 in NGT;

// Confirm our other statements from the paper. First, that the centralizer of the element of order 8 has order exactly 2.

CCE:=ConjugacyClasses(E);
Z1:=[i[3]:i in CCE|i[1] eq 8 and not(i[3]^4 in T)];
S1:=sub<E|Z1[1]>;
S2:=sub<E|Z1[2]>;

// Since we think that the centralizer has order 2, we expect to find two classes of subgroups of order 8

CS1:=Centralizer(T,S1); CS2:=Centralizer(T,S2);
#CS1 eq 2 and #CS2 eq 2;

// Since T contains the subgroup 4^4, we see that the (2-part of the) centralizer has order exactly 2, as verified by another means above.
// Now to prove uniqueness of L in algebraic G.

NT:=sub<NGT|T17,S1>;
subNT:=[i`subgroup:i in Subgroups(NT:OrderEqual:=17*8)|Order(Center(i`subgroup)) eq 1];
#subNT eq 4;
&and[IsConjugate(NGT,i,subNT[1]):i in subNT];


// Now build the matrix of linear relations that the trilinear form must satisfy.
repeat e1:=Random(E); e2:=Random(E); until E eq sub<E|e1,e2>;

mat:=[];
for h in [e1,e2] 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;
end for;

repeat
  h:=rho(Random(GG));
  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.";
