/* Here we let H denote a copy of PSU(3,3) and L denote a copy of 4^2.S3, the normalizer of a Phi2-torus.
   This is unique up to conjugacy in NG(T). First we set up L and H. L requires three matrices, and one
   more for H. We give them as elements of GL(27,9). 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. We give two methods: the
   easiest (but not quickest) is to use Groebner bases, which immediately gives us the answer, but we
   have to compute the whole set of equations first. (It is likely we do not need to do this.) We also
   prove the same result using individual equations directly.

   In QuickCheckDetermination() we also prove that each solution is contained in a copy of G2(3). The
   original element h1 is contained in a G2(3) given by <L,k1>. Conjugating k1 yields two matrices t11
   and t22 that we show lie in E6. This proves that all copies of PSU(3,3) (acting irreducibly on M(E6))
   lie in copies of G2(3), as claimed in the article. */

q:=9;
F<w>:=GF(q);

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

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

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

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

h1:=GL(27,F)![[0,1,0,2,1,w^5,1,w^7,2,1,2,2,1,1,1,0,w^2,2,w^6,0,1,1,w^5,w^3,w,1,w^7],
[1,0,1,0,0,w,0,w^3,0,2,1,2,2,2,1,0,0,0,0,0,0,0,0,0,0,0,0],
[0,1,0,1,2,w^5,2,w^7,1,1,2,2,1,1,1,0,w^6,1,w^2,0,2,1,w,w^7,w^5,1,w^3],
[1,0,2,2,2,0,1,0,1,w,w^3,w^6,w^3,w^5,w^6,w,w^6,1,w^6,w^5,1,2,w^5,0,w^5,1,0],
[2,0,1,2,2,0,1,0,1,w^5,w^7,w^2,w^7,w,w^2,w^5,w^6,1,w^6,w,1,1,w^5,0,w^5,2,0],
[w^3,w^7,w^3,0,0,0,0,w^6,0,w^6,1,w^7,2,w^6,w^3,1,0,0,0,1,0,0,0,0,0,0,0],
[2,0,1,1,1,0,2,0,2,w^7,w^5,w^6,w^5,w^3,w^6,w^5,w^2,2,w^2,w,2,2,0,w^7,0,1,w^7],
[w,w^5,w,0,0,w^2,0,0,0,w^2,1,w^5,2,w^2,w,w^2,0,0,0,w^2,0,0,0,0,0,0,0],
[1,0,2,1,1,0,2,0,2,w^3,w,w^2,w,w^7,w^2,w,w^2,2,w^2,w^5,2,1,0,w^7,0,2,w^7],
[2,1,2,w^3,w^7,w^2,w^5,w^6,w,0,2,1,1,0,2,w^3,w^5,w^5,w,w^3,w,0,2,1,1,0,2],
[1,2,1,w,w^5,1,w^7,1,w^3,2,0,2,0,2,1,w^3,w^3,w^7,w^7,w^3,w^3,0,w^2,w^2,w^6,0,w^6],
[1,1,1,w^2,w^6,w^5,w^2,w^7,w^6,1,2,1,1,1,2,0,0,0,0,0,0,1,w^7,w,w^3,1,w^5],
[2,1,2,w,w^5,2,w^7,2,w^3,1,0,1,0,1,2,w^7,w^3,w^7,w^7,w^7,w^3,0,w^2,w^2,w^6,0,w^6],
[2,1,2,w^7,w^3,w^2,w,w^6,w^5,0,2,1,1,0,2,w^3,w,w,w^5,w^3,w^5,0,1,2,2,0,1],
[2,2,2,w^2,w^6,w,w^2,w^3,w^6,2,1,2,2,2,1,0,0,0,0,0,0,2,w^7,w,w^3,2,w^5],
[0,0,0,w^7,w^3,2,w^3,w^2,w^7,w^5,w^5,0,w,w^5,0,0,0,0,0,0,0,0,2,w^2,1,0,w^6],
[w^6,0,w^2,w^6,w^6,0,w^2,0,w^2,w^3,w^5,0,w^5,w^7,0,0,1,w^6,1,0,w^6,w^2,w^7,w^5,w^7,w^6,w^5],
[2,0,1,2,2,0,1,0,1,w^3,w,0,w,w^7,0,0,w^2,1,w^2,0,1,2,w^5,w^3,w^5,1,w^3],
[w^2,0,w^6,w^6,w^6,0,w^2,0,w^2,w^7,w,0,w,w^3,0,0,1,w^6,1,0,w^6,w^6,w^7,w^5,w^7,w^2,w^5],
[0,0,0,w^3,w^7,2,w^7,w^2,w^3,w^5,w^5,0,w,w^5,0,0,0,0,0,0,0,0,1,w^6,2,0,w^2],
[1,0,2,2,2,0,1,0,1,w^7,w^5,0,w^5,w^3,0,0,w^2,1,w^2,0,1,1,w^5,w^3,w^5,2,w^3],
[1,0,1,1,2,0,1,0,2,0,0,2,0,0,1,0,w^6,2,w^2,0,1,1,w,w^3,w^5,1,w^7],
[w^7,0,w^3,w^3,w^3,0,0,0,0,1,w^2,w,w^2,2,w,2,w^5,w^7,w^5,1,w^7,w^3,2,w^2,2,w^7,w^2],
[w,0,w^5,0,0,0,w,0,w,2,w^2,w^7,w^2,1,w^7,w^6,w^7,w,w^7,w^2,w,w,w^6,2,w^6,w^5,2],
[w^3,0,w^7,w^3,w^3,0,0,0,0,2,w^6,w^5,w^6,1,w^5,1,w^5,w^7,w^5,2,w^7,w^7,2,w^2,2,w^3,w^2],
[1,0,1,2,1,0,2,0,1,0,0,2,0,0,1,0,w^2,1,w^6,0,2,1,w^5,w^7,w,1,w^3],
[w^5,0,w,0,0,0,w,0,w,1,w^6,w^3,w^6,2,w^3,w^2,w^7,w,w^7,w^6,w,w^5,w^6,2,w^6,w,2]];

H:=sub<GL(27,F)|l1,l2,l3,h1>;

k1:=GL(27,9)![[0,2,1,w^2,1,w^2,1,w^6,w^6,2,0,0,0,0,0,0,0,2,w^6,0,0,1,w,w,w^3,1,w^3],
[2,0,1,0,w^5,w^2,w^7,w^6,0,1,1,0,2,0,2,2,w,w^2,1,w^6,w^5,2,w^3,w,w^3,2,w],
[1,1,1,w^7,w^6,w^5,w^2,w^7,w^5,2,1,1,0,1,2,w^5,w^6,w^6,2,w,2,w^6,1,w^7,w^5,w^2,1],
[w^2,0,w,1,0,w,0,w^2,2,1,w^5,0,2,0,w^7,w^7,w^2,w^5,1,1,w^2,0,1,w,1,0,2],
[2,w^3,w^6,0,0,w^7,w^3,w^3,0,w^3,w^2,w^5,w^3,w^7,2,w,0,0,w,w^7,w^5,1,w^7,w^3,2,2,w^7],
[w^2,w^2,w^3,w^3,w^5,2,w^3,w^2,w^2,1,w^5,w^2,0,w^6,w^5,w,w^7,w^6,w^2,w^2,w^2,2,w,w,w^3,w^6,w^2],
[2,w,w^2,0,w,w,0,w^5,0,w,w^6,w^7,w,w^5,2,w^7,w,w^5,0,w^5,0,2,w^5,2,w,1,w^5],
[w^6,w^6,w,w^6,w,w^6,w^7,2,w,1,w^7,w^6,0,w^2,w^7,1,1,1,2,w^5,w^7,w^2,w^6,w,w^3,2,w^3],
[w^6,0,w^3,2,0,w^6,0,w^3,1,1,w^7,0,2,0,w^5,w^2,1,w^2,w,w^7,1,0,2,1,w^3,0,1],
[1,2,1,1,w,1,w^3,1,1,2,0,0,0,2,0,w^7,w^7,w^3,w^3,w^7,w^7,0,w^6,w^7,w^5,0,w^2],
[0,2,2,w^7,w^6,w^7,w^2,w^5,w^5,0,1,2,0,2,2,0,w^7,w^2,1,0,w^7,w^2,w^6,w^3,w,w^6,w^2],
[0,0,2,0,w^7,w^6,w^5,w^2,0,0,2,1,1,2,2,0,1,w^7,w^7,0,w^2,w^3,0,2,2,w,0],
[0,1,0,2,w,0,w^3,0,2,0,0,1,0,0,0,w^6,w^2,2,w^6,2,1,1,w^6,w^3,w,1,w^2],
[0,0,2,0,w^5,w^2,w^7,w^6,0,2,2,2,0,1,2,w^6,w^2,w^2,1,2,1,w^5,0,1,1,w^7,0],
[0,1,1,w^5,2,w^7,2,w^5,w^7,0,2,2,0,2,1,w^5,w^3,w^7,w^7,w,w^3,2,w^2,w^3,w,2,w^6],
[0,2,w^7,w,w^7,w^7,w,2,w^2,w,0,0,w^6,w^6,w^3,2,1,2,w^2,w^3,2,w^2,w^2,w^3,1,2,2],
[0,w^3,w^2,w^2,0,w,w^7,2,2,w,w,2,w^2,w^2,w^5,1,0,w,w^7,2,2,w^6,w^5,w^2,w^5,1,0],
[2,w^6,w^2,w^3,0,w^6,w^3,2,w^2,w^5,w^2,w,1,w^2,w,2,w^3,0,w^6,w^2,w^7,w^7,w^2,2,0,2,w^7],
[w^2,1,2,2,w^7,w^2,0,1,w^7,w^5,2,w,w^6,2,w,w^6,w^5,w^2,0,2,w,w^2,w^3,0,w^2,w^3,2],
[0,w^2,w^3,2,w,w^2,w^3,w^3,w,w,0,0,1,1,w^7,w,2,w^6,2,2,1,w^2,w^2,w^6,w^7,2,2],
[0,w^7,2,w^2,w^3,w^2,0,w,2,w,w,w^2,2,2,w^5,2,2,w^5,w^3,1,0,w^6,0,w^5,2,1,w^5],
[1,2,w^2,0,2,1,1,w^2,0,0,w^2,w^5,2,w^3,1,w^6,w^2,w^5,w^6,w^6,w^2,0,w,w^7,0,w,w],
[w^3,w,1,2,w,w^7,w^3,w^6,1,w^6,w^6,0,w^6,0,w^2,w^6,w^7,w^6,w,w^6,0,w^3,1,1,w^5,w,0],
[w^3,w^3,w^5,w^7,w^5,w^7,1,w^7,2,w,w^5,1,w^5,2,w^5,w,w^6,2,0,w^2,w^7,w^5,1,1,w^3,0,w^5],
[w,w,w^7,2,1,w^5,w^7,w^5,w^5,w^3,w^7,1,w^7,2,w^7,1,w^7,0,w^6,w^5,2,0,w^7,w,1,w^7,1],
[1,2,w^6,0,1,w^6,2,1,0,0,w^6,w^7,2,w,1,2,1,2,w,2,1,w^3,w^3,0,w^5,0,w^3],
[w,w^3,1,1,w,w^2,w^3,w^5,2,w^2,w^2,0,w^2,0,w^6,2,0,w^5,2,2,w^7,w^3,0,w^7,1,w,1]];

K:=sub<GL(27,F)|l1,l2,l3,k1>;

h1 in K;

// Now construct the trilinear form, then the conjugating element into E6

seqs:=[[i,j,k]:i,j,k in [1..NumberOfRows(l1)]|i le j and j le k];

ents:=[ 29, 105, 229, 311, 353, 479, 597, 674, 690, 760, 887, 1007, 1021, 1080, 1208,
1270, 1347, 1501, 1517, 1615, 1719, 1785, 1871, 1906, 1958, 2034, 2052, 2183,
2218, 2273, 2390, 2423, 2440, 2534, 2664, 2800, 2945, 2991, 3035, 3151, 3219,
3304, 3500, 3606, 3649 ];
vals:=[ 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 2, 2, 1, 2, 1, 1, 1, 1, 2, 1,
1, 1, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 1, 2, 1, 2, 2, 1, 1 ];

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

// Now we give the centralizer of L as a 10-parameter matrix scal.

R<a,b,c1,c2,c3,c4,m1,m2,m3,m4>:=PolynomialRing(F,10);
mats1:=MatrixRing(R,9);
mats2:=MatrixRing(R,6);
mats3:=MatrixRing(R,12);

mats:=MatrixRing(R,27);

scal1:=mats1!DiagonalMatrix([a+b,a+b,a+b,c1,c1,c1,c2,c2,c2]);
scal1[2,1]:=b; scal1[3,1]:=b; scal1[3,2]:=b; scal1[1,2]:=b; scal1[1,3]:=b; scal1[2,3]:=b; 
scal2:=mats2![[c3-c4,0,0,c3+c4,0,0],
[0,c3-c4,0,0,2*c3+2*c4,0],
[0,0,c3-c4,0,0,c3+c4],
[c3+c4,0,0,c3-c4,0,0],
[0,2*c3+2*c4,0,0,c3-c4,0],
[0,0,c3+c4,0,0,c3-c4]];

scal3:=mats3!DiagonalMatrix([m1,m1,m1,m1,m1,m1,m4,m4,m4,m4,m4,m4]);
for i in [1..6] do scal3[i+6,i]:=m3; scal3[i,i+6]:=m2; end for;

scal:=mats!DiagonalJoin(<scal1,scal2,scal3>);

// Now the code to find the coefficients for the form
V:=GModule(L);
g:=mats!h1;

// These speed up calculation of the form, to only look at non-zero entries
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:=Matricise(u); v1:=Matricise(v); w1:=Matricise(w);
u2:=Matricise(u)*mats!gg; v2:=Matricise(v)*mats!gg; w2:=Matricise(w)*mats!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;

coeffs0:=[a,b,c1,c2,c3,c4,m1,m2,m3,m4];

// The next two functions allow us to replace one of the coefficients by an expression.

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;

// The fast method first determines all relations, reduces the system, and then repeatedly
// takes Groebner bases and simplification to end up with the answer. We also check
// containment in G2(3).

function QuickCheckDetermination()

rels:=[];
for i in seqs do Append(~rels,ProdRel(i)); end for;
rels2:=Reduce(rels);

// Since a,c1,c2,c3,c4 are all non-zero, we can remove these terms from the expressions.
B2:=GroebnerBasis(rels2); B2:=SimplifyBasis(B2,[a,c1,c2,c3,c4]);
m3^2-m4^2 in B2; // Thus neither m3 nor m4 is zero, since m1m4 - m2m3=0.
m1*m4 + m2*m3 in B2; // Thus neither m1 nor m2 is zero, for otherwise both are.
B2:=SimplifyBasis(B2,[a,c1,c2,c3,c4,m1,m2,m3,m4]);
B2:=SimplifyBasis(GroebnerBasis(B2),[a,c1,c2,c3,c4,m1,m2,m3,m4]);
B2:=SimplifyBasis(GroebnerBasis(B2),[a,c1,c2,c3,c4,m1,m2,m3,m4]);
B2 eq [
    b + w^5*m2,
    m1 + 2*m2,
    c4 + 2*m4,
    c3 + w^2*m4,
    (m2+m4)*(m2+w^2*m4),
    c1 + w^2*m2 + w*m4,
    a + w*m2 + m4,
    m3 + m4,
    c2 + m2 + w^7*m4
];
coeffs1:=ChangeLinearCoefficients(coeffs0,B2);
coeffs2a:=ChangeCoefficient(coeffs1,m2,-m4);
coeffs2b:=ChangeCoefficient(coeffs1,m2,-w^2*m4);

// These are the solutions. We now scale by an element of the normalizer of H and obtain two matrices.
coeffs2a:=ChangeCoefficient(coeffs2a,m4,1);
coeffs2b:=ChangeCoefficient(coeffs2b,m4,1);

coeffs2a eq [w^7,w^5,1,w^6,w^6,1,2,2,2,1];
coeffs2b eq [w^5,w^7,w^6,2,w^6,1,w^6,w^6,2,1];

scala:=GL(27,9)!Evaluate(scal,coeffs2a);
scalb:=GL(27,9)!Evaluate(scal,coeffs2b);

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

// We know that h1 lies in K, so we want to check that t11 and t22 lie in E6.
{CheckRel(i,t11):i in seqs} eq {0};
{CheckRel(i,t22):i in seqs} eq {0};
return "All coefficients are zero, so matrices lie in E6";
end function;

// The next method is faster, but much more difficult as it involves working with the individual equations.
// We just check that we obtain the same answers.

function CheckDetermination()
// We embark on a proof of the result. We check our statements one by one:

Factorization(ProdRel(seqs[2847])*w^3-ProdRel(seqs[2885])) eq
[<c3 + w^2*c4, 1>,<c1*m3 + w^2*c2*m4, 1>];

// Thus either c3 + w^2*c4=0 or c1*m3 + w^2*c2*m4=0. Start with the latter, which will collapse to the former (called (*) ).
// Thus assume that this equation holds.

aa1:=c1*m3+w^2*c2*m4; // So can add and remove copies of aa1 without affecting zero-ness.

Factorization(c1*(c1*(ProdRel(seqs[2885])-c3*aa1)-2*c2*m3*aa1)-w^2*c2^2*m4*aa1) eq [<c1 + w*m4, 1>,<c1 + w^5*m4, 1>,<c1 + w^6*c2, 3>];

// Thus we have one of c1=-w*m4, c1=-w^5*m4, c1=-w^6*c2. Gives us three options.
badcoeffs1a:=ChangeCoefficient(coeffs0,c1,-w*m4);
badcoeffs1b:=ChangeCoefficient(coeffs0,c1,-w^5*m4);
badcoeffs1c:=ChangeCoefficient(coeffs0,c1,-w^6*c2);

Factorization(Evaluate(aa1,badcoeffs1a)) eq [<m4, 1>,<c2 + w^3*m3, 1>];
Factorization(Evaluate(aa1,badcoeffs1b)) eq [<m4, 1>,<c2 + w^7*m3, 1>];
Factorization(Evaluate(aa1,badcoeffs1c)) eq [<m3 + m4, 1>,<c2, 1>];
// In neither cases a nor b can m4 be zero, (as it is up to scalar c1), so second part holds in both cases. In branch c, c2 ne 0, so m3=-m4
badcoeffs2a:=ChangeCoefficient(badcoeffs1a,c2,-w^3*m3);
badcoeffs2b:=ChangeCoefficient(badcoeffs1b,c2,-w^7*m3);
badcoeffs2c:=ChangeCoefficient(badcoeffs1c,m4,-m3);

// Close off branch a first, as this is the easiest.
Factorization(Evaluate(ProdRel(seqs[2842]),badcoeffs2a)) eq [<m4, 1>,<m3, 1>,<c3 + w^2*c4, 1>];
// Since neither m3 nor m4 is zero, we have c3=-w^2*c4. But this is the original equation (*), so we have collapsed to the former.

// Now analyse branch b
Factorization(Evaluate(ProdRel(seqs[652]),badcoeffs2b)) eq [<c4^2 + m3*m4, 1>,<a + 2*b + 2*c4, 1>];
// If the second factor is zero, obtain b=a-c4. If the first factor we obtain:
Factorization(Evaluate(ProdRel(seqs[2842]),badcoeffs2b)+c4*(c4^2 + m3*m4)) eq [<c4 + 2*m3 + m4, 3>];
// So c4=m3-m4. So we obtain two branches:
badcoeffs3bcase1:=ChangeCoefficient(badcoeffs2b,b,a-c4);
badcoeffs3bcase2:=ChangeCoefficient(badcoeffs2b,c4,m3-m4);

Factorization(Evaluate(ProdRel(seqs[3335]),badcoeffs3bcase1)/w^6+Evaluate(ProdRel(seqs[2842]),badcoeffs3bcase1)) eq [<m4, 1>,<m3, 1>,<a + c4, 1>];
// So a=-c4. But at the same time
Factorization(Evaluate(ProdRel(seqs[399]),badcoeffs3bcase1)) eq [<m4, 1>,<m3, 1>,<a + w^6*c3 + 2*c4, 1>];
//Thus w^6*c3=a-c4=c4, so c3=w^6*c4. This is the original equation (*), so this branch collapses into the former as well.

// Now tackle case 2 of branch b
Factorization(Evaluate(ProdRel(seqs[2901]),badcoeffs3bcase2)) eq [<m3 + m4, 2>,<a, 1>];
// Thus m3=-m4
badcoeffs4bcase2:=ChangeCoefficient(badcoeffs3bcase2,m3,-m4);
// We cannot have m4=0, so a=-m4
badcoeffs5bcase2:=ChangeCoefficient(badcoeffs4bcase2,a,-m4);
Evaluate(ProdRel(seqs[2897]),badcoeffs5bcase2)+Evaluate(ProdRel(seqs[2939]),badcoeffs5bcase2) eq m1*m2*m4;
// Since m4 ne 0, either m1 or m2 must be 0.
Evaluate(ProdRel(seqs[656]),badcoeffs5bcase2) eq 2*b*m1*m2 + w^6*m4^3;
// Since m1*m2 =0, we obtain m4=0, a contradiction. This closes off branch b.

// End with branch c.
Factorization(Evaluate(ProdRel(seqs[383]),badcoeffs2c)/w^7-Evaluate(ProdRel(seqs[385]),badcoeffs2c)/w^5) eq [<m1 + 2*m2, 1>,<c3, 1>,<c2, 1>];
// Since c2,c3 ne 0, must have m1=m2.

badcoeffs3c:=ChangeCoefficient(badcoeffs2c,m2,m1);
Factorization(Evaluate(ProdRel(seqs[3324]),badcoeffs3c)-Evaluate(ProdRel(seqs[3303]),badcoeffs3c)+Evaluate(ProdRel(seqs[3307]),badcoeffs3c)/w^6) eq 
[<m3, 2>,<c4 + m3, 1>];
// Thus m3=-c4
badcoeffs4c:=ChangeCoefficient(badcoeffs3c,m3,-c4);
Factorization(Evaluate(ProdRel(seqs[3297]),badcoeffs4c)) eq [<c2 + w^3*c4, 1>,<c2 + w^7*c4, 1>,<a + w^5*c2, 1>];
// Thus c2=-w^3*c4, c2=-w^7*c4, or a=-w^5*c2. Deal with each in turn.
badcoeffs5ccase1:=ChangeCoefficient(badcoeffs4c,c2,-w^3*c4);
badcoeffs5ccase2:=ChangeCoefficient(badcoeffs4c,c2,-w^7*c4);
badcoeffs5ccase3:=ChangeCoefficient(badcoeffs4c,a,-w^5*c2);
// Now kill off each in turn:

Evaluate(ProdRel(seqs[2897]),badcoeffs5ccase1)+Evaluate(ProdRel(seqs[2939]),badcoeffs5ccase1) eq c4*m1^2;//neither of these can be 0
Evaluate(ProdRel(seqs[3301]),badcoeffs5ccase2) eq w^2*c4^3; // c4 cannot be 0
Factorization(Evaluate(ProdRel(seqs[2535]),badcoeffs5ccase3)) eq [<m1, 1>,<c3 + w^2*c4, 1>,<c2, 1>];
// Since m1,c2 cannot be 0, must have c3=-w^2*c4, which is equation (*) above.

// Thus in all cases for this branch we have either ended in contradiction or resulted in (*) holding. Thus we assume (*) holds.

coeffs1:=ChangeCoefficient(coeffs0,c3,-w^2*c4);
// Note that a and the ci cannot be 0.
Factorization(Evaluate(ProdRel(seqs[1676]),coeffs1)) eq [<m2, 1>,<m1, 1>,<a + b + c4, 1>];

// We either have m1=0, m2=0, or b=-a-c4. However,
Evaluate(ProdRel(seqs[383]),coeffs1)/w^7-Evaluate(ProdRel(seqs[385]),coeffs1)/w^5 eq c1*c4*m1 + w^6*c2*c4*m2;
// This means that if one of m1 and m2 is 0, so is the other one. So must have b=-a-c4. Also, we have
Factorization(Evaluate(ProdRel(seqs[246]),coeffs1)) eq [<a + w^7*c1 + w*c2 + c4 + 2*m3 + m4, 3>];
// So we have m4=-(a+w^7*c1+w*c2+c4+2*m3)
coeffs2:=ChangeCoefficients(coeffs1,[b,m4],[-a-c4,-(a+w^7*c1+w*c2+c4+2*m3)]);

// This took a long time to find!
Factorization(w^7*(a+c4)*Evaluate(ProdRel(seqs[1674]),coeffs2)+w^5*(a+c4)*Evaluate(ProdRel(seqs[1633]),coeffs2)
+w^2*(a+w^4*c4)*Evaluate(ProdRel(seqs[487]),coeffs2)) eq [<c4, 2>,<a + w*c4, 1>,<a + w^3*c4, 1>];
// Thus a=-w*c4 or a=-w^3*c4

coeffs3c1:=ChangeCoefficient(coeffs2,a,-w*c4);
coeffs3c2:=ChangeCoefficient(coeffs2,a,-w^3*c4);

Factorization(Evaluate(ProdRel(seqs[1641]),coeffs3c1)/w^7-Evaluate(ProdRel(seqs[1633]),coeffs3c1)/w) eq [<c4, 1>,<c2 + w^2*m1, 1>,<c1, 1>];
Factorization(Evaluate(ProdRel(seqs[1641]),coeffs3c1)-Evaluate(ProdRel(seqs[1674]),coeffs3c1)) eq [<c4, 1>,<c2, 1>,<c1 + 2*m2, 1>];
// Thus m2=c1 and m1=w^2*c2

Factorization(Evaluate(ProdRel(seqs[1641]),coeffs3c2)/w^7-Evaluate(ProdRel(seqs[1633]),coeffs3c2)/w) eq [<c4, 1>,<c2 + w^6*m1, 1>,<c1, 1>];
Factorization(Evaluate(ProdRel(seqs[1641]),coeffs3c2)-Evaluate(ProdRel(seqs[1674]),coeffs3c2)) eq [<c4, 1>,<c2, 1>,<c1 + m2, 1>];
// Thus m2=-c1 and m1=w^6*c2

coeffs4c1:=ChangeCoefficients(coeffs3c1,[m2,m1],[c1,w^2*c2]);
coeffs4c2:=ChangeCoefficients(coeffs3c2,[m2,m1],[-c1,w^6*c2]);

Factorization(Evaluate(ProdRel(seqs[32]),coeffs4c1)+w*Evaluate(ProdRel(seqs[1678]),coeffs4c1)+w*Evaluate(ProdRel(seqs[126]),coeffs4c1)) eq
[<c4, 1>,<c2 + 2*m3, 1>,<c1, 1>];
// Thus m3=c2
Factorization(Evaluate(ProdRel(seqs[32]),coeffs4c2)-w*Evaluate(ProdRel(seqs[1678]),coeffs4c2)+w*Evaluate(ProdRel(seqs[126]),coeffs4c2)) eq
[<c4, 1>,<c2 + w^6*m3, 1>,<c1, 1>];
// Thus m3=w^6*c2
coeffs5c1:=ChangeCoefficient(coeffs4c1,m3,c2);
coeffs5c2:=ChangeCoefficient(coeffs4c2,m3,w^6*c2);

Factorization(Evaluate(ProdRel(seqs[1932]),coeffs5c1)) eq [<c1 + w^2*c4, 3>];
Factorization(Evaluate(ProdRel(seqs[1932]),coeffs5c2)) eq [<c1+2*c4, 3>];
// Thus c1=-w^2*c4 in the first case and c1=c4 in the second case

coeffs6c1:=ChangeCoefficient(coeffs5c1,c1,-w^2*c4);
coeffs6c2:=ChangeCoefficient(coeffs5c2,c1,c4);
Factorization(Evaluate(ProdRel(seqs[126]),coeffs6c1)) eq [<c4, 2>,<c2 + c4, 1>];
Factorization(Evaluate(ProdRel(seqs[126]),coeffs6c2)) eq [<c4, 2>,<c2 + w^2*c4, 1>];
// Finally, we have c2=-c4 or c2=-w^2*c4
coeffs7c1:=ChangeCoefficient(coeffs6c1,c2,-c4);
coeffs7c2:=ChangeCoefficient(coeffs6c2,c2,-w^2*c4);

coeffs7c1 eq [w^5*c4,w^7*c4,w^6*c4,2*c4,w^6*c4,c4,w^6*c4,w^6*c4,2*c4,c4];
coeffs7c2 eq [w^7*c4,w^5*c4,c4,w^6*c4,w^6*c4,c4,2*c4,2*c4,2*c4,c4];
return "checks complete";
end function;


// Code to check that the trilinear form is correct.
function CheckE6Form()

GG:=GroupOfLieType("E6",F);
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,12,16,18,23] do x1[i,i]:=2; end for;
x2:=PermutationMatrix(F,Sym(27)!(1,4)(2,22,26,11,19,21,6,13,7,12,20,23,24,5,10,18,27,3,25,15,9,14,8));

// Confirm that L is a subgroup of E6 via conjugation by (x1*x2)^-1:
x1*x2*l1*x2^-1*x1^-1 in NGT and x1*x2*l2*x2^-1*x1^-1 in NGT and x1*x2*l3*x2^-1*x1^-1 in NGT;

mat:=[];
for h in [l1,l2,l3] 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;

//To go from NGT conjugate as x2*x^-1*MATRIX*x*x2^-1;
for g in [NGT.1,NGT.5,NGT.6,NGT.7] 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..7] 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;
return &and[ftest[i] eq f27[i]:i in [1..#f27]];
end function;

function CheckClassesOfL()
// L is isomorphic to SmallGroup(96,64)
// We have proved that L0 is toral, so may assume that when constructing our classes of L.

GG:=GroupOfLieType("E6",F);
W:=VectorSpace(F,6);
rho:=StandardRepresentation(GG);
rho2:=AdjointRepresentation(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];
Mats2:=[rho2(i):i in [g1,g2,g3,g4,g5,g6]];
Mats2E:=[rho2(i):i in Refs];
NGT:=sub<Over|Mats cat MatsE>;
T:=sub<NGT|Mats>;
CCNGT:=ConjugacyClasses(NGT);
[i:i in [1..#CCNGT]|CCNGT[i,1] eq 3 and Dimension(Socle(GModule(sub<NGT|CCNGT[i,3]>))) eq 9] eq [12,13];
C1:=sub<NGT|T,CCNGT[12,3]>;
C2:=sub<NGT|T,CCNGT[13,3]>;
NC1:=Normalizer(NGT,C1);
NC2:=Normalizer(NGT,C2);
D1,phi1:=PCGroup(NC1);
D2,phi2:=PCGroup(NC2);
ZD1:=Subgroups(D1:OrderEqual:=96);
ZD2:=Subgroups(D2:OrderEqual:=96);
ZD1a:=[i`subgroup:i in ZD1|IsIsomorphic(i`subgroup,SmallGroup(96,64))];
ZD2a:=[i`subgroup:i in ZD2|IsIsomorphic(i`subgroup,SmallGroup(96,64))];
ZC1a:=[i@@phi1:i in ZD1a];
ZC2a:=[i@@phi2:i in ZD2a];

// Now remove those that have the wrong character values on M(E6).
H:=PSU(3,3);
L:=MaximalSubgroups(H)[2]`subgroup;
IrrH:=IrreducibleModules(H,GF(3));
{*i:i in [BrauerCharacter(Restriction(IrrH[5],L))[j]:j in [1..10]]*} eq {* -1^^2, 0, 1^^2, 3^^4, 27 *};
ZC1b:=[i:i in ZC1a|{*ii:ii in [BrauerCharacter(GModule(i))[j]:j in [1..10]]*} eq {*-1^^2,0,1^^2,3^^4,27*}];
ZC2b:=[i:i in ZC2a|{*ii:ii in [BrauerCharacter(GModule(i))[j]:j in [1..10]]*} eq {*-1^^2,0,1^^2,3^^4,27*}];

// All that are left are conjugate.
return &and[IsConjugate(NGT,ZC1b[1],i):i in ZC1b cat ZC2b];
end function;

"The function CheckClassesOfL() checks that L is unique up to conjugacy in G. (This takes about a minute.)";
printf "\n";
"The function CheckE6Form() checks that the trilinear form f27 from E6 is correct. (This takes a couple of minutes.)";
printf "\n";
"The function QuickCheckDetermination() checks that the determination of the possible conjugates into E6 is correct using Groebner bases. It also checks that G2(3) contains H in all cases. (This takes about ten minutes.)";
printf "\n";
"The function CheckDetermination() checks that the determination of the possible conjugates into E6 is correct using individual equations. (This takes a few seconds.)";