#We find all primitive subgroups of $GL(4,31)$ with parameters $e=4$, $a=1$. 
#In the notations of Theorem 5.1, if $H\leq GL(4,31)$ is a primitive solvable subgroup, then it possesses a series

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

#where $U$ is cyclic of order $30$, $F=U\circ E$ and $E$ is extraspecial of order $2^{1+4}$, moreover
#in this case $|H:A|$ equals $1$. Since $b=2$ does not divide $a=1$, Lemma 5.3 implies that there are two nonconjugate subgroups $F_1$ and $F_2$, and
#$N_{GL(4,31)}(F_1)/F_1\simeq Q^+(4,2)$, while  $N_{GL(4,31)}(F_2)/F_2\simeq Q^-(4,2)$


#First we generate the list of all primitive subgroups in $GL(4,31)$, using standart function of the package IRREDSOL

gap> e:=4;; p:=31;; d:=4;; a:=1;;

#Now we generate the list of all primitive subgroups in $GL(4,31)$, using standart function of the package IRREDSOL

gap> l0:=AllIrreducibleSolvableMatrixGroups(Degree,d,Field,GF(p),IsPrimitiveMatrixGroup,true);;

#In this list we choose subgroups, whose order is divisible by the order of $F$ in the notations of Theorem 5.1

gap> l1:=Filtered(l0,x-> (Order(x) mod ((p^a-1)*e*e)) = 0);;

#In view of Lemma 5.3, we choose subgroups satisfying $\vert A/F\vert\in \{12,20,72\}$ and collect these subgroups in the list l2.

gap> l2:=Filtered(l1,x-> (Order(x)/((p^a-1)*e*e)) in  [12,20,72]);;

#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(31),0*Z(31),0*Z(31),Z(31)^6];;
gap> result:=[];;
gap> i:=0;;
gap> for j in l2 do
> n:=Size(j);;
> m:=Size(v^j);
> i:=i+1;
> Append(result,[[i,n,n/m]]);
> Display([i,n,n/m]);
> od;
[ 1, 5760, 1 ]
[ 2, 5760, 1 ]
[ 3, 5760, 1 ]
[ 4, 5760, 2 ]
[ 5, 5760, 2 ]
[ 6, 5760, 1 ]
[ 7, 5760, 1 ]
[ 8, 5760, 1 ]
[ 9, 34560, 2 ]
[ 10, 5760, 1 ]
[ 11, 5760, 2 ]
[ 12, 5760, 1 ]
[ 13, 9600, 1 ]
[ 14, 5760, 1 ]
[ 15, 5760, 1 ]
[ 16, 5760, 1 ]


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


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

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

gap> v:=[Z(31),0*Z(31),Z(31)+Z(31),Z(31)^6];;
gap> result1:=[];;
gap> for j in res1 do
> m:=Size(v^l2[j[1]]);
> Append(result,[[j[1],j[2],j[2]/m]]);
> Display([j[1],j[2],j[2]/m]);
> od;
[ 4, 5760, 1 ]
[ 5, 5760, 1 ]
[ 9, 34560, 1 ]
[ 11, 5760, 1 ]