/* 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,27). 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). This serves as a change-of-basis matrix, 

   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. After a number of
   computations we arrive at exactly two possibilities for the conjugate of h1. (We give the computations
   here, but of course we only give the precise terms that we need. A lot of the work was trial and
   error, looking at small polynomials from the 3654 available.) As a final check, we then test that the
   two matrices we are given really do leave the form invariant. We call these two matrices t1 and t2.

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

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

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

h1:=GL(27,F)![[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]];

// Here is a matrix that generates G2(3), together with H.

k1:=GL(27,F)![[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];

// A basis for the space of trilinear forms on H is given below.

vals1:=[1,2,2,2,2,2,2,2,2,1,1,1,2,1,2,1,2,1,1,2,1,2,1,2,1,1,1,1,2,2,1,2,1,2,2,2,2,1,2,2,2,2,2,1,2,1,2,1,
1,2,2,1,2,1,1,2,2,2,1,2,1,1,1,2,2,2,1,1,2,2,1,1,2,2,1,1,2,2,1,1,2,1,1,2,1,2,2,1,2,1,1,2,1,2,1,1,1,1,1,2,
1,1,2,1,1,1,1,2,2,1,1,1,1,1,1,1,2,1,1,2,1,1,1,2,1,1,2,1,1,1,1,1,2,1,2,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,1,2,
2,2,2,2,2,1,1,1,1,1,1,1];
vals2:=[1,1,1,1,1,1,1,2,2,2,1,1,1,1,2,2,2,1,2,1,2,2,2,2,1,1,1,2,1,2,2,1,2,2,1,1,2,2,1,2,2,1,2,2,1,2,2,2,
1,1,1,2,2,2,1,2,1,1,1,1,1,2,2,2,1,1,2,2,2,1,1,1,1,1,1,2,2,2,1,1,2,2,1,1,1,1,1,2,1,2,2,1,1,2,2,1,2,1,1,1,
2,2,2,2,1,2,1,2,2,1,2,2,2,1,1,2,2,1,2,1,2,1,2,2,1,2,2,1,2,2,1,2,2,2,1,1,2,1,1,1,1,2,2,2,1,2,2,2,1,2,2,1,
2,2,2,2,1,1,1,2,2,1,1,1,2,1,1,1,2,2,2,2,2,2,1,1,1];
vals3:=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,1,1,2,1,2,1,1,2,2,1,2,1,2,1,1,2,2,2,2,2,1,1,1,1,1,1,2,
1,1,2,1,1,1,1,1,2,2,1,1,1,1,2,1,1,1,2,2,1,1,1,1,2,1,1,1,1,2,1,1,2,1,1,2,1,1,1,2,1,1,1,2,1,1,2,1,2,2,2,1,
1,1,1,1,2,1,2,1,2,2,1,1,1,2,1,1,2,2,1,2,1,1,1,2,2,1,2,2,1,2,2,2,1,2,1,1,2,1,1,2,1,2,1,1,1,1,1,2,2,2,1,1,
2,2,2,1,2,2,2,2,1,2,1,1,2,1,1,2,1,2,2,2,2,1];

ents1:=[3,28,29,95,119,139,158,179,196,225,240,257,270,273,282,285,296,299,446,470,509,530,573,576,588,591,605,621,624,633,636,768,
771,792,795,812,815,834,837,855,858,872,875,898,901,916,930,933,946,949,958,961,972,975,1075,1078,1083,1104,1141,1144,1146,1192,
1202,1205,1229,1271,1284,1372,1375,1379,1402,1435,1438,1473,1483,1495,1498,1521,1572,1648,1651,1654,1708,1711,1713,1716,1744,1753,
1789,1833,1901,1904,1955,1962,1965,1989,2031,2074,2129,2132,2167,2180,2189,2192,2215,2293,2339,2342,2344,2347,2375,2387,2420,2488,
2534,2614,2617,2620,2623,2644,2647,2650,2779,2782,2784,2787,2806,2809,2811,2897,2900,2926,2929,2930,2950,2953,3029,3032,3059,3145,
3148,3172,3213,3216,3246,3307,3316,3334,3359,3362,3381,3392,3414,3417,3446,3475,3478,3525,3528,3550,3553,3584,3587,3606];
ents2:=[54,95,119,139,158,179,196,225,240,257,270,282,296,380,443,446,467,470,487,509,512,533,547,550,573,576,588,591,605,621,624,
633,636,650,768,771,792,834,837,855,858,901,913,916,930,933,946,949,958,972,975,1078,1083,1104,1144,1146,1192,1195,1205,1232,
1240,1243,1250,1284,1314,1317,1375,1379,1402,1438,1476,1483,1486,1498,1528,1531,1540,1572,1599,1602,1651,1654,1711,1716,1747,1753,
1756,1764,1795,1798,1833,1857,1860,1952,1955,1962,1965,1989,2000,2010,2031,2034,2044,2047,2074,2077,2167,2177,2180,2189,2192,2215,
2218,2222,2234,2266,2269,2293,2296,2344,2347,2375,2384,2387,2395,2420,2423,2426,2464,2467,2488,2491,2537,2552,2614,2617,2704,2779,
2782,2926,2929,3026,3032,3053,3059,3143,3148,3167,3172,3216,3245,3246,3266,3307,3316,3334,3359,3362,3368,3381,3392,3414,3417,3429,
3446,3475,3478,3487,3528,3534,3553,3565,3587,3596,3609,3612,3626];
ents3:=[92,95,116,119,136,139,158,161,179,182,196,199,222,237,254,273,285,299,405,467,470,487,490,509,530,533,550,573,605,624,650,
771,792,815,834,858,898,901,913,916,933,946,949,958,961,972,1078,1080,1083,1104,1141,1146,1192,1195,1202,1205,1232,1240,1243,
1250,1284,1314,1317,1375,1379,1402,1435,1476,1483,1486,1495,1498,1528,1531,1540,1572,1599,1602,1651,1654,1708,1713,1716,1747,1753,
1756,1764,1795,1798,1833,1857,1860,1906,1955,1962,1965,2000,2010,2031,2034,2047,2074,2077,2180,2189,2192,2215,2218,2222,2234,2269,
2293,2296,2344,2347,2387,2395,2420,2423,2426,2467,2488,2491,2534,2582,2617,2623,2721,2782,2787,2856,2900,2929,2991,3026,3053,3059,
3143,3167,3172,3245,3246,3266,3304,3307,3316,3334,3359,3362,3368,3381,3392,3414,3417,3429,3446,3475,3478,3487,3500,3525,3528,3534,
3550,3553,3565,3584,3587,3596,3606,3618,3631,3640,3649];

f1:=[F3!0:i in [1..#seqs]];
for i in [1..#ents1] do f1[ents1[i]]:=vals1[i]; end for;
f2:=[F3!0:i in [1..#seqs]];
for i in [1..#ents2] do f2[ents2[i]]:=vals2[i]; end for;
f3:=[F3!0:i in [1..#seqs]];
for i in [1..#ents3] do f3[ents3[i]]:=vals3[i]; end for;

function CheckPSL33Forms()
// Build the form space for H=PSL(3,3), which should be 3-dimensional.

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;

h:=h1;
repeat
  for coun in [1..33] 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;
until Dimension(Nullspace(Transpose(Matrix(F,mat)))) eq 3;

Nu:=Nullspace(Transpose(Matrix(F,mat)));
return &and[f1[i] eq (Nu.1)[i]:i in [1..#seqs]],&and[f2[i] eq (Nu.2)[i]:i in [1..#seqs]],&and[f3[i] eq (Nu.3)[i]:i in [1..#seqs]];
end function;

/* 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:=[F3!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,al,be,ga,de>:=PolynomialRing(F,23);
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;

// The general function.
fgen:=[al*R!f1[i]+be*R!f2[i]+ga*R!f3[i]:i in [1..#f1]];


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);

function TestRel(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;
e1:=&+[u1[1,i]*v1[1,j]*w1[1,k]*fgen[Position(seqs,Sort([i,j,k]))]:i,j,k in [1..NumberOfRows(l1)]];
return e1-f27[Position(seqs,seq)];
end function;

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

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;


// We will need to extend the field to test the solutions, although their determination can be done
// over F27. Thus we need functions over F729 as well.

S<az,bz,cz,m1z,m2z,m3z,m4z,m5z,m6z,m7z,m8z,n1z,n2z,n3z,n4z,n5z,n6z,n7z,n8z,alz,bez,gaz,dez>:=PolynomialRing(F729,23);

matsS:=MatrixRing(S,NumberOfRows(l1));
scal2:=matsS!scal;

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

V2:=ChangeRing(GModule(L),F729);
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 CheckRel(seq,gg)
u:=V2.seq[1]; v:=V2.seq[2]; w:=V2.seq[3];
u1:=MatriciseTwo(u); v1:=MatriciseTwo(v); w1:=MatriciseTwo(w);
u2:=MatriciseTwo(u)*GL(27,3^6)!gg; v2:=MatriciseTwo(v)*GL(27,3^6)!gg; w2:=MatriciseTwo(w)*GL(27,3^6)!gg;
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 CheckDetermination()

// Since we may multiply scal by any scalar matrix, we may assume that a=1.

Factorization(TestRel(seqs[3])) eq [<a^3*al + 2, 1>];
// Thus we see that al=1.

coeffs1:=[1,b,c,m1,m2,m3,m4,m5,m6,m7,m8,n1,n2,n3,n4,n5,n6,n7,n8,1,be,ga,de];

[Evaluate(TestRel(seqs[i]),coeffs1):i in [1080,1906,3304,3500]] eq [2*m1^3*ga + m2^3*ga + 2,2*m3^3*ga + m4^3*ga + 2,n1^3*ga + 2*n2^3*ga + 2,n3^3*ga + 2*n4^3*ga + 2];
// The first yields ga(m2^3-m1^3)=1, so m2-m1=de. Note also that from now on we may assume that ga is non-zero!
// The second yields ga(m4^3-m3^3)=1, so m4-m3=de.
// The third yields ga(n1^3-n2^3)=1, so n1-n2=de.
// The fourth yields ga(n3^3-n4^3)=1, so n3-n4=de.

[Evaluate(TestRel(seqs[i]),coeffs1):i in [54,405]] eq [-b^2+b-c+be,-b^3+b^2+c-ga];
// This gives us beta and gamma in terms of b and c, but this makes things more complicated and we prefer to stick with beta and gamma for now.

[Evaluate(TestRel(seqs[i]),coeffs1):i in [2534,2991,3606,3649]] eq
[2*m5^3*ga+m5^3+2*m6^3*ga+1,2*m7^3*ga+m7^3+2*m8^3*ga+1,2*n5^3*ga+n5^3+n6^3*ga+2,2*n7^3*ga+n7^3+n8^3*ga+1];
// The first yields ga(m5^3+m6^3)=1+m5^3, so m5+m6=de(1+m5)
// The second yields ga(m7^3+m8^3)=1+m7^3, so m7+m8=de(1+m7)
// The third yields ga(-n5^3+n6^3)=1-n5^3, so n6-n5=de(1-n5)
// The fourth yields ga(-n7^3+n8^3)=-1-n7^3, so n7-n8=de(1+n7)

// Thus we have eliminated eight of the mi and ni, and expressed be and ga in terms of b and c. This eliminates all of the 'easy' variables.
// We have to work for the others.

coeffs2:=ChangeCoefficients(coeffs1,[m2,m4,m6,m8,n1,n3,n6,n8],[m1+de,m3+de,m5*de+de-m5,de+de*m7-m7,n2+de,n4+de,n5+de-n5*de,n7-de-n7*de]);

aa1:=Evaluate(TestRel(seqs[3334]),coeffs2);
aa2:=Evaluate(TestRel(seqs[3307]),coeffs2);
Factorization(aa1+aa2+ga*de^3) eq [<de, 1>,<ga + 1, 1>,<n2 + 2*n4, 2>];
// This means that ga*de^2, which is 1/de, is equal to (ga+1)(n2-n4)^2, so de(ga+1)(n2-n4)^2=1. In particular, n2 ne n4 and ga ne -1.
Factorization(aa1-aa2) eq [<de, 1>,<n2 + 2*n4, 1>,<n2*ga + n2 + n4*ga + n4 + 2*be*de + 2*ga*de + de, 1>];
// Thus n2*ga + n2 + n4*ga + n4 + 2*be*de + 2*ga*de + de=0.

aa3:=Evaluate(TestRel(seqs[1146]),coeffs2);
aa4:=Evaluate(TestRel(seqs[1083]),coeffs2);
Factorization(aa3+aa4+ga*de^3) eq [<de, 1>,<ga + 1, 1>,<m1 + 2*m3, 2>];
Factorization(aa3-aa4) eq [<de, 1>,<m1 + 2*m3, 1>,<m1*ga + m1 + m3*ga + m3 + 2*be*de + 2*ga*de + de, 1>];
// This means that, as before, 1/de = (ga+1)(m1-m3)^2, so de(ga+1)(m1-m3)^2=1. In particular, m1 ne m3.
// Thus m1*ga + m1 + m3*ga + m3 + 2*be*de + 2*ga*de + de=0.

bb1:=Factorization(aa3-aa4)[3,1];
bb2:=Factorization(aa1-aa2)[3,1];
Factorization(bb1-bb2) eq [<ga + 1, 1>,<m1 + m3 + 2*n2 + 2*n4, 1>];
// Hence m1+m3=n2+n4
// We also have (m1-m3)^2=(n2-n4)^2, since both are 1/de(ga+1). Thus m1-m3=\pm(n2-n4), and together with m1+m3=n2+n4, this yields
// 2*m1=2*n2, or 2*m1=2*n4. Thus either m1=n2, m3=n4, or m1=n4, m3=n2. We will come back to this.

Evaluate(TestRel(seqs[161]),coeffs2)-Evaluate(TestRel(seqs[92]),coeffs2)
-Evaluate(TestRel(seqs[158]),coeffs2)+Evaluate(TestRel(seqs[95]),coeffs2) eq
(be-ga-1)*(n2-n4)*(m1+m3+de);

//We cannot have n2=n4 by the above.

// Suppose that be=ga+1.

badcoeffs3:=ChangeCoefficient(coeffs2,be,ga+1);
badcoeffs3 eq [1,b,c,m1,m1+de,m3,m3+de,m5,m5*de+de-m5,m7,de+de*m7-m7,n2+de,n2,n4+de,n4,n5,n5+de-n5*de,n7,n7-de-n7*de,1,ga+1,ga,de];
Evaluate(TestRel(seqs[92]),badcoeffs3) eq 2*ga*de^2;
// So this definitely cannot work. Thus we always have that m1+m3+de=0. Since m1+m3=n2+n4, this implies n2+n4+de=0. so n4=-n2-de

coeffs3:=ChangeCoefficients(coeffs2,[m3,n4],[-m1-de,-n2-de]);
Evaluate(bb1,coeffs3) eq 2*be*de + ga*de;
//Thus be=ga. 

coeffs4:=ChangeCoefficient(coeffs3,be,ga);
// Go back to our original statement, m1=n2,m3=n4 or m1=n4,m3=n2.
// Thus we obtain two cases: n2=m1, or n2=-m1-de. Test the first case first.
badcoeffs5:=ChangeCoefficient(coeffs4,n2,m1);
// We prove that this cannot be correct.
Evaluate(TestRel(seqs[1284]),badcoeffs5)+Evaluate(TestRel(seqs[2077]),badcoeffs5)-(n5+n7)*Evaluate(TestRel(seqs[95]),badcoeffs5) eq n5+n7;
// Looking at the determinant of scal:
Determinant(Evaluate(scal,coeffs4)) eq de^12*(n5+n7)^3*(n2-de)^3*(m5-m7)^3*(m1-de)^3;
// we see that n5+n7 ne 0, so not both of them can be zero. Hence this yields a contradiction.

// Thus we are left with the following.
coeffs5:=ChangeCoefficient(coeffs4,n2,-m1-de);

// Suppose that de=m1*ga

badcoeffs6:=ChangeCoefficient(coeffs5,de,m1*ga);
Evaluate(TestRel(seqs[1205]),badcoeffs6) eq 1;
// That's very wrong. So m1*ga-2*de cannot equal 0. But there are many coefficients that factorize with m1*ga-2*de in them.

Factorization(Evaluate(TestRel(seqs[2614]),coeffs5)+Evaluate(TestRel(seqs[2617]),coeffs5)) eq [<de, 1>,<n5 + n7, 1>,<m5, 1>,<m1*ga + 2*de, 1>];
// but remember our determinant condition:
Factorization(Determinant(Evaluate(scal,coeffs5))) eq [<de, 12>,<n5 + n7, 3>,<m5 + 2*m7, 3>,<m1 + 2*de, 6>];
// Thus m5=0. We do something similar:
Evaluate(TestRel(seqs[1989]),coeffs5)-Evaluate(TestRel(seqs[2031]),coeffs5) eq -de*n5*(m5-m7)*(m1*ga-de);
// Thus n5=0 as well.
Evaluate(TestRel(seqs[2034]),coeffs5) eq -de*n7*(m7+1)*(m1*ga-de);
Evaluate(TestRel(seqs[3534]),coeffs5) eq de*n7*(n7+1)*(m1*ga-de);
// Thus, since n7 cannot also equal zero, we have m7=n7=-1

coeffs6:=ChangeCoefficients(coeffs5,[m5,n5,m7,n7],[0,0,-1,-1]);
Evaluate(TestRel(seqs[624]),coeffs6) eq b;
// Thus b=0. But we previously computed gamma in terms of b and c, -b^3+b^2+c-ga=0. Thus c=ga as well.
coeffs7:=ChangeCoefficients(coeffs6,[b,c],[0,ga]);

// Now we need to deal with m1, at last, and then we evaluate ga.
Evaluate(TestRel(seqs[792]),coeffs7) eq (m1+ga*de-de)*(m1-ga*de-de);

// If m1=ga*de+de, this is bad.
badcoeffs8:=ChangeCoefficient(coeffs7,m1,de+ga*de);
Evaluate(TestRel(seqs[795]),badcoeffs8) eq 2;
// So we know that m1=de-ga*de.
coeffs8:=ChangeCoefficient(coeffs7,m1,de-ga*de);
Evaluate(TestRel(seqs[92]),coeffs8) eq -de^2*(ga^2+ga-1);
// This means that gamma is a solution to gamma^2+gamma-1=0. These exist in GF(9), and there are of course two of them. We therefore
// make two sets of coefficients, depending on which root we take. They both happen to satisfy de=-ga.
// So our final set of evaluations, before plugging in the value of ga, is
coeffs9:=ChangeCoefficient(coeffs8,de,-ga);
coeffs9 eq [1,0,ga,ga^2-ga,ga^2+ga,-ga^2-ga,-ga^2+ga,0,-ga,2,1,-ga^2+ga,-ga^2-ga,ga^2+ga,ga^2-ga,0,-ga,2,2,1,ga,ga,-ga];

// Now we test that these are solutions, and that they are contained in copies of G2(3). First, it suffices to check that k1
// conjugates into E6.

J:=sub<GL(27,F)|l1,l2,k1>;
h1 in J;

// Write coeffs9 over the field extension.

c9:=[1,0,gaz,gaz^2-gaz,gaz^2+gaz,-gaz^2-gaz,-gaz^2+gaz,0,-gaz,2,1,-gaz^2+gaz,-gaz^2-gaz,gaz^2+gaz,gaz^2-gaz,0,-gaz,2,2,1,gaz,gaz,-gaz];
c101:=[Evaluate(i,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-w2^91,0]):i in c9];
c102:=[Evaluate(i,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-w2^273,0]):i in c9];

scala:=Evaluate(scal2,c101);
scalb:=Evaluate(scal2,c102);

t11:=scala*GL(27,3^6)!k1*scala^-1;
t22:=scalb*GL(27,3^6)!k1*scalb^-1;

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;

// Confirm that L is a subgroup of E6 via conjugation by x:
x1*x2*l1*x2^-1*x1^-1 in NGT and 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 CheckPSL33Forms() checks that the trilinear forms on H are correct. (This takes about five minutes.)";
printf "\n";
"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, and proves that there are copies of G2(3) containing them. (This takes about ten minutes.)";