diff --git a/.gitignore b/.gitignore
index 85d5146..dbdc569 100644
--- a/.gitignore
+++ b/.gitignore
@@ -22,3 +22,5 @@
 /tmp
 /gh-pages
 tmp.json
+profiling
+benchmark.csv
diff --git a/PackageInfo.g b/PackageInfo.g
index 68f57cc..4ae14a1 100644
--- a/PackageInfo.g
+++ b/PackageInfo.g
@@ -99,7 +99,7 @@ PackageDoc := rec(
 
 Dependencies := rec(
   GAP := ">=4.9",
-  NeededOtherPackages := [],
+  NeededOtherPackages := [["nofoma", "1.0"]],
   SuggestedOtherPackages := [],
   ExternalConditions := []),
 
diff --git a/doc/forms.xml b/doc/forms.xml
index 1e2a0c8..ee31f10 100644
--- a/doc/forms.xml
+++ b/doc/forms.xml
@@ -41,6 +41,8 @@
 <#Include SYSTEM "theory.xml">
 <#Include SYSTEM "construct_basic.xml">
 <#Include SYSTEM "morphisms.xml">
+<#Include SYSTEM "_Chapter_Formspace.xml">
+
 </Body>
 
 <Bibliography Databases="forms"/>
diff --git a/init.g b/init.g
index 6f21851..2c82f60 100644
--- a/init.g
+++ b/init.g
@@ -14,3 +14,4 @@
 
 ReadPackage("forms","lib/forms.gd");
 ReadPackage("forms","lib/recognition.gd");
+ReadPackage("forms","lib/formspace.gd");
\ No newline at end of file
diff --git a/lib/forms.gi b/lib/forms.gi
index fa1d6ad..f50bcce 100644
--- a/lib/forms.gi
+++ b/lib/forms.gi
@@ -1699,22 +1699,66 @@ InstallGlobalFunction(Forms_QUAD_EQ,
 
 InstallMethod( FORMS_IsSymplecticMatrix, [IsFFECollColl, IsField],
   function(m,f)
+    local i, j, n;
+    # maybe thorw an error if m is not square
+    n := NrRows(m);
     if Characteristic(f) = 2 then
-       if not ForAll([1..NrRows(m)], i -> IsZero(m[i,i]) ) then
-           return false;
-       fi;
+      # is this needed??? i am not sure how symplectic matrices are definied in char 2... i thought A = A^T is enough!
+      if not ForAll([1..n], i -> IsZero(m[i,i]) ) then
+        return false;
+      fi;
+      for i in [1..n] do
+        for j in [1..i] do
+          if m[i, j] <> m[j, i] then
+            return false;
+          fi;
+        od;
+      od;
+      return true;
     fi;
-    return m = -TransposedMat(m);
+    for i in [1..n] do
+      for j in [1..i] do
+        if m[i, j] <> -m[j, i] then
+          return false;
+        fi;
+      od;
+    od;
+    return true;
   end );
 
 InstallMethod( FORMS_IsSymmetricMatrix, [IsFFECollColl],
   function(m)
-    return m=TransposedMat(m);
+    local n, i, j;
+    n := NrRows(m);
+    for i in [1..n] do
+      # minus one to ignore diagonal entries for symmetric matrices.
+      for j in [1..i - 1] do
+        if m[i, j] <> m[j, i] then
+          return false;
+        fi;
+      od;
+    od;
+    return true;
   end );
 
 InstallMethod( FORMS_IsHermitianMatrix, [IsFFECollColl, IsField],
   function(m,f)
-    return m=Forms_HERM_CONJ(m,Sqrt(Size(f)));
+    local p, q, n, i, j;
+    if DegreeOverPrimeField(f) mod 2 <> 0 then
+      # should this be the behaviour? could also return fail.
+      return false;
+    fi;
+    p := Characteristic(f);
+    q := p^(DegreeOverPrimeField(f)/2);
+    n := NrRows(m);
+    for i in [1..n] do
+      for j in [1..i] do
+        if m[i, j] <> m[j, i]^q then
+          return false;
+        fi;
+      od;
+    od;
+    return true;
   end );
 
 #############################################################################
diff --git a/lib/formspace.gd b/lib/formspace.gd
new file mode 100644
index 0000000..c896295
--- /dev/null
+++ b/lib/formspace.gd
@@ -0,0 +1,8 @@
+DeclareOperation("PreservedFormspace", [IsMatrixGroup, IsVector and IsFFECollection, IsBool]);
+DeclareOperation("PreservedFormspace", [IsMatrixGroup]);
+
+DeclareOperation("PreservedFormsWithScalars", [IsMatrixGroup]);
+DeclareOperation("PreservedSesquilinearFormsWithScalars", [IsMatrixGroup]);
+
+
+DeclareOperation("FilterFormspace", [IsList, IsFinite and IsField, IsBool]);
diff --git a/lib/formspace.gi b/lib/formspace.gi
new file mode 100644
index 0000000..2f1be84
--- /dev/null
+++ b/lib/formspace.gi
@@ -0,0 +1,790 @@
+#! @Chapter Formspace
+#! For a matrix Group $G$ consisting of $n\times n$ matrices over the field $K$. Let $\Lambda : G \to K^\times $ be a group homomorphism. Suppose $h : K \to K$ is a field automorphism of $K$. For $g\in G$ we define $g^* := h(g^T)$ where the automorphism $h$ is applied entrywise.
+#! We define the Formspace $\mathcal{F}_h(G, \Lambda) := \{A\in K^{n\times n} \mid gAg^* = \Lambda(g)A \forall g\in G\}$.
+#! @Section Computing the Formspace
+#! For details on how the form space is computed SEE: somewhere that does not exist yet
+
+# underlying field F, g is in F^{n\times n}, v is in F^n, coeffs are the coefficients of a polynomial p in F[X]. Returns vp(g).
+FORMS_EvaluateMatrixPolynomialWithVector := function(F, n, g, v, coeffs)
+    local res, i, deg;
+    deg := Size(coeffs);
+
+    if deg = 0 then
+        return ZeroVector(F, n);
+    fi;
+    if deg = 1 then
+        return v * coeffs[1];
+    fi;
+    
+    res := v * g * coeffs[deg];
+    for i in [1..deg-2] do
+        res := res + coeffs[deg - i]*v;
+        res := res * g;
+    od;
+    res := res + coeffs[1]*v;
+    return res;
+end;
+
+
+# Given mat in F^{n\times n} F Field, n \in N, mode says whether to apply hom or not. Returns (mat^T)^hom.
+# mode = False \iff mat* := mat^tr, mode = True \iff mat* = komplex_konjugiert(mat^tr)
+FORMS_CalculateAdjoint := function(mat, mode, hom, n, F)
+    local transposed, i, j;
+    transposed := TransposedMat(mat);
+    # ConvertToMatrixRep(transposed, F);
+    if mode = false then
+        return transposed;
+    fi;
+    if mode then
+        return transposed^hom;
+    fi;
+    return fail;
+end;
+
+# tries to find a element g \in <Gens> such that the Frobenius Normal form of g has as few blocks as possible. Lambdas describes a group homomorphism induced by Phi : Gens[i] \mapsto Lambdas[i]
+# Returns [g, Phi(g), FrobeniusNormalForm(g), nrOfTries]
+# nrOfTries contains the number of random elements tested before g was found.
+FORMS_FindCyclicGroupElementAndScalars := function(Gens, Lambdas, n)
+    local known_elements, known_scalars, i, j, run, runs_max, best_known_length, best_known_res, best_known_scalar, best_known_element, g, lamb, s, d, res, interval;
+
+    d := Size(Gens);
+
+    # no point to random group generators in cyclic groups.
+    if d = 1 then 
+        return Concatenation([Gens[1], Lambdas[1]], FrobeniusNormalForm(Gens[1]), [-1]); 
+    fi;
+    runs_max := 3*d + 15; # pretty random idk what to put here
+
+    known_elements := ShallowCopy(Gens);
+    known_scalars := ShallowCopy(Lambdas);
+
+    best_known_length := n + 1;
+
+    for run in [1..runs_max] do
+        i := PseudoRandom([1..d]);
+        interval := [1..d];
+        Remove(interval, i);
+        j := PseudoRandom(interval);
+        # TODO: add inverses? seems fine without
+        g := known_elements[i]*known_elements[j];
+        lamb := known_scalars[i]*known_scalars[j];
+        res := FrobeniusNormalForm(g);
+        s := Size(res[3]);
+        if s = 1 then
+            return Concatenation([g, lamb], res, [run]);
+        fi;
+        if s < best_known_length then
+            best_known_length := s;
+            best_known_element := g;
+            best_known_scalar := lamb;
+            best_known_res := res;
+        fi;
+        known_elements[i] := g;
+        known_scalars[i] := lamb;
+    od;
+   
+    if best_known_length = n + 1 then
+        return Concatenation([Gens[1], Lambdas[1]], FrobeniusNormalForm(Gens[1]), [-1]);
+    fi;
+    return Concatenation([best_known_element, best_known_scalar], best_known_res, [runs_max]);
+end;
+
+# turns the (jn) vector vec in F^{jn} and returns a F^{j times n} matrix
+FORMS_VectorReorganize := function(vec, j, F, n)
+    local A, i;
+    A := NullMat(j, n, F);
+    # A := ZeroMatrix(F, j, n);
+    # ConvertToMatrixRep(A, F);
+    for i in [1..j] do
+        A[i] := vec{[((i - 1) * n + 1)..(i*n)]};
+    od;
+    return A;
+end;
+
+# turns the F^{j times n} matrix into F^{jn} vector
+FORMS_MatrixReorganize := function(mat, j, F, n)
+    local vec, i;
+    vec := ZeroVector(F, j*n);
+    for i in [1..j] do
+        vec{[((i - 1) * n + 1)..(i*n)]} := mat[i];
+    od;
+    return vec;
+end;
+
+# given j times n matrices over field F, this tries to compute a linear combination of the matrices such that their sum is zero.
+# TODO: maybe make this funciton more efficient by only considereing some vectors, a better optimatizion however, would be to not even generate the unused equations
+FORMS_SolveMatrixSystem := function(mats, j, n, F)
+    local eqs, mat, sol, out;
+    eqs := [];
+    for mat in mats do 
+        Add(eqs, FORMS_MatrixReorganize(mat, j, F, n));
+    od;
+    return NullspaceMatDestructive(eqs);
+end;
+
+# Spins a block in the frobenius normal form, with the images.
+FORMS_FrobSpin := function(Images, spin_elem, frob_base_blocks, n, F)
+    local A, j, i, k, end_pos;
+    j := Size(frob_base_blocks);
+    A := NullMat(n, n, F);
+    # A := ZeroMatrix(F, n, n);
+    CopySubMatrix(Images, A, [1..j], frob_base_blocks, [1..n], [1..n]);
+    for i in [1..j] do
+        if i = j then
+            end_pos := n;
+        else
+            end_pos := frob_base_blocks[i + 1] - 1;
+        fi;
+        # causes error if A is matrix obj and Images is plain list... :(( 
+        # A[frob_base_blocks[i]] := Images[i];
+        for k in [(frob_base_blocks[i] + 1)..end_pos] do
+            CopySubMatrix([A[k - 1]*spin_elem], A, [1], [k], [1..n], [1..n]);
+            # A[k] := A[k - 1]*spin_elem;
+        od;
+    od;
+    return A;
+end;
+
+# Evaluates the univariate polynomial p (given as coefficients) in the matrix g \in F^{n\times n}. frob_base = FrobeniusNormalForm(g) must be satisfied and frob_base_inv = Inverse(FrobeniusNormalForm(g)[2]). The reason these two parameters are given, and not computed in the function itself is to not compute FrobeniusNormalForm(g) multiple times when evaluating multiple polynomials in g.
+FORMS_EvaluatePolynomialWithFrobenius := function(p, g, frob_base, frob_base_inv, F, n) 
+    local ws, C, i, end_pos, j, k;
+    ws := [];
+    j := Size(frob_base[3]);
+    for k in [1..j] do
+        Add(ws, FORMS_EvaluateMatrixPolynomialWithVector(F, n, g, frob_base[2][frob_base[3][k]]{[1..n]}, p));
+    od;
+    C := FORMS_FrobSpin(ws, g, frob_base[3], n, F);
+    # TODO: this function is used to build the condition matrices, however since we are only interested in the kernels of these matrices and C is an invertible matrix we can probably eliminate this multiplication and do it at a later stage?
+    return frob_base_inv * C;
+end;
+
+# This computes (and returns) the matrix \mathcal{P}_{h, u} from the bachelors thesis.
+# Here we have u, h \in F^{n\times n}, h_star = h^*, scalar_h = \lambda_h, g_star_inv_scaled = g^{-*}*lambda_g, frob_base = FrobeniusNormalForm(g^{-*}), frob_base_inv = Inverse(FrobeniusNormalForm(g)[2]), frob_base_inv_star = Inverse(frob_base[2]). The reason we to give so many parameters is to avoid computing the same matrices multiple times. TODO: better names!!!!!!
+FORMS_ComputeConditionMatrixFrob := function(u, h, h_star, scalar_h, g_star_inv_scaled, frob_base, frob_base_inv, frob_base_inv_star, F, n)
+    local coeffs_c, coeffs_f, Ps, i, j, b_end, b_start, cpol, fpol;
+    coeffs_c := (u * h) * frob_base_inv;
+    coeffs_f := (u * frob_base_inv) * scalar_h;
+    j := Size(frob_base[3]);
+    Ps := NullMat(n * j, n, F);
+    # Ps := ZeroMatrix(F, n*j, n);
+    for i in [1..j] do
+        if i = j then
+            b_end := n;
+        else
+            b_end := frob_base[3][i + 1] - 1;
+        fi;
+        Ps{[((i - 1)*n + 1)..(i*n)]}{[1..n]} :=
+            FORMS_EvaluatePolynomialWithFrobenius(coeffs_c{[frob_base[3][i]..b_end]}, g_star_inv_scaled, frob_base, frob_base_inv_star, F, n) * h_star - FORMS_EvaluatePolynomialWithFrobenius(coeffs_f{[frob_base[3][i]..b_end]}, g_star_inv_scaled, frob_base, frob_base_inv_star, F, n);
+    od;
+    return Ps;
+end;
+
+## similar to FORMS_FrobSpin but only spins one block namely that of block_index
+FORMS_FrobSpinAtBlock := function(Image, spin_elem, frob_base_blocks, block_index, n, F)
+    local A, j, i, k, end_pos;
+    j := Size(frob_base_blocks);
+    A := NullMat(n, n, F);
+    # A := ZeroMatrix(F, n, n);
+    if block_index = j then
+        end_pos := n;
+    else
+        end_pos := frob_base_blocks[block_index + 1] - 1;
+    fi;
+    A[frob_base_blocks[block_index]] := Image;
+    for k in [(frob_base_blocks[block_index] + 1)..end_pos] do
+        A[k] := A[k - 1]*spin_elem;
+    od;
+    return A;
+end;
+
+# find symplectic and symmetric matrices in Forms
+# returns bases [[symmetric forms], [symplectic forms]] for char <> 2 and [symmetric forms] for char = 2
+FORMS_FilterBilinearForms := function(Forms, F, n)
+    local transposed_equal_result, form, symmetric_base, symplectic_base, transposed_form, TransposedEqual, symmetric_base_vecs, symplectic_base_vecs, char_2_eqs, sol, out, mat, tmp, i, s;
+
+    if Size(Forms) = 0 then
+        return [];
+    fi;
+    
+    if Characteristic(F) <> 2 then
+        if Size(Forms) = 1 then
+            if FORMS_IsSymplecticMatrix(Forms[1], F) then
+                return [[], [Forms[1]]];
+            fi;
+            if FORMS_IsSymmetricMatrix(Forms[1]) then
+                return [[Forms[1]], []];
+            fi;
+            return [];
+        fi;
+
+
+        # maybe not use these mutable bases
+        symmetric_base := MutableBasis(F, [NullMat(n, n, F)]);
+        symplectic_base := MutableBasis(F, [NullMat(n, n, F)]);
+        # symmetric_base := MutableBasis(F, [ZeroMatrix(F, n, n)]);
+        # symplectic_base := MutableBasis(F, [ZeroMatrix(F, n, n)]);
+        # symmetric_base := MutableBasis(F, [], ZeroVector(F, n));
+        # symplectic_base := MutableBasis(F, [], ZeroVector(F, n));
+        for form in Forms do
+            transposed_form := TransposedMat(form);
+            CloseMutableBasis(symmetric_base, transposed_form + form);
+            CloseMutableBasis(symplectic_base, form - transposed_form);
+        od;
+
+        symmetric_base_vecs := BasisVectors(ImmutableBasis(symmetric_base));
+        symplectic_base_vecs := BasisVectors(ImmutableBasis(symplectic_base));
+
+        if Size(symmetric_base_vecs) + Size(symplectic_base_vecs) <> Size(Forms) then 
+            Error("This should not have happend!! there are supposedly ", Size(symmetric_base_vecs), " linearly independent symmetric forms and ", Size(symplectic_base_vecs), " linearly independent symplectic forms, yet the dimension of the formspace is ", Size(Forms), "\n");
+        fi;
+        return [symmetric_base_vecs, symplectic_base_vecs];
+    fi;
+    # TODO: solve this in some efficient way that does not formulate this by solving sets of linear equations of matrices. Instead it would be better to gradually consider entries of the matrices. Then use some heuristic to determine we are done and just check wether the resulting matrices are infact symmetric. this should be faster because now worst case we are solving a system of linear equations that consists of n^2 equations and Size(Forms) indeterminates.
+    # TODO: maybe do these todos for all characteristics?
+
+    if Size(Forms) = 1 then
+        if FORMS_IsSymplecticMatrix(Forms[1], F) then
+            return [Forms[1]];
+        else 
+            return [];
+        fi;
+    fi;
+    
+    char_2_eqs := [];
+    for form in Forms do
+        Add(char_2_eqs, form - TransposedMat(form));
+    od;
+    out := [];
+    sol := FORMS_SolveMatrixSystem(char_2_eqs, n, n, F);
+    for s in sol do
+        tmp := Forms[1]*s[1];
+        for i in [2..Size(s)] do
+            tmp := tmp + s[i]*Forms[i];
+        od;
+        Add(out, tmp);
+    od;
+    return out;
+end;
+
+# tries to filter the F = GF(q^2) vectorspace generated by <Forms> and return the GF(q) vector space A such that A = <Forms> \cap B where B = {A \in F^{n\times n}, A* = A}
+FORMS_FilterUnitaryForms := function(Forms, F, n, hom)
+    local i, j, ent, q, half, O, l, FF, p, tr_form, Base, baseVecs, gf_base, hgf_base, mat, small_field, field_aut_mat, gf_base_vecs, big_aut_mat, to_smaller_field_matrix, transpose_weird_mat, apply_aut_to_mat, eqs, s, form, form_changed, sol, out, tmp, big_field, changed_forms, apply_coefficient_and_get_smaller_matrix, square_base_entry_rep, multiply_with_scalar, form_changed_mult, form_changed_mult_star, form_changed_star, FORMS_ScalarFormIdentifyCheck;
+    if Size(Forms) = 0 then
+        return [];
+    fi;
+    p := Characteristic(F);
+    q := p^(DegreeOverPrimeField(F)/2);
+    
+    # checks if (Form^T)^hom = c Form for some c in F and returns d such that d * Form = ((d Form)^T)^hom if possible or if not possibe fail
+    FORMS_ScalarFormIdentifyCheck := function(Form, F, n, hom, p, q)
+        local lambda, i, j;
+        lambda := fail;
+        for i in [1..n] do
+            for j in [1..n] do
+                if not IsZero(Form[i, j]) then 
+                    if IsZero(Form[j, i]) then
+                        return fail;
+                    fi;
+                    if lambda <> fail and Form[i, j] * lambda <> hom(Form[j, i]) then
+                        return fail;
+                    fi;
+                    if lambda = fail then
+                        lambda := hom(Form[j, i]) * Inverse(Form[i, j]);
+                    fi;
+                fi;
+            od;
+        od;
+        return RootFFE(F, Inverse(lambda), q - 1);
+    end;
+
+    if Size(Forms) = 1 then
+        #checks if A = A* or A = cA* if A = A* return A, if A = cA* we want to return scalar multiples of A, namely lA for l such that c = l^(1-q) iff c^-1 = l^(q-1)
+        # all the solutions then are lA*GF(q) is this correct?? i am not sure if this are indeed all the possible solutions, but it certanly are solutions.
+        l := FORMS_ScalarFormIdentifyCheck(Forms[1], F, n, hom, p, q);
+        if l = fail then
+            return [];
+        fi;
+        
+        return [Forms[1] * l];
+    fi;
+       
+    ## Lemma: If A = A* is a form preserved modulo scalar c in GF(q^2) i.e. gAg* = cA for group elements g and A is not 0. then c in GF(q).
+
+    ## for char = 2 this can not yield all hermitian matrices as it deletes the diagonal.
+    if p <> 2 then
+        Base := MutableBasis(GF(q), [NullMat(n, n, GF(q))]);
+        # Base := MutableBasis(GF(q), [ZeroMatrix(GF(q), n, n)]);
+        
+        # Base := MutableBasis(GF(q), [], ZeroVector(GF(q), n));
+        gf_base := BasisVectors(Basis(GF(GF(q), 2)))[2];
+        hgf_base := hom(gf_base);
+        for FF in Forms do
+            tr_form := TransposedMat(FF^hom);
+            CloseMutableBasis(Base, FF + tr_form);
+            CloseMutableBasis(Base, gf_base * FF + hgf_base*tr_form);
+        od;
+
+        O := [];
+        baseVecs := BasisVectors(ImmutableBasis(Base));
+
+        return baseVecs;
+    fi;
+
+    # the idea for char 2 is to solve the semiliner system of equations. we take a GF(q) basis of GF(q^2) namely <1, delta> and express n times n matrices with this basis. 
+
+    # TODO: this function is potentially really slow, it would be much better two only take a few equations and constain the problem instead of taking the entire n times n matrix. One such example would be the form space preserved by G := Group(SU(200, 2^2).1) then Size(Forms) = 38420 consisting of 200x200 matrices. For cases such as this it might also be good to just have some function that tries to find one random non-deg form instead of generating random ones
+
+    small_field := GF(q);
+    big_field := GF(q^2);
+    gf_base := Basis(GF(small_field, 2));
+    gf_base_vecs :=  BasisVectors(gf_base);
+    # expresses n times n matrices as 2n times n matrices where 2 times 1 collumn vectors contain the coefficients ascoiciated with the basis gf_base
+    to_smaller_field_matrix := function(basis_of_field, small_field_, mat, n, c)
+        local outmat, i, j;
+        outmat := NullMat(2 * n, c, small_field_);
+        # outmat := ZeroMatrix(small_field_, 2 * n, c);
+        for i in [1..n] do
+            for j in [1..c] do
+                outmat{[(2*i - 1)..(2*i)]}[j] := Coefficients(basis_of_field, mat[i][j]);
+            od;
+        od;
+        return outmat;
+    end;
+    # this matrix applies the field automorphism hom to a element expressed in the basis over the smaller field
+    field_aut_mat := to_smaller_field_matrix(gf_base, small_field, [[hom(gf_base_vecs[1]), hom(gf_base_vecs[2])]], 1, 2); 
+    apply_aut_to_mat := function(mat, n, aut_mat)
+        local outmat, i, j;
+        outmat := 1*mat; # ugly hack to force gap to make a new copy of the matrix.
+        for i in [1..n] do
+            for j in [1..n] do
+                outmat{[(2*i - 1)..(2*i)]}[j] := aut_mat * outmat{[(2* i - 1)..(2*i)]}[j];
+            od;
+        od;
+        return outmat;
+    end;
+
+    # transposes 2n times n matrix (as if it were a n times n matrix)
+    transpose_weird_mat := function(mat, n)
+        local outmat, i, j;
+        outmat := 1*mat; # hack to force gap to copy the matrix.
+        for i in [1..n] do
+            for j in [1..n] do
+                outmat{[(2*i-1)..(2*i)]}[j] := mat{[(2*j-1)..(2*j)]}[i];
+            od;
+        od;
+        return outmat;
+    end;
+
+    square_base_entry_rep := Coefficients(gf_base, gf_base_vecs[2]^2);
+    # multiplies a 2n times n matrix over F_q with a scalar from F_q^2
+    multiply_with_scalar := function(scalar,gf_base, square_base_entry_rep, mat, n)
+        local smat, srep;
+        srep := Coefficients(gf_base, scalar);
+        smat := [[srep[1], srep[2] * square_base_entry_rep[1]], 
+        [srep[2], srep[1] + srep[2] * square_base_entry_rep[2]]];
+        return apply_aut_to_mat(mat, n, smat);
+    end;
+
+    eqs := [];
+    # build linear equations that can now be solved over F_q
+    for form in Forms do
+        form_changed := to_smaller_field_matrix(gf_base, small_field, form, n, n);
+        form_changed_star := apply_aut_to_mat(transpose_weird_mat(form_changed, n), n, field_aut_mat);
+        form_changed_mult := multiply_with_scalar(gf_base_vecs[2], gf_base, square_base_entry_rep, form_changed, n);
+        form_changed_mult_star := multiply_with_scalar(hom(gf_base_vecs[2]), gf_base, square_base_entry_rep, form_changed, n);
+
+        Add(eqs, form_changed - form_changed_star);
+        Add(eqs, form_changed_mult - form_changed_mult_star);
+    od;
+
+    sol := FORMS_SolveMatrixSystem(eqs, 2*n, n, small_field);
+    out := [];
+    for s in sol do
+        tmp := (s[1]*gf_base_vecs[1] + s[2]*gf_base_vecs[2])*Forms[1];
+        for i in [2..Size(Forms)] do
+            tmp := tmp + (s[2*i - 1]*gf_base_vecs[1] + s[2*i]*gf_base_vecs[2])*Forms[i];
+        od;
+        Add(out, tmp);
+    od;
+    return out;
+end;
+
+# Compute the formspace for cyclic matrix group. 
+# TODO: optimizations:
+# this function (sometimes) yields a very large formspace 
+# to better recognize forms in this case it would be good to add a function that does not do this (since we only care about non degenerate classical forms). to find a bilinear/symplectic/hermitian form we can just compute a invertibe matrix S such that gS = Sg^{-*} (with frobenius normal form) and hope that S + S^*, S - S^* are also inevertible. Then we have found symmetric/symplectic non degenrate forms This seems like a good idea? idk
+FORMS_CyclicGroupCase := function(Gen, Gen_adjoint_inv_scaled, Lambda, unitary, hom, frob, frob_inv_star_scaled, frob_inv_star_base_change, frob_inv_base_change, F, n)
+    # maybe recoginize the trivial group here as a special case
+    local p, mat, outspace, i, j, w, OutForms;
+
+    outspace := [];
+    for p in frob[1] do
+        mat := FORMS_EvaluatePolynomialWithFrobenius(CoefficientsOfUnivariatePolynomial(p), Gen_adjoint_inv_scaled * Lambda, frob_inv_star_scaled, frob_inv_star_base_change, F, n);
+        Add(outspace, NullspaceMatDestructive(mat));
+    od;
+    OutForms := [];
+    
+    for i in [1..Size(outspace)] do
+        for w in outspace[i] do
+            Add(OutForms, frob_inv_base_change * FORMS_FrobSpinAtBlock(w, Gen_adjoint_inv_scaled, frob[3], i, n, F)); # this can be used to bound the rank very cheaply, it is smaller than the lenght of the ith frobenius block, furthermore we use the different frobenius blocks, to build a non-deg form we should use a form from each block and add them??
+        od;
+    od;
+    return OutForms;
+end;
+
+
+# builds the forms from image vectors and frobenius normal forms and so on, might needs to check if the forms are actually preserved forms
+FORMS_ReturnFormspace := function(needs_checking, W, g_res, g_star_inv_scaled, Lambdas, unitary, hom, Gens, d, F, g_inv_frob, n)
+    local O, w, A, i;
+    # no kernel, return empty
+    O := [];
+    for w in W do
+        A := g_inv_frob * FORMS_FrobSpin(FORMS_VectorReorganize(w, Size(g_res[5]), F, n), g_star_inv_scaled, g_res[5], n, F);
+        if needs_checking then
+            for i in [1..d] do
+                #todo: store adjoints so they do not have to be reused!
+                if Gens[i] * A * FORMS_CalculateAdjoint(Gens[i], unitary, hom, n, F) <> Lambdas[i] * A then
+                    if Size(W) = 1 then
+                        return []; #one dimensional case no further computation needed
+                    fi;
+                    return false;
+                fi;
+            od;
+            Add(O, A);
+        else
+            Add(O, A);
+        fi;
+    od;
+    return O;
+end;
+
+# Returns formspace preserved by the group <Gens> modulo Lambdas. Unitary says wheter to look for unitary forms or not. hom can be the Field Automorphism of order two. g_res = [g, Lambda_g, ## The elements of ## FrobeniusNormalForm(g)]. Where g is randomly determenied. g_inv_frob = Inverse(FrobeniusNormalForm(g)[2]). g_star_inv_scaled = g^{-*} * Lambda_g, g_star_inv_scaled_frob = FrobeniusNormalForm(g^{-*} * Lambda_g), frob_base_inv_star = Inverse(g_star_inv_scaled_frob[2]). d = Size(Gens), F is base field and n is the matrix dimension. 
+FORMS_FormspaceInternal := function(Gens, Lambdas, unitary, hom, g_res, g_inv_frob, g_star_inv_scaled, g_star_inv_scaled_frob, frob_base_inv_star, d, F, n)
+    local k, i, W, first, j, Cond, Conds, h_star, h, w, O, A, failed_check, nspace, vec, stagnatiton, stagnation_max, old_kernel_size, nspace_size;
+    first := true;
+   
+    ## after computing 5 kernels and the formspace not getting smaller, assume that all forms have been found.
+    old_kernel_size := n + 1;
+    stagnatiton := 0;
+    stagnation_max := 2; # test three matrices
+    for i in [1..d] do
+        h := Gens[i];
+        h_star := FORMS_CalculateAdjoint(h, unitary, hom, n, F);
+        for j in [1..n] do
+            # vec := RandomVector(F, n); # does not seem to be better
+            vec := g_res[4][j];
+            Conds := 
+                FORMS_ComputeConditionMatrixFrob(vec, h, h_star, Lambdas[i], g_star_inv_scaled, g_star_inv_scaled_frob, g_inv_frob, frob_base_inv_star, F, n);
+            
+            if not first then
+                nspace := NullspaceMat(W * Conds);
+                nspace_size := Size(nspace);
+                if nspace_size = 0 then 
+                    return [];
+                fi;
+                W := nspace * W;
+            fi;
+            if first then
+                nspace := NullspaceMat(Conds);
+                nspace_size := Size(nspace);
+                if nspace_size = 0 then 
+                    return [];
+                fi;
+                W := nspace;
+                first := false;
+            fi;
+            if nspace_size = 1 then
+                break;
+            fi;
+            if old_kernel_size = nspace_size then
+                stagnatiton := stagnatiton + 1;
+            else 
+                stagnatiton := 0;
+            fi;
+            old_kernel_size := nspace_size;
+            # try exit early
+            if stagnatiton = stagnation_max then
+                stagnatiton := 0;
+                # maybe not do this computation too often..?
+                # stagnation_max := stagnation_max * 2;
+                O := FORMS_ReturnFormspace(true, W, g_res, g_star_inv_scaled, Lambdas, unitary, hom, Gens, d, F, g_inv_frob, n);
+                if O <> false then
+                    return O;
+                fi;
+            fi;
+        od;
+    od;
+    
+    # TODO: only check if W one dimensional or cyclic group element condition is not met, i think this can be determined from frob normal forms if we even need to do this.
+    O := FORMS_ReturnFormspace(true, W, g_res, g_star_inv_scaled, Lambdas, unitary, hom, Gens, d, F, g_inv_frob, n);
+
+    if O = false then
+        # TODO: cyclic matrix conditition and example Group where this is needed if such a group exists. Can only happen for a group that does not preserve non-deg form but preserves deg form i think.
+        # TODO: one could also do this instead of checking if forms are preserved however this seems slower
+    fi;
+    return O;
+end;
+
+# ideas for scalars : 
+## non-degenerate 
+# - do the things already in the forms package with trace/ analyze minimal and characteristic polynomial DONE! :)
+## degenerate : 
+# - do the trick where a group generator with scalar guranteed to be one gets added and compute the eigenvalues of this one these are then all possible choices of scalars, we can even do this a few times to arrow down the possibilities. 
+    ## structure : 
+        # 1. find a random element in the group of form ghg^{-1}h^{-1} (i.e.) it has scalar one with as few frobenius blocks (best case cyclic) in the group.
+        # 2. we compute the eigenvalues of one part of the condition matrix. TODO: specify
+        # 3. repeat the first two steps until scalars are sufficiently constricted
+        # 4. compute formspaces of these scalars.
+
+        ## If there is a frobenius hom wi need to do the above twice. once with the normal conjugation and once with frobenius conjugation?
+## commutator group? idk if this is a good approach, ideally i want it to be able to handle every case, even if the commutator subgroup generated is trivial. the issue with this is that the resulting forms of the commutator group may not be forms of the bigger group. In this case we might need linear combinations and are still missing scalars
+
+#! @Arguments G, L, unitary
+#! @Returns a basis of $\mathcal{F}_h(G, \Lambda)$
+#! @Description
+#!  <A>G</A> is a finitely generated matrix group over a finite field $K$. 
+#!  L is a list of scalars corresponding to the group generators. I.e. $\Lambda(G_i)$ = <A>L[i]</A> for $i = 1, \dots, $<A>Size(GeneratorsOfGroup(G))</A>. 
+#!  If unitary is true, the function uses $h(x) = x$ for all $x \in K$. 
+#!  If unitary is false, the function tries to use a field automorphism of order two for $h$. If no such field automorphism exists, it will return [].
+InstallMethod(PreservedFormspace,
+    "for matrix group over finite field, with given scalars, and search for unitary forms", [IsMatrixGroup, IsVector and IsFFECollection, IsBool],
+    function(G, Lambdas, unitary) 
+        local F, p_exponent, Gens, n, d, first, hom, g_res, g_inv_frob, g_star_inv_unscaled, g_star_inv_scaled_frob, frob_base_inv_star, frob, frob_inv_star, Gen, Gen_adjoint, Gen_adjoint_inv_scaled, frob_inv_star_scaled, frob_inv_star_base_change;
+        F := DefaultFieldOfMatrixGroup(G);
+        p_exponent := DegreeOverPrimeField(F);
+        if unitary and (p_exponent mod 2 <> 0) then
+            # Print("Field does not admit field automorphism of order two!");
+            return [];
+        fi;
+        # check if field is finite?
+        Gens := GeneratorsOfGroup(G);
+        # F := DefaultFieldOfMatrix(Gens[1]);
+        n := NrRows(Gens[1]);
+        d := Size(Gens);
+        hom := fail;
+
+        
+        if unitary then
+            hom := FrobeniusAutomorphism(F)^(p_exponent/2);
+        fi;
+        
+        if d = 1 then
+            Gen := Gens[1];
+            Gen_adjoint := FORMS_CalculateAdjoint(Gen, unitary, hom, n, F);
+            Gen_adjoint_inv_scaled := Lambda[1]*Inverse(Gen_adjoint);
+            frob := FrobeniusNormalForm(Gen);
+            frob_inv_star_scaled := FrobeniusNormalForm(Gen_adjoint_inv_scaled);
+            frob_inv_star_base_change := Inverse(frob_inv_star_scaled[2]);
+
+            return FORMS_CyclicGroupCase(Gen, Gen_adjoint_inv_scaled, Lambda[1], unitary, hom, frob, frob_inv_star_scaled, frob_inv_star_base_change, Inverse(frob[2]), F, n);
+        fi;
+        #contains  group element, scalar, (factors of minopol), Basis change to Frobenius (v, vg, vg^2, ...), Frobenius block lengths, number of iterations to compute
+        g_res := FORMS_FindCyclicGroupElementAndScalars(Gens, Lambdas, n);
+        g_inv_frob := Inverse(g_res[4]);
+        g_star_inv_unscaled := TransposedMat(Inverse(g_res[1]));
+        if unitary then
+            g_star_inv_unscaled := g_star_inv_unscaled^hom;
+        fi;
+        g_star_inv_scaled_frob := FrobeniusNormalForm(g_star_inv_unscaled * g_res[2]);
+        frob_base_inv_star := Inverse(g_star_inv_scaled_frob[2]);
+
+        return FORMS_FormspaceInternal(Gens, Lambdas, unitary, hom, g_res, g_inv_frob, g_star_inv_unscaled * g_res[2], g_star_inv_scaled_frob, frob_base_inv_star, d, F, n);
+    end
+);
+
+#! @Arguments G,
+#! @Returns a basis of $\mathcal{F}_{id}(G, \Lambda_1)$ and $\mathcal{F}_h(G, \Lambda_1)$
+#! @Description
+#!  <A>G</A> is a finitely generated Matrix group over a finite field $K$. 
+#! Here $\Lambda_1(g) := 1$ for all $g\in G$ and $id(x) := x$ for all $x \in K$. Furthermore $h$ is a field automorphism of order two, admitted by $h$. If no such field automorphism exists, the second returned basis is empty.
+InstallMethod(PreservedFormspace,
+    "for matrix group (finds bilinear/unitary forms modulo 1)",
+    [IsMatrixGroup],
+    function(G)
+        local F, p_exponent, Gens, n, d, first, hom, g_res, g_inv_frob, g_star_inv_unscaled, g_star_inv_scaled_frob, frob_base_inv_star, Lambdas, Out, i, Gen, Gen_adjoint, Gen_adjoint_inv_scaled, frob, frob_inv_star_scaled, frob_inv_star_base_change;
+        Out := [];
+        F := DefaultFieldOfMatrixGroup(G);
+        p_exponent := DegreeOverPrimeField(F);
+        # Todo check if F finite?
+        Gens := GeneratorsOfGroup(G);
+        n := NrRows(Gens[1]);
+        d := Size(Gens);
+        hom := fail;
+
+        if d = 1 then
+            Gen := Gens[1];
+            Gen_adjoint := FORMS_CalculateAdjoint(Gen, false, fail, n, F);
+            Gen_adjoint_inv_scaled := One(F)*Inverse(Gen_adjoint); # scaling happens here!!
+            frob := FrobeniusNormalForm(Gen);
+            frob_inv_star_scaled := FrobeniusNormalForm(Gen_adjoint_inv_scaled);
+            frob_inv_star_base_change := Inverse(frob_inv_star_scaled[2]);
+            g_inv_frob := Inverse(frob[2]);
+
+            Add(Out, FORMS_CyclicGroupCase(Gen, Gen_adjoint_inv_scaled, One(F), false, fail, frob, frob_inv_star_scaled, frob_inv_star_base_change, g_inv_frob, F, n));
+
+            if p_exponent mod 2 = 0 then
+                hom := FrobeniusAutomorphism(F)^(p_exponent/2);
+                Gen_adjoint_inv_scaled := Gen_adjoint_inv_scaled^hom;
+                frob_inv_star_base_change := frob_inv_star_base_change^hom;
+                frob_inv_star_scaled[2] := frob_inv_star_scaled[2]^hom;
+
+                Add(Out, FORMS_CyclicGroupCase(Gen, Gen_adjoint_inv_scaled, One(F), true, hom, frob, frob_inv_star_scaled, frob_inv_star_base_change, g_inv_frob, F, n));
+            else 
+                Add(Out, []);
+            fi;
+            return Out;
+        fi;
+    
+        Lambdas := [];
+        for i in [1..d] do
+            Add(Lambdas, One(F));
+        od;
+        #contains  group element, scalar, (factors of minopol), Basis change to Frobenius (v, vg, vg^2, ...), Frobenius block lengths, number of iterations to compute
+        g_res := FORMS_FindCyclicGroupElementAndScalars(Gens, Lambdas, n);
+        g_inv_frob := Inverse(g_res[4]);
+        g_star_inv_unscaled := TransposedMat(Inverse(g_res[1]));
+        g_star_inv_scaled_frob := FrobeniusNormalForm(g_star_inv_unscaled); # maybe this causes a bug with scalars
+        frob_base_inv_star := Inverse(g_star_inv_scaled_frob[2]);
+
+        Add(Out, FORMS_FormspaceInternal(Gens, Lambdas, false, hom, g_res, g_inv_frob, g_star_inv_unscaled * g_res[2], g_star_inv_scaled_frob, frob_base_inv_star, d, F, n));
+        if p_exponent mod 2 = 0 then
+            # is ^hom actually cheaper than computing the frobenius normal form?? investigate! certainly makes the code ugly... oh well
+            hom := FrobeniusAutomorphism(F)^(p_exponent/2);
+            g_star_inv_unscaled := g_star_inv_unscaled^hom;
+            g_star_inv_scaled_frob[2] := g_star_inv_scaled_frob[2]^hom;
+            frob_base_inv_star := frob_base_inv_star^hom;
+            Add(Out, FORMS_FormspaceInternal(Gens, Lambdas, true, hom, g_res, g_inv_frob, g_star_inv_unscaled * g_res[2], g_star_inv_scaled_frob, frob_base_inv_star, d, F, n));
+        else
+            Add(Out, []);
+        fi;
+        return Out;
+    end
+);
+
+#! @Arguments Forms, Field, unitary
+#! @Returns basis of the spaces of symmetric/symplectic matrices or a basis of the hermitian matrices contained in Forms
+#! @Description
+#! In the case where unitary is false, this will return a list that contains to lists of matrices. The first is a basis of the symmetric matrices, the second is a basis of the symplectic matrices. Be carefull: If the characteristic of the field is 2, then symplectic matrices and symmetric matrices are the same. Hence only one basis will be returned. If unitary is false this will return a basis of the hermitian matrices.
+InstallMethod(FilterFormspace, "for list of matrices, F finite field, bool hermitian", [IsList, IsFinite and IsField, IsBool], function(Forms, F, unitary)
+    local n, hom, p_exponent;
+    if Size(Forms) = 0 then
+        return [];
+    fi;
+    n := NrRows(Forms[1]);
+
+    if Size(Forms) = n^2 then
+        # TODO: special case where we should just return a precomputed basis... 
+    fi;
+
+    if not unitary then
+        return FORMS_FilterBilinearForms(Forms, F, n);
+    else
+        p_exponent := DegreeOverPrimeField(F);
+        if p_exponent mod 2 <> 0 then 
+            Error("The given Field ", F, " must admit a field automorphism of order two for unitary = true\n");
+            return [];
+        fi;
+        hom := FrobeniusAutomorphism(F)^(p_exponent/2);
+        return FORMS_FilterUnitaryForms(Forms, F, n, hom);
+    fi;
+end);
+
+
+# computes non degenerate forms and scalars, TODO: random_non_deg should specify whether the function should use random methods to make sure that the returned forms are non deg, this is a todo and not necessary for irreducible groups.
+FORMS_PreservedNonDegFormsWithScalarsOp := function(G, random_non_deg, sesquilinear)
+    local F, Gens, n, p_exponent, hom, newGens, forms_bil, forms_herm, filtered, space, forms_out, filt_bil, filt_herm, Form, quad, x, biglist;
+    F := DefaultFieldOfMatrixGroup(G);
+    Gens := GeneratorsOfGroup(G);
+    n := NrRows(Gens[1]);
+    p_exponent := DegreeOverPrimeField(F);
+    forms_bil := [];
+    forms_herm := [];
+
+    hom := FrobeniusAutomorphism(F)^0; # better way to obtain identity mapping?
+    newGens := ClassicalForms_GeneratorsWithBetterScalarsSesquilinear(G, hom);
+    filt_bil := [[], []];
+    if newGens <> false then
+        # this should be improved, many computations from PreservedFormspace could probably be recycled for example the cyclic group element, and so on... Does not seem to be such a big issue, this is only inefficient for groups that actually have multiple possible scalars which seems to be very rare. Infact TODO: test case with (irreducible) group which preserves forms modulo different sets of scalars. Does such a group even exist?
+        biglist := Cartesian(newGens[2]);
+        # Print(biglist);
+        for x in biglist do
+            forms_bil := PreservedFormspace(Group(newGens[1]), x, false);
+            filtered := FilterFormspace(forms_bil, F, false);
+            if forms_bil = [] then
+                # do nothing .. :)
+            elif Characteristic(F) = 2 then
+                for Form in filtered do
+                    if not sesquilinear then
+                        quad := ClassicalForms_QuadraticForm2(F, Form, newGens[1], x);
+                        if quad = false then
+                            Add(filt_bil[2], Form);
+                        else
+                            Add(filt_bil[1], quad);
+                        fi;
+                    else
+                        Add(filt_bil[2], Form);
+                    fi;
+                od;
+            else 
+                filt_bil[2] := Concatenation(filt_bil[2], filtered[2]);
+                for Form in filtered[1] do
+                    if not sesquilinear then
+                        quad := ClassicalForms_QuadraticForm(F, Form);
+                        Add(filt_bil[1], quad);
+                    else
+                        Add(filt_bil[2], Form);
+                    fi;
+                od;
+            fi;
+        od;
+    fi;
+   
+    filt_herm := [];
+    if p_exponent mod 2 = 0 then
+        hom := FrobeniusAutomorphism(F)^(p_exponent/2);
+        newGens := ClassicalForms_GeneratorsWithBetterScalarsSesquilinear(G, hom);
+        if newGens <> false then
+            biglist := Cartesian(newGens[2]);
+            for x in biglist do
+                Add(forms_herm, PreservedFormspace(Group(newGens[1]), x, true));
+            od;
+        fi;
+        for space in forms_herm do
+            filt_herm := Concatenation(filt_herm, FilterFormspace(space, F, true));
+        od;
+    fi;
+    
+    ## process the forms in forms_bil, forms_herm
+    forms_out := [];
+    for Form in filt_bil[1] do
+        Add(forms_out, QuadraticFormByMatrix(Form, F));             
+    od;
+    for Form in filt_bil[2] do
+        Add(forms_out, BilinearFormByMatrix(Form, F));
+    od;
+    for Form in filt_herm do
+        Add(forms_out, HermitianFormByMatrix(Form, F));
+    od;
+    return forms_out;
+end;
+
+
+#! @Arguments G
+#! @Returns List of preserved Forms by G
+#! @Description
+#! <A>G</A> is a finitely generated matrix group over a finite field $K$. 
+#! This function computes scalars $\lambda$ and forms $\Phi$ such that <A>G</A> preserves $\Phi$ modulo $\lambda$. The function does only find scalars that belong to non-degenerate forms but may return degenerate forms if the group action is reducible. 
+InstallMethod(PreservedFormsWithScalars, "for matrix groups", [IsMatrixGroup], 
+function(G)
+    return FORMS_PreservedNonDegFormsWithScalarsOp(G, true, false);
+end);
+
+#! @Arguments G
+#! @Returns List of preserved sesquilinear Forms by G
+#! @Description 
+#! Very similar to <A>PreservedFormsWithScalars</A> where <A>G</A> is a finitely generated matrix group over a finite field $K$. 
+#! This function computes scalars $\lambda$ and forms $\Phi$ such that <A>G</A> preserves $\Phi$ modulo $\lambda$. The function does only find scalars that belong to non-degenerate forms but may return degenerate forms if the group action is reducible.
+InstallMethod(PreservedSesquilinearFormsWithScalars, "for matrix groups", [IsMatrixGroup], 
+function(G)
+    return FORMS_PreservedNonDegFormsWithScalarsOp(G, true, true);
+end);
diff --git a/lib/recognition_new.gi b/lib/recognition_new.gi
index 5fb035e..c8bb26a 100644
--- a/lib/recognition_new.gi
+++ b/lib/recognition_new.gi
@@ -92,15 +92,6 @@ InstallGlobalFunction( ClassicalForms_GeneratorsWithBetterScalarsSesquilinear,
   function( grp, frob )
     local tries, gens, field, m1, a1, new, i, scalars, root, improvegenerator, res, newgens, champion, len, qq, q;
 
-    # the aim of this function is to replace the matrix m1 by a
-    # matrix that has as few solutions to the scalar equation
-    # lambda^a1[1] = a1[2] as possible. It checks first if a1[1] = 1,
-    # since then lambda is determined. Next we check if a1[2] is a
-    # square. And then we replace m1 by a matrix that leaves the
-    # bilinear form invariant modula the scalar 1. If none of these
-    # are possible, we try to replace m1 by a matrix that has fewer
-    # solutions to the scalar equation.
-
     #field := FieldOfMatrixGroup(grp); #this causes a problem if one has a maximal subgroup of e.g. SU(3,7^2). There are examples of which de FieldOf is GF(7), while DefaultFieldOF is GF(7^2). Then this causes problems.
     field := DefaultFieldOfMatrixGroup(grp);
     qq := Size(field);
@@ -111,10 +102,9 @@ InstallGlobalFunction( ClassicalForms_GeneratorsWithBetterScalarsSesquilinear,
     fi;
 
     # the next function returns a matrix with a list of possible scalars for this matrix.
-    # if it is possible to change the matrix to multiple that preserves up to scalar one, this is achieved.
 
-    improvegenerator := function(m1,i,count,len)
-        local a1, s, j, k, scalars, root; #I think root should be declare locally here!
+    improvegenerator := function(m1)
+        local a1, scalars;
         
         a1 := ClassicalForms_PossibleScalarsSesquilinear(field,m1,frob);
         #Recall that the scalars satisfy the equation lambda^a[1] = a[2]
@@ -123,55 +113,38 @@ InstallGlobalFunction( ClassicalForms_GeneratorsWithBetterScalarsSesquilinear,
         fi;
                 
         if IsList(a1) then
-            root := NthRoot(field,a1[2],a1[1]);
             if a1[1] = 1 then # the matrix m1 has scalar a1[2]
-                #if a1[2] has a square root, we can replace m1 with m1*sqrt{a1[2]};
-                if LogFFE(a1[2],PrimitiveRoot(field)) mod 2 = 0 then
-                    return [m1/NthRoot(field,a1[2],2),[One(field)]];
-                fi;
-                return [m1,[a1[2]]]; #originally, the three lines above this return were not there. Those three lines make sure scalar becomes 1 if possible (basicaly if there is a sqrt).
-            elif LogFFE(root,PrimitiveRoot(field)) mod (q+1) = 0 then
-                return [m1/NthRoot(field,root,q+1),[One(field)]]; #either frob = id, then q+1 = 2, or frob is not trivial, then we take q+1-st root. In both cases, modify m1 to a matrix that has scalar one.
+                return [a1[2]];
             else
                 scalars := AsList(Group(NthRoot(field,a1[2],a1[1]))); # add all possible scalars for m1
-                if count = 0 then
-                    return champion; # add all possible scalars for m1
-                elif Length(scalars) < len then
-                    champion := [m1,scalars];
-                    len := Length(scalars);
-                fi;
-                k := Random(Difference([1..Length(gens)],[i])); #m1 could not be improved, so try to change m1 by multiplying with another generator.
-                m1 := m1*gens[k];
-                return improvegenerator(m1,i,count-1,len);
+                return scalars;
             fi;
         fi;
     end;
 
     # start with 2 random elements,  at most 10 tries
-    tries := 0;
-    gens  := ShallowCopy(GeneratorsOfGroup(grp));
-    if Length(gens) = 1 then
-        Add(gens,PseudoRandom(grp));
-    fi;
-    #We will randomize the generating set in the hope that we obtain generators preserving fewer
-    #scalars.
+    
+    # i have commented this out because my routine for computing preserved forms has extra logic for cyclic groups. However maybe i should investigate this!
+    # tries := 0;
+    # gens  := ShallowCopy(GeneratorsOfGroup(grp));
+    # if Length(gens) = 1 then
+    #     Add(gens,PseudoRandom(grp));
+    # fi;
+    #We will randomize the generating set in the hope that we obtain generators preserving fewer scalars.
      
     scalars := [];
  
-    newgens := ShallowCopy(gens);
+    gens := GeneratorsOfGroup(grp);
         
     for i in [1..Length(gens)] do
-        champion := [gens[i],AsList(Group(PrimitiveRoot(field)))];
-        len := Length(champion[2]);
-        res := improvegenerator(gens[i],i,10,len);
+        res := improvegenerator(gens[i]);
         if res = false then
             return false; #The group does not preserve a bilinear form modulo scalars
         fi;
-        newgens[i] := res[1];
-        scalars[i] := res[2];
+        scalars[i] := res;
     od;
 
-    return [newgens,scalars];
+    return [gens,scalars];
 
   end );
 
@@ -355,7 +328,7 @@ InstallMethod( PreservedFormsOp, [ IsMatrixGroup ],
     #Is this really necessary? Can we not simply delete it?
     if not MTX.IsIrreducible(module)  then
         #Error("Currently the use of MeatAxe requires the module to be absolutely irreducible");
-        Info( InfoForms, 1,  "group is not irreducible and therefore it does not preserve non-degenerate forms\n" );
+        Info( InfoForms, 1,  "group is not irreducible and therefore it does not preserve non-degenerate forms\n" ); #this is wrong. Reducible groups can preserve non-degenerate forms. For example the trivial group G = <1> is reducible and preserves any form.
         return [];
     fi;
 
diff --git a/makedoc.g b/makedoc.g
index 0cacfde..712dc6b 100644
--- a/makedoc.g
+++ b/makedoc.g
@@ -13,7 +13,8 @@ if fail = LoadPackage("AutoDoc", ">= 2019.04.10") then
 fi;
 AutoDoc(rec(
     scaffold := rec( MainPage := false ),
-    gapdoc := rec( main := "forms.xml" ),
+    gapdoc := rec( main := "forms.xml",),
+    autodoc := true,
 ));
 QUIT;
 
diff --git a/read.g b/read.g
index abe27b6..703598f 100644
--- a/read.g
+++ b/read.g
@@ -16,3 +16,4 @@ ReadPackage("forms", "lib/forms.gi");
 ReadPackage("forms", "lib/recognition.gi");
 ReadPackage("forms", "lib/classic.gi");
 ReadPackage("forms", "lib/recognition_new.gi");
+ReadPackage("forms", "lib/formspace.gi");
diff --git a/tst/benchmark_forms.g b/tst/benchmark_forms.g
new file mode 100644
index 0000000..23d2205
--- /dev/null
+++ b/tst/benchmark_forms.g
@@ -0,0 +1,60 @@
+GOMinus := function(nn, qq)
+    return GeneralOrthogonalGroup(-1, nn, qq);
+end;
+
+RandomMatrixInvertible := function(n, q)
+    local F, M, i, j;
+    F := GF(q);
+    while true do
+        M := NullMat(n, n, F);
+        for i in [1..n] do
+            for j in [1..n] do
+                M[i][j] := PseudoRandom(F);
+            od;
+        od;
+        ConvertToMatrixRep(M, F);
+        if RankMat(M) = n then
+            return M;
+        fi;
+    od;
+end;
+
+FileBenchmark := function()
+    local bench_time, number_to_average, Ns, Qs, Gs, Gs_names, n, q, k, g, G, mine_total, theirs_total, i, H_conj, H, start_mine, time_after, time_start, time_after_after, average_mine, average_theirs, waste1, waste2, file_name;
+    bench_time := Runtime();
+    number_to_average := 10;
+    Ns := [100, 300, 500];
+    Qs := [3, 5, 3^2];
+    file_name := "benchmark.csv";
+    Gs := [Sp, SU, GOMinus];
+    Gs_names := ["Sp", "SU", "GOMinus"];
+    PrintTo(file_name, "GroupName;n;q;average_time_ours;average_time_theirs");
+    AppendTo(file_name, "\n");
+    for n in Ns do
+        for q in Qs do
+            for k in [1..Size(Gs)] do
+                g := Gs[k];
+                G := g(n, q);
+                mine_total := 0;
+                theirs_total := 0;
+                for i in [1..number_to_average] do
+                    H_conj := RandomMatrixInvertible(n, q);
+                    H := G^H_conj;
+                    start_mine := Runtime();
+                    waste1 := PreservedFormsWithScalars(H);
+                    time_after := Runtime() - start_mine;
+                    time_start := Runtime();
+                    waste2 := PreservedForms(H);
+                    time_after_after := Runtime() - time_start;
+                    mine_total := mine_total + time_after;
+                    theirs_total := theirs_total + time_after_after;
+                od;
+                # '*1000' to convert to seconds
+                average_mine := Float(mine_total/(number_to_average*1000));
+                average_theirs := Float(theirs_total/(number_to_average*1000));
+                AppendTo(file_name, Gs_names[k],";",n,";",q,";",average_mine,";",average_theirs,"\n");
+            od;
+        od;
+    od;
+    Print("Benchmark took ", Float((Runtime() - bench_time)/1000), " seconds.");
+end;
diff --git a/tst/formspace/custom_test_functions.g b/tst/formspace/custom_test_functions.g
new file mode 100644
index 0000000..b6fa3d2
--- /dev/null
+++ b/tst/formspace/custom_test_functions.g
@@ -0,0 +1,190 @@
+TestPolyEval := function(benchmark)
+    local iters, n, F, mat, coeffs, f, frob, eval, i, time_start, time_average_frob, time_average_normal, normal_eval, x;
+    iters := 50;
+    time_average_frob := 0;
+    time_average_normal := 0;
+    for i in [1..iters] do
+        n := Random([1..200]);
+        F := GF(Random([2, 3, 5])^(Random(1, 3)));
+        mat := RandomMat(n, n, F);
+        coeffs := List([1..Random(1, 300)], x -> Random(F));
+        f := UnivariatePolynomial(F, coeffs);
+        time_start := NanosecondsSinceEpoch();
+        frob := FrobeniusNormalForm(mat);
+        eval := FORMS_EvaluatePolynomialWithFrobenius(coeffs, Matrix(mat), frob, Inverse(frob[2]), F, n);
+        time_average_frob := time_average_frob + (NanosecondsSinceEpoch() - time_start);
+
+        time_start := NanosecondsSinceEpoch();
+        normal_eval := f(mat);
+        time_average_normal := time_average_normal + (NanosecondsSinceEpoch() - time_start);
+        if eval <> normal_eval then
+            Error("Polynomial: ", f, " Matrix: ", mat, " Frobenius Normal Form: ", frob, "\n");
+        fi;
+    od;
+    time_average_frob := Float(time_average_frob) / Float((iters * 1000000)); # convert to ms
+    time_average_normal := Float(time_average_normal) / Float((iters  * 1000000));
+    if benchmark then
+        Print("Frob: ", time_average_frob, " Normal: ", time_average_normal, "\n");
+    fi;
+    return "Ok";
+end;
+
+RandomMatrix := function(n, F)
+    local i, M, j;
+    M := ZeroMatrix(F, n, n);
+    for i in [1..n] do
+        for j in [1..n] do
+            M[i][j] := PseudoRandom(F);
+        od;
+    od;
+    return M;
+end;
+
+
+# only checks whether the returned matrices infact unitary, does not check if all matrices have been found!
+TestChar2Filter := function()
+    local powers, num_tries, matnr_range, ns, i, n, power, F, hom, nr_mats, mats, k, out, o;
+    powers := [2, 4];
+    num_tries := 10;
+    matnr_range := [1..10];
+    ns := [10, 20];
+   
+    for i in [1..num_tries] do
+        for n in ns do
+            for power in powers do
+                F := GF(2^power);
+                hom := FrobeniusAutomorphism(F)^(power/2);
+                nr_mats := PseudoRandom(matnr_range);
+                mats := [];
+                for k in [1..nr_mats] do
+                    Add(mats, RandomMatrix(n, F));
+                od;
+                out := FORMS_FilterUnitaryForms(mats, F, n, hom);
+                for o in out do
+                    if o <> TransposedMat(o^hom) then
+                        Error("Error in computation with mats ", mats);
+                    fi;
+                od;
+            od;
+        od;
+    od;
+    return "Ok";
+end;
+
+TestMatricesAreForms := function(G, Lambdas, unitary, forms)
+    local p_exponent, hom, F, f, g, i, Gens, n, j;
+    F := DefaultFieldOfMatrixGroup(G);
+    Gens := GeneratorsOfGroup(G);
+    n := DimensionsMat(Gens[1])[1];
+    p_exponent := DegreeOverPrimeField(F);
+    hom := fail;
+    if unitary then
+        if p_exponent mod 2 <> 0 then
+            if Size(forms) = 0 then
+                return "Ok";
+            else 
+                Error("Claims to have found unitary form although the field does not admit a field automorphism of order two!");
+            fi;
+        fi;
+        hom := FrobeniusAutomorphism(F)^(p_exponent/2);
+    fi;
+
+    for i in [1..Size(forms)] do
+        f := forms[i];
+        for j in [1..Size(Gens)] do
+            g := Gens[j];
+            if g * f * FORMS_CalculateAdjoint(g, unitary, hom, n, F) <> Lambdas[j] * f then
+                Error("Computed non formspace element ", f, "group " , G, " unitary ", unitary, " Lambdas ", Lambdas);
+            fi;
+        od;
+    od;
+    return "Ok";
+end;
+
+TestMatricesAreForms2 := function(G, forms)
+    local F, n, lambdas, Gens, i;
+    F := DefaultFieldOfMatrixGroup(G);
+    Gens := GeneratorsOfGroup(G);
+    n := DimensionsMat(Gens[1])[1];
+    lambdas := [];
+    for i in [1..Size(Gens)] do
+        Add(lambdas, One(F));
+    od;
+    return [TestMatricesAreForms(G, lambdas, false, forms[1]), TestMatricesAreForms(G, lambdas, true, forms[2])];
+end;
+
+# computes F by solving big system of linear equations.
+# this is a good idea for testing to see that the presered formspace function actually computes the entire formspace not just a subspace. (this is the function that provides the formsace solutions used in other tests)
+# however maybe we should test the test? this seems a little stupid
+TestComputeFormspaceBruteForce := function(G, Lambdas, unitary)
+    local Gens, i, j, F, n, base, b, g, eqs, p_exponent, ComputeMatrixVectorSpaceBasis, VectorToMatrix, MatrixToVector, sol, out, s, hom, v;
+    Gens := GeneratorsOfGroup(G);
+    n := DimensionsMat(Gens[1])[1];
+    F := DefaultFieldOfMatrixGroup(G);
+
+    ComputeMatrixVectorSpaceBasis := function(F, n)
+        local O, i, j, m;
+        O := [];
+        for i in [1..n] do
+            for j in [1..n] do
+                m := NullMat(n, n, F);
+                m[i][j] := One(F);
+                Add(O, m);
+            od;
+        od;
+        return O;
+    end;
+
+    VectorToMatrix := function(vec, F, n)
+        local i, j, m;
+        m := NullMat(n, n, F);
+        for i in [1..n] do
+            for j in [1..n] do
+                m[i][j] := vec[(i-1)*n + j];
+            od;
+        od;
+        return m;
+    end;
+
+    MatrixToVector := function(mat, F, n)
+        local vec, i, j;
+        vec := ZeroVector(F, n^2);
+        for i in [1..n] do
+            for j in [1..n] do
+                vec[(i-1)*n + j] := mat[i][j];
+            od;
+        od;
+        return vec;
+    end;
+
+    p_exponent := DegreeOverPrimeField(F);
+    hom := fail;
+    if unitary then
+        if p_exponent mod 2 <> 0 then
+            return [];
+        fi;
+        hom := FrobeniusAutomorphism(F)^(p_exponent/2);
+    fi;
+    base := ComputeMatrixVectorSpaceBasis(F, n);
+    eqs := [];
+   
+    for j in [1..Size(base)] do  
+        for i in [1..Size(Gens)] do
+            b := base[j];  
+            v :=  MatrixToVector(Gens[i] * b * FORMS_CalculateAdjoint(Gens[i], unitary, hom, n, F) - Lambdas[i] * b, F, n);
+            if i = 1 then
+                Add(eqs, v);
+            else
+                eqs[j] := Concatenation(List(eqs[j]), List(v));
+            fi;
+        od;
+    od;
+    # Print(aaa);
+    sol := NullspaceMat(eqs);
+    out := [];
+    for s in sol do
+        Add(out, VectorToMatrix(s, F, n));
+    od;
+    # Print(aaa);
+    return out;
+end;
diff --git a/tst/formspace/filter_tests/filter_bilinear_forms.tst b/tst/formspace/filter_tests/filter_bilinear_forms.tst
new file mode 100644
index 0000000..7b9b212
--- /dev/null
+++ b/tst/formspace/filter_tests/filter_bilinear_forms.tst
@@ -0,0 +1,15 @@
+gap> START_TEST("Formspace: Filter Bilinear Forms"); # todo test char 2
+gap> F := GF(5^2);;
+gap> Forms := [[[0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0], [1, 0, 0, 0]]*One(GF(5)), [[0, 0, 0, 1], [0, 0, 1, 0], [0, 0, 0, 0], [0, 0, 0, 0]]*One(GF(5))];; # This should give a symmetric form and a syplectic form
+gap> L := FORMS_FilterBilinearForms(Forms, F, 4);;
+gap> Size(L);
+2
+gap> Size(L[1]);
+1
+gap> Size(L[2]);
+1
+gap> FORMS_IsSymmetricMatrix(L[1][1]);
+true
+gap> FORMS_IsSymplecticMatrix(L[2][1], GF(5^2));
+true
+gap> STOP_TEST("Formspace: Filter Bilinear Forms");
diff --git a/tst/formspace/filter_tests/unitary_filter_char_2.tst b/tst/formspace/filter_tests/unitary_filter_char_2.tst
new file mode 100644
index 0000000..5d87a8b
--- /dev/null
+++ b/tst/formspace/filter_tests/unitary_filter_char_2.tst
@@ -0,0 +1,4 @@
+gap> START_TEST("Formspace: unitary_filter_char_2");
+gap> TestChar2Filter();
+"Ok"
+gap> STOP_TEST("Formspace: unitary_filter_char_2");
diff --git a/tst/formspace/generate_formspace_tests.g b/tst/formspace/generate_formspace_tests.g
new file mode 100644
index 0000000..b63bccc
--- /dev/null
+++ b/tst/formspace/generate_formspace_tests.g
@@ -0,0 +1,62 @@
+# you do not need to run this. this creates .tst files that should already exist and be correct!
+
+ReadPackage("forms", "tst/formspace/custom_test_functions.g");
+ReadPackage("forms", "tst/interesting_groups.g");
+
+## TODO: Tests with scalars that are not equal to one are desperately missing!!!
+
+WriteTestFilePreservedFormspaceTest := function(path, name, G, n, F,f_space_expected_normal_d, f_space_expected_unitary_d)
+    local full_name, start_test, end_test, stream, file_name, dir, full_path;
+    full_name := StringFormatted("test_{}.tst", name);
+    full_name := ReplacedString(full_name, ",", "_");
+    full_path := StringFormatted("{}/{}", path, full_name);
+
+    start_test := StringFormatted("gap> START_TEST(\"Formspace: Preserved Formspace {}\");\n", name);
+    end_test := StringFormatted("gap> STOP_TEST(\"Formspace: Preserved Formspace {}\");\n", name);
+    dir := DirectoriesPackageLibrary("forms", path);
+    # create file with this
+    PrintTo(full_path, "");
+
+    file_name := Filename(dir, full_name);
+
+    #write to file with this
+    stream := OutputTextFile(file_name, true);
+    WriteAll(stream, start_test);
+    WriteAll(stream, StringFormatted("gap> G := {};; # Some groups are defined in interesting_groups.g\n", G[1]));
+    WriteAll(stream, StringFormatted("gap> R := PseudoRandom(GL({}, {}));;\n", n, F));
+    WriteAll(stream, "gap> G := G^R;;\n");
+    WriteAll(stream, StringFormatted("gap> f_space_expected_normal_d := {};; # the dimensions of the expected formspaces \n", f_space_expected_normal_d));
+    WriteAll(stream, StringFormatted("gap> f_space_expected_unitary_d := {};;\n", f_space_expected_unitary_d));
+    WriteAll(stream, "gap> L:=PreservedFormspace(G);;\n");
+    WriteAll(stream, "gap> TestMatricesAreForms2(G, L); # found in custom_test_functions.g, Tests whether the matrices returned are forms preserved by G or not\n");
+    WriteAll(stream, "[ \"Ok\", \"Ok\" ]\n");
+    WriteAll(stream, "gap> Size(L[1])=f_space_expected_normal_d;\n");
+    WriteAll(stream, "true\n");
+    WriteAll(stream, "gap> Size(L[2])=f_space_expected_unitary_d;\n");
+    WriteAll(stream, "true\n");
+    WriteAll(stream, end_test);
+    CloseStream(stream);
+end;
+
+GenerateTestsForPreservedFormspace := function()
+    local Groups, R, G, GG, n, F, Gens, path, lambdas, i, conjugated, f_space_expected_normal_d, f_space_expected_unitary_d;
+    path := "tst/formspace/preserved_formspace";
+    Groups := [["GO(5,3)", GO(5, 3)], ["SU(4,5)", SU(4, 5)], ["Sp(4,5)", Sp(4, 5)], ["Gtriv", Gtriv], ["G1", G1], ["G2", G2], ["G3", G3], ["G4", G4], ["G5", G5], ["G6", G6], ["G7", G7], ["G8", G8], ["GP22", GP22], ["Group(SP(4,5).1)", Group(SP(4,5).1)]];
+    for GG in Groups do
+        G := GG[2];
+        Gens := GeneratorsOfGroup(G);
+        n := DimensionsMat(Gens[1])[1];
+        F := DefaultFieldOfMatrixGroup(G);
+        R := PseudoRandom(GL(n, F));
+        lambdas := [];
+        for i in [1..n] do
+            Add(lambdas, One(F));
+        od;
+        conjugated := GG[2]^R;
+        f_space_expected_normal_d := Size(TestComputeFormspaceBruteForce(conjugated, lambdas, false));
+        f_space_expected_unitary_d := Size(TestComputeFormspaceBruteForce(conjugated, lambdas, true));
+        WriteTestFilePreservedFormspaceTest(path, GG[1], GG, n, F, f_space_expected_normal_d, f_space_expected_unitary_d);
+    od;
+end;
+
+
diff --git a/tst/formspace/preserved_formspace/test_G1.tst b/tst/formspace/preserved_formspace/test_G1.tst
new file mode 100644
index 0000000..c847240
--- /dev/null
+++ b/tst/formspace/preserved_formspace/test_G1.tst
@@ -0,0 +1,14 @@
+gap> START_TEST("Formspace: Preserved Formspace G1");
+gap> G := G1;; # Some groups are defined in interesting_groups.g
+gap> R := PseudoRandom(GL(10, GF(5^4)));;
+gap> G := G^R;;
+gap> f_space_expected_normal_d := 2;; # the dimensions of the expected formspaces 
+gap> f_space_expected_unitary_d := 0;;
+gap> L:=PreservedFormspace(G);;
+gap> TestMatricesAreForms2(G, L); # found in custom_test_functions.g, Tests whether the matrices returned are forms preserved by G or not
+[ "Ok", "Ok" ]
+gap> Size(L[1])=f_space_expected_normal_d;
+true
+gap> Size(L[2])=f_space_expected_unitary_d;
+true
+gap> STOP_TEST("Formspace: Preserved Formspace G1");
diff --git a/tst/formspace/preserved_formspace/test_G2.tst b/tst/formspace/preserved_formspace/test_G2.tst
new file mode 100644
index 0000000..e985e5f
--- /dev/null
+++ b/tst/formspace/preserved_formspace/test_G2.tst
@@ -0,0 +1,14 @@
+gap> START_TEST("Formspace: Preserved Formspace G2");
+gap> G := G2;; # Some groups are defined in interesting_groups.g
+gap> R := PseudoRandom(GL(4, GF(5^4)));;
+gap> G := G^R;;
+gap> f_space_expected_normal_d := 0;; # the dimensions of the expected formspaces 
+gap> f_space_expected_unitary_d := 2;;
+gap> L:=PreservedFormspace(G);;
+gap> TestMatricesAreForms2(G, L); # found in custom_test_functions.g, Tests whether the matrices returned are forms preserved by G or not
+[ "Ok", "Ok" ]
+gap> Size(L[1])=f_space_expected_normal_d;
+true
+gap> Size(L[2])=f_space_expected_unitary_d;
+true
+gap> STOP_TEST("Formspace: Preserved Formspace G2");
diff --git a/tst/formspace/preserved_formspace/test_G3.tst b/tst/formspace/preserved_formspace/test_G3.tst
new file mode 100644
index 0000000..ea40d0d
--- /dev/null
+++ b/tst/formspace/preserved_formspace/test_G3.tst
@@ -0,0 +1,14 @@
+gap> START_TEST("Formspace: Preserved Formspace G3");
+gap> G := G3;; # Some groups are defined in interesting_groups.g
+gap> R := PseudoRandom(GL(7, GF(5^2)));;
+gap> G := G^R;;
+gap> f_space_expected_normal_d := 2;; # the dimensions of the expected formspaces 
+gap> f_space_expected_unitary_d := 1;;
+gap> L:=PreservedFormspace(G);;
+gap> TestMatricesAreForms2(G, L); # found in custom_test_functions.g, Tests whether the matrices returned are forms preserved by G or not
+[ "Ok", "Ok" ]
+gap> Size(L[1])=f_space_expected_normal_d;
+true
+gap> Size(L[2])=f_space_expected_unitary_d;
+true
+gap> STOP_TEST("Formspace: Preserved Formspace G3");
diff --git a/tst/formspace/preserved_formspace/test_G4.tst b/tst/formspace/preserved_formspace/test_G4.tst
new file mode 100644
index 0000000..2dcfc74
--- /dev/null
+++ b/tst/formspace/preserved_formspace/test_G4.tst
@@ -0,0 +1,14 @@
+gap> START_TEST("Formspace: Preserved Formspace G4");
+gap> G := G4;; # Some groups are defined in interesting_groups.g
+gap> R := PseudoRandom(GL(14, GF(7)));;
+gap> G := G^R;;
+gap> f_space_expected_normal_d := 2;; # the dimensions of the expected formspaces 
+gap> f_space_expected_unitary_d := 0;;
+gap> L:=PreservedFormspace(G);;
+gap> TestMatricesAreForms2(G, L); # found in custom_test_functions.g, Tests whether the matrices returned are forms preserved by G or not
+[ "Ok", "Ok" ]
+gap> Size(L[1])=f_space_expected_normal_d;
+true
+gap> Size(L[2])=f_space_expected_unitary_d;
+true
+gap> STOP_TEST("Formspace: Preserved Formspace G4");
diff --git a/tst/formspace/preserved_formspace/test_G5.tst b/tst/formspace/preserved_formspace/test_G5.tst
new file mode 100644
index 0000000..c5dff94
--- /dev/null
+++ b/tst/formspace/preserved_formspace/test_G5.tst
@@ -0,0 +1,14 @@
+gap> START_TEST("Formspace: Preserved Formspace G5");
+gap> G := G5;; # Some groups are defined in interesting_groups.g
+gap> R := PseudoRandom(GL(14, GF(11^2)));;
+gap> G := G^R;;
+gap> f_space_expected_normal_d := 2;; # the dimensions of the expected formspaces 
+gap> f_space_expected_unitary_d := 0;;
+gap> L:=PreservedFormspace(G);;
+gap> TestMatricesAreForms2(G, L); # found in custom_test_functions.g, Tests whether the matrices returned are forms preserved by G or not
+[ "Ok", "Ok" ]
+gap> Size(L[1])=f_space_expected_normal_d;
+true
+gap> Size(L[2])=f_space_expected_unitary_d;
+true
+gap> STOP_TEST("Formspace: Preserved Formspace G5");
diff --git a/tst/formspace/preserved_formspace/test_G6.tst b/tst/formspace/preserved_formspace/test_G6.tst
new file mode 100644
index 0000000..cd5bb05
--- /dev/null
+++ b/tst/formspace/preserved_formspace/test_G6.tst
@@ -0,0 +1,14 @@
+gap> START_TEST("Formspace: Preserved Formspace G6");
+gap> G := G6;; # Some groups are defined in interesting_groups.g
+gap> R := PseudoRandom(GL(10, GF(2^8)));;
+gap> G := G^R;;
+gap> f_space_expected_normal_d := 0;; # the dimensions of the expected formspaces 
+gap> f_space_expected_unitary_d := 1;;
+gap> L:=PreservedFormspace(G);;
+gap> TestMatricesAreForms2(G, L); # found in custom_test_functions.g, Tests whether the matrices returned are forms preserved by G or not
+[ "Ok", "Ok" ]
+gap> Size(L[1])=f_space_expected_normal_d;
+true
+gap> Size(L[2])=f_space_expected_unitary_d;
+true
+gap> STOP_TEST("Formspace: Preserved Formspace G6");
diff --git a/tst/formspace/preserved_formspace/test_G7.tst b/tst/formspace/preserved_formspace/test_G7.tst
new file mode 100644
index 0000000..2111c5b
--- /dev/null
+++ b/tst/formspace/preserved_formspace/test_G7.tst
@@ -0,0 +1,14 @@
+gap> START_TEST("Formspace: Preserved Formspace G7");
+gap> G := G7;; # Some groups are defined in interesting_groups.g
+gap> R := PseudoRandom(GL(10, GF(2^8)));;
+gap> G := G^R;;
+gap> f_space_expected_normal_d := 0;; # the dimensions of the expected formspaces 
+gap> f_space_expected_unitary_d := 2;;
+gap> L:=PreservedFormspace(G);;
+gap> TestMatricesAreForms2(G, L); # found in custom_test_functions.g, Tests whether the matrices returned are forms preserved by G or not
+[ "Ok", "Ok" ]
+gap> Size(L[1])=f_space_expected_normal_d;
+true
+gap> Size(L[2])=f_space_expected_unitary_d;
+true
+gap> STOP_TEST("Formspace: Preserved Formspace G7");
diff --git a/tst/formspace/preserved_formspace/test_G8.tst b/tst/formspace/preserved_formspace/test_G8.tst
new file mode 100644
index 0000000..fd46fd2
--- /dev/null
+++ b/tst/formspace/preserved_formspace/test_G8.tst
@@ -0,0 +1,14 @@
+gap> START_TEST("Formspace: Preserved Formspace G8");
+gap> G := G8;; # Some groups are defined in interesting_groups.g
+gap> R := PseudoRandom(GL(10, GF(2^8)));;
+gap> G := G^R;;
+gap> f_space_expected_normal_d := 0;; # the dimensions of the expected formspaces 
+gap> f_space_expected_unitary_d := 2;;
+gap> L:=PreservedFormspace(G);;
+gap> TestMatricesAreForms2(G, L); # found in custom_test_functions.g, Tests whether the matrices returned are forms preserved by G or not
+[ "Ok", "Ok" ]
+gap> Size(L[1])=f_space_expected_normal_d;
+true
+gap> Size(L[2])=f_space_expected_unitary_d;
+true
+gap> STOP_TEST("Formspace: Preserved Formspace G8");
diff --git a/tst/formspace/preserved_formspace/test_GO(5_3).tst b/tst/formspace/preserved_formspace/test_GO(5_3).tst
new file mode 100644
index 0000000..d5aac70
--- /dev/null
+++ b/tst/formspace/preserved_formspace/test_GO(5_3).tst
@@ -0,0 +1,14 @@
+gap> START_TEST("Formspace: Preserved Formspace GO(5,3)");
+gap> G := GO(5,3);; # Some groups are defined in interesting_groups.g
+gap> R := PseudoRandom(GL(5, GF(3)));;
+gap> G := G^R;;
+gap> f_space_expected_normal_d := 1;; # the dimensions of the expected formspaces 
+gap> f_space_expected_unitary_d := 0;;
+gap> L:=PreservedFormspace(G);;
+gap> TestMatricesAreForms2(G, L); # found in custom_test_functions.g, Tests whether the matrices returned are forms preserved by G or not
+[ "Ok", "Ok" ]
+gap> Size(L[1])=f_space_expected_normal_d;
+true
+gap> Size(L[2])=f_space_expected_unitary_d;
+true
+gap> STOP_TEST("Formspace: Preserved Formspace GO(5,3)");
diff --git a/tst/formspace/preserved_formspace/test_GP22.tst b/tst/formspace/preserved_formspace/test_GP22.tst
new file mode 100644
index 0000000..0aaf41a
--- /dev/null
+++ b/tst/formspace/preserved_formspace/test_GP22.tst
@@ -0,0 +1,14 @@
+gap> START_TEST("Formspace: Preserved Formspace GP22");
+gap> G := GP22;; # Some groups are defined in interesting_groups.g
+gap> R := PseudoRandom(GL(3, GF(7^2)));;
+gap> G := G^R;;
+gap> f_space_expected_normal_d := 0;; # the dimensions of the expected formspaces 
+gap> f_space_expected_unitary_d := 1;;
+gap> L:=PreservedFormspace(G);;
+gap> TestMatricesAreForms2(G, L); # found in custom_test_functions.g, Tests whether the matrices returned are forms preserved by G or not
+[ "Ok", "Ok" ]
+gap> Size(L[1])=f_space_expected_normal_d;
+true
+gap> Size(L[2])=f_space_expected_unitary_d;
+true
+gap> STOP_TEST("Formspace: Preserved Formspace GP22");
diff --git a/tst/formspace/preserved_formspace/test_Group(SP(4_5).1).tst b/tst/formspace/preserved_formspace/test_Group(SP(4_5).1).tst
new file mode 100644
index 0000000..d26026c
--- /dev/null
+++ b/tst/formspace/preserved_formspace/test_Group(SP(4_5).1).tst
@@ -0,0 +1,14 @@
+gap> START_TEST("Formspace: Preserved Formspace Group(SP(4,5).1)");
+gap> G := Group(SP(4,5).1);; # Some groups are defined in interesting_groups.g
+gap> R := PseudoRandom(GL(4, GF(5)));;
+gap> G := G^R;;
+gap> f_space_expected_normal_d := 6;; # the dimensions of the expected formspaces 
+gap> f_space_expected_unitary_d := 0;;
+gap> L:=PreservedFormspace(G);;
+gap> TestMatricesAreForms2(G, L); # found in custom_test_functions.g, Tests whether the matrices returned are forms preserved by G or not
+[ "Ok", "Ok" ]
+gap> Size(L[1])=f_space_expected_normal_d;
+true
+gap> Size(L[2])=f_space_expected_unitary_d;
+true
+gap> STOP_TEST("Formspace: Preserved Formspace Group(SP(4,5).1)");
diff --git a/tst/formspace/preserved_formspace/test_Gtriv.tst b/tst/formspace/preserved_formspace/test_Gtriv.tst
new file mode 100644
index 0000000..7c9af42
--- /dev/null
+++ b/tst/formspace/preserved_formspace/test_Gtriv.tst
@@ -0,0 +1,14 @@
+gap> START_TEST("Formspace: Preserved Formspace Gtriv");
+gap> G := Gtriv;; # Some groups are defined in interesting_groups.g
+gap> R := PseudoRandom(GL(10, GF(5)));;
+gap> G := G^R;;
+gap> f_space_expected_normal_d := 100;; # the dimensions of the expected formspaces 
+gap> f_space_expected_unitary_d := 0;;
+gap> L:=PreservedFormspace(G);;
+gap> TestMatricesAreForms2(G, L); # found in custom_test_functions.g, Tests whether the matrices returned are forms preserved by G or not
+[ "Ok", "Ok" ]
+gap> Size(L[1])=f_space_expected_normal_d;
+true
+gap> Size(L[2])=f_space_expected_unitary_d;
+true
+gap> STOP_TEST("Formspace: Preserved Formspace Gtriv");
diff --git a/tst/formspace/preserved_formspace/test_SU(4_5).tst b/tst/formspace/preserved_formspace/test_SU(4_5).tst
new file mode 100644
index 0000000..be18de1
--- /dev/null
+++ b/tst/formspace/preserved_formspace/test_SU(4_5).tst
@@ -0,0 +1,14 @@
+gap> START_TEST("Formspace: Preserved Formspace SU(4,5)");
+gap> G := SU(4,5);; # Some groups are defined in interesting_groups.g
+gap> R := PseudoRandom(GL(4, GF(5^2)));;
+gap> G := G^R;;
+gap> f_space_expected_normal_d := 0;; # the dimensions of the expected formspaces 
+gap> f_space_expected_unitary_d := 1;;
+gap> L:=PreservedFormspace(G);;
+gap> TestMatricesAreForms2(G, L); # found in custom_test_functions.g, Tests whether the matrices returned are forms preserved by G or not
+[ "Ok", "Ok" ]
+gap> Size(L[1])=f_space_expected_normal_d;
+true
+gap> Size(L[2])=f_space_expected_unitary_d;
+true
+gap> STOP_TEST("Formspace: Preserved Formspace SU(4,5)");
diff --git a/tst/formspace/preserved_formspace/test_Sp(4_5).tst b/tst/formspace/preserved_formspace/test_Sp(4_5).tst
new file mode 100644
index 0000000..57749a9
--- /dev/null
+++ b/tst/formspace/preserved_formspace/test_Sp(4_5).tst
@@ -0,0 +1,14 @@
+gap> START_TEST("Formspace: Preserved Formspace Sp(4,5)");
+gap> G := Sp(4,5);; # Some groups are defined in interesting_groups.g
+gap> R := PseudoRandom(GL(4, GF(5)));;
+gap> G := G^R;;
+gap> f_space_expected_normal_d := 1;; # the dimensions of the expected formspaces 
+gap> f_space_expected_unitary_d := 0;;
+gap> L:=PreservedFormspace(G);;
+gap> TestMatricesAreForms2(G, L); # found in custom_test_functions.g, Tests whether the matrices returned are forms preserved by G or not
+[ "Ok", "Ok" ]
+gap> Size(L[1])=f_space_expected_normal_d;
+true
+gap> Size(L[2])=f_space_expected_unitary_d;
+true
+gap> STOP_TEST("Formspace: Preserved Formspace Sp(4,5)");
diff --git a/tst/formspace/recog_preserved_sesqui_with_scalars_test.tst b/tst/formspace/recog_preserved_sesqui_with_scalars_test.tst
new file mode 100644
index 0000000..46331a6
--- /dev/null
+++ b/tst/formspace/recog_preserved_sesqui_with_scalars_test.tst
@@ -0,0 +1,38 @@
+gap> START_TEST("Forms: recog_preserved_sesqui_with_scalars_test.tst"); # Adapted from test_recog.tst for PreservedSesquilinearFormsWithScalars
+gap> g := Sp(6,3);
+Sp(6,3)
+gap> forms := PreservedSesquilinearFormsWithScalars(g);
+[ < bilinear form > ]
+gap> TestPreservedSesquilinearForms(g,forms);
+true
+gap> g := SU(4,3);
+SU(4,3)
+gap> forms := PreservedSesquilinearFormsWithScalars(g);
+[ < hermitian form > ]
+gap> TestPreservedSesquilinearForms(g,forms);
+true
+gap> g := SO(1,4,3);
+SO(+1,4,3)
+gap> forms := PreservedSesquilinearFormsWithScalars(g);
+[ < bilinear form > ]
+gap> TestPreservedSesquilinearForms(g,forms);
+true
+gap> g := SO(1,4,4);
+GO(+1,4,4)
+gap> forms := PreservedSesquilinearFormsWithScalars(g);
+[ < bilinear form > ]
+gap> TestPreservedSesquilinearForms(g,forms);
+true
+gap> g := SO(-1,4,3);
+SO(-1,4,3)
+gap> forms := PreservedSesquilinearFormsWithScalars(g);
+[ < bilinear form > ]
+gap> TestPreservedSesquilinearForms(g,forms);
+true
+gap> g := SO(5,3);
+SO(0,5,3)
+gap> forms := PreservedSesquilinearFormsWithScalars(g);
+[ < bilinear form > ]
+gap> TestPreservedSesquilinearForms(g,forms);
+true
+gap> STOP_TEST("recog_preserved_sesqui_with_scalars_test.tst", 10000 );
diff --git a/tst/formspace/recog_preserved_with_scalars_test.tst b/tst/formspace/recog_preserved_with_scalars_test.tst
new file mode 100644
index 0000000..6883e5b
--- /dev/null
+++ b/tst/formspace/recog_preserved_with_scalars_test.tst
@@ -0,0 +1,36 @@
+gap> START_TEST("Forms: recog_preserved_with_scalars_test.tst"); # Adapted from test_recog.tst for PreservedFormsWithScalars
+gap> g := Sp(6,3);
+Sp(6,3)
+gap> forms := PreservedFormsWithScalars(g);
+[ < bilinear form > ]
+gap> TestPreservedSesquilinearForms(g,forms);
+true
+gap> g := SU(4,3);
+SU(4,3)
+gap> forms := PreservedFormsWithScalars(g);
+[ < hermitian form > ]
+gap> TestPreservedSesquilinearForms(g,forms);
+true
+gap> g := SO(1,4,3);
+SO(+1,4,3)
+gap> forms := PreservedFormsWithScalars(g);
+[ < quadratic form > ]
+gap> TestPreservedSesquilinearForms(g,[BilinearFormByQuadraticForm(forms[1])]);
+true
+gap> g := SO(1,4,4);
+GO(+1,4,4)
+gap> forms := PreservedFormsWithScalars(g); # Because of Char = 2 BilinearFormByQuadraticForm(forms[1]) will not work
+[ < quadratic form > ]
+gap> g := SO(-1,4,3);
+SO(-1,4,3)
+gap> forms := PreservedFormsWithScalars(g);
+[ < quadratic form > ]
+gap> TestPreservedSesquilinearForms(g,[BilinearFormByQuadraticForm(forms[1])]);
+true
+gap> g := SO(5,3);
+SO(0,5,3)
+gap> forms := PreservedFormsWithScalars(g);
+[ < quadratic form > ]
+gap> TestPreservedSesquilinearForms(g,[BilinearFormByQuadraticForm(forms[1])]);
+true
+gap> STOP_TEST("recog_preserved_with_scalars_test.tst", 10000 );
diff --git a/tst/formspace/test_poly_eval.tst b/tst/formspace/test_poly_eval.tst
new file mode 100644
index 0000000..76b182e
--- /dev/null
+++ b/tst/formspace/test_poly_eval.tst
@@ -0,0 +1,4 @@
+gap> START_TEST("Formspace: Poly Eval with Frobenius Normal Form");
+gap> TestPolyEval(false); # false to specify that we are not benchmarking..
+"Ok"
+gap> STOP_TEST("Formspace: Poly Eval with Frobenius Normal Form");
diff --git a/tst/interesting_groups.g b/tst/interesting_groups.g
new file mode 100644
index 0000000..a144a9d
--- /dev/null
+++ b/tst/interesting_groups.g
@@ -0,0 +1,404 @@
+## obtained through the classicalmaximals package
+
+# FindIrredNotAbs := function(Gs)
+#     local i;
+#     for i in [1..Length(Gs)] do
+#         if IsIrreducible(Gs[i]) and not IsAbsolutelyIrreducible(Gs[i]) then
+#             return i;
+#         fi;
+#     od;
+# end;
+
+#         gap --packagedirs $PWD
+#         Read("tst/formspace/custom_test_functions.g");
+
+## various micselaneous functions..
+
+Disp := function(Fspaces) 
+    local F, i;
+    Print("bilinear:\n");
+    for i in [1..Size(Fspaces[1])] do
+        F := Fspaces[1][i];
+        Print("--------------------------------------\n");
+        Display(F);
+    od;
+    Print("unitary:\n");
+    for i in [1..Size(Fspaces[2])] do
+        F := Fspaces[2][i];
+        Print("--------------------------------------\n");
+        Display(F);
+    od;
+end;
+
+PolyEval := function(f, M)
+     local F, n, frob;
+     n := DimensionsMat(M)[1];
+     F := DefaultFieldOfMatrix(M);
+     frob := FrobeniusNormalForm(M);
+     return FORMS_EvaluatePolynomialWithFrobenius(CoefficientsOfUnivariatePolynomial(f), M, frob, Inverse(frob[2]), F, n);
+end;
+
+
+## groups
+
+## groups from issue https://github.com/gap-packages/forms/issues/83
+G83_1 := Group([ [ [ Z(3^2)^3, Z(3^2)^6, Z(3^2)^7 ], [ Z(3^2)^7, Z(3^2)^2, Z(3)^0 ], [ Z(3^2)^7, Z(3^2)^7, 0*Z(3) ] ], [ [ 0*Z(3), 0*Z(3), Z(3) ], [ Z(3^2)^7, Z(3)^0, Z(3^2)^7 ], [ Z(3), Z(3^2), Z(3) ] ], [ [ Z(3^2), 0*Z(3), 0*Z(3) ], [ 0*Z(3), Z(3^2), 0*Z(3) ], [ 0*Z(3), 0*Z(3), Z(3^2) ] ] ]);
+
+G83_2 := Group(Z(5)^0*[
+[ [ 0, 1, 0, 0, 0, 0, 0, 0 ],
+  [ 0, 0, 0, 0, 1, 0, 0, 0 ],
+  [ 0, 0, 0, 0, 0, 0, 0, 1 ],
+  [ 2, 3, 4, 3, 1, 3, 2, 4 ],
+  [ 1, 4, 0, 0, 1, 0, 0, 0 ],
+  [ 0, 2, 2, 2, 4, 4, 2, 1 ],
+  [ 2, 3, 4, 2, 1, 3, 3, 4 ],
+  [ 2, 2, 3, 0, 1, 0, 0, 3 ] ],
+[ [ 0, 0, 1, 0, 0, 0, 0, 0 ],
+  [ 0, 0, 0, 0, 0, 1, 0, 0 ],
+  [ 2, 4, 0, 3, 2, 1, 2, 2 ],
+  [ 2, 4, 0, 2, 4, 4, 4, 4 ],
+  [ 4, 4, 1, 0, 1, 0, 3, 2 ],
+  [ 1, 1, 2, 3, 2, 4, 1, 3 ],
+  [ 3, 2, 0, 3, 1, 4, 1, 4 ],
+  [ 3, 3, 1, 4, 0, 3, 4, 3 ] ],
+[ [ 0, 0, 0, 1, 0, 0, 0, 0 ],
+  [ 0, 0, 0, 0, 0, 0, 1, 0 ],
+  [ 2, 2, 3, 1, 3, 3, 2, 4 ],
+  [ 3, 0, 0, 0, 0, 0, 0, 0 ],
+  [ 1, 0, 2, 1, 0, 3, 2, 1 ],
+  [ 1, 3, 2, 4, 3, 1, 2, 2 ],
+  [ 0, 3, 0, 0, 0, 0, 0, 0 ],
+  [ 0, 1, 3, 0, 3, 1, 0, 1 ] ]
+]);
+
+
+# trivial group over GF(5)
+Gtriv := Group([[ Z(5^4), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5^4)^623 ]]^0);
+
+## subgroup of SP(10, 5^4), which is not irreducible (2d formspace)
+G1 := Group([ [ Z(5^4), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ], 
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5^4)^623 ] ], [ [ Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5) ],
+  [ Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0 ],
+  [ 0*Z(5), Z(5)^2, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ] ], [ [ Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), Z(5^4), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5^4)^623, 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0 ] ],
+  [ [ Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^2, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0 ] ]
+  );
+
+## reducible subgroup of U(4, 5^4) with 2d formspace
+G2 := Group([ [ Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5) ], [ 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5) ], [ 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0 ] ], [ [ Z(5^4), 0*Z(5), 0*Z(5), 0*Z(5) ], [ 0*Z(5), Z(5^4)^467, Z(5^4)^194, 0*Z(5) ],
+  [ 0*Z(5), Z(5^4)^532, Z(5^4)^441, 0*Z(5) ], [ 0*Z(5), 0*Z(5), 0*Z(5), Z(5^4)^599 ] ], [ [ Z(5^4)^547, Z(5^4)^217, Z(5^4)^324, Z(5)^0 ], [ Z(5^4)^612, Z(5^4)^65, Z(5), 0*Z(5) ],
+  [ Z(5^4)^433, Z(5^2)^19, Z(5^4)^377, 0*Z(5) ], [ Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5) ] ], [ [ Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5) ], [ 0*Z(5), Z(5^2)^21, Z(5^4)^247, 0*Z(5) ],
+  [ 0*Z(5), Z(5^4)^585, Z(5^2)^7, 0*Z(5) ], [ 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0 ] ]);
+
+## subgroup of omega(7, 5^2) which is reducible and preserves both unitary and bilinear forms
+G3 := Group([ [ Z(5^2)^5, Z(5^2)^9, 0*Z(5), 0*Z(5), 0*Z(5), Z(5^2)^21, Z(5^2)^13 ],
+  [ Z(5^2)^21, Z(5^2)^21, 0*Z(5), Z(5^2)^21, 0*Z(5), 0*Z(5), Z(5^2)^21 ],
+  [ 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5) ],
+  [ 0*Z(5), Z(5^2)^15, 0*Z(5), Z(5^2)^9, 0*Z(5), Z(5^2)^15, 0*Z(5) ],
+  [ 0*Z(5), 0*Z(5), 0*Z(5), 0*Z(5), Z(5)^0, 0*Z(5), 0*Z(5) ],
+  [ Z(5^2)^9, 0*Z(5), 0*Z(5), Z(5^2)^21, 0*Z(5), Z(5^2)^21, Z(5^2)^9 ],
+  [ Z(5^2)^13, Z(5^2)^9, 0*Z(5), 0*Z(5), 0*Z(5), Z(5^2)^21, Z(5^2)^5 ] ],
+  [ [ Z(5^2)^7, 0*Z(5), Z(5^2)^2, 0*Z(5), Z(5^2)^2, 0*Z(5), Z(5^2)^9 ],
+  [ 0*Z(5), Z(5^2)^10, Z(5^2)^10, Z(5^2)^23, Z(5^2)^22, Z(5^2)^3, 0*Z(5) ],
+  [ Z(5^2)^2, Z(5^2)^10, Z(5^2)^16, Z(5), Z(5^2)^22, Z(5^2)^10, Z(5^2)^14 ],
+  [ 0*Z(5), Z(5^2)^17, Z(5)^0, Z(5^2)^5, Z(5)^2, Z(5^2)^17, 0*Z(5) ],
+  [ Z(5^2)^2, Z(5^2)^22, Z(5^2)^22, Z(5)^3, Z(5^2)^16, Z(5^2)^22, Z(5^2)^14 ],
+  [ 0*Z(5), Z(5^2)^3, Z(5^2)^10, Z(5^2)^23, Z(5^2)^22, Z(5^2)^10, 0*Z(5) ],
+  [ Z(5^2)^9, 0*Z(5), Z(5^2)^14, 0*Z(5), Z(5^2)^14, 0*Z(5), Z(5^2)^7 ] ],
+  [ [ Z(5^2)^15, Z(5^2)^4, Z(5^2)^3, Z(5^2)^22, Z(5^2)^3, Z(5)^2, Z(5^2)^21 ],
+  [ Z(5)^2, Z(5)^0, Z(5^2)^9, Z(5^2)^9, Z(5^2)^9, Z(5^2)^19, Z(5^2)^19 ],
+  [ Z(5^2)^13, Z(5^2)^7, Z(5)^3, Z(5^2)^22, Z(5), Z(5)^0, Z(5^2)^13 ],
+  [ Z(5^2)^4, Z(5^2)^8, 0*Z(5), Z(5^2)^21, 0*Z(5), Z(5^2), Z(5^2)^4 ],
+  [ Z(5^2)^13, Z(5^2)^7, Z(5), Z(5^2)^22, Z(5)^3, Z(5)^0, Z(5^2)^13 ],
+  [ Z(5^2)^19, Z(5)^0, Z(5^2)^21, Z(5^2)^9, Z(5^2)^21, Z(5^2)^19, Z(5)^2 ],
+  [ Z(5^2)^9, Z(5^2)^16, Z(5^2)^3, Z(5^2)^10, Z(5^2)^3, Z(5)^0, Z(5^2)^3 ] ]
+  );
+
+## Subgroup of O+(14, 7) which is irreducible, but not absolutely irreducible 
+## Intersting is that the foms package seems to reliably find a hermitian and quadratic form, It seems that it is hard to find a cyclic group element in this particular group! 
+G4 := Group([ [ Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), Z(7)^4, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^3, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^4, Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^3, Z(7)^4, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^4, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^3, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^5, 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^3, Z(7)^2, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^5, Z(7)^3, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^0, 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^2, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^3, 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^0 ]
+ ], [ [ Z(7)^3, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^2, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), Z(7)^3, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^3, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^3, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^3, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), Z(7)^5, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^0, 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^3, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^3, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^3, 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^3, 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^3, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^3 ]
+ ], [ [ Z(7)^3, Z(7)^2, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^5, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ Z(7)^4, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^0, Z(7)^5, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^4 ],
+  [ Z(7), Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^2, Z(7)^4 ],
+  [ 0*Z(7), 0*Z(7), Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^0, 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), Z(7)^4, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^0, 0*Z(7), 0*Z(7), 0*Z(7) ],
+  [ 0*Z(7), Z(7)^5, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^5, 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), 0*Z(7), Z(7)^3 ]
+ ]);
+
+
+
+# Subgroup of O-(14, 11^2) another group which is irreducible but not absolutely irreducible. it is intersting that once again my program fails to find a cyclic element. 
+G5 := Group([ [ 0*Z(11), 0*Z(11), Z(11)^0, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11) ], [ 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^0, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11) ],
+  [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^0, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11) ], [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^0, 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11) ],
+  [ Z(11)^0, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^3, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), Z(11)^6 ], [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11), 0*Z(11), Z(11^2)^81, 0*Z(11) ],
+  [ Z(11)^3, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^5, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11) ], [ Z(11^2)^6, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^5, 0*Z(11),
+       0*Z(11), 0*Z(11), 0*Z(11), Z(11^2)^6, 0*Z(11) ],
+  [ 0*Z(11), Z(11^2)^39, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^8, Z(11^2)^66, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+      Z(11)^0, 0*Z(11) ], [ Z(11)^4, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+      0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11) ], [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+      Z(11)^0, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11) ],
+  [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^0, 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11) ], [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), Z(11)^0, 0*Z(11), 0*Z(11), 0*Z(11) ],
+  [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^0,
+       0*Z(11), 0*Z(11) ] ], [ [ Z(11^2)^29, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       Z(11^2)^119, 0*Z(11) ], [ 0*Z(11), Z(11^2)^30, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^5 ],
+  [ 0*Z(11), 0*Z(11), Z(11)^5, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11) ], [ 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^5, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11) ],
+  [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^5, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11) ], [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^5, 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11) ],
+  [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^5, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11) ], [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^5, 0*Z(11),
+       0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11) ],
+  [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^5, 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11) ], [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       Z(11)^5, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11) ],
+  [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^5, 0*Z(11),
+       0*Z(11), 0*Z(11) ], [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11), Z(11)^5, 0*Z(11), 0*Z(11) ],
+  [ Z(11^2)^86, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       Z(11^2)^29, 0*Z(11) ], [ 0*Z(11), Z(11^2)^27, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11^2)^30 ] ],[ [ Z(11^2)^90, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       Z(11)^0, 0*Z(11) ], [ 0*Z(11), Z(11^2)^89, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11^2)^59 ],
+  [ 0*Z(11), 0*Z(11), Z(11)^0, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11) ], [ 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^0, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11) ],
+  [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11^2)^90, 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^0, 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11) ], [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11^2)^89, 0*Z(11), 0*Z(11), 0*Z(11),
+       Z(11^2)^59, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11) ],
+  [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^0, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11) ], [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^0, 0*Z(11),
+       0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11) ],
+  [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11^2)^87, 0*Z(11), 0*Z(11), 0*Z(11), Z(11^2)^90, 0*Z(11), 0*Z(11), 0*Z(11),
+      0*Z(11), 0*Z(11) ], [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11^2)^26, 0*Z(11), 0*Z(11), 0*Z(11),
+       Z(11^2)^89, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11) ],
+  [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11)^0, 0*Z(11),
+       0*Z(11), 0*Z(11) ], [ 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11), Z(11)^0, 0*Z(11), 0*Z(11) ],
+  [ Z(11^2)^87, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       Z(11^2)^90, 0*Z(11) ], [ 0*Z(11), Z(11^2)^26, 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11),
+       0*Z(11), 0*Z(11), 0*Z(11), 0*Z(11), Z(11^2)^89 ] ]);
+
+# subgroup of U(10, 2^8) which sometimes lands in the scalar case c A = A* (for hermitian forms) (also non irreducible)
+G6 := Group([ [ Z(2^8), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), Z(2^8)^254, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^16, 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^239 ] ], [ [ Z(2)^0, 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), Z(2)^0 ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2) ] ], [ [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^4), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^4)^14, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0 ] ], [ [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, Z(2^4)^2, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^4)^13, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0 ] ], [ [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2) ],
+  [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0 ] ], [ [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0 ] ], [ [ Z(2^8), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^254, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^16, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^239 ] ]);
+
+# Subgroup of U(10, 2^8) (non irreducible) that lands in the annoying case of multiple unitary forms 
+G7 := Group([ [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2),
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2),
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0,
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0 ] ], [ [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2),
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^67, Z(2^8)^33, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^84, Z(2^8)^169, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^239, 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2),
+       Z(2)^0, 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0 ] ], [ [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2),
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ Z(2^8)^222, 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^183, Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0 ],
+  [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ Z(2^8)^123, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2),
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2) ] ], [ [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2),
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^4)^12, Z(2^4)^3, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^4)^6, Z(2^2)^2, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0,
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0 ] ]);
+
+# Subgroup of U(10, 2^8) (non irreducible) that lands in the annoying case of multiple unitary forms (similar as G7)
+G8 := Group([ [ Z(2^8), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2),
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^169, Z(2^8)^33, 0*Z(2),
+       0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^84, Z(2^8)^67, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0,
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^239 ] ], [ [ Z(2^8)^222, 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^123, Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0 ],
+  [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2),
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ Z(2^8)^183, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0,
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2) ],
+  [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ] ], [ [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2),
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^67, Z(2^8)^33, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^84, Z(2^8)^169, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^8)^239, 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2),
+       Z(2)^0, 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0 ] ], [ [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2),
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), Z(2^8)^222, 0*Z(2), 0*Z(2), Z(2^8)^183, Z(2)^0, 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2) ],
+  [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), Z(2^8)^123, 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2),
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0 ] ], [ [ Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2),
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^4)^12, Z(2^4)^3, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2^4)^6, Z(2^2)^2, 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2), 0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0,
+       0*Z(2), 0*Z(2) ], [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0, 0*Z(2) ],
+  [ 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), 0*Z(2), Z(2)^0 ] ]);
+
+# The group from the forms package pull request https://github.com/gap-packages/forms/issues/22
+GP22 := Group([ [ Z(7)^0, 0*Z(7), 0*Z(7) ], [ Z(7^2)^33, Z(7^2)^14, Z(7^2)^26 ], [ Z(7^2)^19, Z(7^2)^31, Z(7^2)^5 ] ], [ [ Z(7^2)^39, Z(7^2)^9, Z(7)^3 ], [ Z(7^2)^25, Z(7)^2, Z(7^2)^6 ], [ Z(7^2)^7, Z(7)^4, Z(7^2)^28 ] ]);
diff --git a/tst/testall.g b/tst/testall.g
index 8667153..7deec07 100644
--- a/tst/testall.g
+++ b/tst/testall.g
@@ -6,6 +6,9 @@ if not IsBound(DescribesInvariantQuadraticForm) then
   Add( exclude, "adv/classic.tst" );
 fi;
 
+ReadPackage("forms", "tst/interesting_groups.g");
+ReadPackage("forms", "tst/formspace/custom_test_functions.g");
+
 TestDirectory(DirectoriesPackageLibrary("forms", "tst"),
     rec(
       exitGAP := true,