/* Here we let H denote a copy of PSL(3,3) and L denote a copy of 13.3, the normalizer of a Phi3-torus.
   This is unique up to conjugacy in NG(T). First we set up L and H. L requires two matrices, and one
   more for H. We give them as elements of GL(27,729). We construct L so it is very close to belonging to
   the Magma version of E6, with only a slight tweak required. This means that the trilinear form is
   very compact, and will help us with our calculations. We offer code to confirm that the trilinear form
   is as suggested, using the function CheckE6Form().

   We then construct the centralizer of L in GL(27,F), and check that it indeed centralizes L. Since it
   has the correct number of parameters, this means that it must be the whole centralizer. (And it is
   clear from the module structure of L on the 27-space.)

   We can then launch into working out the coefficients necessary for the conjugate of h1 (the matrix for
   H) to lie in our copy of E6, i.e., for the trilinear form to be h1-invariant. This version of the file
   uses Groebner bases to do this, rather than working by hand. We obtain exactly two solutions.

   The final step is to give two matrices, t11 and t22, which together with L generate a copy of G2(3).
   We check that t11 and t22 also leave the form invariant, proving first that the two solutions above are
   genuine solutions, and second that all copies of PSL(3,3) (acting irreducibly on M(E6)) lie in copies
   of G2(3), as claimed in the article. */

F729<w2>:=GF(729);
F<w>:=GF(27);

l1:=GL(27,F729)![[w^2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,w^24,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,w^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,w^14,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,w^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,w^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,w^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,w^8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,w^20,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,w^18,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,w^22,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,w^14,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,w^20,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,w^18,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,w^16,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,w^16,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,w^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,w^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,w^22,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,w^24,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,w^2,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,w^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,w^8,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,w^6]];

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

x2:=PermutationMatrix(F729,Sym(27)!(1,4,12,27,5,22,3)(2,17,26,16,23,20,10,21,24,7,13,25,14,18)(6,8,15)(9,19,11));

l1:=x2^-1*l1*x2;
l2:=x2^-1*l2*x2;
L:=sub<GL(27,F729)|l1,l2>;

h1:=GL(27,F729)![[0,0,0,2,2,2,0,0,0,1,1,1,0,0,0,1,1,1,0,0,0,1,1,1,0,0,0],
[0,0,0,2,0,1,1,1,1,0,2,1,1,1,1,1,0,2,1,1,1,2,1,0,0,0,0],
[2,0,0,0,0,2,2,1,0,0,1,0,0,2,1,0,0,1,2,1,0,2,1,1,1,1,1],
[1,2,0,0,0,2,0,1,1,0,2,0,0,2,2,0,1,0,0,0,1,0,1,0,1,1,1],
[0,2,0,2,0,0,1,0,1,0,0,2,2,0,2,0,0,1,1,0,0,0,0,1,1,1,1],
[2,2,0,0,2,0,1,1,0,2,0,0,2,2,0,1,0,0,0,1,0,1,0,0,1,1,1],
[1,2,2,1,2,2,0,2,0,0,2,0,0,1,2,2,1,0,0,1,0,2,2,2,0,2,2],
[0,1,2,2,1,2,0,0,2,0,0,2,2,0,1,0,2,1,0,0,1,2,2,2,2,0,2],
[1,0,2,2,2,1,2,0,0,2,0,0,1,2,0,1,0,2,1,0,0,2,2,2,2,2,0],
[1,1,0,0,1,2,0,1,1,2,1,0,2,0,0,1,0,2,1,2,2,0,2,2,1,0,1],
[2,1,0,2,0,1,1,0,1,0,2,1,0,2,0,2,1,0,2,1,2,2,0,2,1,1,0],
[0,1,0,1,2,0,1,1,0,1,0,2,0,0,2,0,2,1,2,2,1,2,2,0,0,1,1],
[1,2,2,0,1,0,0,1,2,2,1,0,1,0,1,0,0,2,0,2,2,0,2,0,2,2,2],
[0,1,2,0,0,1,2,0,1,0,2,1,1,1,0,2,0,0,2,0,2,0,0,2,2,2,2],
[1,0,2,1,0,0,1,2,0,1,0,2,0,1,1,0,2,0,2,2,0,2,0,0,2,2,2],
[1,2,0,0,0,1,0,1,0,1,0,0,2,1,1,0,2,0,0,1,1,0,2,0,1,2,2],
[0,2,0,1,0,0,0,0,1,0,1,0,1,2,1,0,0,2,1,0,1,0,0,2,2,1,2],
[2,2,0,0,1,0,1,0,0,0,0,1,1,1,2,2,0,0,1,1,0,2,0,0,2,2,1],
[2,1,1,2,0,1,0,0,1,1,0,0,2,1,0,1,2,2,0,0,2,1,1,1,1,1,0],
[0,2,1,1,2,0,1,0,0,0,1,0,0,2,1,2,1,2,2,0,0,1,1,1,0,1,1],
[2,0,1,0,1,2,0,1,0,0,0,1,1,0,2,2,2,1,0,2,0,1,1,1,1,0,1],
[1,0,0,0,2,2,1,1,1,2,2,1,2,2,0,1,0,1,1,2,2,0,1,1,2,0,0],
[1,0,0,2,0,2,1,1,1,1,2,2,0,2,2,1,1,0,2,1,2,1,0,1,0,2,0],
[1,0,0,2,2,0,1,1,1,2,1,2,2,0,2,0,1,1,2,2,1,1,1,0,0,0,2],
[2,0,1,2,0,0,2,2,0,1,0,2,1,2,1,0,1,0,2,1,0,1,0,0,1,1,1],
[2,1,1,0,2,0,0,2,2,2,1,0,1,1,2,0,0,1,0,2,1,0,1,0,1,1,1],
[0,2,1,0,0,2,2,0,2,0,2,1,2,1,1,1,0,0,1,0,2,0,0,1,1,1,1]];

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

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

/* We next give the symmetric trilinear form that comes from E6. This was found by choosing elements of E6, 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.
*/
vals:=[1,2,2,1,1,1,1,1,1,2,2,2,1,1,1,2,1,1,2,1,2,2,1,1,1,2,1,1,2,2,1,2,2,1,2,2,1,1,2,1,2,1,2,1,2,1,1,1,2,1,1,1,
1,1,1,1,1,1,2];
ents:=[3,28,29,95,119,139,158,179,196,222,237,254,273,285,299,470,490,509,547,573,588,636,650,795,872,898,961,1080,1141,1205,1287,
1435,1498,1575,1708,1716,1836,1906,1992,2003,2170,2225,2378,2429,2534,2623,2787,2900,2991,3023,3140,3242,3304,3500,3528,3553,3587,
3606,3649];

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

// Now set up the centralizer of L in GL(27,F), which is given by all possible matrices scal below.
// Since the centralizer is 19-dimensional, we need 19 parameters in a polynomial ring. We also need three parameters for a general
// trilinear form on H, namely fgen=al*f1+be*f2+ga*f3. Finally, in the proof it will become clear that we need the cube root of ga^-1.
// We add one more parameter, de, under the standing assumption that de^3*ga=1. Note that we do not take the quotient ring, as then
// you cannot factorize.

R<a,b,c,m1,m2,m3,m4,m5,m6,m7,m8,n1,n2,n3,n4,n5,n6,n7,n8>:=PolynomialRing(F729,19);
mats:=MatrixRing(R,NumberOfRows(l1));
scal:=mats!DiagonalMatrix([a,a,a,m1,m1,m1,m4,m4,m4,m5,m5,m5,m8,m8,m8,n1,n1,n1,n4,n4,n4,n5,n5,n5,n8,n8,n8]);
scal[2,1]:=b; scal[3,2]:=b; scal[3,1]:=c;
for i in [1..3] do
  scal[i+3,i+6]:=m2; scal[i+6,i+3]:=m3;
  scal[i+9,i+12]:=m6; scal[i+12,i+9]:=m7;
  scal[i+15,i+18]:=n2; scal[i+18,i+15]:=n3;
  scal[i+21,i+24]:=n6; scal[i+24,i+21]:=n7;
end for;

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

// Make a vector into a matrix.

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

// This computes the difference between the form values on three vectors u,v,w, under multiplication by scal and by g*scal.
// This is zero if g^scal lies in E6.

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];
u1:=u; v1:=v; w1:=w;
u2:=u^gg; v2:=v^gg; w2:=w^gg;
e1:=&+[u1[i[1]]*v1[i[2]]*w1[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;

coeffs0:=[a,b,c,m1,m2,m3,m4,m5,m6,m7,m8,n1,n2,n3,n4,n5,n6,n7,n8];

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;

// This function automatically replaces all variables that appear as the leading term
// in a linear equation by the rest of the expression.

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;

// This function removes all powers from all polynomials in B, and removes all factors from
// polynomials in B that lie in X. (Here X should be a set of polynomials known to be
// non-zero.)

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 CheckDetermination()

rels:=[ProdRel(i):i in seqs];

B:=GroebnerBasis(ChangeOrder(ideal<R|rels>,"grevlex"));
B:=SimplifyBasis(B,[a]); B:=GroebnerBasis(B);
B:=SimplifyBasis(B,[a]); B:=GroebnerBasis(B);
coeffs1:=ChangeLinearCoefficients(coeffs0,[R!i:i in B]);

Determinant(Evaluate(scal,coeffs1)) eq -a^3*(m2+m4)^3*(n7-n8)^6*(n2+n4)^3*(n5*n8-n6*n7)^3*(m5*m8+m6*m8-m6*n7+m6*n8)^3;

// Thus each of these factors is non-zero.

B2:=GroebnerBasis(SimplifyBasis([Evaluate(i,coeffs1):i in B],[a,m2+m4,n7-n8,n2+n4]));
n8 in B2; // Thus this is zero, and we need to alter our non-zero terms to include n7.

coeffs2:=ChangeCoefficient(coeffs1,n8,0);
B2:=GroebnerBasis(SimplifyBasis([Evaluate(i,coeffs2):i in B2],[a,m2+m4,n7,n2+n4]));
// This now has a single non-linear equation.
coeffs2:=ChangeLinearCoefficients(coeffs2,B2);
B2:=GroebnerBasis([Evaluate(R!i,coeffs2):i in B2]);

B2 eq [(n6 + w2^91*n7)*(n6 + w2^273*n7)];

// In the direct version of this we set a=1. Let's do that, so n6=2. For n7, we use the first equation above.

coeffs3a:=ChangeCoefficients(coeffs2,[n6,n7],[2,1/w2^91]);
coeffs3b:=ChangeCoefficients(coeffs2,[n6,n7],[2,1/w2^273]);
scala:=GL(27,F729)!Evaluate(scal,coeffs3a);
scalb:=GL(27,F729)!Evaluate(scal,coeffs3b);

J:=sub<GL(27,F729)|l1,l2,k1>;
h1 in J;
// So we only need to check k1.

t11:=k1^scala;
t22:=k1^scalb;

return {CheckRel(i,t11):i in seqs} eq {0} and {CheckRel(i,t22):i in seqs} eq {0};

end function;

function CheckE6Form()
GG:=GroupOfLieType("E6",27);
W:=VectorSpace(F,6);
rho:=StandardRepresentation(GG);
Over:=GL(27,F);
g1:=elt<GG|W![w^2,1,1,1,1,1]>;
g2:=elt<GG|W![1,w^2,1,1,1,1]>;
g3:=elt<GG|W![1,1,w^2,1,1,1]>;
g4:=elt<GG|W![1,1,1,w^2,1,1]>;
g5:=elt<GG|W![1,1,1,1,w^2,1]>;
g6:=elt<GG|W![1,1,1,1,1,w^2]>;
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>;

x1:=IdentityMatrix(F,27); for i in [4,10,13,14,15,24] do x1[i,i]:=2; end for;
x1[3,3]:=0; x1[3,18]:=0; x1[3,22]:=1;
x1[18,3]:=0; x1[18,18]:=1; x1[18,22]:=1;
x1[22,3]:=2; x1[22,18]:=1; x1[22,22]:=2;
x:=GL(27,F)!x1*GL(27,F)!x2;
// Confirm that L is a subgroup of E6 via conjugation by x^-1:
GL(27,F)!(x1*x2*l1*x2^-1*x1^-1) in NGT and GL(27,F)!(x1*x2*l2*x2^-1*x1^-1) in NGT;

mat:=[];
for h in [l1,l2] do
  for nn in [1..#seqs] do aa:=seqs[nn,1]; bb:=seqs[nn,2]; cc:=seqs[nn,3];
    val:=[F!0:i in [1..#seqs]];
    for i in [1..NumberOfRows(h)] do if(h[aa,i] ne 0) then for j in [1..NumberOfRows(h)] do if(h[bb,j] ne 0) then for k in [1..NumberOfRows(h)] do
      val[Position(seqs,Sort([i,j,k]))]+:=h[aa,i]*h[bb,j]*h[cc,k];
    end for; end if; end for; end if; end for;
    val[nn]-:=1; Append(~mat,val); delete val;
  end for;
end for;

for g in [NGT.1,NGT.6] do
  h:=x2^-1*(x1^-1*g*x1)*x2;
  for nn in [1..#seqs] do aa:=seqs[nn,1]; bb:=seqs[nn,2]; cc:=seqs[nn,3];
    val:=[F!0:i in [1..#seqs]];
    for i in [1..NumberOfRows(h)] do if(h[aa,i] ne 0) then for j in [1..NumberOfRows(h)] do if(h[bb,j] ne 0) then for k in [1..NumberOfRows(h)] do
      val[Position(seqs,Sort([i,j,k]))]+:=h[aa,i]*h[bb,j]*h[cc,k];
    end for; end if; end for; end if; end for;
    val[nn]-:=1; Append(~mat,val); delete val;
  end for;
end for;

repeat
  g:=rho(Random(GG));
  h:=x2^-1*(x1^-1*g*x1)*x2;
  for coun in [1..32] do nn:=Random([1..#seqs]); aa:=seqs[nn,1]; bb:=seqs[nn,2]; cc:=seqs[nn,3];
    val:=[F!0:i in [1..#seqs]];
    for i in [1..NumberOfRows(h)] do if(h[aa,i] ne 0) then for j in [1..NumberOfRows(h)] do if(h[bb,j] ne 0) then for k in [1..NumberOfRows(h)] do
      val[Position(seqs,Sort([i,j,k]))]+:=h[aa,i]*h[bb,j]*h[cc,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 f27[i]:i in [1..#f27]];
end function;


"The function CheckE6Form() checks that the trilinear form f27 from E6 is correct. (This takes about three minutes.)";
printf "\n";
"The function CheckDetermination() checks that the determination of the possible conjugates into E6 is correct, using Groebner bases, and that there are copies of G2(3) containing them. (This takes about ten minutes.)";