#We find all primitive subgroups of $GL(16,3)$ with parameters $e=8$, $a=2$. In this case $b=2$ divides $a$, so we use item (i) of Lemma 5.5 to find primitive 
#solvable subgroups in $GL(16,3)$. We start with subgroups in $GL(8,3^b)=GL(8,9)$ and then obtain subgroups in $GL(16,3)$
#In the notations of Theorem 5.1, if $H$ is a primitive solvable subgroup, then it possesses a series

#$$1<U<F<A<H,$$

#where $U$ is cyclic of order $3^2-1=8$, $F=U\circ E$ and $E$ is extraspecial of order $2^{1+8}$, moreover
#in this case $|H:A|$ divides $2$. 

gap> e:=8;; p:=3;; d:=16;; a:=2;;

#Since $b=2$, we take generators of primitive subgroups in $GL(8,9)$ from the file generators-GL(8,9).txt. as a list gens for every matrix 
#g in gens we substitute each entry $g_{i,j}\in GF(3^2)$ by a 2\times2-matrix over $GF(3)$ 




gap> bas:= Basis( AsVectorSpace( GF(p),GF(p^a) ) );;
gap> genS:=[];;                               
gap> for i in [1..Size(gens)] do
> genS[i]:=[];;
> for j in gens[i] do             
> Append(genS[i],[BlownUpMat(bas,j)]);;           
> od;
> od;


#Now we produce Singer cycle in $GL(2,3)$ that agrees with generators obtained above and element from normalizer inducing the field automorphism of order $2$. 
#We find a central element in subgroup H with generators genS[63], its center is cyclic of order $8$.

gap> H:=Subgroup(GL(16,3),genS[63]);;
gap> Zz:=Center(H);;
gap> Order(Zz);
8
gap> gn:=GeneratorsOfGroup(Zz);
[ < immutable compressed matrix 16x16 over GF(3) >, < immutable compressed matrix 16x
    16 over GF(3) >, < immutable compressed matrix 16x16 over GF(3) > ]
gap> Order(gn[1]);
4
gap> Order(gn[2]);
4
gap> Order(gn[3]);
8
gap> Display(gn[3]);
 . 2 . . . . . . . . . . . . . .
 2 2 . . . . . . . . . . . . . .
 . . . 2 . . . . . . . . . . . .
 . . 2 2 . . . . . . . . . . . .
 . . . . . 2 . . . . . . . . . .
 . . . . 2 2 . . . . . . . . . .
 . . . . . . . 2 . . . . . . . .
 . . . . . . 2 2 . . . . . . . .
 . . . . . . . . . 2 . . . . . .
 . . . . . . . . 2 2 . . . . . .
 . . . . . . . . . . . 2 . . . .
 . . . . . . . . . . 2 2 . . . .
 . . . . . . . . . . . . . 2 . .
 . . . . . . . . . . . . 2 2 . .
 . . . . . . . . . . . . . . . 2
 . . . . . . . . . . . . . . 2 2


#Now we produce t so that $I_9\otimes t=x$ and find $s$ such that $t^s=t^5$

gap> z:=Z(3);;
gap> t:=[[0*z,z],[z,z]];;
gap> Display(t);                
. 2
2 2
gap> s:=RepresentativeAction(GL(a,p),t,t^p);;
gap> s16:=KroneckerProduct(IdentityMat(8,Z(3)),s);;
gap> Display(s16);   
 2 2 . . . . . . . . . . . . . .
 . 1 . . . . . . . . . . . . . .
 . . 2 2 . . . . . . . . . . . .
 . . . 1 . . . . . . . . . . . .
 . . . . 2 2 . . . . . . . . . .
 . . . . . 1 . . . . . . . . . .
 . . . . . . 2 2 . . . . . . . .
 . . . . . . . 1 . . . . . . . .
 . . . . . . . . 2 2 . . . . . .
 . . . . . . . . . 1 . . . . . .
 . . . . . . . . . . 2 2 . . . .
 . . . . . . . . . . . 1 . . . .
 . . . . . . . . . . . . 2 2 . .
 . . . . . . . . . . . . . 1 . .
 . . . . . . . . . . . . . . 2 2
 . . . . . . . . . . . . . . . 1

#We check, that s16 normalizes Fitting subgroup of H

gap> F:=FittingSubgroup(H);;
gap> Order(F);
1944
gap> F^s16=F;
true


#We checked that s16 normalizes all subgroups with generators from the list genS, so we need not to generate the full normalizer

gap> for j in genS do
> l:=Subgroup(GL(16,3),j);
> Display(l^s16=l);
> od;
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true
true



#We add s16 to all generatos

gap> for j in genS do
> Append(j,[s16]);
> od;

#Now we run the following loop. Each entry of the list result has three items: the first is the
#number of the group in the list genS, the second is the order of the group, and the third is the order of the stabilizer of v in H. 
#If the order of the stabilizer is 1, then |v^H=|H| and Corollary 2.5 can be applied.

gap> v:=[z,z^2,0*z,z^2,z,z^2,z,0*z,z,z,z^2,z^2,z,0*z,z,z^2];;
gap> result:=[];;
gap> for k in [1..Size(genS)] do
> H:=Subgroup(GL(16,3),genS[k]);;
> n:=Size(H);;
> m:=Size(v^H);
> result[k]:=[k,n,n/m];
> Display(result[k]);
> od;
[ 1, 9216, 1 ]
[ 2, 9216, 1 ]
[ 3, 14336, 1 ]
[ 4, 15360, 1 ]
[ 5, 18432, 1 ]
[ 6, 18432, 1 ]
[ 7, 18432, 1 ]
[ 8, 18432, 1 ]
[ 9, 27648, 1 ]
[ 10, 27648, 1 ]
[ 11, 27648, 1 ]
[ 12, 30720, 2 ]
[ 13, 30720, 2 ]
[ 14, 30720, 1 ]
[ 15, 43008, 1 ]
[ 16, 36864, 1 ]
[ 17, 55296, 2 ]
[ 18, 55296, 1 ]
[ 19, 55296, 1 ]
[ 20, 55296, 1 ]
[ 21, 55296, 1 ]
[ 22, 55296, 1 ]
[ 23, 55296, 1 ]
[ 24, 55296, 1 ]
[ 25, 61440, 4 ]
[ 26, 61440, 1 ]
[ 27, 61440, 4 ]
[ 28, 82944, 1 ]
[ 29, 110592, 2 ]
[ 30, 110592, 1 ]
[ 31, 110592, 1 ]
[ 32, 110592, 1 ]
[ 33, 110592, 1 ]
[ 34, 110592, 1 ]
[ 36, 110592, 1 ]
[ 37, 110592, 1 ]
[ 38, 110592, 2 ]
[ 39, 110592, 2 ]
[ 40, 122880, 8 ]
[ 41, 165888, 1 ]
[ 42, 165888, 1 ]
[ 43, 165888, 2 ]
[ 44, 221184, 1 ]
[ 45, 221184, 1 ]
[ 46, 221184, 1 ]
[ 47, 221184, 1 ]
[ 48, 221184, 2 ]
[ 49, 221184, 1 ]
[ 50, 221184, 2 ]
[ 51, 221184, 1 ]
[ 52, 221184, 1 ]
[ 53, 221184, 1 ]
[ 54, 331776, 2 ]
[ 55, 331776, 2 ]
[ 56, 442368, 2 ]
[ 57, 442368, 1 ]
[ 58, 663552, 1 ]
[ 59, 663552, 1 ]
[ 60, 663552, 1 ]
[ 61, 663552, 2 ]
[ 62, 1327104, 1 ]
[ 63, 1327104, 2 ]

#Now we choose groups such that |v^H|<>|H| in the list res1


gap> res1:=Filtered(result,x->x[3]<>1);;
gap> Size(res1);
17

#We choose another vector v and repeat the cycle for the remaining groups 

gap> v:=[z,z^2,0*z,z^2,z,z^2,z,0*z,z,z,z^2,0*z^2,z,0*z,0*z,0*z^2];;
gap> result1:=[];;        
gap> for j in res1 do           
> H:=Subgroup(GL(16,3),genS[j[1]]);;
> n:=j[2];;
> m:=Size(v^H);;
> Append(result1,[[j[1],j[2],n/m]]);
> Display([j[1],j[2],n/m]);
> od;
[ 12, 30720, 1 ]
[ 13, 30720, 2 ]
[ 17, 55296, 1 ]
[ 25, 61440, 2 ]
[ 27, 61440, 1 ]
[ 29, 110592, 1 ]
[ 38, 110592, 1 ]
[ 39, 110592, 1 ]
[ 40, 122880, 1 ]
[ 43, 165888, 1 ]
[ 48, 221184, 1 ]
[ 50, 221184, 1 ]
[ 54, 331776, 1 ]
[ 55, 331776, 1 ]
[ 56, 442368, 1 ]
[ 61, 663552, 2 ]
[ 63, 1327104, 1 ]


#Now we choose groups such that |v^H|<>|H| in the list res2

gap> res2:=Filtered(result1,x->x[3]<>1);;
gap> Size(res2);

#We choose another vector v and repeat the cycle for the remaining groups 

gap> v:=[z,0*z^2,0*z,0*z^2,z,z^2,z,0*z,z,0*z,z^2,0*z^2,z,0*z,0*z,0*z^2];;
gap> result2:=[];;        
gap> for j in res2 do           
> H:=Subgroup(GL(16,3),genS[j[1]]);;
> n:=j[2];;
> m:=Size(v^H);;
> Append(result2,[[j[1],j[2],n/m]]);
> Display([j[1],j[2],n/m]);
> od;
[ 13, 30720, 1 ]
[ 25, 61440, 1 ]
[ 61, 663552, 1 ]