/* In this case we let H denote a copy of PSL(2,7), L denote the subgroup 7.3, and let G denote E6(29^2).
   We give matrices defining L, an element of order 7 and an element of order 3. 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 work form there, but we
   prefer a less synthetic proof. Indeed, it seems much harder to do it this way, as the equations are
   incredibly opaque.
*/

q:=29^2;
F<ww>:=GF(q);

l1:=GL(27,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,28,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28],
[0,0,0,0,0,0,28,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28,0],
[0,0,0,28,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,28,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,28,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28,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,28,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,28,0,0,0,0,0,0,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,28,0,0,0,0,0,0,0,0,0,0,0,0,0,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,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,28,0,0,0,0,0,0,0],
[28,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28,0,0,0,0,0],
[0,28,0,0,0,0,0,0,0,0,0,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]];


l2:=GL(27,F)!DiagonalMatrix([25,7,20,24,20,25,23,16,24,25,23,1,1,16,20,7,1,20,24,7,23,24,16,23,16,7,25]);

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


// Here is the centralizer in GL(27,F) for L, even as a subgroup of NGT. This makes computation incredibly easy

a1:=[1,22,23];
a2:=[27,4,8];
a3:=[6,9,25];
a4:=[10,19,14];

b1:=[2,18,24];
b2:=[20,3,21];
b3:=[26,5,7];
b4:=[16,15,11];

R<a,b,c,m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13,m14,m15,m16,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n13,n14,n15,n16,
t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15>:=PolynomialRing(F,50);
mats:=MatrixRing(R,27);
scal:=mats!DiagonalMatrix([a:i in [1..27]]);
scal[12,13]:=b; scal[13,17]:=-b; scal[17,12]:=-b;
scal[12,17]:=c; scal[13,12]:=-c; scal[17,13]:=c;

// m12 should be negated to make describing the action of the torus S easier.

for i in [1..#a1] do
  scal[a1[i],a1[i]]:=m1; scal[a1[i],a2[i]]:=m2; scal[a1[i],a3[i]]:=m3; scal[a1[i],a4[i]]:=m4;
  scal[a2[i],a1[i]]:=m5; scal[a2[i],a2[i]]:=m6; scal[a2[i],a3[i]]:=m7; scal[a2[i],a4[i]]:=m8;
  scal[a3[i],a1[i]]:=m9; scal[a3[i],a2[i]]:=m10; scal[a3[i],a3[i]]:=m11; scal[a3[i],a4[i]]:=-m12;
  scal[a4[i],a1[i]]:=m13; scal[a4[i],a2[i]]:=m14; scal[a4[i],a3[i]]:=m15; scal[a4[i],a4[i]]:=m16;
end for;
for i in [1..#b1] do
  scal[b1[i],b1[i]]:=n1; scal[b1[i],b2[i]]:=n2; scal[b1[i],b3[i]]:=n3; scal[b1[i],b4[i]]:=n4;
  scal[b2[i],b1[i]]:=n5; scal[b2[i],b2[i]]:=n6; scal[b2[i],b3[i]]:=n7; scal[b2[i],b4[i]]:=n8;
  scal[b3[i],b1[i]]:=n9; scal[b3[i],b2[i]]:=n10; scal[b3[i],b3[i]]:=n11; scal[b3[i],b4[i]]:=n12;
  scal[b4[i],b1[i]]:=n13; scal[b4[i],b2[i]]:=n14; scal[b4[i],b3[i]]:=n15; scal[b4[i],b4[i]]:=n16;
end for;

// But some of them are minus (or in m12's case, plus)

scal[23,25]:=-m3;
scal[22,19]:=-m4;
scal[23,14]:=-m4;
scal[8,25]:=-m7;
scal[27,10]:=-m8;
scal[25,23]:=-m9;
scal[25,8]:=-m10;
scal[9,19]:=m12;
scal[14,23]:=-m13;
scal[19,22]:=-m13;
scal[14,8]:=-m14;
scal[19,4]:=-m14;
scal[19,9]:=-m15;


scal[18,3]:=-n2;
scal[18,5]:=-n3;
scal[18,15]:=-n4;
scal[20,2]:=-n5;
scal[21,24]:=-n5;
scal[7,24]:=-n9;
scal[26,2]:=-n9;
scal[15,18]:=-n13;

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

h1:=GL(27,F)![[20,0,0,0,15,0,13,0,0,0,0,0,0,0,0,0,0,0,0,0,0,23,4,0,0,27,0],
[0,6,0,0,0,4,0,0,5,0,0,18,18,0,0,0,11,11,0,0,0,0,0,12,27,0,0],
[0,0,14,0,0,0,0,0,0,23,0,5,3,9,0,0,8,0,28,2,16,0,0,0,0,0,0],
[0,0,0,4,0,0,0,9,0,0,15,0,0,0,2,13,0,0,0,0,0,0,0,0,0,0,23],
[1,0,0,0,6,0,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,3,0,0,4,0],
[0,9,0,0,0,6,0,0,12,0,0,2,2,0,0,0,27,10,0,0,0,0,0,4,11,0,0],
[26,0,0,0,9,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,4,0,0,23,0],
[0,0,0,9,0,0,0,6,0,0,2,0,0,0,16,15,0,0,0,0,0,0,0,0,0,0,4],
[0,4,0,0,0,12,0,0,18,0,0,27,27,0,0,0,2,9,0,0,0,0,0,19,6,0,0],
[0,0,7,0,0,0,0,0,0,13,0,27,15,27,0,0,13,0,14,6,4,0,0,0,0,0,0],
[0,0,0,1,0,0,0,4,0,0,4,0,0,0,9,23,0,0,0,0,0,0,0,0,0,0,26],
[0,26,28,0,0,2,0,0,27,27,0,14,26,13,0,0,7,3,15,11,10,0,0,3,27,0,0],
[0,26,11,0,0,2,0,0,27,15,0,26,22,2,0,0,15,3,16,19,1,0,0,3,27,0,0],
[0,0,4,0,0,0,0,0,0,27,0,13,2,14,0,0,15,0,16,22,6,0,0,0,0,0,0],
[0,0,0,4,0,0,0,3,0,0,9,0,0,0,6,4,0,0,0,0,0,0,0,0,0,0,1],
[0,0,0,26,0,0,0,1,0,0,23,0,0,0,4,20,0,0,0,0,0,0,0,0,0,0,25],
[0,3,10,0,0,27,0,0,2,13,0,7,15,15,0,0,26,26,2,1,11,0,0,26,2,0,0],
[0,11,0,0,0,27,0,0,4,0,0,11,11,0,0,0,18,17,0,0,0,0,0,6,24,0,0],
[0,0,6,0,0,0,0,0,0,14,0,15,16,16,0,0,2,0,2,25,22,0,0,0,0,0,0],
[0,0,2,0,0,0,0,0,0,28,0,3,21,6,0,0,24,0,20,13,15,0,0,0,0,0,0],
[0,0,16,0,0,0,0,0,0,9,0,8,24,28,0,0,3,0,6,15,2,0,0,0,0,0,0],
[23,0,0,0,2,0,15,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,9,0,0,13,0],
[4,0,0,0,16,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,6,0,0,15,0],
[0,12,0,0,0,5,0,0,2,0,0,11,11,0,0,0,18,6,0,0,0,0,0,18,4,0,0],
[0,10,0,0,0,11,0,0,6,0,0,27,27,0,0,0,2,25,0,0,0,0,0,9,17,0,0],
[25,0,0,0,4,0,23,0,0,0,0,0,0,0,0,0,0,0,0,0,0,26,1,0,0,20,0],
[0,0,0,23,0,0,0,4,0,0,13,0,0,0,15,27,0,0,0,0,0,0,0,0,0,0,20]];

/* This is how h1 was constructed.

c1:=[27,4,8,16,15,11];
c2:=[1,22,23,26,5,7];

cc2:=Matrix([[20,23,4,27,15,13],
[23,4,9,13,2,15],
[4,9,6,15,16,2],
[25,26,1,20,4,23],
[1,4,3,4,6,9],
[26,1,4,23,9,4]]);

d1:=[6,9,25,2,18,24,10,19,14,20,3,21,12,13,17];

dd2:=Matrix([[6,12,11,9,10,4,0,0,0,0,0,0,2,2,27],
[12,18,6,4,9,19,0,0,0,0,0,0,27,27,2],
[11,6,17,10,25,9,0,0,0,0,0,0,27,27,2],
[4,5,27,6,11,12,0,0,0,0,0,0,18,18,11],
[27,4,24,11,17,6,0,0,0,0,0,0,11,11,18],
[5,2,4,12,6,18,0,0,0,0,0,0,11,11,18],
[0,0,0,0,0,0,13,14,27,6,7,4,27,15,13],
[0,0,0,0,0,0,14,2,16,25,6,22,15,16,2],
[0,0,0,0,0,0,27,16,14,22,4,6,13,2,15],
[0,0,0,0,0,0,28,20,6,13,2,15,3,21,24],
[0,0,0,0,0,0,23,28,9,2,14,16,5,3,8],
[0,0,0,0,0,0,9,6,28,15,16,2,8,24,3],
[2,27,27,26,3,3,27,15,13,11,28,10,14,26,7],
[2,27,27,26,3,3,15,16,2,19,11,1,26,22,15],
[27,2,2,3,26,26,13,2,15,1,10,11,7,15,26]]);

h1:=ZeroMatrix(F,27);
for i in [1..6] do for j in [1..6] do
  h1[c1[i],c1[j]]:=cc2[i,j];
  h1[c2[i],c2[j]]:=cc2[i,j];
end for; end for;

for i in [1..15] do for j in [1..15] do
  h1[d1[i],d1[j]]:=dd2[i,j];
end for; end for;

// The square of a,b,c is this.
//[ 1 10  9]
//[28  9  9]
//[ 0 20 19]

*/

// A few commands for substituting in varibles and simplify Groebner bases.
coeffs1:=[R.i:i in [1..50]];

function ChangeCoefficient(coeffs,coeff,target)
coeffs2:=coeffs;
nn:=Position(coeffs1,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(coeffs1,coeff[i]);
  coeffs2[nn]:=target[i];
end for;
return [Evaluate(i,coeffs2):i in coeffs];
end function;

function ChangeLinearCoefficients(coeffs,rels)

relsd:=[i:i in rels|Degree(LeadingTerm(i)) eq 1];
return ChangeCoefficients(coeffs,[LeadingTerm(i):i in relsd],[-i+LeadingTerm(i):i in relsd]);

end function;

function SimplifyBasis(B,X)

if B eq [] or B eq [1] then return B; end if;
Bfact:={Factorization(i):i in B|i ne 0};
Bfact2:={[i[1]:i in j|not(i[1] in X)]:j in Bfact};
// If any of the polynomials are now 1 then the basis dies: 
if(Minimum([#i:i in Bfact2]) eq 0) then return [1]; end if;
B2:=[&*i:i in Bfact2];
return B2;
end function;

function StripBasis(B,X)
if B eq [] then return []; end if;
repeat B:=SimplifyBasis(B,X); B:=GroebnerBasis(B); until #{{j[2]:j in Factorization(i)}:i in B} eq 1;
return B;
end function;


// Commands for generating the relations.

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

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;

function CentRel(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); v2:=Matricise(v); w2:=Matricise(w);
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;

// Here is a reduced version of the relations. We can check it is correct with the code below.

relsR:=[
    n1*n13^2 + 28*n2*n14^2 + n2*n15*n16 + n3*n14*n16 + n3*n15^2 + n4*n14*n15 + n4*n16^2,
    n1*n9*n13 + n2*n10*n14 + 14*n2*n11*n16 + 14*n2*n12*n15 + 14*n3*n10*n16 + 28*n3*n11*n15 + 14*n3*n12*n14 + 14*n4*n10*n15 + 14*n4*n11*n14 + 28*n4*n12*n16,
    n1*n9^2 + 28*n2*n10^2 + n2*n11*n12 + n3*n10*n12 + n3*n11^2 + n4*n10*n11 + n4*n12^2,
    m13*n13*n14 + 28*m15*n13*n16 + m16*n13*n15 + 20*n5*n13^2 + 20*n6*n14^2 + 9*n6*n15*n16 + 9*n7*n14*n16 + 9*n7*n15^2 + 9*n8*n14*n15 + 9*n8*n16^2,
    m13*n9*n14 + 28*m13*n10*n13 + 28*m15*n9*n16 + m15*n12*n13 + m16*n9*n15 + 28*m16*n11*n13 + 11*n5*n9*n13 + 18*n6*n10*n14 + 20*n6*n11*n16 + 20*n6*n12*n15 + 20*n7*n10*n16 + 11*n7*n11*n15 + 20*n7*n12*n14 + 20*n8*n10*n15 + 20*n8*n11*n14 + 11*n8*n12*n16,
    m13*n9*n10 + 28*m15*n9*n12 + m16*n9*n11 + 9*n5*n9^2 + 9*n6*n10^2 + 20*n6*n11*n12 + 20*n7*n10*n12 + 20*n7*n11^2 + 20*n8*n10*n11 + 20*n8*n12^2,
    m13*n1*n14 + m13*n2*n13 + 28*m15*n1*n16 + 28*m15*n4*n13 + m16*n1*n15 + m16*n3*n13 + 26*n2*n7*n16 + 3*n2*n8*n15 + 3*n3*n6*n16 + 26*n3*n8*n14 + 26*n4*n6*n15 + 3*n4*n7*n14,
    m13*n1*n10 + 28*m13*n2*n9 + 28*m15*n1*n12 + m15*n4*n9 + m16*n1*n11 + 28*m16*n3*n9 + 26*n2*n7*n12 + 3*n2*n8*n11 + 3*n3*n6*n12 + 26*n3*n8*n10 + 26*n4*n6*n11 + 3*n4*n7*n10,
    m13*n1*n6 + 28*m13*n2*n5 + 28*m15*n1*n8 + m15*n4*n5 + m16*n1*n7 + 28*m16*n3*n5 + 9*n1*n5^2 + 20*n2*n6^2 + 9*n2*n7*n8 + 9*n3*n6*n8 + 9*n3*n7^2 + 9*n4*n6*n7 + 9*n4*n8^2,
    m13*n1*n2 + 28*m15*n1*n4 + m16*n1*n3 + 9*n1^2*n5 + 9*n2^2*n6 + 20*n2*n3*n8 + 20*n2*n4*n7 + 20*n3^2*n7 + 20*n3*n4*n6 + 20*n4^2*n8,
    m13*m14*n16 + m14*m15*n15 + m14*m16*n14 + 6*n5^2*n13 + 23*n6^2*n14 + 6*n6*n7*n16 + 6*n6*n8*n15 + 6*n7^2*n15 + 6*n7*n8*n14 + 6*n8^2*n16,
    m13*m14*n12 + m14*m15*n11 + m14*m16*n10 + 23*n5^2*n9 + 23*n6^2*n10 + 6*n6*n7*n12 + 6*n6*n8*n11 + 6*n7^2*n11 + 6*n7*n8*n10 + 6*n8^2*n12,
    m13*m14*n8 + 20*m13*n5*n6 + m14*m15*n7 + m14*m16*n6 + 9*m15*n5*n8 + 20*m16*n5*n7,
    m13*m14*n4 + m14*m15*n3 + m14*m16*n2 + 23*n1*n5^2 + 6*n2*n6^2 + 23*n2*n7*n8 + 23*n3*n6*n8 + 23*n3*n7^2 + 23*n4*n6*n7 + 23*n4*n8^2,
    m13^3 + 3*m13*m15*m16 + 28*m14^3 + m15^3 + 28*m16^3 + 25*n5^3 + 25*n6^3 + 12*n6*n7*n8 + 4*n7^3 + 4*n8^3,
    m9*n13*n14 + 28*m11*n13*n16 + 28*m12*n13*n15,
    m9*n9*n14 + 28*m9*n10*n13 + 28*m11*n9*n16 + m11*n12*n13 + 28*m12*n9*n15 + m12*n11*n13,
    m9*n9*n10 + 28*m11*n9*n12 + 28*m12*n9*n11,
    m9*n5*n14 + 28*m9*n6*n13 + 28*m11*n5*n16 + m11*n8*n13 + 28*m12*n5*n15 + m12*n7*n13 + 5*n1*n5*n13 + 5*n2*n6*n14 + 22*n2*n7*n16 + 2*n2*n8*n15 + 2*n3*n6*n16 + 24*n3*n7*n15 + 22*n3*n8*n14 + 22*n4*n6*n15 + 2*n4*n7*n14 + 24*n4*n8*n16,
    m9*n5*n10 + m9*n6*n9 + 28*m11*n5*n12 + 28*m11*n8*n9 + 28*m12*n5*n11 + 28*m12*n7*n9 + 24*n1*n5*n9 + 5*n2*n6*n10 + 22*n2*n7*n12 + 2*n2*n8*n11 + 2*n3*n6*n12 + 24*n3*n7*n11 + 22*n3*n8*n10 + 22*n4*n6*n11 + 2*n4*n7*n10 + 24*n4*n8*n12,
    m9*n5*n6 + 28*m11*n5*n8 + 28*m12*n5*n7 + 13*n1*n5^2 + 16*n2*n6^2 + 13*n2*n7*n8 + 13*n3*n6*n8 + 13*n3*n7^2 + 13*n4*n6*n7 + 13*n4*n8^2,
    m9*n1*n6 + 28*m9*n2*n5 + 28*m11*n1*n8 + m11*n4*n5 + 28*m12*n1*n7 + m12*n3*n5 + 21*n1^2*n5 + 21*n2^2*n6 + 8*n2*n3*n8 + 8*n2*n4*n7 + 8*n3^2*n7 + 8*n3*n4*n6 + 8*n4^2*n8,
    m9*n1*n2 + 28*m11*n1*n4 + 28*m12*n1*n3 + 21*n1^3 + 8*n2^3 + 5*n2*n3*n4 + 21*n3^3 + 21*n4^3,
    m9*m14*n16 + m10*m13*n16 + m10*m15*n15 + m10*m16*n14 + m11*m14*n15 + 28*m12*m14*n14 + 15*n1*n5*n13 + 15*n2*n6*n14 + 14*n2*n7*n16 + 14*n3*n7*n15 + 14*n3*n8*n14 + 14*n4*n6*n15 + 14*n4*n8*n16,
    m9*m14*n12 + m10*m13*n12 + m10*m15*n11 + m10*m16*n10 + m11*m14*n11 + 28*m12*m14*n10 + 14*n1*n5*n9 + 15*n2*n6*n10 + 14*n2*n7*n12 + 14*n3*n7*n11 + 14*n3*n8*n10 + 14*n4*n6*n11 + 14*n4*n8*n12,
    m9*m14*n8 + m10*m13*n8 + m10*m15*n7 + m10*m16*n6 + m11*m14*n7 + 28*m12*m14*n6 + 28*n1*n5^2 + n2*n6^2 + 28*n2*n7*n8 + 28*n3*n6*n8 + 28*n3*n7^2 + 28*n4*n6*n7 + 28*n4*n8^2,
    m9*m14*n4 + m10*m13*n4 + m10*m15*n3 + m10*m16*n2 + m11*m14*n3 + 28*m12*m14*n2 + 14*n1^2*n5 + 14*n2^2*n6 + 15*n2*n3*n8 + 15*n2*n4*n7 + 15*n3^2*n7 + 15*n3*n4*n6 + 15*n4^2*n8,
    m9*m13^2 + m9*m15*m16 + 28*m10*m14^2 + m11*m13*m16 + m11*m15^2 + 28*m12*m13*m15 + m12*m16^2 + 9*n1*n5^2 + 20*n2*n6^2 + 9*n2*n7*n8 + 9*n3*n6*n8 + 9*n3*n7^2 + 9*n4*n6*n7 + 9*n4*n8^2,
    m9*m10*n16 + 12*m9*n1*n14 + 12*m9*n2*n13 + m10*m11*n15 + 28*m10*m12*n14 + 17*m11*n1*n16 + 17*m11*n4*n13 + 17*m12*n1*n15 + 17*m12*n3*n13 + 12*n1^2*n13 + 17*n2^2*n14 + 12*n2*n3*n16 + 12*n2*n4*n15 + 12*n3^2*n15 + 12*n3*n4*n14 + 12*n4^2*n16,
    m9*m10*n12 + 12*m9*n1*n10 + 17*m9*n2*n9 + m10*m11*n11 + 28*m10*m12*n10 + 17*m11*n1*n12 + 12*m11*n4*n9 + 17*m12*n1*n11 + 12*m12*n3*n9 + 17*n1^2*n9 + 17*n2^2*n10 + 12*n2*n3*n12 + 12*n2*n4*n11 + 12*n3^2*n11 + 12*n3*n4*n10 + 12*n4^2*n12,
    m9*m10*n8 + m10*m11*n7 + 28*m10*m12*n6 + 5*n1^2*n5 + 5*n2^2*n6 + 24*n2*n3*n8 + 24*n2*n4*n7 + 24*n3^2*n7 + 24*n3*n4*n6 + 24*n4^2*n8,
    m9*m10*n4 + m10*m11*n3 + 28*m10*m12*n2 + 17*n1^3 + 12*n2^3 + 22*n2*n3*n4 + 17*n3^3 + 17*n4^3,
    m9^2*m13 + m9*m11*m16 + 28*m9*m12*m15 + 28*m10^2*m14 + m11^2*m15 + 28*m11*m12*m13 + 28*m12^2*m16 + 16*n1^2*n5 + 16*n2^2*n6 + 13*n2*n3*n8 + 13*n2*n4*n7 + 13*n3^2*n7 + 13*n3*n4*n6 + 13*n4^2*n8,
    m9^3 + 26*m9*m11*m12 + 28*m10^3 + m11^3 + m12^3 + 22*n1^3 + 7*n2^3 + 8*n2*n3*n4 + 22*n3^3 + 22*n4^3,
    m5*n5*n14 + 28*m5*n6*n13 + 28*m7*n5*n16 + m7*n8*n13 + 28*m8*n5*n15 + m8*n7*n13 + 14*n5*n13^2 + 14*n6*n14^2 + 15*n6*n15*n16 + 15*n7*n14*n16 + 15*n7*n15^2 + 15*n8*n14*n15 + 15*n8*n16^2,
    m5*n5*n10 + m5*n6*n9 + 28*m7*n5*n12 + 28*m7*n8*n9 + 28*m8*n5*n11 + 28*m8*n7*n9 + 15*n5*n9*n13 + 14*n6*n10*n14 + 15*n6*n11*n16 + 15*n7*n11*n15 + 15*n7*n12*n14 + 15*n8*n10*n15 + 15*n8*n12*n16,
    m5*n5*n6 + 28*m7*n5*n8 + 28*m8*n5*n7 + 8*m13*n5*n14 + 21*m13*n6*n13 + 21*m15*n5*n16 + 8*m15*n8*n13 + 8*m16*n5*n15 + 21*m16*n7*n13 + 14*n5^2*n13 + 15*n6^2*n14 + 14*n6*n7*n16 + 14*n6*n8*n15 + 14*n7^2*n15 + 14*n7*n8*n14 + 14*n8^2*n16,
    m5*n1*n14 + m5*n2*n13 + 28*m7*n1*n16 + 28*m7*n4*n13 + 28*m8*n1*n15 + 28*m8*n3*n13,
    m5*n1*n10 + 28*m5*n2*n9 + 28*m7*n1*n12 + m7*n4*n9 + 28*m8*n1*n11 + m8*n3*n9 + 15*n2*n11*n16 + 14*n2*n12*n15 + 14*n3*n10*n16 + 15*n3*n12*n14 + 15*n4*n10*n15 + 14*n4*n11*n14,
    m5*n1*n6 + 28*m5*n2*n5 + 28*m7*n1*n8 + m7*n4*n5 + 28*m8*n1*n7 + m8*n3*n5 + 15*n1*n5*n13 + 15*n2*n6*n14 + 19*n2*n7*n16 + 24*n2*n8*n15 + 24*n3*n6*n16 + 14*n3*n7*n15 + 19*n3*n8*n14 + 19*n4*n6*n15 + 24*n4*n7*n14 + 14*n4*n8*n16,
    m5*n1*n2 + 28*m7*n1*n4 + 28*m8*n1*n3 + 10*m9*n1*n14 + 10*m9*n2*n13 + 19*m11*n1*n16 + 19*m11*n4*n13 + 19*m12*n1*n15 + 19*m12*n3*n13 + 15*n1^2*n13 + 14*n2^2*n14 + 15*n2*n3*n16 + 15*n2*n4*n15 + 15*n3^2*n15 + 15*n3*n4*n14 + 15*n4^2*n16,
    m5*m14*n16 + m6*m13*n16 + m6*m15*n15 + m6*m16*n14 + m7*m14*n15 + 28*m8*m14*n14 + 10*n5*n13^2 + 10*n6*n14^2 + 19*n6*n15*n16 + 19*n7*n14*n16 + 19*n7*n15^2 + 19*n8*n14*n15 + 19*n8*n16^2,
    m5*m14*n12 + m6*m13*n12 + m6*m15*n11 + m6*m16*n10 + m7*m14*n11 + 28*m8*m14*n10 + 19*n5*n9*n13 + 10*n6*n10*n14 + 19*n6*n11*n16 + 19*n7*n11*n15 + 19*n7*n12*n14 + 19*n8*n10*n15 + 19*n8*n12*n16,
    m5*m14*n8 + m6*m13*n8 + m6*m15*n7 + m6*m16*n6 + m7*m14*n7 + 28*m8*m14*n6 + 15*m13*n5*n14 + 14*m13*n6*n13 + 14*m15*n5*n16 + 15*m15*n8*n13 + 15*m16*n5*n15 + 14*m16*n7*n13 + 10*n5^2*n13 + 19*n6^2*n14 + 10*n6*n7*n16 + 10*n6*n8*n15 + 10*n7^2*n15 + 10*n7*n8*n14 + 10*n8^2*n16,
    m5*m14*n4 + m6*m13*n4 + m6*m15*n3 + m6*m16*n2 + m7*m14*n3 + 28*m8*m14*n2 + 19*n1*n5*n13 + 19*n2*n6*n14 + 3*n2*n7*n16 + 7*n2*n8*n15 + 7*n3*n6*n16 + 10*n3*n7*n15 + 3*n3*n8*n14 + 3*n4*n6*n15 + 7*n4*n7*n14 + 10*n4*n8*n16,
    m5*m13^2 + m5*m15*m16 + 28*m6*m14^2 + m7*m13*m16 + m7*m15^2 + 28*m8*m13*m15 + m8*m16^2 + 19*m13*n5*n14 + 10*m13*n6*n13 + 10*m15*n5*n16 + 19*m15*n8*n13 + 19*m16*n5*n15 + 10*m16*n7*n13 + 26*n5^2*n13 + 3*n6^2*n14 + 26*n6*n7*n16 + 26*n6*n8*n15 + 26*n7^2*n15 + 26*n7*n8*n14 + 26*n8^2*n16,
    m5*m11*m16 + m5*m12*m15 + 28*m7*m9*m16 + 28*m7*m12*m13 + 28*m8*m9*m15 + m8*m11*m13 + 13*n1*n5*n13 + 13*n2*n6*n14 + 15*n2*n7*n16 + n2*n8*n15 + n3*n6*n16 + 16*n3*n7*n15 + 15*n3*n8*n14 + 15*n4*n6*n15 + n4*n7*n14 + 16*n4*n8*n16,
    m5*m10*n16 + m6*m9*n16 + m6*m11*n15 + 28*m6*m12*n14 + m7*m10*n15 + 28*m8*m10*n14,
    m5*m10*n12 + m6*m9*n12 + m6*m11*n11 + 28*m6*m12*n10 + m7*m10*n11 + 28*m8*m10*n10 + 8*n2*n11*n16 + 21*n2*n12*n15 + 21*n3*n10*n16 + 8*n3*n12*n14 + 8*n4*n10*n15 + 21*n4*n11*n14,
    m5*m10*n8 + m6*m9*n8 + m6*m11*n7 + 28*m6*m12*n6 + m7*m10*n7 + 28*m8*m10*n6 + 21*n1*n5*n13 + 21*n2*n6*n14 + 22*n2*n7*n16 + 15*n2*n8*n15 + 15*n3*n6*n16 + 8*n3*n7*n15 + 22*n3*n8*n14 + 22*n4*n6*n15 + 15*n4*n7*n14 + 8*n4*n8*n16,
    m5*m10*n4 + m6*m9*n4 + m6*m11*n3 + 28*m6*m12*n2 + m7*m10*n3 + 28*m8*m10*n2 + 14*m9*n1*n14 + 14*m9*n2*n13 + 15*m11*n1*n16 + 15*m11*n4*n13 + 15*m12*n1*n15 + 15*m12*n3*n13 + 13*n1^2*n13 + 16*n2^2*n14 + 13*n2*n3*n16 + 13*n2*n4*n15 + 13*n3^2*n15 + 13*n3*n4*n14 + 13*n4^2*n16,
    m5*m9*m13 + 28*m5*m12*m15 + 28*m6*m10*m14 + m7*m9*m16 + m7*m11*m15 + 28*m8*m11*m13 + 28*m8*m12*m16 + n1*n5*n13 + n2*n6*n14 + 9*n2*n7*n16 + 19*n2*n8*n15 + 19*n3*n6*n16 + 28*n3*n7*n15 + 9*n3*n8*n14 + 9*n4*n6*n15 + 19*n4*n7*n14 + 28*n4*n8*n16,
    m5*m9^2 + 28*m5*m11*m12 + 28*m6*m10^2 + 28*m7*m9*m12 + m7*m11^2 + 28*m8*m9*m11 + m8*m12^2 + 25*m9*n1*n14 + 25*m9*n2*n13 + 4*m11*n1*n16 + 4*m11*n4*n13 + 4*m12*n1*n15 + 4*m12*n3*n13 + 12*n1^2*n13 + 17*n2^2*n14 + 12*n2*n3*n16 + 12*n2*n4*n15 + 12*n3^2*n15 + 12*n3*n4*n14 + 12*n4^2*n16,
    m5*m6*n16 + 15*m5*n13*n14 + m6*m7*n15 + 28*m6*m8*n14 + 14*m7*n13*n16 + 14*m8*n13*n15 + 22*n13^3 + 7*n14^3 + 8*n14*n15*n16 + 22*n15^3 + 22*n16^3,
    m5*m6*n12 + 7*m5*n9*n14 + 22*m5*n10*n13 + m6*m7*n11 + 28*m6*m8*n10 + 22*m7*n9*n16 + 7*m7*n12*n13 + 22*m8*n9*n15 + 7*m8*n11*n13 + 7*n9*n13^2 + 7*n10*n14^2 + 22*n10*n15*n16 + 22*n11*n14*n16 + 22*n11*n15^2 + 22*n12*n14*n15 + 22*n12*n16^2,
    m5*m6*n8 + m6*m7*n7 + 28*m6*m8*n6 + 15*n5*n13^2 + 15*n6*n14^2 + 14*n6*n15*n16 + 14*n7*n14*n16 + 14*n7*n15^2 + 14*n8*n14*n15 + 14*n8*n16^2,
    m5*m6*n4 + m6*m7*n3 + 28*m6*m8*n2,
    m5^2*m13 + m5*m7*m16 + 28*m5*m8*m15 + 28*m6^2*m14 + m7^2*m15 + 28*m7*m8*m13 + 28*m8^2*m16 + 10*n5*n13^2 + 10*n6*n14^2 + 19*n6*n15*n16 + 19*n7*n14*n16 + 19*n7*n15^2 + 19*n8*n14*n15 + 19*n8*n16^2,
    m5^2*m9 + 28*m5*m7*m12 + 28*m5*m8*m11 + 28*m6^2*m10 + m7^2*m11 + 28*m7*m8*m9 + m8^2*m12,
    m5^3 + 26*m5*m7*m8 + 15*m5*n13*n14 + 28*m6^3 + m7^3 + 14*m7*n13*n16 + m8^3 + 14*m8*n13*n15 + 7*n13^3 + 22*n14^3 + 21*n14*n15*n16 + 7*n15^3 + 7*n16^3,
    m1*n13*n14 + 28*m3*n13*n16 + m4*n13*n15 + 15*m5*n9*n14 + 14*m5*n10*n13 + 14*m7*n9*n16 + 15*m7*n12*n13 + 14*m8*n9*n15 + 15*m8*n11*n13,
    m1*n9*n14 + 28*m1*n10*n13 + 28*m3*n9*n16 + m3*n12*n13 + m4*n9*n15 + 28*m4*n11*n13 + 27*m5*n9*n10 + 2*m7*n9*n12 + 2*m8*n9*n11,
    m1*n5*n14 + 28*m1*n6*n13 + 28*m3*n5*n16 + m3*n8*n13 + m4*n5*n15 + 28*m4*n7*n13 + 15*n5*n9*n13 + 14*n6*n10*n14 + 15*n6*n12*n15 + 15*n7*n10*n16 + 15*n7*n11*n15 + 15*n8*n11*n14 + 15*n8*n12*n16,
    m1*n5*n10 + m1*n6*n9 + 28*m3*n5*n12 + 28*m3*n8*n9 + m4*n5*n11 + m4*n7*n9 + 14*n5*n9^2 + 14*n6*n10^2 + 15*n6*n11*n12 + 15*n7*n10*n12 + 15*n7*n11^2 + 15*n8*n10*n11 + 15*n8*n12^2,
    m1*n5*n6 + 28*m3*n5*n8 + m4*n5*n7 + 8*m13*n5*n10 + 8*m13*n6*n9 + 21*m15*n5*n12 + 21*m15*n8*n9 + 8*m16*n5*n11 + 8*m16*n7*n9 + 15*n5^2*n9 + 15*n6^2*n10 + 14*n6*n7*n12 + 14*n6*n8*n11 + 14*n7^2*n11 + 14*n7*n8*n10 + 14*n8^2*n12,
    m1*n1*n14 + m1*n2*n13 + 28*m3*n1*n16 + 28*m3*n4*n13 + m4*n1*n15 + m4*n3*n13 + 14*n2*n11*n16 + 15*n2*n12*n15 + 15*n3*n10*n16 + 14*n3*n12*n14 + 14*n4*n10*n15 + 15*n4*n11*n14,
    m1*n1*n10 + 28*m1*n2*n9 + 28*m3*n1*n12 + m3*n4*n9 + m4*n1*n11 + 28*m4*n3*n9,
    m1*n1*n6 + 28*m1*n2*n5 + 28*m3*n1*n8 + m3*n4*n5 + m4*n1*n7 + 28*m4*n3*n5 + 14*n1*n5*n9 + 15*n2*n6*n10 + 19*n2*n7*n12 + 24*n2*n8*n11 + 24*n3*n6*n12 + 14*n3*n7*n11 + 19*n3*n8*n10 + 19*n4*n6*n11 + 24*n4*n7*n10 + 14*n4*n8*n12,
    m1*n1*n2 + 28*m3*n1*n4 + m4*n1*n3 + 10*m9*n1*n10 + 19*m9*n2*n9 + 19*m11*n1*n12 + 10*m11*n4*n9 + 19*m12*n1*n11 + 10*m12*n3*n9 + 14*n1^2*n9 + 14*n2^2*n10 + 15*n2*n3*n12 + 15*n2*n4*n11 + 15*n3^2*n11 + 15*n3*n4*n10 + 15*n4^2*n12,
    m1*m14*n16 + m2*m13*n16 + m2*m15*n15 + m2*m16*n14 + m3*m14*n15 + m4*m14*n14 + 19*n5*n9*n13 + 10*n6*n10*n14 + 19*n6*n12*n15 + 19*n7*n10*n16 + 19*n7*n11*n15 + 19*n8*n11*n14 + 19*n8*n12*n16,
    m1*m14*n12 + m2*m13*n12 + m2*m15*n11 + m2*m16*n10 + m3*m14*n11 + m4*m14*n10 + 10*n5*n9^2 + 10*n6*n10^2 + 19*n6*n11*n12 + 19*n7*n10*n12 + 19*n7*n11^2 + 19*n8*n10*n11 + 19*n8*n12^2,
    m1*m14*n8 + m2*m13*n8 + m2*m15*n7 + m2*m16*n6 + m3*m14*n7 + m4*m14*n6 + 15*m13*n5*n10 + 15*m13*n6*n9 + 14*m15*n5*n12 + 14*m15*n8*n9 + 15*m16*n5*n11 + 15*m16*n7*n9 + 19*n5^2*n9 + 19*n6^2*n10 + 10*n6*n7*n12 + 10*n6*n8*n11 + 10*n7^2*n11 + 10*n7*n8*n10 + 10*n8^2*n12,
    m1*m14*n4 + m2*m13*n4 + m2*m15*n3 + m2*m16*n2 + m3*m14*n3 + m4*m14*n2 + 10*n1*n5*n9 + 19*n2*n6*n10 + 3*n2*n7*n12 + 7*n2*n8*n11 + 7*n3*n6*n12 + 10*n3*n7*n11 + 3*n3*n8*n10 + 3*n4*n6*n11 + 7*n4*n7*n10 + 10*n4*n8*n12,
    m1*m13^2 + m1*m15*m16 + 28*m2*m14^2 + m3*m13*m16 + m3*m15^2 + m4*m13*m15 + 28*m4*m16^2 + 19*m13*n5*n10 + 19*m13*n6*n9 + 10*m15*n5*n12 + 10*m15*n8*n9 + 19*m16*n5*n11 + 19*m16*n7*n9 + 3*n5^2*n9 + 3*n6^2*n10 + 26*n6*n7*n12 + 26*n6*n8*n11 + 26*n7^2*n11 + 26*n7*n8*n10 + 26*n8^2*n12,
    m1*m11*m16 + m1*m12*m15 + 28*m3*m9*m16 + 28*m3*m12*m13 + m4*m9*m15 + 28*m4*m11*m13 + 16*n1*n5*n9 + 13*n2*n6*n10 + 15*n2*n7*n12 + n2*n8*n11 + n3*n6*n12 + 16*n3*n7*n11 + 15*n3*n8*n10 + 15*n4*n6*n11 + n4*n7*n10 + 16*n4*n8*n12,
    m1*m10*n16 + m2*m9*n16 + m2*m11*n15 + 28*m2*m12*n14 + m3*m10*n15 + m4*m10*n14 + 21*n2*n11*n16 + 8*n2*n12*n15 + 8*n3*n10*n16 + 21*n3*n12*n14 + 21*n4*n10*n15 + 8*n4*n11*n14,
    m1*m10*n12 + m2*m9*n12 + m2*m11*n11 + 28*m2*m12*n10 + m3*m10*n11 + m4*m10*n10,
    m1*m10*n8 + m2*m9*n8 + m2*m11*n7 + 28*m2*m12*n6 + m3*m10*n7 + m4*m10*n6 + 8*n1*n5*n9 + 21*n2*n6*n10 + 22*n2*n7*n12 + 15*n2*n8*n11 + 15*n3*n6*n12 + 8*n3*n7*n11 + 22*n3*n8*n10 + 22*n4*n6*n11 + 15*n4*n7*n10 + 8*n4*n8*n12,
    m1*m10*n4 + m2*m9*n4 + m2*m11*n3 + 28*m2*m12*n2 + m3*m10*n3 + m4*m10*n2 + 14*m9*n1*n10 + 15*m9*n2*n9 + 15*m11*n1*n12 + 14*m11*n4*n9 + 15*m12*n1*n11 + 14*m12*n3*n9 + 16*n1^2*n9 + 16*n2^2*n10 + 13*n2*n3*n12 + 13*n2*n4*n11 + 13*n3^2*n11 + 13*n3*n4*n10 + 13*n4^2*n12,
    m1*m9*m13 + 28*m1*m12*m15 + 28*m2*m10*m14 + m3*m9*m16 + m3*m11*m15 + m4*m11*m13 + m4*m12*m16 + 28*n1*n5*n9 + n2*n6*n10 + 9*n2*n7*n12 + 19*n2*n8*n11 + 19*n3*n6*n12 + 28*n3*n7*n11 + 9*n3*n8*n10 + 9*n4*n6*n11 + 19*n4*n7*n10 + 28*n4*n8*n12,
    m1*m9^2 + 28*m1*m11*m12 + 28*m2*m10^2 + 28*m3*m9*m12 + m3*m11^2 + m4*m9*m11 + 28*m4*m12^2 + 25*m9*n1*n10 + 4*m9*n2*n9 + 4*m11*n1*n12 + 25*m11*n4*n9 + 4*m12*n1*n11 + 25*m12*n3*n9 + 17*n1^2*n9 + 17*n2^2*n10 + 12*n2*n3*n12 + 12*n2*n4*n11 + 12*n3^2*n11 + 12*n3*n4*n10 + 12*n4^2*n12,
    m1*m7*m16 + m1*m8*m15 + 28*m3*m5*m16 + 28*m3*m8*m13 + m4*m5*m15 + 28*m4*m7*m13 + 19*n6*n11*n16 + 10*n6*n12*n15 + 10*n7*n10*n16 + 19*n7*n12*n14 + 19*n8*n10*n15 + 10*n8*n11*n14,
    m1*m7*m12 + 28*m1*m8*m11 + 28*m3*m5*m12 + m3*m8*m9 + 28*m4*m5*m11 + m4*m7*m9 + 8*n2*n11*n16 + 21*n2*n12*n15 + 21*n3*n10*n16 + 8*n3*n12*n14 + 8*n4*n10*n15 + 21*n4*n11*n14,
    m1*m6*n16 + m2*m5*n16 + m2*m7*n15 + 28*m2*m8*n14 + m3*m6*n15 + m4*m6*n14 + 14*m5*n9*n14 + 15*m5*n10*n13 + 15*m7*n9*n16 + 14*m7*n12*n13 + 15*m8*n9*n15 + 14*m8*n11*n13 + 14*n9*n13^2 + 14*n10*n14^2 + 15*n10*n15*n16 + 15*n11*n14*n16 + 15*n11*n15^2 + 15*n12*n14*n15 + 15*n12*n16^2,
    m1*m6*n12 + m2*m5*n12 + m2*m7*n11 + 28*m2*m8*n10 + m3*m6*n11 + m4*m6*n10 + 28*m5*n9*n10 + m7*n9*n12 + m8*n9*n11 + 15*n9^2*n13 + 14*n10^2*n14 + 15*n10*n11*n16 + 15*n10*n12*n15 + 15*n11^2*n15 + 15*n11*n12*n14 + 15*n12^2*n16,
    m1*m6*n8 + m2*m5*n8 + m2*m7*n7 + 28*m2*m8*n6 + m3*m6*n7 + m4*m6*n6 + 28*n5*n9*n13 + n6*n10*n14 + 14*n6*n11*n16 + 14*n6*n12*n15 + 14*n7*n10*n16 + 28*n7*n11*n15 + 14*n7*n12*n14 + 14*n8*n10*n15 + 14*n8*n11*n14 + 28*n8*n12*n16,
    m1*m6*n4 + m2*m5*n4 + m2*m7*n3 + 28*m2*m8*n2 + m3*m6*n3 + m4*m6*n2,
    m1*m5*m13 + 28*m1*m8*m15 + 28*m2*m6*m14 + m3*m5*m16 + m3*m7*m15 + m4*m7*m13 + m4*m8*m16 + 19*n5*n9*n13 + 10*n6*n10*n14 + 19*n6*n12*n15 + 19*n7*n10*n16 + 19*n7*n11*n15 + 19*n8*n11*n14 + 19*n8*n12*n16,
    m1*m5*m9 + 28*m1*m8*m11 + 28*m2*m6*m10 + 28*m3*m5*m12 + m3*m7*m11 + m4*m7*m9 + 28*m4*m8*m12 + 4*n2*n11*n16 + 25*n2*n12*n15 + 25*n3*n10*n16 + 4*n3*n12*n14 + 4*n4*n10*n15 + 25*n4*n11*n14,
    m1*m5^2 + 28*m1*m7*m8 + 28*m2*m6^2 + 28*m3*m5*m8 + m3*m7^2 + m4*m5*m7 + 28*m4*m8^2 + 7*m5*n9*n14 + 22*m5*n10*n13 + 22*m7*n9*n16 + 7*m7*n12*n13 + 22*m8*n9*n15 + 7*m8*n11*n13 + 22*n9*n13^2 + 22*n10*n14^2 + 7*n10*n15*n16 + 7*n11*n14*n16 + 7*n11*n15^2 + 7*n12*n14*n15 + 7*n12*n16^2,
    m1*m2*n16 + m2*m3*n15 + m2*m4*n14 + 14*m5*n9*n10 + 15*m7*n9*n12 + 15*m8*n9*n11 + 22*n9^2*n13 + 7*n10^2*n14 + 22*n10*n11*n16 + 22*n10*n12*n15 + 22*n11^2*n15 + 22*n11*n12*n14 + 22*n12^2*n16,
    m1*m2*n12 + 14*m1*n9*n10 + m2*m3*n11 + m2*m4*n10 + 15*m3*n9*n12 + 14*m4*n9*n11 + 7*n9^3 + 7*n10^3 + 8*n10*n11*n12 + 22*n11^3 + 22*n12^3,
    m1*m2*n8 + m2*m3*n7 + m2*m4*n6 + 15*n5*n9^2 + 15*n6*n10^2 + 14*n6*n11*n12 + 14*n7*n10*n12 + 14*n7*n11^2 + 14*n8*n10*n11 + 14*n8*n12^2,
    m1*m2*n4 + m2*m3*n3 + m2*m4*n2,
    m1^2*m13 + m1*m3*m16 + m1*m4*m15 + 28*m2^2*m14 + m3^2*m15 + m3*m4*m13 + 28*m4^2*m16 + 10*n5*n9^2 + 10*n6*n10^2 + 19*n6*n11*n12 + 19*n7*n10*n12 + 19*n7*n11^2 + 19*n8*n10*n11 + 19*n8*n12^2,
    m1^2*m9 + 28*m1*m3*m12 + m1*m4*m11 + 28*m2^2*m10 + m3^2*m11 + m3*m4*m9 + m4^2*m12,
    m1^2*m5 + 28*m1*m3*m8 + m1*m4*m7 + 28*m2^2*m6 + m3^2*m7 + m3*m4*m5 + m4^2*m8 + 14*m5*n9*n10 + 15*m7*n9*n12 + 15*m8*n9*n11 + 7*n9^2*n13 + 22*n10^2*n14 + 7*n10*n11*n16 + 7*n10*n12*n15 + 7*n11^2*n15 + 7*n11*n12*n14 + 7*n12^2*n16,
    m1^3 + 3*m1*m3*m4 + 14*m1*n9*n10 + 28*m2^3 + m3^3 + 15*m3*n9*n12 + 28*m4^3 + 14*m4*n9*n11 + 22*n9^3 + 22*n10^3 + 21*n10*n11*n12 + 7*n11^3 + 7*n12^3,
    b*m13*n15 + 28*b*m15*n14 + 28*b*m16*n16 + c*m14*n13 + 10*m13*n5*n14 + 19*m13*n6*n13 + 19*m15*n5*n16 + 10*m15*n8*n13 + 10*m16*n5*n15 + 19*m16*n7*n13 + 7*n1*n5*n13 + 7*n2*n6*n14 + 17*n2*n7*n16 + 5*n2*n8*n15 + 5*n3*n6*n16 + 22*n3*n7*n15 + 17*n3*n8*n14 + 17*n4*n6*n15 + 5*n4*n7*n14 + 22*n4*n8*n16 + 2*n5^2*n13 + 27*n6^2*n14 +
        2*n6*n7*n16 + 2*n6*n8*n15 + 2*n7^2*n15 + 2*n7*n8*n14 + 2*n8^2*n16,
    b*m13*n11 + 28*b*m15*n10 + 28*b*m16*n12 + 28*c*m14*n9 + 10*m13*n5*n10 + 10*m13*n6*n9 + 19*m15*n5*n12 + 19*m15*n8*n9 + 10*m16*n5*n11 + 10*m16*n7*n9 + 22*n1*n5*n9 + 7*n2*n6*n10 + 17*n2*n7*n12 + 5*n2*n8*n11 + 5*n3*n6*n12 + 22*n3*n7*n11 + 17*n3*n8*n10 + 17*n4*n6*n11 + 5*n4*n7*n10 + 22*n4*n8*n12 + 27*n5^2*n9 + 27*n6^2*n10 +
        2*n6*n7*n12 + 2*n6*n8*n11 + 2*n7^2*n11 + 2*n7*n8*n10 + 2*n8^2*n12,
    b*m13*n7 + 28*b*m15*n6 + 28*b*m16*n8 + 28*c*m14*n5 + 5*m13*n5*n6 + 24*m15*n5*n8 + 5*m16*n5*n7 + 11*n1*n5^2 + 18*n2*n6^2 + 11*n2*n7*n8 + 11*n3*n6*n8 + 11*n3*n7^2 + 11*n4*n6*n7 + 11*n4*n8^2 + 10*n5^3 + 10*n6^3 + 28*n6*n7*n8 + 19*n7^3 + 19*n8^3,
    b*m13*n3 + 28*b*m15*n2 + 28*b*m16*n4 + c*m14*n1 + 22*n1^2*n5 + 28*n1*n5^2 + 22*n2^2*n6 + 7*n2*n3*n8 + 7*n2*n4*n7 + n2*n6^2 + 28*n2*n7*n8 + 7*n3^2*n7 + 7*n3*n4*n6 + 28*n3*n6*n8 + 28*n3*n7^2 + 7*n4^2*n8 + 28*n4*n6*n7 + 28*n4*n8^2,
    b*m9*n15 + 28*b*m11*n14 + b*m12*n16 + c*m10*n13 + 2*m9*n1*n14 + 2*m9*n2*n13 + 27*m11*n1*n16 + 27*m11*n4*n13 + 27*m12*n1*n15 + 27*m12*n3*n13 + 6*n1^2*n13 + 24*n1*n5*n13 + 23*n2^2*n14 + 6*n2*n3*n16 + 6*n2*n4*n15 + 24*n2*n6*n14 + 24*n2*n7*n16 + 10*n2*n8*n15 + 6*n3^2*n15 + 6*n3*n4*n14 + 10*n3*n6*n16 + 5*n3*n7*n15 + 24*n3*n8*n14 +
        6*n4^2*n16 + 24*n4*n6*n15 + 10*n4*n7*n14 + 5*n4*n8*n16,
    b*m9*n11 + 28*b*m11*n10 + b*m12*n12 + 28*c*m10*n9 + 2*m9*n1*n10 + 27*m9*n2*n9 + 27*m11*n1*n12 + 2*m11*n4*n9 + 27*m12*n1*n11 + 2*m12*n3*n9 + 23*n1^2*n9 + 5*n1*n5*n9 + 23*n2^2*n10 + 6*n2*n3*n12 + 6*n2*n4*n11 + 24*n2*n6*n10 + 24*n2*n7*n12 + 10*n2*n8*n11 + 6*n3^2*n11 + 6*n3*n4*n10 + 10*n3*n6*n12 + 5*n3*n7*n11 + 24*n3*n8*n10 +
        6*n4^2*n12 + 24*n4*n6*n11 + 10*n4*n7*n10 + 5*n4*n8*n12,
    b*m9*n7 + 28*b*m11*n6 + b*m12*n8 + 28*c*m10*n5 + 6*n1^2*n5 + 5*n1*n5^2 + 6*n2^2*n6 + 23*n2*n3*n8 + 23*n2*n4*n7 + 24*n2*n6^2 + 5*n2*n7*n8 + 23*n3^2*n7 + 23*n3*n4*n6 + 5*n3*n6*n8 + 5*n3*n7^2 + 23*n4^2*n8 + 5*n4*n6*n7 + 5*n4*n8^2,
    b*m9*n3 + 28*b*m11*n2 + b*m12*n4 + c*m10*n1 + 13*n1^3 + 9*n1^2*n5 + 16*n2^3 + 9*n2^2*n6 + 10*n2*n3*n4 + 20*n2*n3*n8 + 20*n2*n4*n7 + 13*n3^3 + 20*n3^2*n7 + 20*n3*n4*n6 + 13*n4^3 + 20*n4^2*n8,
    b*m5*n15 + 28*b*m7*n14 + b*m8*n16 + c*m6*n13 + 15*n5*n13^2 + 15*n6*n14^2 + 14*n6*n15*n16 + 14*n7*n14*n16 + 14*n7*n15^2 + 14*n8*n14*n15 + 14*n8*n16^2,
    b*m5*n11 + 28*b*m7*n10 + b*m8*n12 + 28*c*m6*n9 + 28*n2*n11*n16 + n2*n12*n15 + n3*n10*n16 + 28*n3*n12*n14 + 28*n4*n10*n15 + n4*n11*n14 + 14*n5*n9*n13 + 15*n6*n10*n14 + 28*n6*n11*n16 + 15*n6*n12*n15 + 15*n7*n10*n16 + 14*n7*n11*n15 + 28*n7*n12*n14 + 28*n8*n10*n15 + 15*n8*n11*n14 + 14*n8*n12*n16,
    b*m5*n7 + 28*b*m7*n6 + b*m8*n8 + 28*c*m6*n5 + 17*m13*n5*n14 + 12*m13*n6*n13 + 12*m15*n5*n16 + 17*m15*n8*n13 + 17*m16*n5*n15 + 12*m16*n7*n13 + 11*n2*n7*n16 + 18*n2*n8*n15 + 18*n3*n6*n16 + 11*n3*n8*n14 + 11*n4*n6*n15 + 18*n4*n7*n14 + 18*n5^2*n13 + 11*n6^2*n14 + 18*n6*n7*n16 + 18*n6*n8*n15 + 18*n7^2*n15 + 18*n7*n8*n14 +
        18*n8^2*n16,
    b*m5*n3 + 28*b*m7*n2 + b*m8*n4 + c*m6*n1 + 11*m9*n1*n14 + 11*m9*n2*n13 + 18*m11*n1*n16 + 18*m11*n4*n13 + 18*m12*n1*n15 + 18*m12*n3*n13 + 26*n1*n5*n13 + 26*n2*n6*n14 + 15*n2*n7*n16 + 17*n2*n8*n15 + 17*n3*n6*n16 + 3*n3*n7*n15 + 15*n3*n8*n14 + 15*n4*n6*n15 + 17*n4*n7*n14 + 3*n4*n8*n16,
    b*m1*n15 + 28*b*m3*n14 + 28*b*m4*n16 + c*m2*n13 + n2*n11*n16 + 28*n2*n12*n15 + 28*n3*n10*n16 + n3*n12*n14 + n4*n10*n15 + 28*n4*n11*n14 + 14*n5*n9*n13 + 15*n6*n10*n14 + 15*n6*n11*n16 + 28*n6*n12*n15 + 28*n7*n10*n16 + 14*n7*n11*n15 + 15*n7*n12*n14 + 15*n8*n10*n15 + 28*n8*n11*n14 + 14*n8*n12*n16,
    b*m1*n11 + 28*b*m3*n10 + 28*b*m4*n12 + 28*c*m2*n9 + 15*n5*n9^2 + 15*n6*n10^2 + 14*n6*n11*n12 + 14*n7*n10*n12 + 14*n7*n11^2 + 14*n8*n10*n11 + 14*n8*n12^2,
    b*m1*n7 + 28*b*m3*n6 + 28*b*m4*n8 + 28*c*m2*n5 + 17*m13*n5*n10 + 17*m13*n6*n9 + 12*m15*n5*n12 + 12*m15*n8*n9 + 17*m16*n5*n11 + 17*m16*n7*n9 + 11*n2*n7*n12 + 18*n2*n8*n11 + 18*n3*n6*n12 + 11*n3*n8*n10 + 11*n4*n6*n11 + 18*n4*n7*n10 + 11*n5^2*n9 + 11*n6^2*n10 + 18*n6*n7*n12 + 18*n6*n8*n11 + 18*n7^2*n11 + 18*n7*n8*n10 +
        18*n8^2*n12,
    b*m1*n3 + 28*b*m3*n2 + 28*b*m4*n4 + c*m2*n1 + 11*m9*n1*n10 + 18*m9*n2*n9 + 18*m11*n1*n12 + 11*m11*n4*n9 + 18*m12*n1*n11 + 11*m12*n3*n9 + 3*n1*n5*n9 + 26*n2*n6*n10 + 15*n2*n7*n12 + 17*n2*n8*n11 + 17*n3*n6*n12 + 3*n3*n7*n11 + 15*n3*n8*n10 + 15*n4*n6*n11 + 17*n4*n7*n10 + 3*n4*n8*n12,
    a*m14*n13 + c*m13*n15 + 28*c*m15*n14 + 28*c*m16*n16 + 13*m13*n5*n14 + 16*m13*n6*n13 + 16*m15*n5*n16 + 13*m15*n8*n13 + 13*m16*n5*n15 + 16*m16*n7*n13 + 22*n1*n5*n13 + 22*n2*n6*n14 + 12*n2*n7*n16 + 24*n2*n8*n15 + 24*n3*n6*n16 + 7*n3*n7*n15 + 12*n3*n8*n14 + 12*n4*n6*n15 + 24*n4*n7*n14 + 7*n4*n8*n16 + 3*n5^2*n13 + 26*n6^2*n14 +
        3*n6*n7*n16 + 3*n6*n8*n15 + 3*n7^2*n15 + 3*n7*n8*n14 + 3*n8^2*n16,
    a*m14*n9 + 28*c*m13*n11 + c*m15*n10 + c*m16*n12 + 16*m13*n5*n10 + 16*m13*n6*n9 + 13*m15*n5*n12 + 13*m15*n8*n9 + 16*m16*n5*n11 + 16*m16*n7*n9 + 22*n1*n5*n9 + 7*n2*n6*n10 + 17*n2*n7*n12 + 5*n2*n8*n11 + 5*n3*n6*n12 + 22*n3*n7*n11 + 17*n3*n8*n10 + 17*n4*n6*n11 + 5*n4*n7*n10 + 22*n4*n8*n12 + 3*n5^2*n9 + 3*n6^2*n10 + 26*n6*n7*n12 +
        26*n6*n8*n11 + 26*n7^2*n11 + 26*n7*n8*n10 + 26*n8^2*n12,
    a*m14*n5 + 28*c*m13*n7 + c*m15*n6 + c*m16*n8 + 25*m13*n5*n6 + 4*m15*n5*n8 + 25*m16*n5*n7 + 11*n1*n5^2 + 18*n2*n6^2 + 11*n2*n7*n8 + 11*n3*n6*n8 + 11*n3*n7^2 + 11*n4*n6*n7 + 11*n4*n8^2 + 8*n5^3 + 8*n6^3 + 5*n6*n7*n8 + 21*n7^3 + 21*n8^3,
    a*m14*n1 + c*m13*n3 + 28*c*m15*n2 + 28*c*m16*n4 + 7*n1^2*n5 + 2*n1*n5^2 + 7*n2^2*n6 + 22*n2*n3*n8 + 22*n2*n4*n7 + 27*n2*n6^2 + 2*n2*n7*n8 + 22*n3^2*n7 + 22*n3*n4*n6 + 2*n3*n6*n8 + 2*n3*n7^2 + 22*n4^2*n8 + 2*n4*n6*n7 + 2*n4*n8^2,
    a*m13*n15 + 28*a*m15*n14 + 28*a*m16*n16 + 28*b*m14*n13 + 3*m13*n5*n14 + 26*m13*n6*n13 + 26*m15*n5*n16 + 3*m15*n8*n13 + 3*m16*n5*n15 + 26*m16*n7*n13 + 7*n1*n5*n13 + 7*n2*n6*n14 + 17*n2*n7*n16 + 5*n2*n8*n15 + 5*n3*n6*n16 + 22*n3*n7*n15 + 17*n3*n8*n14 + 17*n4*n6*n15 + 5*n4*n7*n14 + 22*n4*n8*n16 + n5^2*n13 + 28*n6^2*n14 +
        n6*n7*n16 + n6*n8*n15 + n7^2*n15 + n7*n8*n14 + n8^2*n16,
    a*m13*n11 + 28*a*m15*n10 + 28*a*m16*n12 + b*m14*n9 + 3*m13*n5*n10 + 3*m13*n6*n9 + 26*m15*n5*n12 + 26*m15*n8*n9 + 3*m16*n5*n11 + 3*m16*n7*n9 + 22*n1*n5*n9 + 7*n2*n6*n10 + 17*n2*n7*n12 + 5*n2*n8*n11 + 5*n3*n6*n12 + 22*n3*n7*n11 + 17*n3*n8*n10 + 17*n4*n6*n11 + 5*n4*n7*n10 + 22*n4*n8*n12 + 28*n5^2*n9 + 28*n6^2*n10 + n6*n7*n12 +
        n6*n8*n11 + n7^2*n11 + n7*n8*n10 + n8^2*n12,
    a*m13*n7 + 28*a*m15*n6 + 28*a*m16*n8 + b*m14*n5 + 28*m13*n5*n6 + m15*n5*n8 + 28*m16*n5*n7 + 11*n1*n5^2 + 18*n2*n6^2 + 11*n2*n7*n8 + 11*n3*n6*n8 + 11*n3*n7^2 + 11*n4*n6*n7 + 11*n4*n8^2 + 11*n5^3 + 11*n6^3 + 25*n6*n7*n8 + 18*n7^3 + 18*n8^3,
    a*m13*n3 + 28*a*m15*n2 + 28*a*m16*n4 + 28*b*m14*n1 + 22*n1^2*n5 + 3*n1*n5^2 + 22*n2^2*n6 + 7*n2*n3*n8 + 7*n2*n4*n7 + 26*n2*n6^2 + 3*n2*n7*n8 + 7*n3^2*n7 + 7*n3*n4*n6 + 3*n3*n6*n8 + 3*n3*n7^2 + 7*n4^2*n8 + 3*n4*n6*n7 + 3*n4*n8^2,
    a*m10*n13 + c*m9*n15 + 28*c*m11*n14 + c*m12*n16 + 27*m9*n1*n14 + 27*m9*n2*n13 + 2*m11*n1*n16 + 2*m11*n4*n13 + 2*m12*n1*n15 + 2*m12*n3*n13 + 23*n1^2*n13 + 7*n1*n5*n13 + 6*n2^2*n14 + 23*n2*n3*n16 + 23*n2*n4*n15 + 7*n2*n6*n14 + 9*n2*n7*n16 + 13*n2*n8*n15 + 23*n3^2*n15 + 23*n3*n4*n14 + 13*n3*n6*n16 + 22*n3*n7*n15 + 9*n3*n8*n14 +
        23*n4^2*n16 + 9*n4*n6*n15 + 13*n4*n7*n14 + 22*n4*n8*n16,
    a*m10*n9 + 28*c*m9*n11 + c*m11*n10 + 28*c*m12*n12 + 2*m9*n1*n10 + 27*m9*n2*n9 + 27*m11*n1*n12 + 2*m11*n4*n9 + 27*m12*n1*n11 + 2*m12*n3*n9 + 23*n1^2*n9 + 7*n1*n5*n9 + 23*n2^2*n10 + 6*n2*n3*n12 + 6*n2*n4*n11 + 22*n2*n6*n10 + 20*n2*n7*n12 + 16*n2*n8*n11 + 6*n3^2*n11 + 6*n3*n4*n10 + 16*n3*n6*n12 + 7*n3*n7*n11 + 20*n3*n8*n10 +
        6*n4^2*n12 + 20*n4*n6*n11 + 16*n4*n7*n10 + 7*n4*n8*n12,
    a*m10*n5 + 28*c*m9*n7 + c*m11*n6 + 28*c*m12*n8 + 6*n1^2*n5 + 10*n1*n5^2 + 6*n2^2*n6 + 23*n2*n3*n8 + 23*n2*n4*n7 + 19*n2*n6^2 + 10*n2*n7*n8 + 23*n3^2*n7 + 23*n3*n4*n6 + 10*n3*n6*n8 + 10*n3*n7^2 + 23*n4^2*n8 + 10*n4*n6*n7 + 10*n4*n8^2,
    a*m10*n1 + c*m9*n3 + 28*c*m11*n2 + c*m12*n4 + 16*n1^3 + 2*n1^2*n5 + 13*n2^3 + 2*n2^2*n6 + 19*n2*n3*n4 + 27*n2*n3*n8 + 27*n2*n4*n7 + 16*n3^3 + 27*n3^2*n7 + 27*n3*n4*n6 + 16*n4^3 + 27*n4^2*n8,
    a*m9*n15 + 28*a*m11*n14 + a*m12*n16 + 28*b*m10*n13 + 2*m9*n1*n14 + 2*m9*n2*n13 + 27*m11*n1*n16 + 27*m11*n4*n13 + 27*m12*n1*n15 + 27*m12*n3*n13 + 6*n1^2*n13 + 12*n1*n5*n13 + 23*n2^2*n14 + 6*n2*n3*n16 + 6*n2*n4*n15 + 12*n2*n6*n14 + 14*n2*n7*n16 + 3*n2*n8*n15 + 6*n3^2*n15 + 6*n3*n4*n14 + 3*n3*n6*n16 + 17*n3*n7*n15 + 14*n3*n8*n14
        + 6*n4^2*n16 + 14*n4*n6*n15 + 3*n4*n7*n14 + 17*n4*n8*n16,
    a*m9*n11 + 28*a*m11*n10 + a*m12*n12 + b*m10*n9 + 2*m9*n1*n10 + 27*m9*n2*n9 + 27*m11*n1*n12 + 2*m11*n4*n9 + 27*m12*n1*n11 + 2*m12*n3*n9 + 23*n1^2*n9 + 17*n1*n5*n9 + 23*n2^2*n10 + 6*n2*n3*n12 + 6*n2*n4*n11 + 12*n2*n6*n10 + 14*n2*n7*n12 + 3*n2*n8*n11 + 6*n3^2*n11 + 6*n3*n4*n10 + 3*n3*n6*n12 + 17*n3*n7*n11 + 14*n3*n8*n10 +
        6*n4^2*n12 + 14*n4*n6*n11 + 3*n4*n7*n10 + 17*n4*n8*n12,
    a*m9*n7 + 28*a*m11*n6 + a*m12*n8 + b*m10*n5 + 6*n1^2*n5 + 14*n1*n5^2 + 6*n2^2*n6 + 23*n2*n3*n8 + 23*n2*n4*n7 + 15*n2*n6^2 + 14*n2*n7*n8 + 23*n3^2*n7 + 23*n3*n4*n6 + 14*n3*n6*n8 + 14*n3*n7^2 + 23*n4^2*n8 + 14*n4*n6*n7 + 14*n4*n8^2,
    a*m9*n3 + 28*a*m11*n2 + a*m12*n4 + 28*b*m10*n1 + 13*n1^3 + 22*n1^2*n5 + 16*n2^3 + 22*n2^2*n6 + 10*n2*n3*n4 + 7*n2*n3*n8 + 7*n2*n4*n7 + 13*n3^3 + 7*n3^2*n7 + 7*n3*n4*n6 + 13*n4^3 + 7*n4^2*n8,
    a*m6*n13 + c*m5*n15 + 28*c*m7*n14 + c*m8*n16 + 19*n5*n13^2 + 19*n6*n14^2 + 10*n6*n15*n16 + 10*n7*n14*n16 + 10*n7*n15^2 + 10*n8*n14*n15 + 10*n8*n16^2,
    a*m6*n9 + 28*c*m5*n11 + c*m7*n10 + 28*c*m8*n12 + 28*n2*n11*n16 + n2*n12*n15 + n3*n10*n16 + 28*n3*n12*n14 + 28*n4*n10*n15 + n4*n11*n14 + 19*n5*n9*n13 + 10*n6*n10*n14 + 7*n6*n11*n16 + 12*n6*n12*n15 + 12*n7*n10*n16 + 19*n7*n11*n15 + 7*n7*n12*n14 + 7*n8*n10*n15 + 12*n8*n11*n14 + 19*n8*n12*n16,
    a*m6*n5 + 28*c*m5*n7 + c*m7*n6 + 28*c*m8*n8 + 5*m13*n5*n14 + 24*m13*n6*n13 + 24*m15*n5*n16 + 5*m15*n8*n13 + 5*m16*n5*n15 + 24*m16*n7*n13 + 11*n2*n7*n16 + 18*n2*n8*n15 + 18*n3*n6*n16 + 11*n3*n8*n14 + 11*n4*n6*n15 + 18*n4*n7*n14 + 3*n5^2*n13 + 26*n6^2*n14 + 3*n6*n7*n16 + 3*n6*n8*n15 + 3*n7^2*n15 + 3*n7*n8*n14 + 3*n8^2*n16,
    a*m6*n1 + c*m5*n3 + 28*c*m7*n2 + c*m8*n4 + 18*m9*n1*n14 + 18*m9*n2*n13 + 11*m11*n1*n16 + 11*m11*n4*n13 + 11*m12*n1*n15 + 11*m12*n3*n13 + 8*n1*n5*n13 + 8*n2*n6*n14 + 28*n2*n7*n16 + 22*n2*n8*n15 + 22*n3*n6*n16 + 21*n3*n7*n15 + 28*n3*n8*n14 + 28*n4*n6*n15 + 22*n4*n7*n14 + 21*n4*n8*n16,
    a*m5*n15 + 28*a*m7*n14 + a*m8*n16 + 28*b*m6*n13 + 4*n5*n13^2 + 4*n6*n14^2 + 25*n6*n15*n16 + 25*n7*n14*n16 + 25*n7*n15^2 + 25*n8*n14*n15 + 25*n8*n16^2,
    a*m5*n11 + 28*a*m7*n10 + a*m8*n12 + b*m6*n9 + 28*n2*n11*n16 + n2*n12*n15 + n3*n10*n16 + 28*n3*n12*n14 + 28*n4*n10*n15 + n4*n11*n14 + 25*n5*n9*n13 + 4*n6*n10*n14 + 23*n6*n11*n16 + 2*n6*n12*n15 + 2*n7*n10*n16 + 25*n7*n11*n15 + 23*n7*n12*n14 + 23*n8*n10*n15 + 2*n8*n11*n14 + 25*n8*n12*n16,
    a*m5*n7 + 28*a*m7*n6 + a*m8*n8 + b*m6*n5 + 7*m13*n5*n14 + 22*m13*n6*n13 + 22*m15*n5*n16 + 7*m15*n8*n13 + 7*m16*n5*n15 + 22*m16*n7*n13 + 11*n2*n7*n16 + 18*n2*n8*n15 + 18*n3*n6*n16 + 11*n3*n8*n14 + 11*n4*n6*n15 + 18*n4*n7*n14 + 8*n5^2*n13 + 21*n6^2*n14 + 8*n6*n7*n16 + 8*n6*n8*n15 + 8*n7^2*n15 + 8*n7*n8*n14 + 8*n8^2*n16,
    a*m5*n3 + 28*a*m7*n2 + a*m8*n4 + 28*b*m6*n1 + 11*m9*n1*n14 + 11*m9*n2*n13 + 18*m11*n1*n16 + 18*m11*n4*n13 + 18*m12*n1*n15 + 18*m12*n3*n13 + 11*n1*n5*n13 + 11*n2*n6*n14 + 13*n2*n7*n16 + 5*n2*n8*n15 + 5*n3*n6*n16 + 18*n3*n7*n15 + 13*n3*n8*n14 + 13*n4*n6*n15 + 5*n4*n7*n14 + 18*n4*n8*n16,
    a*m2*n13 + c*m1*n15 + 28*c*m3*n14 + 28*c*m4*n16 + 28*n2*n11*n16 + n2*n12*n15 + n3*n10*n16 + 28*n3*n12*n14 + 28*n4*n10*n15 + n4*n11*n14 + 10*n5*n9*n13 + 19*n6*n10*n14 + 17*n6*n11*n16 + 22*n6*n12*n15 + 22*n7*n10*n16 + 10*n7*n11*n15 + 17*n7*n12*n14 + 17*n8*n10*n15 + 22*n8*n11*n14 + 10*n8*n12*n16,
    a*m2*n9 + 28*c*m1*n11 + c*m3*n10 + c*m4*n12 + 10*n5*n9^2 + 10*n6*n10^2 + 19*n6*n11*n12 + 19*n7*n10*n12 + 19*n7*n11^2 + 19*n8*n10*n11 + 19*n8*n12^2,
    a*m2*n5 + 28*c*m1*n7 + c*m3*n6 + c*m4*n8 + 5*m13*n5*n10 + 5*m13*n6*n9 + 24*m15*n5*n12 + 24*m15*n8*n9 + 5*m16*n5*n11 + 5*m16*n7*n9 + 11*n2*n7*n12 + 18*n2*n8*n11 + 18*n3*n6*n12 + 11*n3*n8*n10 + 11*n4*n6*n11 + 18*n4*n7*n10 + 26*n5^2*n9 + 26*n6^2*n10 + 3*n6*n7*n12 + 3*n6*n8*n11 + 3*n7^2*n11 + 3*n7*n8*n10 + 3*n8^2*n12,
    a*m2*n1 + c*m1*n3 + 28*c*m3*n2 + 28*c*m4*n4 + 18*m9*n1*n10 + 11*m9*n2*n9 + 11*m11*n1*n12 + 18*m11*n4*n9 + 11*m12*n1*n11 + 18*m12*n3*n9 + 21*n1*n5*n9 + 8*n2*n6*n10 + 28*n2*n7*n12 + 22*n2*n8*n11 + 22*n3*n6*n12 + 21*n3*n7*n11 + 28*n3*n8*n10 + 28*n4*n6*n11 + 22*n4*n7*n10 + 21*n4*n8*n12,
    a*m1*n15 + 28*a*m3*n14 + 28*a*m4*n16 + 28*b*m2*n13 + n2*n11*n16 + 28*n2*n12*n15 + 28*n3*n10*n16 + n3*n12*n14 + n4*n10*n15 + 28*n4*n11*n14 + 25*n5*n9*n13 + 4*n6*n10*n14 + 2*n6*n11*n16 + 23*n6*n12*n15 + 23*n7*n10*n16 + 25*n7*n11*n15 + 2*n7*n12*n14 + 2*n8*n10*n15 + 23*n8*n11*n14 + 25*n8*n12*n16,
    a*m1*n11 + 28*a*m3*n10 + 28*a*m4*n12 + b*m2*n9 + 4*n5*n9^2 + 4*n6*n10^2 + 25*n6*n11*n12 + 25*n7*n10*n12 + 25*n7*n11^2 + 25*n8*n10*n11 + 25*n8*n12^2,
    a*m1*n7 + 28*a*m3*n6 + 28*a*m4*n8 + b*m2*n5 + 7*m13*n5*n10 + 7*m13*n6*n9 + 22*m15*n5*n12 + 22*m15*n8*n9 + 7*m16*n5*n11 + 7*m16*n7*n9 + 11*n2*n7*n12 + 18*n2*n8*n11 + 18*n3*n6*n12 + 11*n3*n8*n10 + 11*n4*n6*n11 + 18*n4*n7*n10 + 21*n5^2*n9 + 21*n6^2*n10 + 8*n6*n7*n12 + 8*n6*n8*n11 + 8*n7^2*n11 + 8*n7*n8*n10 + 8*n8^2*n12,
    a*m1*n3 + 28*a*m3*n2 + 28*a*m4*n4 + 28*b*m2*n1 + 11*m9*n1*n10 + 18*m9*n2*n9 + 18*m11*n1*n12 + 11*m11*n4*n9 + 18*m12*n1*n11 + 11*m12*n3*n9 + 18*n1*n5*n9 + 11*n2*n6*n10 + 13*n2*n7*n12 + 5*n2*n8*n11 + 5*n3*n6*n12 + 18*n3*n7*n11 + 13*n3*n8*n10 + 13*n4*n6*n11 + 5*n4*n7*n10 + 18*n4*n8*n12,
    a*b*c + 21*m13*n5*n6 + 8*m15*n5*n8 + 21*m16*n5*n7 + 8*n1^3 + 12*n1*n5^2 + 21*n2^3 + 24*n2*n3*n4 + 17*n2*n6^2 + 12*n2*n7*n8 + 8*n3^3 + 12*n3*n6*n8 + 12*n3*n7^2 + 8*n4^3 + 12*n4*n6*n7 + 12*n4*n8^2 + 17*n5^3 + 17*n6^3 + 7*n6*n7*n8 + 12*n7^3 + 12*n8^3,
    a^2*c + 28*a*b^2 + 28*b*c^2 + 27*m13*n5*n6 + 2*m15*n5*n8 + 27*m16*n5*n7 + 24*n1^3 + 5*n2^3 + 14*n2*n3*n4 + 24*n3^3 + 24*n4^3 + 10*n5^3 + 10*n6^3 + 28*n6*n7*n8 + 19*n7^3 + 19*n8^3,
    a^2*b + a*c^2 + 28*b^2*c + 3*m13*n5*n6 + 26*m15*n5*n8 + 3*m16*n5*n7 + 5*n1^3 + 24*n2^3 + 15*n2*n3*n4 + 5*n3^3 + 5*n4^3 + 3*n5^3 + 3*n6^3 + 20*n6*n7*n8 + 26*n7^3 + 26*n8^3,
    a^3 + b^3 + 28*c^3 + 24*m13*n5*n6 + 5*m15*n5*n8 + 24*m16*n5*n7 + 5*n1^3 + 14*n1*n5^2 + 24*n2^3 + 15*n2*n3*n4 + 15*n2*n6^2 + 14*n2*n7*n8 + 5*n3^3 + 14*n3*n6*n8 + 14*n3*n7^2 + 5*n4^3 + 14*n4*n6*n7 + 14*n4*n8^2 + 7*n5^3 + 7*n6^3 + 8*n6*n7*n8 + 22*n7^3 + 22*n8^3
];

function CheckRelations()

rels:=[];
for i in [#rels+1..#seqs] do Append(~rels,ProdRel(seqs[i])); end for;

Sort(~rels);

return Sort(Reduce(rels)) eq relsR;
end function;

function CheckDetermination()

relsC:=[];
for i in [#relsC+1..#seqs] do Append(~relsC,CentRel(seqs[i])); end for;

// We first analyse the centralizer in E6 of L.

BCent:=GroebnerBasis(ChangeOrder(ideal<R|relsC>,"grevlex")); BCent:=[R!i:i in BCent];
n1^3-1 in GroebnerBasis(BCent);

// Thus we replace n1 by t11 to alert us to this being a cube root of unity.

cCent:=ChangeCoefficient(coeffs1,n1,t11);
// All of the equations are homogeneous in BCent, but we can move t11 around. 

[n1*n6 + 28*m6*n16,
n1*n7 + 28*m6*n14,
n1*n8 + m6*n15,
n1*n10 + a*n15,
n1*n11 + 28*a*n16,
n1*n12 + a*n14] subset BCent;

// Move n1 to the right (and relabel t11) and we obtain our first substitutions.

cCent:=ChangeCoefficients(cCent,[n6,n7,n8,n10,n11,n12],[t11^2*i:i in [m6*n16,m6*n14,-m6*n15,-a*n15,a*n16,-a*n14]]);
cCent:=ChangeLinearCoefficients(cCent,BCent);

BCent:=GroebnerBasis(ChangeOrder(ideal<R|[Evaluate(i,cCent):i in BCent]>,"grevlex")); BCent:=[R!i:i in BCent];

[a^2 + 28*m6*t11,a*m6 + 28*t11^2] subset BCent;

// The first yields t11^2*a^2=m6. Thus the second yields a^3=1.
// Hence we may replace a by t12, another cube root of unity.
// And replace m6 by t11^2*t12^2, while we are at it.

cCent:=ChangeCoefficients(cCent,[a,m6],[t12,t11^2*t12^2]);
BCent:=GroebnerBasis(ChangeOrder(ideal<R|[Evaluate(i,cCent):i in BCent]>,"grevlex")); BCent:=[R!i:i in BCent];

[m4*t11 + m9*t12,
m12*t11 + m13*t12,
28*m1*t12 + m16*t11,
28*m4*t12 + m15*t11,
m11*t11 + 28*m16*t12] subset BCent;

cCent:=ChangeCoefficients(cCent,[m9,m13,m1,m15,m11],[-t11*t12^2*m4,-t11*t12^2*m12,t11*t12^2*m16,t12*t11^2*m4,t12*t11^2*m16]);
BCent:=GroebnerBasis(ChangeOrder(ideal<R|[Evaluate(i,cCent):i in BCent]>,"grevlex")); BCent:=[R!i:i in BCent];

// We will be left with m4,m12,m16 as variables, subject to a single cubic equation.
// So let's eliminate variables that can be expressed in this form.

[m3*t11 - m12*t12,
m4*m12 + m16^2 - n16*t12,
m4*m16 + m12^2 - n15*t11,
m4^2*t11 - m12*m16*t11 + n14*t12^2] subset BCent;

cCent:=ChangeCoefficients(cCent,[m3,n16,n15],[t12*t11^2*m12,t12^2*(m4*m12 + m16^2),t11^2*(m4*m16 + m12^2)]);
cCent:=ChangeCoefficients(cCent,[n14],[-t12*(m4^2*t11 - m12*m16*t11)]);

// We now relabel the elements m4,m12 and m16 by t1,t2,t3, with a sign to make things nicer.

cCent:=ChangeCoefficients(cCent,[m4,m12,m16],[-t1,t2,t3]);

BCent:=GroebnerBasis([Evaluate(i,cCent):i in BCent]);

BCent eq [t1^3+t2^3+t3^3-3*t1*t2*t3-1,t11^3-1,t12^3-1];

// Identify the torus S, the normal subgroup.

cCentS:=ChangeCoefficients(cCent,[t11,t12],[1,1]);

scalS:=Evaluate(scal,cCentS);

// Check the action of S on (m9,m12,m11) given in the paper.

[scal[9,j]:j in [22,9,19]] eq [m9,m11,m12];
[(scal*scalS)[9,j]:j in [22,9,19]] eq 
[m9*t3 + m11*t1 + m12*t2,
 m9*t2 + m11*t3 + m12*t1,
 m9*t1 + m11*t2 + m12*t3];

// We also need an element of order 3. Notice that [0,0,t11^2] is an element of S of order 3.
// So we may choose [t1,t2,t3,t11,t12] to be [0,0,t11^2,t11,1].

cCentw:=ChangeCoefficients(cCent,[t1,t2,t3,t12],[0,0,t11^2,1]);
scalw:=Evaluate(scal,cCentw);
[(scal*scalw)[9,j]:j in [22,9,19]] eq [m9*t11^3,m11*t11^4,m12*t11^2];

// Since t11^3=1, this is [m9,m9*t11,m12*t11^2]. This obviously permutes the three
// orbits of dimension 2 and the three orbits of dimension 1.

// Now we need nine elements that stabilize [m9,m11,m12]=[0,1,0].

scalCent:=Evaluate(scal,cCent);

[(scal*scalCent)[9,j]:j in [22,9,19]] eq 
[
    m9*t3*t11*t12^2 + m11*t1*t11*t12^2 + m12*t2*t11*t12^2,
    m9*t2*t11^2*t12 + m11*t3*t11^2*t12 + m12*t1*t11^2*t12,
    m9*t1 + m11*t2 + m12*t3
];

// Thus m9=0,m11=1,m12=0, and we need to stabilize this.

// thus t1*t11*t12^2=0, t2=0 (first and second equations) and

// t3*t11^2*t12=1.

cCentx:=ChangeCoefficients(cCent,[t1,t2,t3,t11,t12],[0,0,t11^2,t11,t11^2]);
cCenty:=ChangeCoefficients(cCent,[t1,t2,t3,t11,t12],[0,0,t12^2,1,t12]);
scalx:=Evaluate(scal,cCentx);
scaly:=Evaluate(scal,cCenty);


// We want to check that the centralizer in GL of H acts as we expect.
// We also check that the scalar that affects m9 to m12 also affects n1 to n4.

H:=sub<GL(27,F)|l1,l2,h1>;
MH:=GModule(H);
homsH:=AHom(MH,MH);
matH:=t1*mats!(homsH.1)+t2*mats!(homsH.2)+t3*mats!(homsH.5)+t4*mats!(homsH.6)+t5*mats!(homsH.3)+t6*mats!(homsH.4);

[scal[i,j]:i,j in [1,27]] eq 
[m1, m2,
 m5, m6];

[(matH*scal)[i,j]:i,j in [1,27]] eq
[m1*t1 + m5*t2, m2*t1 + m6*t2,
 m1*t4 + m5*t3, m2*t4 + m6*t3];

// So matH acts like a 2x2 matrix on m1,m2,m5,m6. It does the same on m3,m4,m7,m8.

[scal[i,j]:j in [6,10],i in [1,27]] eq 
[m3,  m4,
 m7, -m8];

[(matH*scal)[i,j]:j in [6,10],i in [1,27]] eq 
[m3*t1 + m7*t2, m4*t1 - m8*t2,
 m3*t4 + m7*t3, m4*t4 - m8*t3];

// The action on m9,m10,m11,m12 is multiplication by t5:
[scal[9,j]:j in [22,4,9,19]] eq [m9,m10,m11,m12];
[(matH*scal)[9,j]:j in [22,4,9,19]] eq [i*t5:i in [m9,m10,m11,m12]];

// and same for n1 to n4.
[scal[2,j]:j in [2,20,26,16]] eq [n1,n2,n3,n4];
[(matH*scal)[2,j]:j in [2,20,26,16]] eq [i*t5:i in [n1,n2,n3,n4]];

// t6, the last variable, affects n5 to n8 and m13 to m16.
[scal[3,j]:j in [18,3,5,15]] eq [n5,n6,n7,n8];
[(matH*scal)[3,j]:j in [18,3,5,15]] eq [i*t6:i in [n5,n6,n7,n8]];

[scal[10,j]:j in [1,27,6,10]] eq [m13,m14,m15,m16];
[(matH*scal)[10,j]:j in [1,27,6,10]] eq [i*t6:i in [m13,m14,m15,m16]];


// We see from this that we may choose [m1,m2,m5,m6] as a 2x2 matrix, so one of
// [0,0,0,0], [1,i,0,0], [0,1,0,0] and [1,0,0,1].

// We can choose [m9,m11,m12] from one of the orbits given in the paper.
// These are [0,1,0], [-1,1,0], [1,1,1] and [0,0,0]. We do [0,0,0] first as it is easiest, then
// [0,1,0], where the solution is found. Then [-1,1,0] is a bit easier and yields no solutions.
// Finally, [1,1,1] is fairly easy to eliminate.

// The last centralizer element allows use to normalize one of [n5,n6,n7,n8,m13,m14,m15,m16],
// and is only used when obtaining the solutions in the [0,1,0] case.

// We also have the elements x and y of the centralizer constructed above. These will act
// regularly on the set of solutions we obtain, thus yielding a unique G-class.

///////////////////
// Option 1: [m9,m11,m12]=[0,0,0]
///////////////////

// We start by eliminating the highly degenerate cases of all zeros. First, [m9,m11,m12]=[0,0,0].

coef1:=ChangeCoefficients(coeffs1,[m9,m11,m12],[0,0,0]);
rel1:=Reduce([Evaluate(i,coef1):i in relsR]);
m10^3 in rel1;

// Thus the determinant is 0, a contradiction.

///////////////////
// Option 2: [m1,m2,m5,m6]=[0,0,0,0]
///////////////////

coef2:=ChangeCoefficients(coeffs1,[m1,m2,m5,m6],[0,0,0,0]);
// We can now reduce [m3,m4,m7,m8] instead. Notice that this has to be rank 2, so the identity.

coef2:=ChangeCoefficients(coef2,[m4,m7],[0,0]);
rel2:=Reduce([Evaluate(i,coef2):i in relsR]);
coef2:=ChangeCoefficients(coef2,[m3,m8],[1,1]);
rel2:=[Evaluate(i,coef2):i in rel2];

GroebnerBasis(rel2) eq [1];

// This yields a contradiction.

///////////////////
// Option 3a: [m9,m11,m12]=[0,1,0] and [m1,m2,m5,m6]=[1,0,0,1].
///////////////////

// We now consider the easier option for [m9,m11,m12], namely [0,1,0].

// We first assume that [m1,m2,m5,m6] is invertible, and find a unique solution.

// This fixes seven values and is a 'generic point'. We will use the last specialization right at the end.

co2a:=ChangeCoefficients(coeffs1,[m1,m2,m5,m6,m9,m11,m12],[1,0,0,1,0,1,0]);
rels2a:=[Evaluate(i,co2a):i in relsR];
[m7^2 + 28*m10,m3^2 + m4] subset rels2a;
co2a:=ChangeCoefficients(co2a,[m10,m4],[m7^2,-m3^2]);
rels2a:=[Evaluate(i,co2a):i in rels2a];
BB2a:=StripBasis(GroebnerBasis(ChangeOrder(ideal<R|[i:i in rels2a|#Monomials(i) le 4]>,"grevlex")),[]);
m3^2*n2 - m3*n3 - n4 in BB2a;

co2a:=ChangeCoefficient(co2a,n4,m3^2*n2 - m3*n3);
BB2a:=GroebnerBasis(ChangeOrder(ideal<R|[Evaluate(i,co2a):i in BB2a]>,"grevlex"));

// Suppose that m3=1.

co2a1:=ChangeCoefficient(co2a,m3,1);
rels2a1:=[Evaluate(i,co2a1):i in rels2a]; rels2a1:=[R!i:i in rels2a1];

[m7^2*n3 + 17*n1^3,28*m7^6 + 22*n1^3 + 1] subset rels2a1;

// By the second equation, either m7 ne 0 or n1 ne 0, and if n1 ne 0 then m7,n3 ne 0 by the first equation. Thus m7 ne 0.

rels2a1:=SimplifyBasis(rels2a1,[m7]);
BB2a1:=StripBasis(GroebnerBasis(ChangeOrder(ideal<R|[i:i in rels2a1|#Monomials(i) le 9]>,"grevlex")),[m7]);

// As n4=m3^2*n2 - m3*n3, cannot have n1,n2,n3=0. If n2=n3=0 then n1=0:
n1^3 + 11*n1*n2 + 18*n1*n3 in BB2a1;
// Thus n13=0.
n2*n13 in BB2a1 and n3*n13 in BB2a1;
co2a1:=ChangeCoefficient(co2a1,n13,0);
co2a1:=ChangeLinearCoefficients(co2a1,BB2a1);
rels2a1:=[Evaluate(i,co2a1):i in rels2a1]; rels2a1:=[R!i:i in rels2a1];
GroebnerBasis([i:i in rels2a1|#Monomials(i) le 2]) eq [1];

// Thus m3 is not equal to 1. By scal1 this means that m3^3 ne 1.

[i*j:i in [n9,n13],j in [n12,n16]] subset [R!i:i in BB2a];

// One of these doesn't work.

co2a2:=ChangeCoefficients(co2a,[n9,n13],[0,0]); co2a2:=[R!i:i in co2a2];
rels2a2:=[Evaluate(i,co2a2):i in rels2a];
BB2a2:=StripBasis(GroebnerBasis(ChangeOrder(ideal<R|[Evaluate(R!i,co2a2):i in BB2a] cat [i:i in rels2a2|#Monomials(i) le 4]>,"grevlex")),[]);
BB2a2:=StripBasis(GroebnerBasis(ChangeOrder(ideal<R|[R!i:i in BB2a2] cat [i:i in rels2a2|#Monomials(i) eq 7]>,"grevlex")),[]);
m3^3-1 in BB2a2;

// So n12=n16=0

co2a:=ChangeCoefficients(co2a,[n12,n16],[0,0]); co2a:=[R!i:i in co2a];
rels2a:=[Evaluate(i,co2a):i in rels2a];
rels2a:=SimplifyBasis(rels2a,[]);

// Further reduce relsR. To do this we take account of the zeros in the basis.

zz:=[i:i in [1..35]|co2a[i] eq 0];
coo:=coeffs1; for i in zz do coo[i]:=0; end for;

relsRo:=SimplifyBasis([Evaluate(i,coo):i in relsR],[m3-1,m3 + ww^140,m3 + ww^700]); // This is the factorization of m3^3-1.
relsRo:=Reduce(relsRo);
rel2a:=SimplifyBasis([Evaluate(i,co2a):i in relsRo],[m3-1,m3 + ww^140,m3 + ww^700]);

m3*m16 + m13 in rel2a;
co2a:=ChangeCoefficient(co2a,m13,-m3*m16);
rel2a:=SimplifyBasis([Evaluate(i,co2a):i in rel2a],[m3-1,m3 + ww^140,m3 + ww^700]);

// Suppose that m3=0.
co2a3:=ChangeCoefficient(co2a,m3,0);
rel2a3:=[Evaluate(i,co2a3):i in rel2a];
BB2a3:=StripBasis(GroebnerBasis(rel2a3),[]);
[a,b,c] subset BB2a3;
// Thus m3 ne 0.

rel2a:=SimplifyBasis(rel2a,[m3]);

// Suppose that m7=0.
co2a4:=ChangeCoefficient(co2a,m7,0);
rel2a4:=[Evaluate(i,co2a4):i in rel2a];
StripBasis(GroebnerBasis(ChangeOrder(ideal<R|[R!i:i in rel2a4|#Monomials(i) le 6]>,"grevlex")),[]) eq [1];

// Thus m7 ne 0 as well.
rel2a:=SimplifyBasis(rel2a,[m3,m7,m3-1]);

m3*n10 - n11 in rel2a;

co2a:=ChangeCoefficient(co2a,n11,m3*n10);
rel2a:=[Normalize(Evaluate(i,co2a)):i in rel2a];
rel2a:=SimplifyBasis(rel2a,[m3,m7,m3-1,m3 + ww^140,m3 + ww^700]);

// This cracks it.
BB2a:=StripBasis(GroebnerBasis(ChangeOrder(ideal<R|[i:i in rel2a|#Monomials(i) le 6]>,"grevlex")),[m3,m7,m3 + ww^140,m3 + ww^700]);
BB2a:=GroebnerBasis([R!i:i in BB2a]);

co2a:=ChangeLinearCoefficients(co2a,BB2a);
rel2a:=[Normalize(Evaluate(i,co2a)):i in rel2a];
BB2a:=GroebnerBasis(rel2a);

co2a:=ChangeLinearCoefficients(co2a,BB2a);
rel2a:=[Normalize(Evaluate(i,co2a)):i in rel2a];
BB2a:=GroebnerBasis(rel2a);

BB2a eq [n14^3-24,n15^3-5];

// We have our last variable that we can set to 0 or 1. Also, n8 is still free
// It also divides the entire row of the n-matrix containing it.

&and[IsDivisibleBy(co2a[i],n8):i in [24..27]];

// Thus it cannot be zero.

ctall:=ChangeCoefficient(co2a,n8,1);

// This yields nine solutions. We also have nine elements of the centralizer
// in E6 of L that are outside of S.

// Let us tidy up ctall by removing all of the powers of n14 and n15 above 2.
II:=ideal<R|BB2a>;
RR:=quo<R|II>;
for i in [1..35] do
  ctall[i]:=R!(RR!ctall[i]);
end for;

[ctall[i]:i in [1..35]] eq 
[
    7*n14^2*n15^2 + 12*n15^2,
    7*n14^2*n15^2 + 7*n15^2,
    22*n14^2*n15^2 + 19*n15^2,
    1,
    0,
    17*n14^2*n15,
    24*n14*n15^2,
    0,
    1,
    17*n15^2,
    19*n14^2,
    0,
    24*n15,
    1,
    0,
    19*n14^2*n15^2,
    18*n14*n15,
    0,
    4*n15,
    16*n14,
    7*n14*n15,
    28*n15^2,
    28*n14^2,
    2*n14^2,
    9*n14^2*n15,
    10*n14*n15^2,
    1,
    16*n15,
    7*n15^2,
    15*n14^2,
    0,
    12*n14*n15^2,
    n14,
    n15,
0];

scalall:=Evaluate(scal,ctall);

// We see that the nine elements are determined by their values on [n1,n2,n3,n4].

[scal[2,j]:j in [2,20,26,16]] eq [n1,n2,n3,n4];

[scalall[2,j]:j in [2,20,26,16]] eq [16*n14,7*n14*n15,-n15^2,-n14^2];


// Multiply by scalx
[(scalall*scalx)[2,j]:j in [2,20,26,16]] eq [16*n14*t11,7*n14*n15*t11^16,-n15^2*t11^12,-n14^2*t11^8];

// Since t11^3=1, this has the effect of multiplying n14 by t11 and leaving n15 alone.


// Multiply by scaly
[(scalall*scaly)[2,j]:j in [2,20,26,16]] eq [16*n14,7*n14*n15*t12^8,-n15^2*t12^7,-n14^2*t12^6];

// Since t12^3=1, this has the effect of multiplying n15 by t12^2 and leaving n14 alone.

// Thus <x,y> acts regularly on the nine solutions, and therefore there is a unique G-class,
// as claimed. To confirm, we construct two copies of PSU(3,3) containing H and lying in E6.

ct1:=ChangeCoefficients(ctall,[n14,n15],[7,-7]);
scal1:=GL(27,F)!Evaluate(scal,ct1);

h11:=h1^scal1;

// Double-check that h11 does indeed lie in E6.

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

// Now we need two matrices that, together with h11, generate PSU(3,3);

j1a:=GL(27,F)![[7,8,7,22,3,16,22,5,17,10,19,25,26,10,15,14,24,22,9,24,1,18,1,10,0,16,27],
[6,27,1,19,20,12,6,15,27,14,22,11,11,24,16,0,20,20,13,10,3,10,21,17,0,3,21],
[26,23,15,0,18,14,1,14,14,4,25,26,9,28,11,0,13,12,23,23,15,2,21,11,0,17,26],
[18,10,21,18,19,23,7,23,10,17,21,19,26,7,12,13,11,16,21,10,16,10,4,6,13,11,21],
[17,1,27,4,8,24,7,13,6,18,14,15,28,24,5,11,8,3,19,0,23,27,6,13,27,0,6],
[5,17,4,3,13,6,20,3,21,9,25,22,2,26,19,26,0,28,13,25,1,28,8,9,24,28,15],
[16,9,26,5,23,5,15,28,16,10,9,26,9,4,9,5,20,15,5,27,19,8,16,5,11,16,0],
[28,14,11,10,25,27,21,11,25,10,24,23,22,21,5,27,13,5,0,0,2,14,15,4,3,1,26],
[6,9,19,25,4,9,16,8,5,12,3,8,3,22,11,22,24,6,0,3,27,11,26,3,5,22,3],
[22,3,14,15,26,15,0,4,24,24,16,7,20,27,1,0,27,23,1,24,7,17,8,2,3,19,11],
[26,13,16,23,3,9,21,9,16,8,27,3,0,22,18,20,5,17,25,17,27,4,8,5,0,24,14],
[16,18,25,0,7,25,8,27,25,7,26,1,7,20,1,17,11,3,12,11,26,6,1,22,4,1,7],
[22,28,18,26,17,28,16,25,10,2,27,0,26,28,1,24,13,7,17,11,19,25,13,17,22,28,17],
[27,7,16,19,26,1,20,4,24,28,23,28,1,17,14,16,26,15,9,1,8,24,28,9,23,8,9],
[25,24,16,5,24,1,17,14,17,1,8,16,17,3,27,15,15,19,26,13,11,22,10,23,13,26,1],
[12,26,23,17,28,26,2,10,17,7,0,25,17,20,16,28,12,11,27,25,1,7,23,17,19,9,11],
[15,0,20,3,13,11,18,3,12,21,0,2,6,14,23,2,10,8,3,18,20,7,22,8,24,2,27],
[19,4,22,11,28,11,7,25,12,18,20,4,15,23,21,3,18,18,10,26,26,9,5,0,0,17,0],
[18,7,1,0,7,28,4,25,27,4,12,15,12,1,6,25,3,17,2,3,26,12,4,8,11,4,25],
[19,16,18,2,4,20,19,20,22,13,27,25,8,9,28,17,20,13,16,17,24,26,24,19,17,16,4],
[15,10,14,13,28,12,8,1,24,5,6,18,8,24,9,19,24,8,20,3,16,22,18,2,2,16,1],
[19,1,28,26,20,16,3,20,24,17,17,15,11,5,0,28,11,19,13,21,22,18,19,0,12,11,3],
[10,28,27,19,11,21,19,26,6,28,0,25,26,17,12,8,12,7,14,8,4,6,18,13,0,14,19],
[28,13,19,24,2,6,28,13,23,18,4,26,28,24,23,7,22,24,0,27,28,9,6,19,0,27,5],
[24,11,2,6,24,10,12,28,2,2,25,22,6,11,19,14,12,23,0,16,10,11,21,8,17,19,1],
[12,10,3,24,16,4,14,4,12,12,13,23,23,5,14,7,28,23,26,28,7,14,27,15,3,25,27],
[20,24,2,7,19,22,15,12,2,3,9,8,19,19,10,16,4,9,25,15,6,21,16,23,9,3,14]];

j1b:=GL(27,F)![[9,20,21,18,27,17,13,20,24,3,13,18,7,22,22,26,5,23,16,3,9,13,14,4,15,4,15],
[17,28,19,0,7,5,23,1,22,7,28,17,15,11,13,16,20,17,13,1,8,0,18,9,19,25,25],
[26,19,3,4,1,4,11,3,12,2,23,14,9,2,17,25,27,18,10,5,10,24,6,9,28,27,13],
[9,27,5,20,3,10,20,14,26,20,11,9,22,2,5,25,22,5,9,21,8,15,10,3,3,4,16],
[15,19,4,21,8,19,4,0,21,21,26,23,23,0,25,18,1,3,22,2,26,15,12,9,0,6,14],
[2,6,16,0,17,3,19,24,23,12,27,5,23,17,24,27,7,0,8,3,28,14,7,28,20,22,9],
[23,22,22,20,10,8,2,10,3,20,27,15,23,10,8,21,17,16,26,2,8,25,3,4,1,26,14],
[22,3,26,25,27,10,19,3,20,26,6,25,27,17,6,6,21,24,11,24,28,4,9,23,19,17,10],
[2,3,10,4,24,12,27,23,4,8,23,15,7,27,15,28,1,18,25,17,14,16,8,15,0,27,15],
[2,0,16,19,23,21,27,11,16,23,4,28,6,11,27,4,27,4,15,15,8,1,7,12,22,25,15],
[23,1,7,24,6,18,9,9,28,6,21,15,21,27,1,18,4,3,9,2,15,22,10,21,6,4,28],
[4,24,23,27,9,26,28,26,5,10,2,9,10,0,12,13,10,27,16,19,23,5,19,18,5,16,26],
[25,1,20,15,7,5,0,20,19,27,17,3,24,16,10,24,17,27,25,27,4,6,22,8,3,8,28],
[21,4,12,16,8,11,23,10,15,7,27,7,21,17,27,5,25,8,19,26,21,20,2,5,1,10,26],
[14,16,12,1,22,24,1,9,9,20,13,17,23,19,17,18,20,17,26,12,20,0,26,4,0,24,21],
[8,26,6,22,20,13,10,13,25,10,27,10,24,11,21,24,23,10,1,5,16,16,23,21,5,10,23],
[11,11,26,8,8,21,16,8,22,3,11,1,10,23,24,16,4,18,28,7,1,7,18,17,21,9,22],
[10,24,11,2,25,25,10,9,25,15,25,20,12,28,23,12,11,11,14,26,28,0,22,26,19,23,21],
[16,20,23,14,24,12,23,2,22,9,21,13,17,4,7,7,23,20,7,6,22,23,16,17,13,3,6],
[8,14,5,25,19,12,3,16,20,20,23,26,7,4,11,3,19,18,14,14,4,22,9,28,16,6,27],
[10,5,20,17,16,25,20,28,3,6,19,24,12,6,25,2,28,2,1,5,19,26,14,28,1,17,6],
[26,13,28,10,26,21,7,4,26,3,7,15,7,12,6,23,26,6,25,24,23,12,20,12,24,26,9],
[25,6,28,4,9,5,9,13,4,11,22,5,1,16,17,1,23,28,8,12,7,18,21,9,20,2,14],
[28,3,9,9,27,2,24,3,1,16,9,0,26,28,26,8,27,0,19,14,20,11,16,6,6,18,1],
[18,21,13,16,21,1,21,5,10,4,8,15,12,26,5,12,17,18,21,24,19,14,4,28,22,18,6],
[26,11,6,6,21,18,15,28,5,14,4,8,21,12,1,11,1,9,15,8,26,7,1,1,21,2,25],
[19,10,3,24,28,19,0,8,12,2,23,17,11,20,6,16,19,16,8,18,16,20,13,9,28,26,16]];

J1a:=sub<GL(27,F)|[l1,l2,h11,j1a]>;
J1b:=sub<GL(27,F)|[l1,l2,h11,j1b]>;

&and[IsIsomorphic(i,PSU(3,3)):i in [J1a,J1b]];

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

// If you want, you can verify that J1a and J1b generate a group E6(29).
// This code is commented out as it's completely unnecessary.

/*
GG:=sub<GL(27,F)|[j1a,j1b,h11]>;
GG_Tree:=CompositionTree(GG:Verify:=true);
CompositionTreeNonAbelianFactors(GG) eq [* <"E",6,29>*];
CompositionTreeOrder(GG) eq 29^36*(29^12-1)*(29^9-1)*(29^8-1)*(29^6-1)*(29^5-1)*(29^2-1);
*/

///////////////////
// Option 3b: [m9,m11,m12]=[0,1,0] and [m1,m2,m5,m6]=[1,m2,0,0]
///////////////////

// The next case is where [m1,m2,m5,m6]=[1,m2,0,0].

co2b:=ChangeCoefficients(coeffs1,[m1,m5,m6,m9,m11,m12],[1,0,0,0,1,0]);
rels2b:=[Evaluate(i,co2b):i in relsR];

m7^2 in rels2b;

co2b:=ChangeCoefficient(co2b,m7,0);
rels2b:=[Evaluate(i,co2b):i in relsR];

// m8 is the only variable from row 2 of the m-matrix, so m8 ne 0.

rels2b:=SimplifyBasis(rels2b,[m8]);

x1:=m8 + 21*n2*n11*n16 + 8*n2*n12*n15 + 8*n3*n10*n16 + 21*n3*n12*n14 + 21*n4*n10*n15 + 8*n4*n11*n14;
x2:=m8 + 25*n2*n11*n16 + 4*n2*n12*n15 + 4*n3*n10*n16 + 25*n3*n12*n14 + 25*n4*n10*n15 + 4*n4*n11*n14;

x1 in rels2b and x2 in rels2b;
x1/21-x2/25 eq 11*m8; // But m8 ne 0. This brings this branch to a close.

///////////////////
// Option 3c: [m9,m11,m12]=[0,1,0] and [m1,m2,m5,m6]=[0,1,0,0]
///////////////////

// Or we can have that [m1,m2,m5,m6]=[0,1,0,0].

co2c:=ChangeCoefficients(coeffs1,[m1,m2,m5,m6,m9,m11,m12],[0,1,0,0,0,1,0]);
rels2c:=[Evaluate(i,co2c):i in relsR];

m7^2 in rels2c;

co2c:=ChangeCoefficient(co2c,m7,0);
rels2c:=[Evaluate(i,co2c):i in relsR];

// m8 is the only variable from row 2 of the m-matrix, so m8 ne 0.

rels2c:=SimplifyBasis(rels2c,[m8]);

BB2c:=StripBasis(GroebnerBasis(ChangeOrder(ideal<R|rels2c>,"grevlex")),[]);
[n13,n14,n15,n16] subset [R!i:i in BB2c];

// This contradiction finishes this branch, where [m1,m2,m5,m6] has rank 1.
// We have already considered [m1,m2,m5,m6]=[0,0,0,0].

// Thus this completes the case [m9,m11,m12] = [0,1,0].


///////////////////
// Option 4: [m9,m11,m12]=[-1,1,0]
///////////////////

// Now, we have [m9,m11,m12]=[-1,1,0]. We start by eliminating the degenerate cases
// for [m1,m2,m5,m6]. Thus m5=m6=0. We still have an element left from the centralizer,
// so may choose m7=1 or m7=0, m8=1.

coe1:=ChangeCoefficients(coeffs1,[m12,m5,m6],[0,0,0]);
rel1:=Reduce([Evaluate(i,coe1):i in relsR]);

coe1:=ChangeCoefficients(coe1,[m9,m11],[-1,1]);
rel1:=[Evaluate(i,coe1):i in rel1];

m7*(m7+m8) in rel1;
// If m7=0 then m8=1.
coe1a:=ChangeCoefficients(coe1,[m7,m8],[0,1]);
// If m7=-m8, then m7=1.
coe1b:=ChangeCoefficients(coe1,[m7,m8],[1,-1]);

rel1a:=[Evaluate(i,coe1a):i in rel1];
rel1b:=[Evaluate(i,coe1b):i in rel1];

GroebnerBasis(rel1a) eq [1];

// Now branch b.

&and[i*(n14+n16) in rel1b:i in [a,b,c]]; // Not all of a,b,c can be zero, so n14=-n16.

coe1b:=ChangeCoefficient(coe1b,n14,-n16);
rel1b:=[Evaluate(i,coe1b):i in rel1];

BB1b:=StripBasis(GroebnerBasis(ChangeOrder(ideal<R|rel1b>,"grevlex")),[]);

&and[i*j in BB1b:i in [a,b,c],j in [n2+n4,n3-n4,n15-n16]]; // Thus each of these is zero
n13 in BB1b;

coe1b:=ChangeCoefficients(coe1b,[n2,n3,n15,n13],[-n4,n4,n16,0]);
BB1b:=StripBasis(GroebnerBasis(ChangeOrder(ideal<R|[Evaluate(i,coe1b):i in BB1b]>,"grevlex")),[]);

n1 in BB1b; // Thus n1=n13=0.

[-n4,n4,n4,-n16,n16,n16] eq [coe1b[i]:i in [21,22,23,33,34,35]];

// Thus these two rows are scalar multiples, the determinant is zero, and this case ends.

// The final option is where [m1,m2,m5,m6] is non-degenerate, so [1,0,0,1].

coe2:=ChangeCoefficients(coeffs1,[m1,m2,m5,m6,m9,m11,m12],[1,0,0,1,-1,1,0]);
rel2:=([Evaluate(i,coe2):i in relsR]);

// It turns out that this can be easily solved using the Groebner Basis algorithm.

BB2:=GroebnerBasis(rel2);
[a,b,c] subset BB2; // These cannot all be zero, so this is a contradiction.


///////////////////
// Option 5: [m9,m11,m12]=[1,1,1]
///////////////////

// Now, we have [m9,m11,m12]=[1,1,1]. We start by eliminating the degenerate cases
// for [m1,m2,m5,m6]. Thus m5=m6=0. We still have an element left from the centralizer,
// so may choose m7=1 or m7=0, m8=1. First, make a cube root of unity.

ze:=ww^280; ze^3 eq 1;

coe1:=ChangeCoefficients(coeffs1,[m9,m11,m12,m5,m6],[t1,t1,t1,0,0]);
rel1:=Reduce([Evaluate(i,coe1):i in relsR]);

t1*(m7+ze*m8)*(m7+ze^2*m8) in rel1;
// Up to field automorphism, we may set m7=1 and m8=-ze.

coe1:=ChangeCoefficients(coeffs1,[m7,m8],[t1,-ze*t1]);
rel1:=[Evaluate(i,coe1):i in rel1];

&and[i*t1*(n14+ze*n16) in rel1:i in [a,b,c]];

// Since a,b,c cannot all be zero, we have n14=-ze*n16.

coe1:=ChangeCoefficient(coe1,n14,-ze*n16);

// That's enough playing around. Let's do this.
coe1:=ChangeCoefficient(coe1,t1,1);
rel1:=[Evaluate(i,coe1):i in rel1];

BB1:=StripBasis(GroebnerBasis(ChangeOrder(ideal<R|rel1>,"grevlex")),[]);

&and[i*(n2+ze*n4) in BB1:i in [a,b,c]];
&and[i*(n10+ze*n12) in BB1:i in [a,b,c]];

&and[i*(n3-ze^2*n4) in BB1:i in [a,b,c]];
&and[i*(n11-ze^2*n12) in BB1:i in [a,b,c]];

// Thus each of the bracketed terms is zero.
coe1:=ChangeCoefficients(coe1,[n2,n10,n3,n11],[-ze*n4,-ze*n12,ze^2*n4,ze^2*n12]);
BB1:=StripBasis(GroebnerBasis(ChangeOrder(ideal<R|[Evaluate(R!i,coe1):i in BB1]>,"grevlex")),[]);
n1 in BB1 and n9 in BB1;
coe1:=ChangeCoefficients(coe1,[n1,n9],[0,0]);
// The row conaining n1 is a scalar multiple of that containing n9.

[coe1[i]:i in [20..23]] eq [0,-ze*n4,ze^2*n4,n4];
[coe1[i]:i in [28..31]] eq [0,-ze*n12,ze^2*n12,n12];

// This collapses the case.

// Now assume that [m1,m2,m5,m6] is non-degenerate, so [1,0,0,1].

coe2:=ChangeCoefficients(coeffs1,[m9,m11,m12,m1,m2,m5,m6],[t1,t1,t1,t1,0,0,t1]);
rel2:=Reduce([Evaluate(i,coe2):i in relsR]);
coe2:=ChangeCoefficient(coe2,t1,1);
rel2:=[Evaluate(i,coe2):i in rel2];

BB2:=StripBasis(GroebnerBasis(ChangeOrder(ideal<R|rel2>,"grevlex")),[]);

(m7+ze*m8)*(m7+ze^2*m8) in BB2; // As before, we set m7=-ze*m8

coe2:=ChangeCoefficient(coe2,m7,-ze*m8);

BB2:=StripBasis(GroebnerBasis(ChangeOrder(ideal<R|[Evaluate(R!i,coe2):i in BB2]>,"grevlex")),[]);
BB2:=GroebnerBasis([R!i:i in BB2]);
coe2:=ChangeLinearCoefficients(coe2,BB2);
// The row conaining n1 is a scalar multiple of that containing n9.

[coe2[i]:i in [20..23]] eq [0,-ze^2*n4,ze*n4,n4];
[coe2[i]:i in [28..31]] eq [0,-ze^2*n12,ze*n12,n12];

// This collapses the case, and completes the proof.

return "";

end function;

// If you want to check the trilinear form, which has been used before, then here is some code.

function CheckE6Form()
qq:=29; FF:=GF(qq);
z:=PrimitiveElement(FF);
GG:=GroupOfLieType("E6",qq:Isogeny:="SC");
W:=VectorSpace(FF,6);
rho:=StandardRepresentation(GG);
Over:=GL(27,qq);
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>;

GL(27,FF)!l1 in NGT and GL(27,FF)!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:=[FF!0:i in [1..#seqs]];
    for i in [1..NumberOfRows(l2)] do if(h[aa,i] ne 0) then for j in [1..NumberOfRows(l2)] do if(h[bb,j] ne 0) then for k in [1..NumberOfRows(l2)] 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(FF,mat))).1;
return &and[ftest[i] eq f27[i]:i in [1..#seqs]];
end function;

"The function CheckRelations() checks that the reduced form of the relations is correct.";
"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.";
