diff --git a/gap/ClassicalMaximals.gi b/gap/ClassicalMaximals.gi
index 2c8cf4db..d4405488 100644
--- a/gap/ClassicalMaximals.gi
+++ b/gap/ClassicalMaximals.gi
@@ -198,13 +198,17 @@ C7SubgroupsSpecialLinearGroupGeneric := function(n, q)
     highestPowern := Gcd(factorisationOfnExponents);
     
     divisorsHighestPowern := DivisorsInt(highestPowern);
+
     for t in divisorsHighestPowern{[2..Length(divisorsHighestPowern)]} do
         m := RootInt(n, t);
+        if m < 3 then
+            continue;
+        fi;
         tensorInducedSubgroup := TensorInducedDecompositionStabilizerInSL(m, t, q);
         # Cf. Tables 3.5.A and 3.5.G in [3]
         numberOfConjugates := Gcd(q - 1, m ^ (t - 1));
         if m mod 4 = 2 and t = 2 and q mod 4 = 3 then
-            numberOfConjugates := QuoInt(numberOfConjugates, 2);
+            numberOfConjugates := Gcd(q - 1, m) / 2;
         fi;
         result := Concatenation(result, 
                                 ConjugatesInGeneralGroup(tensorInducedSubgroup,
@@ -455,7 +459,7 @@ C2SubgroupsSpecialUnitaryGroupGeneric := function(n, q)
     result := List(divisorListOfn, t -> SUNonDegenerateImprimitives(n, q, t));
     # type GL(n / 2, q ^ 2).2 subgroups
     if IsEvenInt(n) then
-        Add(result, SUIsotropicImprimitives);
+        Add(result, SUIsotropicImprimitives(n, q));
     fi;
 
     return result;
@@ -488,6 +492,74 @@ C4SubgroupsSpecialUnitaryGroupGeneric := function(n, q)
     return result;
 end;
 
+C3SubgroupsSpecialUnitaryGroupGeneric := function(n, q)
+    return List(Filtered(PrimeDivisors(n), IsOddInt), 
+                s -> GammaLMeetSU(n, q, s));
+end;
+
+C5SubgroupsSpecialUnitaryGroupGeneric := function(n, q)
+    local factorisation, p, e, generatorGUMinusSU, primeDivisorsOfe,
+    degreeOfExtension, f, subfieldGroup, numberOfConjugates, result, epsilon;
+    
+    factorisation := PrimePowersInt(q);
+    p := factorisation[1];
+    e := factorisation[2];
+    generatorGUMinusSU := GU(n, q).1;
+    primeDivisorsOfe := PrimeDivisors(e);
+
+    result := [];
+    # type GU subgroups
+    for degreeOfExtension in primeDivisorsOfe do
+        if IsEvenInt(degreeOfExtension) then
+            continue;
+        fi;
+        f := QuoInt(e, degreeOfExtension);
+        subfieldGroup := SubfieldSL(n, p, e, f);
+        # Cf. Tables 3.5.B and 3.5.G in [3]
+        numberOfConjugates := Gcd(n, QuoInt(q + 1, p ^ f + 1));
+        result := Concatenation(result,
+                                ConjugatesInGeneralGroup(subfieldGroup, 
+                                                         generatorGUMinusSU, 
+                                                         numberOfConjugates));
+    od;
+
+    # type GO subgroups
+    if IsOddInt(q) then
+        if IsOddInt(n) then 
+            subfieldGroup := OrthogonalSubfieldSU(0, n, q);
+            # Cf. Tables 3.5.B and 3.5.G in [3]
+            numberOfConjugates := Gcd(n, q + 1);
+            result := Concatenation(result, 
+                                    ConjugatesInGeneralGroup(subfieldGroup,
+                                                             generatorGUMinusSU,
+                                                             numberOfConjugates));
+        else 
+            for epsilon in [-1, 1] do
+                subfieldGroup := OrthogonalSubfieldSU(epsilon, n, q);
+                # Cf. Tables 3.5.B and 3.5.G in [3]
+                numberOfConjugates := QuoInt(Gcd(q + 1, n), 2);
+                result := Concatenation(result,
+                                        ConjugatesInGeneralGroup(subfieldGroup,
+                                                                 generatorGUMinusSU,
+                                                                 numberOfConjugates));
+            od;
+        fi;
+    fi;
+
+    # type Sp subgroups
+    if IsEvenInt(n) then
+        subfieldGroup := SymplecticSubfieldSU(n, q);
+        # Cf. Tables 3.5.B and 3.5.G in [3]
+        numberOfConjugates := Gcd(QuoInt(n, 2), q + 1);
+        result := Concatenation(result,
+                                ConjugatesInGeneralGroup(subfieldGroup,
+                                                         generatorGUMinusSU,
+                                                         numberOfConjugates));
+    fi;
+
+    return result;
+end;
+
 C6SubgroupsSpecialUnitaryGroupGeneric := function(n, q)
     local factorisationOfq, p, e, factorisationOfn, r, m, result,
     generatorGUMinusSU, numberOfConjugates, extraspecialNormalizerSubgroup;
@@ -521,7 +593,7 @@ C6SubgroupsSpecialUnitaryGroupGeneric := function(n, q)
         fi;
     elif m >= 2 then
         # n = 2 ^ m >= 4
-        if e = 1 and (q - 1) mod 4 <> 0 then
+        if e = 1 and 2 * e = OrderMod(p, 4) then
             extraspecialNormalizerSubgroup := ExtraspecialNormalizerInSU(2, m, q);
             # Cf. Tables 3.5.A and 3.5.G in [3]
             numberOfConjugates := Gcd(n, q + 1);
@@ -538,10 +610,46 @@ C6SubgroupsSpecialUnitaryGroupGeneric := function(n, q)
     return result;
 end;
 
+C7SubgroupsSpecialUnitaryGroupGeneric := function(n, q)
+    local m, t, factorisationOfn, factorisationOfnExponents, highestPowern,
+    result, divisorsHighestPowern, numberOfConjugates, tensorInducedSubgroup, 
+    generatorGUMinusSU;
+
+    result := [];
+    generatorGUMinusSU := GU(n, q).1;
+    factorisationOfn := PrimePowersInt(n);
+    # get all exponents of prime factorisation of n
+    factorisationOfnExponents := factorisationOfn{Filtered([1..Length(factorisationOfn)], 
+                                                  IsEvenInt)};
+    # n can be written as k ^ highestPowern with k an integer and highestPowern
+    # is maximal with this property
+    highestPowern := Gcd(factorisationOfnExponents);
+    
+    divisorsHighestPowern := DivisorsInt(highestPowern);
+    for t in divisorsHighestPowern{[2..Length(divisorsHighestPowern)]} do
+        m := RootInt(n, t);
+        if (m = 3 and q = 2) or m < 3 then
+            continue;
+        fi;
+        tensorInducedSubgroup := TensorInducedDecompositionStabilizerInSU(m, t, q);
+        # Cf. Tables 3.5.B and 3.5.G in [3]
+        numberOfConjugates := Gcd(q + 1, m ^ (t - 1));
+        if m mod 4 = 2 and t = 2 and q mod 4 = 1 then
+            numberOfConjugates := Gcd(q + 1, m) / 2;
+        fi;
+        result := Concatenation(result, 
+                                ConjugatesInGeneralGroup(tensorInducedSubgroup,
+                                                         generatorGUMinusSU, 
+                                                         numberOfConjugates));
+    od;
+
+    return result;
+end;
 
 InstallGlobalFunction(MaximalSubgroupClassRepsSpecialUnitaryGroup,
 function(n, q, classes...)
-    local maximalSubgroups;
+    local maximalSubgroups, subfieldGroup, numberOfConjugates,
+    generatorGUMinusSU;
 
     if Length(classes) = 0 then
         classes := [1..9];
@@ -562,6 +670,7 @@ function(n, q, classes...)
         Error("PSU(3, 2) is soluble");
     fi;
 
+    generatorGUMinusSU := GU(n, q).1;
 
     maximalSubgroups := [];
 
@@ -595,15 +704,29 @@ function(n, q, classes...)
                 # q = 2
                 Add(maximalSubgroups, SUNonDegenerateImprimitives(n, q, 4));
             fi;
-        else 
-            # n = 6 and q = 2
+        elif n = 6 and q = 2 then 
             # Cf. Theorem 6.3.10 in [1]
-            Add(maximalSubgroups, SUNonDegenerateImprimitives(n, q, 3));
             Add(maximalSubgroups, SUNonDegenerateImprimitives(n, q, 2));
             Add(maximalSubgroups, SUIsotropicImprimitives(n, q));
         fi;
     fi;
 
+    if 3 in classes then
+        # Class C3 subgroups ######################################################
+        # Cf. Propositions 3.2.3 (n = 3), 3.3.4 (n = 4), 3.4.3 (n = 5), 
+        #                  3.5.5 (n = 6), 3.6.3 (n = 7), 3.7.5 (n = 8), 
+        #                  3.8.3 (n = 9), 3.9.5 (n = 10), 3.10.3 (n = 11), 
+        #                  3.11.7 (n = 12) in [1]
+        if not (n = 6 and q = 2) and not (n = 3 and q = 5)
+                                 and not (n = 3 and q = 3)
+                                 and not (n = 5 and q = 2) then
+            maximalSubgroups := Concatenation(maximalSubgroups, 
+                                              C3SubgroupsSpecialUnitaryGroupGeneric(n, q));
+        fi;
+        # There are no maximal C3 subgroups in the cases excluded above, cf.
+        # Proposition 3.5.5 and Theorem 6.3.10 in [1]
+    fi;
+
     if 4 in classes then
         # Class C4 subgroups ######################################################
         # Cf. Propositions 3.5.6 (n = 6), 3.7.7 (n = 8), 3.9.6 (n = 10), 
@@ -612,6 +735,37 @@ function(n, q, classes...)
                                           C4SubgroupsSpecialUnitaryGroupGeneric(n, q));
     fi;
 
+    if 5 in classes then
+        # Class C5 subgroups ######################################################
+        # Cf. Propositions 3.2.4 (n = 3), 3.3.5 (n = 4), 3.4.3 (n = 5), 
+        #                  3.5.7 (n = 6), 3.6.3 (n = 7), 3.7.8 (n = 8),
+        #                  3.8.4 (n = 9), 3.9.7 (n = 10), 3.10.3 (n = 11),
+        #                  3.11.9 (n = 12) in [1]
+        if not (n = 3 and q = 3) and not (n = 3 and q = 5) and not (n = 4 and q = 3) then
+            maximalSubgroups := Concatenation(maximalSubgroups,
+                                              C5SubgroupsSpecialUnitaryGroupGeneric(n, q));
+        # There are no maximal C5 subgroups for n = 3 and q = 3 or n = 3 and q = 5, 
+        # cf. Proposition 3.2.4 and Theorem 6.3.10 in [1]
+        elif n = 4 and q = 3 then
+            # type Sp
+            subfieldGroup := SymplecticSubfieldSU(n, q);
+            # Cf. Tables 3.5.B and 3.5.G in [3]
+            numberOfConjugates := 2;
+            maximalSubgroups := Concatenation(maximalSubgroups, 
+                                              ConjugatesInGeneralGroup(subfieldGroup,
+                                                                       generatorGUMinusSU,
+                                                                       numberOfConjugates));
+            # type GO-
+            subfieldGroup := OrthogonalSubfieldSU(-1, n, q);
+            # Cf. Tables 3.5.B and 3.5.G in [3]
+            numberOfConjugates := 2;
+            maximalSubgroups := Concatenation(maximalSubgroups,
+                                              ConjugatesInGeneralGroup(subfieldGroup,
+                                                                       generatorGUMinusSU,
+                                                                       numberOfConjugates));
+        fi;
+    fi;
+
     if 6 in classes then
         # Class C6 subgroups ######################################################
         # Cf. Lemma 3.1.6 (n = 2) and Propositions 3.2.5 (n = 3), 3.3.6 (n = 4),
@@ -627,6 +781,13 @@ function(n, q, classes...)
         fi;
     fi;
 
+    if 7 in classes then
+        # Class C7 subgroups ######################################################
+        # Cf. Proposition 3.8.6 (n = 9) in [1]
+        # For all other n, class C7 is empty.
+        maximalSubgroups := Concatenation(maximalSubgroups,
+                                          C7SubgroupsSpecialUnitaryGroupGeneric(n, q));
+    fi;
 
     return maximalSubgroups;
 end);
diff --git a/gap/ExtraspecialNormalizerMatrixGroups.gi b/gap/ExtraspecialNormalizerMatrixGroups.gi
index cc95279a..ed53b5f5 100644
--- a/gap/ExtraspecialNormalizerMatrixGroups.gi
+++ b/gap/ExtraspecialNormalizerMatrixGroups.gi
@@ -217,11 +217,12 @@ end;
 # If this function is called in the course of the computation of subgroups of
 # SU(d, q) then the argument q of the function is actually q ^ 2.
 OddExtraspecialNormalizerInSL := function(r, m, q, type...)
-    local d, listOfUi, listOfVi, V, generatorsOfNormalizerInGL, scalarMultiplierUi, 
+    local F, d, listOfUi, listOfVi, V, generatorsOfNormalizerInGL, scalarMultiplierUi, 
     scalarMultiplierVi, generators, generatingScalar, result, zeta, rootOfq;
 
+    F := GF(q);
     d := r ^ m;
-    zeta := PrimitiveElement(GF(q));
+    zeta := PrimitiveElement(F);
 
     if Length(type) = 0 then
         type := "L";
@@ -238,14 +239,14 @@ OddExtraspecialNormalizerInSL := function(r, m, q, type...)
     # We always need a generating element of Z(SL(d, q))
     if type = "L" then
         generatingScalar := zeta ^ (QuoInt(q - 1, Gcd(q - 1, d))) 
-                                * IdentityMat(d, GF(q));
+                                * IdentityMat(d, F);
     elif type = "U" then
         # remember that the variable q is actually q ^ 2 in the unitary case so
         # this is actually just q
         rootOfq := RootInt(q);
         generatingScalar := (zeta ^ (rootOfq - 1)) ^ QuoInt(rootOfq + 1,
                                                             Gcd(rootOfq + 1, d)) 
-                                * IdentityMat(d, GF(q));
+                                * IdentityMat(d, F);
     fi;
 
     # Note that not only det(Xi) = det(Yi) = 1, but as d is odd we
@@ -280,7 +281,7 @@ OddExtraspecialNormalizerInSL := function(r, m, q, type...)
     # after rescaling Vi - which is exactly what we needed for Vi to preserve
     # the unitary form given by the identity matrix (see above).
 
-    scalarMultiplierVi := ScalarToNormalizeDeterminant(V, r, GF(q));
+    scalarMultiplierVi := ScalarToNormalizeDeterminant(V, r, F);
     listOfVi := List(listOfVi, Vi -> scalarMultiplierVi * Vi);
 
     if d = 3 then
@@ -305,7 +306,7 @@ OddExtraspecialNormalizerInSL := function(r, m, q, type...)
             # Ui * I_d * HermitianConjugate(Ui).
 
             scalarMultiplierUi := ScalarToNormalizeDeterminant(listOfUi[1],
-                                                               d, GF(q));
+                                                               d, F);
             listOfUi := List(listOfUi, Ui -> scalarMultiplierUi * Ui);
         else
             # Note that Length(listOfUi) = m = 1 here and use 
@@ -320,6 +321,7 @@ OddExtraspecialNormalizerInSL := function(r, m, q, type...)
                                 generatorsOfNormalizerInGL.listOfYi,
                                 listOfUi, listOfVi,
                                 generatorsOfNormalizerInGL.listOfWi);
+    generators := List(generators, M -> ImmutableMatrix(F, M));
     result := Group(generators);
     # Size according to Table 2.9 of [1]
     if d = 3 and ((q - 4) mod 9 = 0 or (q - 7) mod 9 = 0) then
@@ -338,7 +340,7 @@ end;
 # If this function is called in the course of the computation of subgroups of
 # SU(d, q) then the argument q of the function is actually q ^ 2.
 SymplecticTypeNormalizerInSL := function(m, q, type...)
-    local generatorsOfNormalizerInGL, d, listOfUi, listOfVi, listOfWi,
+    local F, generatorsOfNormalizerInGL, d, listOfUi, listOfVi, listOfWi,
     generatingScalar, scalarMultiplierVi, i, scalarMultiplierUiAndWi, p, e, 
     factorization, generators, result, zeta, U1InGL, rootOfq;
     
@@ -353,12 +355,13 @@ SymplecticTypeNormalizerInSL := function(m, q, type...)
         type := type[1];
     fi;
 
+    F := GF(q);
     d := 2 ^ m;
     # q = p ^ e with p prime
     factorization := PrimePowersInt(q);
     p := factorization[1];
     e := factorization[2];
-    zeta := PrimitiveElement(GF(q));
+    zeta := PrimitiveElement(F);
 
     generatorsOfNormalizerInGL := SymplecticTypeNormalizerInGL(m, q, type);
     listOfUi := generatorsOfNormalizerInGL.listOfUi;
@@ -368,14 +371,14 @@ SymplecticTypeNormalizerInSL := function(m, q, type...)
     # We always need a generating element of Z(SL(d, q))
     if type = "L" then
         generatingScalar := zeta ^ (QuoInt(q - 1, Gcd(q - 1, d))) 
-                                * IdentityMat(d, GF(q));
+                                * IdentityMat(d, F);
     elif type = "U" then
         # remember that the variable q is actually q ^ 2 in the unitary case so
         # this is actually just q
         rootOfq := RootInt(q);
         generatingScalar := (zeta ^ (rootOfq - 1)) ^ QuoInt(rootOfq + 1,
                                                             Gcd(rootOfq + 1, d)) 
-                                * IdentityMat(d, GF(q));
+                                * IdentityMat(d, F);
     fi;
 
     # Note that det(Xi) = det(Yi) = 1, so we do not need to rescale these to
@@ -389,7 +392,7 @@ SymplecticTypeNormalizerInSL := function(m, q, type...)
 
         if type = "L" then
             scalarMultiplierVi := ScalarToNormalizeDeterminant(listOfVi[1], 
-                                                               d, GF(q));
+                                                               d, F);
         elif type = "U" then
             # We have to be a bit more sophisticated in case U for the
             # generators to preserve the unitary form given by the identity
@@ -398,7 +401,7 @@ SymplecticTypeNormalizerInSL := function(m, q, type...)
                 # This will give the expression Vi * I_d * HermitianConjugate(Vi)
                 # and additional factor of 2 ^ (-1), which is what we need; see
                 # the elif-case below, which is analogous, for more details.
-                scalarMultiplierVi := RootFFE(GF(q), 4 * zeta ^ 0, 4) ^ (-1);
+                scalarMultiplierVi := RootFFE(F, 4 * zeta ^ 0, 4) ^ (-1);
             else
                 # Remember that the variable q is actually q ^ 2 in the unitary
                 # case so p ^ QuoInt(e, 2) is actually q. Since 4 | q + 1 in
@@ -411,8 +414,8 @@ SymplecticTypeNormalizerInSL := function(m, q, type...)
                 #                        = 2 * (2 / q) * (-1)
                 #                        = 2
                 # by quadratic reciprocity, which is what we want.
-                i := RootFFE(GF(q), -1 * zeta ^ 0, 2);
-                scalarMultiplierVi := RootFFE(GF(q), 2 * i, 2) ^ (-1);
+                i := RootFFE(F, -1 * zeta ^ 0, 2);
+                scalarMultiplierVi := RootFFE(F, 2 * i, 2) ^ (-1);
             fi;
         fi;
 
@@ -450,9 +453,9 @@ SymplecticTypeNormalizerInSL := function(m, q, type...)
         # out wrong if only 4 | q + 1 and not 8 | q + 1.
 
         scalarMultiplierUiAndWi := ScalarToNormalizeDeterminant(listOfUi[1],
-                                                                d, GF(q));
+                                                                d, F);
         scalarMultiplierVi := ScalarToNormalizeDeterminant(listOfVi[1],
-                                                           d, GF(q));
+                                                           d, F);
         listOfUi := List(listOfUi, Ui -> scalarMultiplierUiAndWi * Ui);
         listOfWi := List(listOfWi, Wi -> scalarMultiplierUiAndWi * Wi);
         listOfVi := List(listOfVi, Vi -> scalarMultiplierVi * Vi);
@@ -477,6 +480,7 @@ SymplecticTypeNormalizerInSL := function(m, q, type...)
                                 generatorsOfNormalizerInGL.listOfYi,
                                 listOfUi, listOfVi, listOfWi);
 
+    generators := List(generators, M -> ImmutableMatrix(F, M));
     result := Group(generators);
 
     # Size according to Table 2.9 of [1]
@@ -568,7 +572,7 @@ end);
 # Construction as in Proposition 9.5 of [2]
 BindGlobal("ExtraspecialNormalizerInSU",
 function(r, m, q)
-    local result, F;
+    local F, result;
     if not r ^ m > 2 then
         ErrorNoReturn("<r> ^ <m> must be at least 2 in the unitary case, but",
                       " <r> = ", r, " and <m> = ", m);
@@ -577,6 +581,7 @@ function(r, m, q)
                       "unitary case, but <q> =", q, " and <r> =", r);
     fi;
 
+    F := GF(q ^ 2);
     if IsOddInt(r) then 
         result := OddExtraspecialNormalizerInSL(r, m, q ^ 2, "U");
     else
@@ -587,11 +592,8 @@ function(r, m, q)
     # (see the constructor functions above for more info).
     # We conjugate the group so that it preserves the standard GAP form
     # Antidiag(1, ..., 1). 
-    F := GF(q ^ 2);
     SetInvariantSesquilinearForm(result, rec(matrix := IdentityMat(r ^ m, F)));
-    result := ChangeFixedSesquilinearForm(result,
-                                          "U",
-                                          AntidiagonalMat(r ^ m, F));
+    result := ConjugateToStandardForm(result, "U");
 
     return result;
 end);
diff --git a/gap/Forms.gd b/gap/Forms.gd
new file mode 100644
index 00000000..bac4138f
--- /dev/null
+++ b/gap/Forms.gd
@@ -0,0 +1,3 @@
+DeclareGlobalFunction("ChangeFixedSesquilinearForm");
+DeclareGlobalFunction("ConjugateToStandardForm");
+DeclareGlobalFunction("UnitaryForm");
diff --git a/gap/Forms.gi b/gap/Forms.gi
new file mode 100644
index 00000000..f6c32a0c
--- /dev/null
+++ b/gap/Forms.gi
@@ -0,0 +1,204 @@
+# If <type> = "S" or type = "O" then <group> must have the attribute
+# InvariantBilinearForm.
+# Also, one need to ensure that the attribute DefaultFieldOfMatrixGroup is set
+# correctly for <group>; this can be done, for example, by making the
+# generators used during construction of the group immutable matrices over the
+# appropriate field.
+InstallGlobalFunction("ChangeFixedSesquilinearForm",
+function(group, type, gramMatrix)
+    local gapForm, newForm, gapToCanonical, canonicalToNew, field, formMatrix;
+    if not type in ["S", "O", "U"] then
+        ErrorNoReturn("<type> must be one of 'S', 'U', 'O', but <type> = ",
+                      type);
+    fi;
+    field := DefaultFieldOfMatrixGroup(group);
+    if type = "S" or type = "O" then
+        gapForm := BilinearFormByMatrix(InvariantBilinearForm(group).matrix, 
+                                        field);
+        newForm := BilinearFormByMatrix(gramMatrix, field);
+    else
+        if HasInvariantSesquilinearForm(group) then
+            gapForm := HermitianFormByMatrix(InvariantSesquilinearForm(group).matrix,
+                                             field);
+        else    
+            formMatrix := UnitaryForm(group);
+            if formMatrix = fail then
+                Error("No preserved unitary form found for <group>");
+            fi;
+            gapForm := HermitianFormByMatrix(formMatrix, field);
+        fi;
+        newForm := HermitianFormByMatrix(gramMatrix, field);
+    fi;
+    # the following if condition can only ever be fulfilled if <group> is an
+    # orthogonal group; there the case of even dimension is problematic since,
+    # in that case, there are two similarity classes of bilinear forms
+    if not WittIndex(gapForm) = WittIndex(newForm) then
+       ErrorNoReturn("The form preserved by <group> must be similar to the form ", 
+                     "described by the Gram matrix <gramMatrix>.");
+    fi;
+    gapToCanonical := BaseChangeHomomorphism(BaseChangeToCanonical(gapForm), 
+                                             field);
+    canonicalToNew := BaseChangeHomomorphism(BaseChangeToCanonical(newForm) ^ (-1), 
+                                             field);
+    return Group(canonicalToNew(gapToCanonical(GeneratorsOfGroup(group))));
+end);
+
+# Can only deal with sesquilinear forms, not with quadratic forms as of yet.
+# If <type> is one of "S", "O+", "O-" or "O" then <group> must have the
+# attribute InvariantBilinearForm.
+# Also, one need to ensure that the attribute DefaultFieldOfMatrixGroup is set
+# correctly for <group>; this can be done, for example, by making the
+# generators used during construction of the group immutable matrices over the
+# appropriate field.
+InstallGlobalFunction("ConjugateToStandardForm",
+function(group, type)
+    local d, F, q, gapForm, broadType;
+
+    # determining d (dimension of matrix group), F (base field) and q (order of
+    # F) plus some sanity checks
+    if not type in ["S", "O+", "O-", "O", "U"] then
+        ErrorNoReturn("<type> must be one of 'S', 'U', 'O+', 'O-', 'O'",
+                      " but <type> =", type);
+    fi;
+    F := DefaultFieldOfMatrixGroup(group);
+    d := DimensionOfMatrixGroup(group);
+    if type = "O" and IsEvenInt(d) then
+        ErrorNoReturn("<type> cannot be 'O' if the dimension of <group> is even");
+    elif type in ["O+", "O-"] and IsOddInt(d) then
+        ErrorNoReturn("<type> cannot be 'O+' or 'O-' if the dimension of",
+                      " <group> is odd");
+    fi;
+    if type in ["S", "O", "O+", "O-"] then
+        q := Size(F);
+    else
+        if IsSquareInt(Size(F)) then
+            q := RootInt(Size(F));
+        else
+            ErrorNoReturn("<group> must have a base field of square order", 
+                          " when <type> = 'U'");
+        fi;
+    fi;
+    if type in ["O", "O+", "O-"] and IsEvenInt(q) then
+        ErrorNoReturn("ConjugateToGAPForm cannot deal with orthogonal groups in",
+                      " even characteristic yet");
+    fi;
+    
+    # get standard GAP form
+    if type = "S" then
+        gapForm := InvariantBilinearForm(Sp(d, q)).matrix;
+    elif type = "U" then
+        gapForm := InvariantSesquilinearForm(GU(d, q)).matrix;
+    elif type = "O" then
+        gapForm := InvariantBilinearForm(GO(d, q)).matrix;
+    elif type = "O+" then
+        gapForm := InvariantBilinearForm(GO(1, d, q)).matrix;
+    elif type = "O-" then
+        gapForm := InvariantBilinearForm(GO(-1, d, q)).matrix;
+    fi;
+
+    if type in ["O", "O+", "O-"] then
+        broadType := "O";
+    else
+        broadType := type;
+    fi;
+
+    return ChangeFixedSesquilinearForm(group, broadType, gapForm);
+end);
+
+
+ConjugateModule := function(M, q)
+  return GModuleByMats(List(MTX.Generators(M), A -> ApplyFunctionToEntries(A, x -> x ^ q)), 
+                       MTX.Field(M));
+end;
+
+# Assuming that the group G acts absolutely irreducibly, try to find a unitary
+# form which is G-invariant or prove that no such form exists.
+#
+# We use this function instead of PreservedSesquilinearForms from the Forms
+# package since PreservedSesquilinearForms seems to be buggy and unreliable. 
+# As an example, take the group generated by
+#     [   1    0    0  ]         [ z^39 z^9  z^24 ]
+#     [ z^33 z^14 z^26 ]   and   [ z^25 z^16 z^6  ]
+#     [ z^19 z^31 z^5  ]         [ z^7  z^32 z^28 ]
+# where z = Z(49), which does preserve a unitary form, but this is not
+# recognised by PreservedSesquilinearForms, even after some 1000 calls of the
+# function.
+#
+# In general, this function should only be used if one can be sure that <G>
+# preserves a unitary form (but one does not know which one). 
+InstallGlobalFunction("UnitaryForm",
+function(G)
+    local d, F, q, M, inverseHermitianConjugateM, formMatrix, row, col, x,
+    scalar, counter;
+
+    d := DimensionOfMatrixGroup(G);
+    F := DefaultFieldOfMatrixGroup(G);
+    if not IsFinite(F) then
+        ErrorNoReturn("The base field of <G> must be finite");
+    fi;
+    if not IsEvenInt(DegreeOverPrimeField(F)) then
+        return fail;
+    fi;
+    q := RootInt(Size(F));
+
+    M := GModuleByMats(GeneratorsOfGroup(G), F);
+    # An element A of G acts as A ^ (-T) in MTX.DualModule(M) and hence as 
+    # HermitianConjugate(A, q) ^ (-1) in inverseHermitianConjugateM
+    inverseHermitianConjugateM := ConjugateModule(MTX.DualModule(M), q);
+
+    counter := 0;
+    scalar := fail;
+    # As the MeatAxe is randomised, we might have to make some more trials to
+    # find a preserved unitary form if there is one; breaking after 1000 trials
+    # is just a "safety net" in case a group <G> that does not preserve a
+    # unitary form is input.
+    while scalar = fail and counter < 1000 do
+        counter := counter + 1;
+
+        # If f: M -> inverseHermitianConjugateM is an isomorphism, it must respect
+        # multiplication by group elements, i.e. for A in G
+        #       f(x * A) = f(x) * HermitianConjugate(A, q) ^ (-1).
+        # Let f be given by the matrix F, i.e. f: x -> x * F. Then we have
+        #       (x * A) * F = x * F * HermitianConjugate(A, q) ^ (-1).
+        # Putting these results together for all vectors x gives
+        #       A * F = F * HermitianConjugate(A, q) ^ (-1)
+        # <==>  A * F * HermitianConjugate(A, q) = F,
+        # which is what we need.
+        formMatrix := MTX.IsomorphismModules(M, inverseHermitianConjugateM);
+
+        # We now need to ensure that formMatrix is actually the matrix of a
+        # unitary form, which can be achieved by multiplying it by a scalar
+        if formMatrix <> fail then
+            if formMatrix <> HermitianConjugate(formMatrix, q) then
+                # find a non-zero entry of formMatrix
+                row := First([1..NrRows(formMatrix)], x -> not IsZero(formMatrix[x]));
+                col := First([1..NrCols(formMatrix)], x -> not IsZero(formMatrix[row][x]));
+                if not IsZero(formMatrix[col, row]) then
+                    # this must be 1 for formMatrix to be hermitian
+                    x := formMatrix[row, col] * formMatrix[col, row] ^ (-q);
+                    # multiplying formMatrix by scalar will ensure that x = 1, i.e. that
+                    # formMatrix is hermitian
+                    scalar := RootFFE(x, q - 1);
+                fi;
+
+                if IsZero(formMatrix[col, row]) or scalar = fail then
+                    if not MTX.IsAbsolutelyIrreducible(M) then
+                        Error("UnitaryForm failed - group is not absolutely irreducible");
+                    fi;
+                    continue;
+                fi;
+
+                # make formMatrix hermitian
+                formMatrix := scalar * formMatrix;
+            fi;
+
+            if formMatrix <> HermitianConjugate(formMatrix, q) and not MTX.IsAbsolutelyIrreducible(M) then
+                Error("UnitaryForm failed - group is not absolutely irreducible");
+            fi;
+
+            return formMatrix;
+        fi;
+    od;
+
+    return fail;
+end);
diff --git a/gap/ImprimitiveMatrixGroups.gi b/gap/ImprimitiveMatrixGroups.gi
index 3e46ff11..338134d4 100644
--- a/gap/ImprimitiveMatrixGroups.gi
+++ b/gap/ImprimitiveMatrixGroups.gi
@@ -104,12 +104,11 @@ function(d, q, t)
                                    List([m + 2..d], i -> zeta ^ 0)));
     Add(generators, E);
 
+    generators := List(generators, M -> ImmutableMatrix(F, M));
     result := Group(generators);
     # change back fixed form into standard GAP form Antidiag(1, ..., 1)
     SetInvariantSesquilinearForm(result, rec(matrix := IdentityMat(d, F)));
-    result := ChangeFixedSesquilinearForm(result,
-                                          "U",
-                                          AntidiagonalMat(d, F));
+    result := ConjugateToStandardForm(result, "U");
     # Size according to Table 2.5 of [1]
     SetSize(result, Size(SU(m, q)) ^ t * (q + 1) ^ (t - 1) * Factorial(t));
     
@@ -169,6 +168,7 @@ function(d, q)
                                    [zeta ^ (-q - 1)]));
     Add(generators, D);
 
+    generators := List(generators, M -> ImmutableMatrix(F, M));
     result := Group(generators);
     # Size according to Table 2.5 of [1]
     SetSize(result, Size(SL(d / 2, q ^ 2)) * (q - 1) * 2);
diff --git a/gap/ReducibleMatrixGroups.gi b/gap/ReducibleMatrixGroups.gi
index 32dc9ba5..6705613e 100644
--- a/gap/ReducibleMatrixGroups.gi
+++ b/gap/ReducibleMatrixGroups.gi
@@ -106,6 +106,7 @@ function(d, q, k)
         Add(generators, D);
     fi;
 
+    generators := List(generators, M -> ImmutableMatrix(F, M));
     result := Group(generators);
     # Size according to Table 2.3 of [1]
     if d - 2 * k > 0 then
@@ -355,6 +356,7 @@ function(d, q, k)
         Add(generators, determinantShiftMatrix);
     fi;
 
+    generators := List(generators, M -> ImmutableMatrix(F, M));
     result := Group(generators);
     # Size according to Table 2.3 of [1]
     SetSize(result, Size(SU(k, q)) * Size(SU(d - k, q)) * (q + 1));
diff --git a/gap/SemilinearMatrixGroups.gi b/gap/SemilinearMatrixGroups.gi
index 2c0968e9..c298df0f 100644
--- a/gap/SemilinearMatrixGroups.gi
+++ b/gap/SemilinearMatrixGroups.gi
@@ -1,9 +1,13 @@
-# Return a subgroup of <M>GL(s, q)</M> isomorphic to the group <M>\Gamma L(1,
-# q ^ s)</M> of semilinear transformations of the vector space <M>F'^1</M>
-# over the field <C>F' := GF(q ^ s)</C>. See <Ref Func="GammaLMeetSL"/> for
-# further details.
+# TODO
+# This does not seem to construct the full GammaL(1, q ^ s): to get the full
+# GammaL(1, q ^ s), we would need the Frobenius x -> x ^ p as a generator, but
+# we only take x -> x ^ q. However, the results of the other functions are as
+# expected -- so it seems that we do not construct the full GammaL here, but
+# that we don't need it anyway?
+# TODO
+#
 # Construction as in Lemma 6.1 of [2]
-BindGlobal("CLASSICALMAXIMALS_GammaLDimension1",
+BindGlobal("GammaLDimension1",
 function(s, q)
     local primitivePolynomial, A, x, xq, B, row, i;
     # Let w be a primitive element of GF(q ^ s) over GF(q).
@@ -43,13 +47,13 @@ end);
 # Construction as in Proposition 6.3 of [2]
 BindGlobal("GammaLMeetSL",
 function(n, q, s)
-    local As, rootAs, Bs, Cs, Fs, m, gammaL1, Y, A, B, C, D, DBlock, ZBlock, i,
+    local As, Bs, Cs, Fs, m, gammaL1, Y, A, B, C, D, DBlock, ZBlock, i,
     range, result;
     if n mod s <> 0 or not IsPrime(s) then
         ErrorNoReturn("<s> must be prime and a divisor of <n> but <s> = ", s,
                       " and <n> = ", n);
     fi;
-    gammaL1 := CLASSICALMAXIMALS_GammaLDimension1(s, q);
+    gammaL1 := GammaLDimension1(s, q);
     # Let w be a primitive element of GF(q ^ s) over GF(q). Since As is the
     # companion matrix of the minimal polynomial of w over GF(q), its
     # determinant is (-1) ^ s times the constant term of said minimal
@@ -78,7 +82,6 @@ function(n, q, s)
 
     A := IdentityMat(n, GF(q));
     A{[1..s]}{[1..s]} := As;
-    A{[1..s]}{[1..s]} := As;
     A{[s + 1..2 * s]}{[s + 1..2 * s]} := As ^ -1;
     Y := SL(m, q ^ s).2;
     B := KroneckerProduct(Y, IdentityMat(s, GF(q)));
@@ -102,3 +105,94 @@ function(n, q, s)
     SetSize(result, Size(SL(n / s, q ^ s)) * (q ^ s - 1) / (q - 1) * s);
     return result;
 end);
+
+# Return a block matrix where the block in place (i, j) is A ^ k if and only if
+# the entry M[i, j] is omega ^ k (if M[i, j] = 0 then the corresponding block
+# is zero as well).
+Theta := function(M, A, omega)
+    local result, i, j, exponent, dimensionOfA;
+
+    if not NumberRows(A) = NumberColumns(A) then
+        ErrorNoReturn("<A> must be a square matrix but <A> = ", A);
+    fi;
+
+    dimensionOfA := NumberRows(A);
+    result := NullMat(NumberRows(M) * dimensionOfA, 
+                      NumberColumns(M) * dimensionOfA, 
+                      DefaultFieldOfMatrix(A));
+
+    for i in [1..NumberRows(M)] do
+        for j in [1..NumberColumns(M)] do
+            if not IsZero(M[i, j]) then
+                exponent := LogFFE(M[i, j], omega);
+                result{[(i - 1) * dimensionOfA + 1..i * dimensionOfA]}
+                      {[(j - 1) * dimensionOfA + 1..j * dimensionOfA]} 
+                          := A ^ exponent;
+            fi;
+        od;
+    od;
+
+    return result;
+end;
+
+# Construction as in Proposition 6.6 of [2]
+BindGlobal("GammaLMeetSU",
+function(d, q, s)
+    local F, As, Bs, Cs, Fs, m, gammaL1, Y, A, B, C, D, i,
+    range, result, AsInGU, omega, generators;
+    if d mod s <> 0 or not IsPrime(s) or not IsOddInt(s) then
+        ErrorNoReturn("<s> must be an odd prime and a divisor of <d> but <s> = ", s,
+                      " and <n> = ", d);
+    fi;
+    F := GF(q ^ 2);
+    gammaL1 := GammaLDimension1(s, q ^ 2);
+    # Let w be a primitive element of GF(q ^ (2 * s)) over GF(q ^ 2). Since As is the
+    # companion matrix of the minimal polynomial of w over GF(q ^ 2), its
+    # determinant is (-1) ^ s times the constant term of said minimal
+    # polynomial. By Vieta, this constant term is (-1) ^ s * the product of
+    # all Galois conjugates of w. Hence, det(As) = w ^ ((q ^ (2 * s) - 1) / (q ^ 2 - 1)).
+    As := gammaL1.A;
+    # By Lemma 6.2 det(Bs) = (-1) ^ (s - 1). Bs has order s as it is the
+    # Frobenius of the extension GF(q ^ 2) <= GF(q ^ (2 * s)).
+    Bs := gammaL1.B;
+    # Raise As to the (q ^ s - 1)-th power so that the result is in GammaU(1, q ^ s).
+    # AsInGU now has order q ^ s + 1 (as As had order q ^ (2 * s) - 1).
+    # det(AsInGU) = w ^ ((q ^ (2 * s) - 1) / (q + 1) * (q ^ s - 1) / (q - 1))
+    AsInGU := As ^ (q ^ s - 1);
+    # s is odd so (q ^ s - 1) / (q - 1) = 1 + q + ... + q ^ (s - 1) and q + 1
+    # are coprime. Hence, we have to raise AsInGU to the (q + 1)-th power to
+    # make the determinant 1. Cs has order (q ^ s + 1) / (q + 1).
+    Cs := AsInGU ^ (q + 1);
+    m := QuoInt(d, s);
+    if m = 1 then
+        # note that we require s to be odd
+        generators := [Bs, Cs];
+        generators := List(generators, M -> ImmutableMatrix(F, M));
+        result := Group(generators);
+        # Size according to Table 2.6 of [1]
+        SetSize(result, Size(SU(d / s, q ^ s)) * (q ^ s + 1) / (q + 1) * s);
+        # conjugate the result so that it preserves the standard unitary form
+        return ConjugateToStandardForm(result, "U");
+    fi;
+
+    omega := PrimitiveElement(GF(q ^ (2 * s)));
+    # The following two matrices generate SU(m, q ^ s) as a subgroup of SU(d, q)
+    A := Theta(SU(m, q ^ s).1, As, omega);
+    B := Theta(SU(m, q ^ s).2, As, omega);
+    # Note that GU(m, q ^ s).1 ^ (q + 1) has determinant 1.
+    C := Theta(GU(m, q ^ s).1 ^ (q + 1), As, omega);
+    # det(D) = 1
+    D := IdentityMat(d, GF(q));
+    for i in [0..m - 1] do
+        range := [i * s + 1..(i + 1) * s];
+        D{range}{range} := Bs;
+    od;
+
+    generators := [A, B, C, D];
+    generators := List(generators, M -> ImmutableMatrix(F, M));
+    result := Group(generators);
+    # Size according to Table 2.6 of [1]
+    SetSize(result, Size(SU(d / s, q ^ s)) * (q ^ s + 1) / (q + 1) * s);
+    # conjugate the result so that it preserves the standard unitary form 
+    return ConjugateToStandardForm(result, "U");
+end);
diff --git a/gap/SubfieldMatrixGroups.gi b/gap/SubfieldMatrixGroups.gi
index 0d180f4a..6a8e963c 100644
--- a/gap/SubfieldMatrixGroups.gi
+++ b/gap/SubfieldMatrixGroups.gi
@@ -49,8 +49,8 @@ end);
 
 BindGlobal("UnitarySubfieldSU",
 function(d, p, e, f)
-    local A, B, C, D, c, k, q, matrixForCongruence, lambda, zeta, omega, z, X,
-        result;
+    local F, A, B, C, D, c, k, q, matrixForCongruence, lambda, zeta, omega, z, X,
+        result, generators;
 
     if e mod f <> 0 or not IsPrimeInt(QuoInt(e, f)) or not IsOddInt(QuoInt(e, f)) then
         ErrorNoReturn("<f> must be a divisor of <e> and their quotient must be",
@@ -58,16 +58,19 @@ function(d, p, e, f)
     fi;
 
     q := p ^ e;
+    F := GF(q ^ 2);
     A := SU(d, p ^ f).1;
     B := SU(d, p ^ f).2;
-    zeta := PrimitiveElement(GF(q ^ 2));
+    zeta := PrimitiveElement(F);
     k := Gcd(q + 1, d);
     c := QuoInt(k * Lcm(p ^ f + 1, QuoInt(q + 1, k)), q + 1);
     # generates the center of SU(d, q)
-    C := zeta ^ QuoInt(q ^ 2 - 1, k) * IdentityMat(d, GF(q ^ 2));
+    C := zeta ^ QuoInt(q ^ 2 - 1, k) * IdentityMat(d, F);
 
     if c = Gcd(p ^ f + 1, d) then
-        result := Group(A, B, C);
+        generators := [A, B, C];
+        # generators := List(generators, M -> ImmutableMatrix(F, M));
+        result := Group(generators);
         # Size according to Table 2.8 of [1]
         SetSize(result, Size(SU(d, p ^ f)) * Gcd(QuoInt(q + 1, p ^ f + 1), d));
         return result;
@@ -85,9 +88,11 @@ function(d, p, e, f)
     # by calculating (d / k) ^ (-1) (mod (q + 1) / k).
     lambda := - QuoInt(z, (q - 1) * k) * ((d / k) ^ (-1) mod ((q + 1) / k));
     # det(X) = 1 by construction of lambda
-    X := zeta ^ (lambda * (q - 1)) * IdentityMat(d, GF(q ^ 2));
+    X := zeta ^ (lambda * (q - 1)) * IdentityMat(d, F);
 
-    result := Group(A, B, C, X * D);
+    generators := [A, B, C, X * D];
+    generators := List(generators, M -> ImmutableMatrix(F, M));
+    result := Group(generators);
     # Size according to Table 2.8 of [1]
     SetSize(result, Size(SU(d, p ^ f)) * Gcd(QuoInt(q + 1, p ^ f + 1), d)); 
 
@@ -96,20 +101,21 @@ end);
 
 BindGlobal("SymplecticSubfieldSU",
 function(d, q)
-    local generators, zeta, k, C, c, result, D, form;
+    local F, generators, zeta, k, C, c, result, D, form;
 
     if IsOddInt(d) then
         ErrorNoReturn("<d> must be even but <d> = ", d);
     fi;
 
+    F := GF(q ^ 2);
     form := AntidiagonalMat(Concatenation(List([1..d / 2], i -> 1),
                                           List([1..d / 2], i -> -1)) * Z(q ^ 2)^0,
-                            GF(q ^ 2));
+                            F);
     generators := ShallowCopy(GeneratorsOfGroup(Sp(d, q)));
-    zeta := PrimitiveElement(GF(q ^ 2));
+    zeta := PrimitiveElement(F);
     k := Gcd(q + 1, d);
     # generates the center of SU(d, q)
-    C := zeta ^ QuoInt(q ^ 2 - 1, k) * IdentityMat(d, GF(q ^ 2));
+    C := zeta ^ QuoInt(q ^ 2 - 1, k) * IdentityMat(d, F);
     Add(generators, C);
     c := QuoInt(Gcd(2, q - 1) * Gcd(q + 1, d / 2), Gcd(q + 1, d));
 
@@ -122,7 +128,8 @@ function(d, q)
                                        List([1..d / 2], i -> - zeta ^ (- (q + 1) / 2))));
         Add(generators, D);
     fi;
-     
+    
+    generators := List(generators, M -> ImmutableMatrix(F, M));
     result := Group(generators);
     if IsOddInt(q) then
         # The result preserves the unitary form given by 
@@ -131,9 +138,7 @@ function(d, q)
         # Antidiag(1, ..., 1). (If q is even, this is not necessary.)
         SetInvariantSesquilinearForm(result, 
                                      rec(matrix := - zeta ^ QuoInt(q + 1, 2) * form));
-        result := ChangeFixedSesquilinearForm(result, 
-                                              "U",
-                                              AntidiagonalMat(d, GF(q ^ 2)));
+        result := ConjugateToStandardForm(result, "U");
     fi;
     # Size according to Table 2.8 of [1]
     SetSize(result, Size(Sp(d, q)) * Gcd(q + 1, d / 2));
@@ -143,7 +148,7 @@ end);
 
 BindGlobal("OrthogonalSubfieldSU",
 function(epsilon, d, q)
-    local zeta, k, C, generators, SOChangedForm, result,
+    local F, zeta, k, C, generators, SOChangedForm, result,
     generatorsOfOrthogonalGroup, D, E, i, W, n, form;
 
     if IsEvenInt(q) then
@@ -159,10 +164,11 @@ function(epsilon, d, q)
                       "<epsilon> = ", epsilon);
     fi;
 
-    zeta := PrimitiveElement(GF(q ^ 2));
+    F := GF(q ^ 2);
+    zeta := PrimitiveElement(F);
     k := Gcd(q + 1, d);
     # generates the center of SU(d, q)
-    C := zeta ^ QuoInt(q ^ 2 - 1, k) * IdentityMat(d, GF(q ^ 2));
+    C := zeta ^ QuoInt(q ^ 2 - 1, k) * IdentityMat(d, F);
     generators := [C];
 
     if IsOddInt(d) then
@@ -170,7 +176,7 @@ function(epsilon, d, q)
                                                      "O",
                                                      AntidiagonalMat(d, GF(q)));
         generators := Concatenation(generators, GeneratorsOfGroup(SOChangedForm));
-
+        generators := List(generators, M -> ImmutableMatrix(F, M));
         result := Group(generators);
     else
         generatorsOfOrthogonalGroup := GeneratorsOfOrthogonalGroup(epsilon, d, q);
@@ -231,6 +237,7 @@ function(epsilon, d, q)
             # the standard orthogonal form in this case is Antidiag(1, ..., 1),
             # i.e. has the same form matrix as the unitary form we want, so we
             # do not need to conjugate the result
+            generators := List(generators, M -> ImmutableMatrix(F, M));
             result := Group(generators);
 
         elif epsilon = -1 then
@@ -259,6 +266,7 @@ function(epsilon, d, q)
             fi;
             Add(generators, W);
 
+            generators := List(generators, M -> ImmutableMatrix(F, M));
             result := Group(generators);
 
             # We still have to change the preserved unitary form to the
@@ -269,19 +277,13 @@ function(epsilon, d, q)
                 form := DiagonalMat(Concatenation([zeta ^ (q + 1)],
                                                   List([2..d], i -> zeta ^ 0)));
                 SetInvariantSesquilinearForm(result, rec(matrix := form));
-                result := ChangeFixedSesquilinearForm(result,
-                                                      "U",
-                                                      AntidiagonalMat(d, 
-                                                                      GF(q ^ 2)));
+                result := ConjugateToStandardForm(result, "U");
             else
                 # The form preserved by the constructed subgroup is given by
                 # the identity matrix.
-                form := IdentityMat(d, GF(q ^ 2));
+                form := IdentityMat(d, F);
                 SetInvariantSesquilinearForm(result, rec(matrix := form));
-                result := ChangeFixedSesquilinearForm(result,
-                                                      "U",
-                                                      AntidiagonalMat(d, 
-                                                                      GF(q ^ 2)));
+                result := ConjugateToStandardForm(result, "U");
             fi;
         fi;
     fi;
diff --git a/gap/TensorInducedMatrixGroups.gi b/gap/TensorInducedMatrixGroups.gi
index 9e59213d..51597581 100644
--- a/gap/TensorInducedMatrixGroups.gi
+++ b/gap/TensorInducedMatrixGroups.gi
@@ -1,7 +1,8 @@
 ######################## WORKING OUT FOR ELSE BRANCH BELOW ####################
 #
 #     Let m be even and assume that a generator of the
-#     TensorWreathProduct(SL(m, q), Sym(t)) has determinant not equal to 1.
+#     TensorWreathProduct(SL(m, q), Sym(t)) (or the same with SU(m, q) instead
+#     of SL(m, q), respectively) has determinant not equal to 1.
 #
 #     We generate TensorWreathProduct(SL(m, q), Sym(t)) from 
 #       * (m - 1)-fold Kronecker products of the generators of SL(m, q) with 
@@ -49,6 +50,10 @@
 #     filtered out in an if-statement before and the case q even will always
 #     give determinant 1 as -1 = 1 in characteristic 2.
 #
+#     In the unitary case, the analysis above still holds true; however, now
+#     the case t = 2, m = 2 mod 4 and q = 1 mod 4 was filtered out before so 
+#     t = 2, m = 2 mod 4 and q = 3 mod 4 remains.
+#
 ###############################################################################     
 
 
@@ -109,10 +114,13 @@ function(m, t, q)
     for i in [3..t] do
         U := KroneckerProduct(U, IdentityMat(m, GF(q)));
     od;
+    # det(U) = 1
     E := D ^ QuoInt(Gcd(q - 1, d), Gcd(q - 1, m ^ (t - 1)));
     for i in [2..t] do
         E := KroneckerProduct(E, IdentityMat(m, GF(q)));
     od;
+    # det(E) = zeta ^ (Gcd(q - 1, d) / Gcd(q - 1, m ^ (t - 1)) * m ^ (t - 1))
+    #        = zeta ^ (Gcd(q - 1, d) / Gcd(q - 1, m ^ (t - 1)) * d / m)
 
     # Write mu = zeta ^ k for some k. We want 
     # zeta ^ Gcd(q - 1, d) = det(mu * I_d) = mu ^ d = zeta ^ (kd), thus 
@@ -122,6 +130,10 @@ function(m, t, q)
     # k = 1 / (d / Gcd(q - 1, d)) mod ((q - 1) / Gcd(q - 1, d)).
     mu := zeta ^ (1 / (d / Gcd(q - 1, d)) mod ((q - 1) / Gcd(q - 1, d)));
     S := mu ^ (- d / (Gcd(q - 1, d / m) * m)) * IdentityMat(d, GF(q));
+    # det(S) = det(mu * I_d) ^ (- d / (Gcd(q - 1, d / m) * m))
+    #        = (zeta ^ Gcd(q - 1, d)) ^ (- d / (Gcd(q - 1, d / m) * m))
+    #        = zeta ^ (- Gcd(q - 1, d) / Gcd(q - 1, m ^ (t - 1)) * d / m)
+    #        = det(E) ^ (-1)
 
     gens := Concatenation(generatorsOfHInSL, [C, U, S * E]);
     result := Group(gens);
@@ -136,3 +148,102 @@ function(m, t, q)
     return result;
 end);
 
+# Construction as in Proposition 10.4 of [2]
+# Note, though, that the structure of G / Z(G) given there is incorrect and
+# that one should rather consult Table 2.10 of [1] on that (which, however, 
+# gives the structure of G, not G / Z(G)!).
+BindGlobal("TensorInducedDecompositionStabilizerInSU",
+function(m, t, q)
+    local F, gensOfSUm, I, D, C, generatorsOfHInSU, gens, i, H, E, U, S, zeta, mu,
+    result, scalingMatrix, d, generator, k;
+    if not t > 1 or not m > 1 then
+        ErrorNoReturn("<t> must be greater than 1 and <m> must be greater than 1 but <t> = ", 
+                      t, " and <m> = ", m);
+    fi;
+
+    F := GF(q ^ 2);
+    d := m ^ t;
+    zeta := PrimitiveElement(F);
+    D := DiagonalMat(Concatenation([zeta], 
+                                   List([1..m - 2], i -> zeta ^ 0),
+                                   [zeta ^ (- q)]));
+    # generates the center of SU(d, q)
+    C := zeta ^ ((q ^ 2 - 1) / Gcd(q + 1, d)) * IdentityMat(d, F);
+
+    if t = 2 and m mod 4 = 2 and q mod 4 = 1 then
+        gensOfSUm := GeneratorsOfGroup(SU(m, q));
+        I := IdentityMat(m, F);
+        # Let Z = Z(SU(d, q)). Then these generate the group 
+        # Z.(SU(m, q) o SU(m, q)) (to see this, realize the first factor of the
+        # central product as all Kronecker Products I * M with M in SU(m, q)
+        # and, similarly, the second factor as the Kronecker Products M * I).
+        generatorsOfHInSU := [KroneckerProduct(gensOfSUm[1], I),
+                              KroneckerProduct(gensOfSUm[2], I),
+                              KroneckerProduct(I, gensOfSUm[1]),
+                              KroneckerProduct(I, gensOfSUm[2])];
+    else
+        H := TensorWreathProduct(SU(m, q), SymmetricGroup(t));
+        generatorsOfHInSU := [];
+        for generator in GeneratorsOfGroup(H) do
+            if DeterminantMat(generator) = zeta ^ 0 then
+                Add(generatorsOfHInSU, generator);
+            else
+                # det = -1 for odd permutation
+                if IsOddInt(m) then
+                    Add(generatorsOfHInSU, -1 * generator);
+                else
+                    # In this case, we have t = 2, m = 2 mod 4 and q = 3 mod 4
+                    # (see working out above).
+
+                    # This has determinant ((det D) ^ ((q + 1) / 4)) ^ m 
+                    # = ((zeta ^ (1 - q)) ^ ((q + 1) / 4)) ^ m, which, using m even,
+                    # becomes (zeta ^ (- (q ^ 2 - 1) / 2)) ^ (m / 2) = (-1) ^ (m / 2)
+                    # and this is -1 due to m = 2 mod 4.
+                    scalingMatrix := KroneckerProduct(D ^ QuoInt(q + 1, 4), 
+                                                      IdentityMat(m, F));
+                    # det(generator * scalingMatrix) = -1 * (-1) = 1
+                    Add(generatorsOfHInSU,(generator * scalingMatrix));
+                fi;
+            fi;
+        od;
+    fi;
+
+    U := KroneckerProduct(D, D ^ (-1));
+    for i in [3..t] do
+        U := KroneckerProduct(U, IdentityMat(m, F));
+    od;
+    # det(U) = 1
+    E := D ^ QuoInt(Gcd(q + 1, d), Gcd(q + 1, m ^ (t - 1)));
+    for i in [2..t] do
+        E := KroneckerProduct(E, IdentityMat(m, F));
+    od;
+    # det(E) = zeta ^ ((1 - q) * Gcd(q + 1, d) / Gcd(q + 1, m ^ (t - 1)) * m ^ (t - 1))
+    #        = zeta ^ ((1 - q) * Gcd(q + 1, d) / Gcd(q + 1, d / m) * d / m))
+
+    # Write mu = zeta ^ ((q - 1) * k) for some k. We want 
+    # det(mu * I_d) = zeta ^ ((q - 1) * Gcd(q + 1, d)), hence
+    # (q - 1) * k * d = (q - 1) * Gcd(q + 1, d) mod (q ^ 2 - 1).
+    # Divide through by (q - 1) and Gcd(q + 1, d) to obtain
+    # k * d / Gcd(q + 1, d) = 1 mod ((q + 1) / Gcd(q + 1, d)). 
+    # Now d / Gcd(q + 1, d) is invertible and we can take 
+    # k = 1 / (d / Gcd(q + 1, d)) mod ((q + 1) / Gcd(q + 1, d)).
+    k := 1 / (d / Gcd(q + 1, d)) mod ((q + 1) / Gcd(q + 1, d));
+    mu := zeta ^ ((q - 1) * k);
+    S := mu ^ (d / (Gcd(q + 1, d / m) * m)) * IdentityMat(d, F);
+    # det(S) = det(mu * I_d) ^ (d / (Gcd(q + 1, d / m) * m))
+    #        = zeta ^ ((q - 1) * Gcd(q + 1, d) / Gcd(q + 1, d / m) * d / m)
+    #        = det(E) ^ (-1)
+
+    gens := Concatenation(generatorsOfHInSU, [C, U, S * E]);
+    gens := List(gens, M -> ImmutableMatrix(F, M));
+    result := Group(gens);
+    # Size according to Table 2.10 of [1]
+    if t = 2 and m mod 4 = 2 and q mod 4 = 1 then
+        SetSize(result, Gcd(q + 1, m) * Size(PSU(m, q)) ^ 2 * Gcd(q + 1, m) ^ 2);
+    else
+        SetSize(result, Gcd(q + 1, m) * Size(PSU(m, q)) ^ t 
+                                      * Gcd(q + 1, m ^ (t - 1)) * Gcd(q + 1, m) ^ (t - 1) 
+                                      * Factorial(t));
+    fi;
+    return result;
+end);
diff --git a/gap/TensorProductMatrixGroups.gi b/gap/TensorProductMatrixGroups.gi
index 93e77b30..4da93177 100644
--- a/gap/TensorProductMatrixGroups.gi
+++ b/gap/TensorProductMatrixGroups.gi
@@ -91,12 +91,11 @@ function(d1, d2, q)
         od;
     fi;
 
+    generators := List(generators, M -> ImmutableMatrix(F, M));
     result := Group(generators);
     # change back fixed form into standard GAP form Antidiag(1, ..., 1)
     SetInvariantSesquilinearForm(result, rec(matrix := IdentityMat(d, F)));
-    result := ChangeFixedSesquilinearForm(result,
-                                          "U",
-                                          AntidiagonalMat(d, F));
+    result := ConjugateToStandardForm(result, "U");
     # Size according to Table 2.7 in [1]
     SetSize(result, Size(SU(d1, q)) * Size(SU(d2, q)) * Gcd(q + 1, d1, d2));
     return result;
diff --git a/gap/Utils.gd b/gap/Utils.gd
index 7a1afa80..f3bd819a 100644
--- a/gap/Utils.gd
+++ b/gap/Utils.gd
@@ -12,7 +12,6 @@
 #! written into the matrix row by row.
 DeclareGlobalFunction("MatrixByEntries");
 
-DeclareGlobalFunction("ChangeFixedSesquilinearForm");
 DeclareGlobalFunction("AntidiagonalMat");
 DeclareGlobalFunction("SolveQuadraticCongruence");
 DeclareGlobalFunction("ApplyFunctionToEntries");
diff --git a/gap/Utils.gi b/gap/Utils.gi
index 74459c56..2f4a3e93 100644
--- a/gap/Utils.gi
+++ b/gap/Utils.gi
@@ -9,39 +9,6 @@ function(field, nrRows, nrCols, entries)
     return ImmutableMatrix(field, m);
 end);
 
-# Note that <group> must have the attribute InvariantBilinearForm or
-# InvariantSesquilinearForm for this to work.
-InstallGlobalFunction("ChangeFixedSesquilinearForm",
-function(group, type, gramMatrix)
-    local gapForm, newForm, gapToCanonical, canonicalToNew, field;
-    if not type in ["S", "O", "U"] then
-        ErrorNoReturn("<type> must be one of 'S', 'U', 'O', but <type> = ",
-                      type);
-    fi;
-    field := DefaultFieldOfMatrixGroup(group);
-    if type = "S" or type = "O" then
-        gapForm := BilinearFormByMatrix(InvariantBilinearForm(group).matrix, 
-                                        field);
-        newForm := BilinearFormByMatrix(gramMatrix, field);
-    else
-        gapForm := HermitianFormByMatrix(InvariantSesquilinearForm(group).matrix,
-                                         field);
-        newForm := HermitianFormByMatrix(gramMatrix, field);
-    fi;
-    # the following if condition can only ever be fulfilled if <group> is an
-    # orthogonal group; there the case of even dimension is problematic since,
-    # in that case, there are two similarity classes of bilinear forms
-    if not WittIndex(gapForm) = WittIndex(newForm) then
-       ErrorNoReturn("The form preserved by <group> must be similar to the form ", 
-                     "described by the Gram matrix <gramMatrix>.");
-    fi;
-    gapToCanonical := BaseChangeHomomorphism(BaseChangeToCanonical(gapForm), 
-                                             field);
-    canonicalToNew := BaseChangeHomomorphism(BaseChangeToCanonical(newForm) ^ (-1), 
-                                             field);
-    return Group(canonicalToNew(gapToCanonical(GeneratorsOfGroup(group))));
-end);
-
 InstallGlobalFunction("AntidiagonalMat",
 function(entries, field)
     local d, m, i;
@@ -184,6 +151,10 @@ ReflectionMatrix := function(n, q, gramMatrix, v)
     local reflectionMatrix, i, basisVector, reflectBasisVector, beta;
     reflectionMatrix := NullMat(n, n, GF(q));
     beta := BilinearFormByMatrix(gramMatrix);
+    if IsZero(EvaluateForm(beta, v, v)) then
+        ErrorNoReturn("The vector <v> must have non-zero norm with respect to",
+                      " the bilinear form given by <gramMatrix>");
+    fi;
     for i in [1..n] do
         basisVector := List([1..n], j -> 0 * Z(q));
         basisVector[i] := Z(q) ^ 0;
diff --git a/init.g b/init.g
index 569539a7..24a39513 100644
--- a/init.g
+++ b/init.g
@@ -10,3 +10,4 @@ ReadPackage( "ClassicalMaximals", "gap/ReducibleMatrixGroups.gd");
 ReadPackage( "ClassicalMaximals", "gap/ImprimitiveMatrixGroups.gd");
 ReadPackage( "ClassicalMaximals", "gap/SemilinearMatrixGroups.gd");
 ReadPackage( "ClassicalMaximals", "gap/SubfieldMatrixGroups.gd");
+ReadPackage( "ClassicalMaximals", "gap/Forms.gd");
diff --git a/read.g b/read.g
index 7159511f..629624d4 100644
--- a/read.g
+++ b/read.g
@@ -14,3 +14,4 @@ ReadPackage( "ClassicalMaximals", "gap/TensorInducedMatrixGroups.gi");
 ReadPackage( "ClassicalMaximals", "gap/ClassicalNormalizerMatrixGroups.gi");
 ReadPackage( "ClassicalMaximals", "gap/ClassicalMaximals.gi");
 ReadPackage( "ClassicalMaximals", "gap/AlmostSimpleDefiningCharacteristic.gi");
+ReadPackage( "ClassicalMaximals", "gap/Forms.gi");
diff --git a/tst/quick/ClassicalMaximals.tst b/tst/quick/ClassicalMaximals.tst
deleted file mode 100644
index 9897d2d8..00000000
--- a/tst/quick/ClassicalMaximals.tst
+++ /dev/null
@@ -1,2 +0,0 @@
-gap> ClassicalMaximalsGeneric("L", 2, 3);
-Error, SL(2, 2) and SL(2, 3) are soluble
diff --git a/tst/quick/ClassicalNormalizerMatrixGroups.tst b/tst/quick/ClassicalNormalizerMatrixGroups.tst
deleted file mode 100644
index 7aa48955..00000000
--- a/tst/quick/ClassicalNormalizerMatrixGroups.tst
+++ /dev/null
@@ -1,56 +0,0 @@
-gap> n := 4;; q := 3;;
-gap> G := SymplecticNormalizerInSL(n, q);;
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> n := 4;; q := 5;;
-gap> G := SymplecticNormalizerInSL(n, q);;
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> n := 6;; q := 4;;
-gap> G := SymplecticNormalizerInSL(n, q);;
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> n := 4;; q := 9;;
-gap> G := UnitaryNormalizerInSL(n, q);;
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> n := 3;; q := 9;;
-gap> G := UnitaryNormalizerInSL(n, q);;
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> n := 4;; q := 4;;
-gap> G := UnitaryNormalizerInSL(n, q);;
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> epsilon := 0;; n := 3;; q := 5;;
-gap> G := OrthogonalNormalizerInSL(epsilon, n, q);;
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> epsilon := -1;; n := 6;; q := 5;;
-gap> G := OrthogonalNormalizerInSL(epsilon, n, q);;
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> epsilon := 1;; n := 6;; q := 5;;
-gap> G := OrthogonalNormalizerInSL(epsilon, n, q);;
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> epsilon := -1;; n := 4;; q := 3;;
-gap> G := OrthogonalNormalizerInSL(epsilon, n, q);;
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> epsilon := 1;; n := 4;; q := 3;;
-gap> G := OrthogonalNormalizerInSL(epsilon, n, q);;
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> epsilon := -1;; n := 4;; q := 5;;
-gap> G := OrthogonalNormalizerInSL(epsilon, n, q);;
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> epsilon := 1;; n := 4;; q := 5;;
-gap> G := OrthogonalNormalizerInSL(epsilon, n, q);;
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> epsilon := -1;; n := 6;; q := 3;;
-gap> G := OrthogonalNormalizerInSL(epsilon, n, q);;
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
diff --git a/tst/quick/Forms.tst b/tst/quick/Forms.tst
new file mode 100644
index 00000000..dd626df1
--- /dev/null
+++ b/tst/quick/Forms.tst
@@ -0,0 +1,54 @@
+gap> UnitaryForm(SU(4, 3)) = InvariantSesquilinearForm(SU(4, 3)).matrix;
+true
+gap> TestFormChangingFunctions := function(args)
+>   local n, q, type, gramMatrix, standardGroup, conjugatedGroup, broadType,
+>   standardGramMatrix, twiceConjugatedGroup;
+>   n := args[1];
+>   q := args[2];
+>   type := args[3];
+>   gramMatrix := args[4];
+>   if type = "U" then
+>       standardGroup := SU(n, q);
+>   elif type = "S" then
+>       standardGroup := Sp(n, q);
+>   elif type = "O" then
+>       standardGroup := SO(n, q);
+>   elif type = "O+" then
+>       standardGroup := SO(1, n, q);
+>   elif type = "O-" then
+>       standardGroup := SO(-1, n, q);
+>   fi;
+>   if type in ["O", "O+", "O-"] then
+>       broadType := "O";
+>   else
+>       broadType := type;
+>   fi;
+>   conjugatedGroup := ChangeFixedSesquilinearForm(standardGroup, broadType, gramMatrix);
+>   if type = "U" then
+>       standardGramMatrix := InvariantSesquilinearForm(standardGroup).matrix;
+>       SetInvariantSesquilinearForm(conjugatedGroup, rec(matrix := gramMatrix));
+>   else
+>       standardGramMatrix := InvariantBilinearForm(standardGroup).matrix;
+>       SetInvariantBilinearForm(conjugatedGroup, rec(matrix := gramMatrix));
+>   fi;
+>   twiceConjugatedGroup := ConjugateToStandardForm(conjugatedGroup, type);
+>   if type = "U" then
+>       return ForAll(GeneratorsOfGroup(conjugatedGroup), 
+>                     g -> g * gramMatrix * HermitianConjugate(g, q) = gramMatrix)
+>              and ForAll(GeneratorsOfGroup(twiceConjugatedGroup), 
+>                         g -> g * standardGramMatrix * HermitianConjugate(g, q) = standardGramMatrix);
+>   else
+>       return ForAll(GeneratorsOfGroup(conjugatedGroup),
+>                     g -> g * gramMatrix * TransposedMat(g) = gramMatrix)
+>              and ForAll(GeneratorsOfGroup(twiceConjugatedGroup),
+>                         g -> g * standardGramMatrix * TransposedMat(g) = standardGramMatrix);
+>   fi;
+> end;;
+gap> testsFormChangingFunctions := [[3, 7, "U", IdentityMat(3, GF(7))],
+>                                   [6, 3, "S", AntidiagonalMat(Z(3) ^ 0 * [1, -1, 1, -1, 1, -1], GF(3))],
+>                                   [5, 5, "O", IdentityMat(5, GF(5))],
+>                                   [4, 7, "O+", AntidiagonalMat(4, GF(7))],
+>                                   [4, 7, "O-", Z(7) ^ 0 * DiagonalMat([Z(7), 1, 1, 1])],
+>                                   [6, 7, "O-", IdentityMat(6, GF(7))]];;
+gap> ForAll(testsFormChangingFunctions, TestFormChangingFunctions);
+true
diff --git a/tst/quick/TensorInducedMatrixGroups.tst b/tst/quick/TensorInducedMatrixGroups.tst
deleted file mode 100644
index 28901837..00000000
--- a/tst/quick/TensorInducedMatrixGroups.tst
+++ /dev/null
@@ -1,24 +0,0 @@
-gap> m := 3;; t := 2;; q := 5;;
-gap> G := TensorInducedDecompositionStabilizerInSL(m, t, q);;
-gap> IsSubset(SL(m ^ t, q), GeneratorsOfGroup(G));
-true
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> m := 2;; t := 2;; q := 7;;
-gap> G := TensorInducedDecompositionStabilizerInSL(m, t, q);;
-gap> IsSubset(SL(m ^ t, q), GeneratorsOfGroup(G));
-true
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> m := 2;; t := 2;; q := 5;;
-gap> G := TensorInducedDecompositionStabilizerInSL(m, t, q);;
-gap> IsSubset(SL(m ^ t, q), GeneratorsOfGroup(G));
-true
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
-gap> m := 3;; t := 3;; q := 3;;
-gap> G := TensorInducedDecompositionStabilizerInSL(m, t, q);;
-gap> IsSubset(SL(m ^ t, q), GeneratorsOfGroup(G));
-true
-gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
-true
diff --git a/tst/quick/Utils.tst b/tst/quick/Utils.tst
index 2cfcf81c..e90a2667 100644
--- a/tst/quick/Utils.tst
+++ b/tst/quick/Utils.tst
@@ -1,3 +1,84 @@
-gap> M := MatrixByEntries(GF(2), 2, 2, [[1,2,1],[2,1,1]]);;
-gap> IsOne(M^2);
+gap> M := MatrixByEntries(GF(2), 2, 2, [[1, 2, 1], [2, 1, 1]]);;
+gap> IsOne(M ^ 2);
+true
+gap> M := AntidiagonalMat(7, GF(9));;
+gap> IsOne(M ^ 2);
+true
+gap> M := GL(7, 5 ^ 2).1;;
+gap> N := ApplyFunctionToEntries(M, x -> x ^ 5);;
+gap> M = ApplyFunctionToEntries(N, x -> x ^ 5);
+true
+gap> M := GU(5, 7).1;;
+gap> M = HermitianConjugate(HermitianConjugate(M, 7), 7);
+true
+gap> sol := SolveQuadraticCongruence(Z(7), 7);;
+gap> sol.a ^ 2 + sol.b ^ 2 = Z(7);
+true
+gap> M := SquareSingleEntryMatrix(GF(9), 3, 1, 3);;
+gap> IsZero(M ^ 2);
+true
+gap> QuoCeil(5, 3) = QuoCeil(6, 3);
+true
+gap> M := ReflectionMatrix(5, 9, AntidiagonalMat(5, GF(9)), Z(9) ^ 0 * [1, 1, 1, 1, 1]);; 
+gap> IsOne(M ^ 2);
+true
+gap> x := SolveFrobeniusEquation("S", - Z(7) ^ 0, 7);;
+gap> x + x ^ 7 = - Z(7) ^ 0;
+true
+gap> x := SolveFrobeniusEquation("P", - Z(7) ^ 0, 7);;
+gap> x * x ^ 7 = - Z(7) ^ 0;
+true
+gap> TestGeneratorsOfOrthogonalGroup := function(args)
+>   local epsilon, n, q, F, zeta, gen, gens, rightDets, gramMatrix, rightForm;
+>   epsilon := args[1];
+>   n := args[2];
+>   q := args[3];
+>   F := GF(q);
+>   zeta := PrimitiveElement(F);
+>   gens := GeneratorsOfOrthogonalGroup(epsilon, n, q);
+>   rightDets := true;
+>   rightForm := true;
+>   for gen in gens.generatorsOfSO do
+>       if not IsOne(Determinant(gen)) then
+>           rightDets := false;
+>       fi;
+>   od;
+>   if not IsOne(- Determinant(gens.D)) then
+>       rightDets := false;
+>   fi;
+>   if epsilon <> 0 then
+>       if not Determinant(gens.E) = (epsilon * zeta) ^ (n / 2) then
+>           rightDets := false;
+>       fi;
+>   fi;
+>   if epsilon = 0 then
+>       gramMatrix := IdentityMat(n, F);
+>   elif epsilon = 1 then
+>       gramMatrix := AntidiagonalMat(n, F);
+>   elif epsilon = -1 then
+>       if IsOddInt(n * (q - 1) / 4) then
+>           gramMatrix := IdentityMat(n, F);
+>       else
+>           gramMatrix := IdentityMat(n, F);
+>           gramMatrix[1, 1] := zeta;
+>       fi;
+>   fi;
+>   for gen in Concatenation(gens.generatorsOfSO, [gens.D]) do
+>       if not gen * gramMatrix * TransposedMat(gen) = gramMatrix then
+>           rightForm := false;
+>       fi;
+>   od;
+>   if epsilon = 0 then
+>       if not gens.E * gramMatrix * TransposedMat(gens.E) = zeta ^ 2 * gramMatrix then
+>           rightForm := false;
+>       fi;
+>   else
+>       if not gens.E * gramMatrix * TransposedMat(gens.E) = zeta * gramMatrix then
+>           rightForm := false;
+>       fi;
+>   fi;
+>   return (rightDets and rightForm);
+> end;;
+gap> testsGeneratorsOfOrthogonalGroup := [[-1, 6, 7], [-1, 4, 9], [1, 4, 7], [0, 5, 5]];;
+gap> ForAll(testsGeneratorsOfOrthogonalGroup, TestGeneratorsOfOrthogonalGroup);
 true
diff --git a/tst/standard/ClassicalMaximals.tst b/tst/standard/ClassicalMaximals.tst
new file mode 100644
index 00000000..b1ed8c4c
--- /dev/null
+++ b/tst/standard/ClassicalMaximals.tst
@@ -0,0 +1,28 @@
+gap> TestClassicalMaximalsUnitary := function(args)
+>   return Length(ClassicalMaximalsGeneric("U", args[1], args[2])) = args[3];
+> end;;
+gap> testsClassicalMaximalsUnitary := [[3, 3, 3], [3, 4, 4], [3, 5, 2], 
+>                                      [3, 7, 5], [3, 8, 7], [3, 9, 5], 
+>                                      [3, 11, 8], [3, 13, 5], [3, 16, 4],
+>                                      [3, 17, 10], [3, 19, 5],
+>                                      [4, 2, 5], [4, 3, 10], [4, 4, 7],
+>                                      [4, 5, 10], [4, 7, 16], [4, 8, 8],
+>                                      [4, 9, 10], [4, 11, 14], [4, 13, 10],
+>                                      [4, 16, 7], [4, 17, 10], [4, 19, 14],
+>                                      [5, 2, 5], [5, 3, 7], [5, 4, 11], 
+>                                      [5, 5, 7], [5, 7, 7], [5, 8, 7],
+>                                      [5, 9, 16], [5, 11, 7], [5, 13, 7], 
+>                                      [5, 16, 6], [5, 17, 7], [5, 19, 16],
+>                                      [6, 2, 10], [6, 3, 14], [6, 4, 12],
+>                                      [6, 5, 20], [6, 7, 14], [6, 8, 17],
+>                                      [7, 2, 8], [7, 3, 9], [7, 8, 9],
+>                                      [7, 13, 22],
+>                                      [8, 2, 11], [8, 3, 25], [8, 4, 13],
+>                                      [8, 5, 17], [8, 8, 14],
+>                                      [9, 2, 14], [9, 3, 13], [9, 4, 12],
+>                                      [9, 5, 20], [9, 8, 17],
+>                                      [10, 2, 14], [10, 3, 18],
+>                                      [11, 2, 12], [11, 3, 13],
+>                                      [12, 2, 21], [12, 3, 27]];;
+gap> ForAll(testsClassicalMaximalsUnitary, TestClassicalMaximalsUnitary);
+true
diff --git a/tst/standard/ExtraspecialNormalizerMatrixGroups.tst b/tst/standard/ExtraspecialNormalizerMatrixGroups.tst
index 7bebf0d2..82115992 100644
--- a/tst/standard/ExtraspecialNormalizerMatrixGroups.tst
+++ b/tst/standard/ExtraspecialNormalizerMatrixGroups.tst
@@ -94,3 +94,80 @@ gap> IsSubset(SU(r ^ m, q), GeneratorsOfGroup(G));
 true
 gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
 true
+gap> TestOddExtraspecialGroup := function(args)
+>   local r, m, q, gens;
+>   r := args[1];
+>   m := args[2];
+>   q := args[3];
+>   gens := OddExtraspecialGroup(r, m, q);
+>   return Size(Group(Concatenation(gens.listOfXi, gens.listOfYi))) = r ^ (2 * m + 1);
+> end;;
+gap> testsOddExtraspecialGroup := [[5, 1, 11], [3, 1, 7], [3, 2, 13]];;
+gap> ForAll(testsOddExtraspecialGroup, TestOddExtraspecialGroup);
+true
+gap> TestOddExtraspecialNormalizerInGL := function(args)
+>   local r, m, q, gensNormalizer, gensExtraspecial, extraspecialGroup,
+>   normalizer;
+>   r := args[1];
+>   m := args[2];
+>   q := args[3];
+>   gensNormalizer := OddExtraspecialNormalizerInGL(r, m, q);
+>   gensExtraspecial := OddExtraspecialGroup(r, m, q);
+>   extraspecialGroup := Group(Concatenation(gensExtraspecial.listOfXi,
+>                                            gensExtraspecial.listOfYi));
+>   normalizer := Group(Concatenation(gensNormalizer.listOfXi,
+>                                     gensNormalizer.listOfYi,
+>                                     gensNormalizer.listOfUi,
+>                                     gensNormalizer.listOfVi,
+>                                     gensNormalizer.listOfWi,
+>                                     [gensNormalizer.generatingScalar]));
+>   return ForAll(GeneratorsOfGroup(normalizer), g -> extraspecialGroup ^ g = extraspecialGroup);
+> end;;
+gap> testsOddExtraspecialNormalizerInGL := [[5, 1, 11], [3, 1, 7], [3, 2, 13]];;
+gap> ForAll(testsOddExtraspecialNormalizerInGL, TestOddExtraspecialNormalizerInGL);
+true
+gap> TestSymplecticTypeNormalizerInGL := function(args)
+>   local r, m, q, gensNormalizer, gensExtraspecial, extraspecialGroup,
+>   normalizer, zeta, F;
+>   r := 2;
+>   m := args[1];
+>   q := args[2];
+>   F := GF(q);
+>   zeta := PrimitiveElement(F);
+>   gensNormalizer := SymplecticTypeNormalizerInGL(m, q);
+>   gensExtraspecial := OddExtraspecialGroup(r, m, q);
+>   extraspecialGroup := Group(Concatenation(gensExtraspecial.listOfXi,
+>                                            gensExtraspecial.listOfYi,
+>                                            [zeta ^ ((q - 1) / 4) * IdentityMat(r ^ m, F)]));
+>   normalizer := Group(Concatenation(gensNormalizer.listOfXi,
+>                                     gensNormalizer.listOfYi,
+>                                     gensNormalizer.listOfUi,
+>                                     gensNormalizer.listOfVi,
+>                                     gensNormalizer.listOfWi,
+>                                     [gensNormalizer.generatingScalar]));
+>   return ForAll(GeneratorsOfGroup(normalizer), g -> extraspecialGroup ^ g = extraspecialGroup);
+> end;;
+gap> testsSymplecticTypeNormalizerInGL := [[3, 5], [2, 5], [2, 9]];;
+gap> ForAll(testsSymplecticTypeNormalizerInGL, TestSymplecticTypeNormalizerInGL);
+true
+gap> TestExtraspecial2MinusTypeNormalizerInGL := function(args)
+>   local r, m, q, gensNormalizer, extraspecialGroup, normalizer, zeta, F;
+>   r := 2;
+>   m := args[1];
+>   q := args[2];
+>   F := GF(q);
+>   zeta := PrimitiveElement(F);
+>   gensNormalizer := Extraspecial2MinusTypeNormalizerInGL(m, q);
+>   extraspecialGroup := Group(Concatenation(gensNormalizer.listOfXi,
+>                                            gensNormalizer.listOfYi));
+>   normalizer := Group(Concatenation(gensNormalizer.listOfXi,
+>                                     gensNormalizer.listOfYi,
+>                                     gensNormalizer.listOfUi,
+>                                     gensNormalizer.listOfVi,
+>                                     gensNormalizer.listOfWi,
+>                                     [gensNormalizer.generatingScalar]));
+>   return ForAll(GeneratorsOfGroup(normalizer), g -> extraspecialGroup ^ g = extraspecialGroup);
+> end;;
+gap> testsExtraspecial2MinusTypeNormalizerInGL := [[1, 9], [1, 5], [1, 7], [2, 5]];;
+gap> ForAll(testsExtraspecial2MinusTypeNormalizerInGL, TestExtraspecial2MinusTypeNormalizerInGL);
+true
diff --git a/tst/quick/ImprimitiveMatrixGroups.tst b/tst/standard/ImprimitiveMatrixGroups.tst
similarity index 100%
rename from tst/quick/ImprimitiveMatrixGroups.tst
rename to tst/standard/ImprimitiveMatrixGroups.tst
diff --git a/tst/standard/SemilinearMatrixGroups.tst b/tst/standard/SemilinearMatrixGroups.tst
index e23463aa..b12c2d81 100644
--- a/tst/standard/SemilinearMatrixGroups.tst
+++ b/tst/standard/SemilinearMatrixGroups.tst
@@ -22,3 +22,31 @@ gap> IsSubset(SL(n, q), GeneratorsOfGroup(G));
 true
 gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
 true
+gap> n := 3;; q := 5;; s := 3;;
+gap> G := GammaLMeetSU(n, q, s);;
+gap> IsSubset(SU(n, q), GeneratorsOfGroup(G));
+true
+gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
+true
+gap> n := 6;; q := 3;; s := 3;;
+gap> G := GammaLMeetSU(n, q, s);;
+gap> IsSubset(SU(n, q), GeneratorsOfGroup(G));
+true
+gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
+true
+gap> n := 3;; q := 7;; s := 3;;
+gap> G := GammaLMeetSU(n, q, s);;
+gap> IsSubset(SU(n, q), GeneratorsOfGroup(G));
+true
+gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
+true
+gap> TestGammaLDimension1 := function(args)
+>   local q, s, gens;
+>   q := args[1];
+>   s := args[2];
+>   gens := GammaLDimension1(s, q);
+>   return Order(gens.A) = (q ^ s - 1) and Order(gens.B) = s;
+> end;;
+gap> testsGammaLDimension1 := [[3, 2], [2, 2], [5, 3], [4, 3]];;
+gap> ForAll(testsGammaLDimension1, TestGammaLDimension1);
+true
diff --git a/tst/standard/TensorInducedMatrixGroups.tst b/tst/standard/TensorInducedMatrixGroups.tst
new file mode 100644
index 00000000..b5a101bf
--- /dev/null
+++ b/tst/standard/TensorInducedMatrixGroups.tst
@@ -0,0 +1,50 @@
+gap> m := 3;; t := 2;; q := 5;;
+gap> G := TensorInducedDecompositionStabilizerInSL(m, t, q);;
+gap> IsSubset(SL(m ^ t, q), GeneratorsOfGroup(G));
+true
+gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
+true
+gap> m := 2;; t := 2;; q := 7;;
+gap> G := TensorInducedDecompositionStabilizerInSL(m, t, q);;
+gap> IsSubset(SL(m ^ t, q), GeneratorsOfGroup(G));
+true
+gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
+true
+gap> m := 2;; t := 2;; q := 5;;
+gap> G := TensorInducedDecompositionStabilizerInSL(m, t, q);;
+gap> IsSubset(SL(m ^ t, q), GeneratorsOfGroup(G));
+true
+gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
+true
+gap> m := 3;; t := 3;; q := 3;;
+gap> G := TensorInducedDecompositionStabilizerInSL(m, t, q);;
+gap> IsSubset(SL(m ^ t, q), GeneratorsOfGroup(G));
+true
+gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
+true
+gap> m := 2;; t := 2;; q := 7;;
+gap> G := TensorInducedDecompositionStabilizerInSU(m, t, q);;
+gap> IsSubset(SU(m ^ t, q), GeneratorsOfGroup(G));
+true
+gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
+true
+gap> m := 2;; t := 2;; q := 5;;
+gap> G := TensorInducedDecompositionStabilizerInSU(m, t, q);;
+gap> IsSubset(SU(m ^ t, q), GeneratorsOfGroup(G));
+true
+gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
+true
+gap> m := 3;; t := 2;; q := 3;;
+gap> G := TensorInducedDecompositionStabilizerInSU(m, t, q);;
+gap> IsSubset(SU(m ^ t, q), GeneratorsOfGroup(G));
+true
+gap> Size(Group(GeneratorsOfGroup(G))) = Size(G);
+true
+gap> m := 3;; t := 3;; q := 3;;
+gap> G := TensorInducedDecompositionStabilizerInSU(m, t, q);;
+gap> IsSubset(SU(m ^ t, q), GeneratorsOfGroup(G));
+true
+gap> m := 3;; t := 2;; q := 5;;
+gap> G := TensorInducedDecompositionStabilizerInSU(m, t, q);;
+gap> IsSubset(SU(m ^ t, q), GeneratorsOfGroup(G));
+true