/* In this case we let H denote a copy of PSL(2,13), L denote the subgroup 13.6, and let G denote F4(7^12).
We give matrices defining L, an element of order 13 and an element of order 6. We give an element x that
conjugates our version of L into the normalizer of a torus of F4(7^12), as constructed using the GroupOfLieType
command. We give H as a pair of matrices, one of which lies in L and the other chosen at random. A posteriori,
one could produce H as a subgroup of F4 and show that no other overgroup of L preserves the F4 trilinear form,
but we prefer a less synthetic proof.
*/
F<w>:=GF(7^12);
z:=w^(Order(w) div 13);

l1:=DiagonalMatrix(F,[z^i:i in [5,7,10,6,9,9,11,12,1,8,10,11,0,0,3,2,12,5,2,1,4,4,7,3,6,8]]);

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



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

x:=IdentityMatrix(F,26); for i in [1,12,25,8,15,4,7] do x[i,i]:=6; end for;
x[13,13]:=3; x[13,14]:=6; x[14,13]:=5; x[14,14]:=6;

// We need the indices of the trilinear form

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 F4. This was found by choosing elements of F4, 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.
*/

ents:=[260,273,295,313,327,332,571,584,597,608,628,644,840,855,883,896,926,942,1081,1098,1158,1171,1181,1209,1295,
1314,1410,1423,1433,1452,1581,1613,1640,1653,1674,1773,1807,1836,1848,1861,2012,2039,2052,2072,2167,2196,2207,
2220,2361,2374,2384,2483,2494,2507,2615,2628,2717,2822];

vals:=[1,1,6,6,1,6,6,2,4,6,6,1,6,6,4,2,1,6,1,6,4,2,6,6,6,1,4,2,6,6,1,6,2,4,6,6,6,6,1,1,6,2,4,6,6,6,1,1,2,4,1,1,1,
1,4,2,1,1];

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



// Now set up the centralizer of L in GL(26,F), which is given by all possible matrices scal below.

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

R<a,b,m1,m2,m3,m4,m5,m6,m7,m8>:=PolynomialRing(F,10);
mats:=MatrixRing(R,NumberOfRows(l1));
scal:=mats!1;
for i in [1..6] do
  scal[aa1[i],aa1[i]]:=m1; scal[aa1[i],aa2[i]]:=m2;
  scal[aa2[i],aa1[i]]:=m3; scal[aa2[i],aa2[i]]:=m4;
end for;
for i in [1..6] do
  scal[bb1[i],bb1[i]]:=m5; scal[bb1[i],bb2[i]]:=m6;
  scal[bb2[i],bb1[i]]:=m7; scal[bb2[i],bb2[i]]:=m8;
end for;
scal[13,13]:=a; scal[14,14]:=b;

// 12-dimensional for H:

ttlogs:=[8634333500,5075185700,7843725500,13382216900,10627795100,5188129700,5630964300,3523852800,2615359500,10097193600,7104883500,0,
5075185700,7843725500,13382216900,10627795100,5188129700,8634333500,3523852800,2615359500,10097193600,7104883500,0,5630964300,
7843725500,13382216900,10627795100,5188129700,8634333500,5075185700,2615359500,10097193600,7104883500,0,5630964300,3523852800,
13382216900,10627795100,5188129700,8634333500,5075185700,7843725500,10097193600,7104883500,0,5630964300,3523852800,2615359500,
10627795100,5188129700,8634333500,5075185700,7843725500,13382216900,7104883500,0,5630964300,3523852800,2615359500,10097193600,
5188129700,8634333500,5075185700,7843725500,13382216900,10627795100,0,5630964300,3523852800,2615359500,10097193600,7104883500,
7507952400,5400840900,4492347600,11974181700,8981871600,1876988100,5188129700,8634333500,5075185700,7843725500,13382216900,10627795100,
5400840900,4492347600,11974181700,8981871600,1876988100,7507952400,8634333500,5075185700,7843725500,13382216900,10627795100,5188129700,
4492347600,11974181700,8981871600,1876988100,7507952400,5400840900,5075185700,7843725500,13382216900,10627795100,5188129700,8634333500,
11974181700,8981871600,1876988100,7507952400,5400840900,4492347600,7843725500,13382216900,10627795100,5188129700,8634333500,5075185700,
8981871600,1876988100,7507952400,5400840900,4492347600,11974181700,13382216900,10627795100,5188129700,8634333500,5075185700,7843725500,
1876988100,7507952400,5400840900,4492347600,11974181700,8981871600,10627795100,5188129700,8634333500,5075185700,7843725500,13382216900];

tt:=GL(12,F)![w^i:i in ttlogs];

// 14-dimensional for H:

uulogs:=[0,13072326800,4613762400,9227524800,0,4613762400,9227524800,0,4613762400,9227524800,0,4613762400,9227524800,0,
7689604000,0,4613762400,0,9227524800,4613762400,0,9227524800,9227524800,4613762400,0,9227524800,4613762400,0,
2306881200,2306881200,5264955150,9172111650,8839632750,6512280450,4062101550,752136450,11182514850,9071167950,8133026850,1566039150,10962274050,7529482350,
11534406000,6920643600,9172111650,8839632750,6512280450,4062101550,752136450,5264955150,9071167950,8133026850,1566039150,10962274050,7529482350,11182514850,
6920643600,11534406000,8839632750,6512280450,4062101550,752136450,5264955150,9172111650,8133026850,1566039150,10962274050,7529482350,11182514850,9071167950,
2306881200,2306881200,6512280450,4062101550,752136450,5264955150,9172111650,8839632750,1566039150,10962274050,7529482350,11182514850,9071167950,8133026850,
11534406000,6920643600,4062101550,752136450,5264955150,9172111650,8839632750,6512280450,10962274050,7529482350,11182514850,9071167950,8133026850,1566039150,
6920643600,11534406000,752136450,5264955150,9172111650,8839632750,6512280450,4062101550,7529482350,11182514850,9071167950,8133026850,1566039150,10962274050,
2306881200,11534406000,6568752450,4457405550,3519264450,10793563950,6348511650,2915719950,752136450,5264955150,9172111650,8839632750,6512280450,4062101550,
11534406000,2306881200,4457405550,3519264450,10793563950,6348511650,2915719950,6568752450,5264955150,9172111650,8839632750,6512280450,4062101550,752136450,
6920643600,6920643600,3519264450,10793563950,6348511650,2915719950,6568752450,4457405550,9172111650,8839632750,6512280450,4062101550,752136450,5264955150,
2306881200,11534406000,10793563950,6348511650,2915719950,6568752450,4457405550,3519264450,8839632750,6512280450,4062101550,752136450,5264955150,9172111650,
11534406000,2306881200,6348511650,2915719950,6568752450,4457405550,3519264450,10793563950,6512280450,4062101550,752136450,5264955150,9172111650,8839632750,
6920643600,6920643600,2915719950,6568752450,4457405550,3519264450,10793563950,6348511650,4062101550,752136450,5264955150,9172111650,8839632750,6512280450];
uu1:=[w^i:i in uulogs]; uu1[1]:=0; uu1[16]:=0;
uu:=GL(14,F)!uu1; 

h1:=Matrix(F,26,26,[0:i in [1..26^2]]);
ab1:=aa1 cat bb1;
ab2:=[13,14] cat aa2 cat bb2;
for i in [1..12] do for j in [1..12] do h1[ab1[i],ab1[j]]:=tt[i,j]; end for; end for;
for i in [1..14] do for j in [1..14] do h1[ab2[i],ab2[j]]:=uu[i,j]; end for; end for;
h1:=GL(26,F)!h1;
H:=sub<GL(26,F)|[l1,l2,h1]>;

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

// We now show that g, which is the second generator of H (the first lies in L, so is fine) satisfies the equations mentioned
V:=GModule(L);
g:=mats!h1;

intseqs1:=[seqs[i]:i in ents];
intseqs:=&join{{[i[1],i[2],i[3]],[i[1],i[3],i[2]],[i[2],i[1],i[3]],[i[2],i[3],i[1]],[i[3],i[1],i[2]],[i[3],i[2],i[1]]}:i in intseqs1};

function ProdRel(seq)

u:=V.seq[1]; v:=V.seq[2]; w:=V.seq[3];
u1:=Matricise(u)*scal; v1:=Matricise(v)*scal; w1:=Matricise(w)*scal;
u2:=Matricise(u)*g*scal; v2:=Matricise(v)*g*scal; w2:=Matricise(w)*g*scal;
e1:=&+[u1[1,i[1]]*v1[1,i[2]]*w1[1,i[3]]*f26[Position(seqs,Sort([i[1],i[2],i[3]]))]:i in intseqs];
e2:=&+[u2[1,i[1]]*v2[1,i[2]]*w2[1,i[3]]*f26[Position(seqs,Sort([i[1],i[2],i[3]]))]:i in intseqs];

return e1-e2;
end function;

function CheckRel(seq,gg)
u:=V.seq[1]; v:=V.seq[2]; w:=V.seq[3];
u2:=u^gg; v2:=v^gg; w2:=w^gg;
e1:=&+[u[i[1]]*v[i[2]]*w[i[3]]*f26[Position(seqs,Sort([i[1],i[2],i[3]]))]:i in intseqs];
e2:=&+[u2[i[1]]*v2[i[2]]*w2[i[3]]*f26[Position(seqs,Sort([i[1],i[2],i[3]]))]:i in intseqs];
return e1-e2;
end function;

coeffs0:=[a,b,m1,m2,m3,m4,m5,m6,m7,m8];

function ChangeCoefficient(coeffs,coeff,target)
coeffs2:=coeffs;
nn:=Position(coeffs0,coeff);
coeffs2[nn]:=target;
return [Evaluate(i,coeffs2):i in coeffs];
end function;

function ChangeCoefficients(coeffs,coeff,target)
coeffs2:=coeffs;
for i in [1..#coeff] do
  nn:=Position(coeffs0,coeff[i]);
  coeffs2[nn]:=target[i];
end for;
return [Evaluate(i,coeffs2):i in coeffs];
end function;

function CheckDetermination()

// Note that m1 cannot equal 0. If it did, then 
Evaluate(ProdRel([4,4,19]),[a,b,0,m2,m3,m4,m5,m6,m7,m8])*w^11025805050+Evaluate(ProdRel([4,7,10]),[a,b,0,m2,m3,m4,m5,m6,m7,m8])*w^1798280250 eq m2^3;
// Note that m2 cannot equal 0. If it did, then 
Evaluate(ProdRel([4,4,19]),[a,b,m1,0,m3,m4,m5,m6,m7,m8])*w^6412042650+Evaluate(ProdRel([4,7,10]),[a,b,m1,0,m3,m4,m5,m6,m7,m8])*w^1798280250 eq m1^3;

//Thus we can express m5 in terms of m6:
(ProdRel([5,8,20])*w^11534406000-ProdRel([4,7,20])*w^7350536700)/w^3675268350 eq m1*m2*(m5-4*m6);

// Thus m6=2*m5
coeffs1:=ChangeCoefficient(coeffs0,m6,2*m5);
// We see from this that m5=0 if and only if m6=0 (as m1,m2 ne 0).
(w^8718923850*ProdRel([4,4,19])+2*ProdRel([5,8,20]))/w^1368387150 eq m1*(m2*m5-m2*m6+w^5982149550*m5*m6);

Evaluate((w^8718923850*ProdRel([4,4,19])+2*ProdRel([5,8,20]))/w^1368387150,coeffs1) eq -m1*m5*(m2-w^10595911950*m5);
// Therefore m2=w^10595911950*m5
coeffs2:=ChangeCoefficient(coeffs1,m2,w^10595911950*m5);

// Now produce m1.

5*Evaluate(ProdRel([4,5,19]),coeffs2) eq m5^2*(m1-w^12902793150*m5);
coeffs3:=ChangeCoefficient(coeffs2,m1,w^12902793150*m5);

Evaluate(ProdRel([1,7,8]),coeffs3)*w^9397058450-Evaluate(ProdRel([1,5,18]),coeffs3)*w^4783296050 eq m5^2*(a-w^13072326800*m3);
coeffs4:=ChangeCoefficient(coeffs3,m3,w^768960400*a);
coeffs4 eq [a,b,w^12902793150*m5,w^10595911950*m5,w^768960400*a,m4,m5,2*m5,m7,m8];

Evaluate(ProdRel([4,5,14]),coeffs4)/w^6751109950 eq m5^2*(m8-w^7689604000*a);
Evaluate(ProdRel([1,5,18]),coeffs4)/w^2906307950 eq m5^2*(b-w^3844802000*a);
coeffs5:=ChangeCoefficients(coeffs4,[m8,b],[w^7689604000*a,w^3844802000*a]);

Evaluate(ProdRel([1,4,4]),coeffs5) eq m5^2*(m4-2*m7);
Evaluate(ProdRel([1,10,10]),coeffs5)/2 eq m5^2*(m7-w^5382722800*a);
// Thus m4=2*m7=2*w^5382722800*a=w^9996485200*a

coeffs6:=ChangeCoefficients(coeffs5,[m4,m7],[w^9996485200*a,w^5382722800*a]);

// This is a two-parameter family, and of course the centralizer of H is two-parameter, so H is unique covering L.

// Check that this is a solution. First, that all matrices with the specialization above yield the same conjugate of h1.

scal6:=Evaluate(scal,coeffs6);
sca:=GL(26,F)!Evaluate(scal6,[1:i in [1..10]]);

gg:=h1^sca;
mats!sca^-1*scal6*mats!gg eq mats!gg*mats!sca^-1*scal6;

// Then that gg does lie in F4.

{*CheckRel(i,gg):i in seqs*} eq {*0^^3276*};

return "";
end function;

// We now give code that can prove that the trilinear form we gave is the correct one. The normalizer of a torus
// has a small space of trilinear forms for F4. We use two random generators to quickly reduce to that space, and
// and then choose a random element of F4 to reduce the space of trilinear forms down to 1 dimension.

function CheckF4Form()

q:=7^12;
F<w>:=GF(q);
GG:=GroupOfLieType("F4",q);
z:=w^(Order(w) div 13);
W:=VectorSpace(GF(q),4);
rho:=StandardRepresentation(GG);
Over:=GL(26,q);
g1:=elt<GG|W![z,1,1,1]>;
g2:=elt<GG|W![1,z,1,1]>;
g3:=elt<GG|W![1,1,z,1]>;
g4:=elt<GG|W![1,1,1,z]>;
Refs:=Reflections(GG);
Mats:=[rho(i):i in [g1,g2,g3,g4]];
MatsE:=[rho(i):i in Refs];
NGT:=sub<Over|Mats cat MatsE>;

// Confirm that L is a subgroup of F4 via conjugation by x:

x^-1*l1*x in NGT and x^-1*l2*x in NGT;

// And that the centralizer of the 6 is trivial, as claimed in the paper.

T:=sub<NGT|Mats>;
l2a:=NGT!(x^-1*l2*x);
Order(Centralizer(T,l2a)) eq 1;

// Now build the matrix of linear relations that the trilinear form must satisfy.

mat:=[];
repeat y1:=Random(NGT); y2:=Random(NGT); until NGT eq sub<NGT|y1,y2>;
for h in [x*y1*x^-1,x*y2*x^-1] do
  for nn in [1..#seqs] do a:=seqs[nn,1]; b:=seqs[nn,2]; c:=seqs[nn,3];
    val:=[F!0:i in [1..#seqs]];
    for i in [1..NumberOfRows(h)] do if(h[a,i] ne 0) then for j in [1..NumberOfRows(h)] do if(h[b,j] ne 0) then for k in [1..NumberOfRows(h)] do
      val[Position(seqs,Sort([i,j,k]))]+:=h[a,i]*h[b,j]*h[c,k];
    end for; end if; end for; end if; end for;
    val[nn]-:=1; Append(~mat,val); delete val;
  end for;
  delete h;
end for;

repeat
  h:=x*rho(Random(GG))*x^-1;
  for iii in [1..5] do nn:=Random([1..#seqs]); a:=seqs[nn,1]; b:=seqs[nn,2]; c:=seqs[nn,3];
    val:=[F!0:i in [1..#seqs]];
    for i in [1..NumberOfRows(h)] do if(h[a,i] ne 0) then for j in [1..NumberOfRows(h)] do if(h[b,j] ne 0) then for k in [1..NumberOfRows(h)] do
      val[Position(seqs,Sort([i,j,k]))]+:=h[a,i]*h[b,j]*h[c,k];
    end for; end if; end for; end if; end for;
    val[nn]-:=1; Append(~mat,val); delete val;
  end for;
  delete h;
until Dimension(Nullspace(Transpose(Matrix(F,mat)))) eq 1;

ftest:=Nullspace(Transpose(Matrix(F,mat))).1;
"Is the trilinear form given in this document correct?";
return &and[ftest[i] eq f26[i]:i in [1..#f26]];

end function;

"The function CheckF4Form() checks that the trilinear form f26 from F4 is correct, and the statement about elements of orders 2 and 3 made in the paper";
"The function CheckDetermination() checks that the determination of the possible conjugates into F4 is correct.";
