diff --git a/.gitignore b/.gitignore
index 85d5146..ccd5fde 100644
--- a/.gitignore
+++ b/.gitignore
@@ -22,3 +22,4 @@
 /tmp
 /gh-pages
 tmp.json
+profiling
\ No newline at end of file
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/formspace.gd b/lib/formspace.gd
new file mode 100644
index 0000000..817f6c0
--- /dev/null
+++ b/lib/formspace.gd
@@ -0,0 +1,4 @@
+DeclareOperation("PreservedFormspace", [IsMatrixGroup, IsVector and IsFFECollection, IsBool]);
+DeclareOperation("PreservedFormspace", [IsMatrixGroup]);
+
+DeclareOperation("FilterFormspace", [IsList, IsFinite and IsField, IsBool]);
diff --git a/lib/formspace.gi b/lib/formspace.gi
new file mode 100644
index 0000000..5de3e6e
--- /dev/null
+++ b/lib/formspace.gi
@@ -0,0 +1,681 @@
+#! @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;
+    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)
+    local cur_group_element, cur_scalar, known_elements, i, mode, n, res, j, known_scalars, best_known_element_index, best_known_res, best_known_length, mod_elem, g, e;
+
+    known_elements := ShallowCopy(Gens);
+    known_scalars := ShallowCopy(Lambdas);
+
+    n := NrRows(Gens[1]);
+    i := Size(known_elements);
+    mod_elem := 1;
+    best_known_length := n + 1;
+    while i < 25 do # 25 is a magic number. maybe investigate a good number here
+        # no accidental identity mat
+        if i < 10 then # dont start with inverting to not accidentally make the identity
+            mode := 1;
+        else
+            mode := Random(1, 2);
+        fi;
+        j := Random([1..Size(known_elements)]);
+        cur_group_element := known_elements[j];
+        cur_scalar := known_scalars[j];
+
+        if mode = 1 then
+            j := Random([1..Size(known_elements)]);
+            cur_group_element := cur_group_element * known_elements[j];
+            cur_scalar := cur_scalar * known_scalars[j];
+        elif mode = 2 then
+            j := Random([1..Size(known_elements)]);
+            cur_group_element := cur_group_element / known_elements[j];
+            cur_scalar := cur_scalar / Inverse(known_scalars[j]);
+        fi;
+        if i mod mod_elem = 0 then
+            res := FrobeniusNormalForm(cur_group_element);
+            if Size(res[3]) = 1 then
+                return Concatenation([cur_group_element, cur_scalar], res, [i]); 
+            fi;
+            if Size(res[3]) <= best_known_length and not (cur_group_element in Gens) then
+                best_known_element_index := i + 1;
+                best_known_res := res;
+                best_known_length := Size(res[3]);
+            fi;
+        fi;
+        Add(known_elements, cur_group_element);
+        Add(known_scalars, cur_scalar);
+        i := i + 1;
+    od;
+    # Print("only found element of length ", best_known_length, "\n");
+    if best_known_length = n + 1 then
+        return Concatenation([Gens[1], Lambdas[1]], FrobeniusNormalForm(Gens[1]), [-1]);
+    fi;
+    return Concatenation([known_elements[best_known_element_index], known_scalars[best_known_element_index]], best_known_res, [i]);
+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);
+    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;
+
+# Spins the stuff
+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);
+    for i in [1..j] do
+        if i = j then
+            end_pos := n;
+        else
+            end_pos := frob_base_blocks[i + 1] - 1;
+        fi;
+        A[frob_base_blocks[i]] := Images[i];
+        for k in [(frob_base_blocks[i] + 1)..end_pos] do
+            A[k] := A[k - 1]*spin_elem;
+        od;
+    od;
+    ConvertToMatrixRep(A, F);
+    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);
+    ## FrobBasSpin needed?
+    # this mulitiplication might not be needed here 
+    # just multiply the found basis in the end when all the condition matrices null spaces got found already (i.e. just return C here)#
+    # this could potentially save up to n*d*B matrix matrix multiplications so quite significant!!! definetly investigate!
+    # ConvertToMatrixRep(C, F);
+    return frob_base_inv * C; #to actually evaluate the polynomial
+    
+    # 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 ten billion 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;
+    # Display(coeffs_c);
+    j := Size(frob_base[3]);
+    Ps := NullMat(n * j, n, F);
+    ConvertToMatrixRep(Ps, F);
+    # Print(frob_base[3], "\n");
+    for i in [1..j] do
+        if i = j then
+            b_end := n;
+        else
+            b_end := frob_base[3][i + 1] - 1;
+        fi;
+        # Print("[", ((i - 1)*n + 1), ",", (i*n), "\n");
+        #Print("->", frob_base[3][i],",",b_end, "\n");
+        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);
+        # cpol := UnivariatePolynomial(F, coeffs_c);
+        # fpol := UnivariatePolynomial(F, coeffs_f);
+        # Ps{[((i - 1)*n + 1)..(i*n)]}{[1..n]} :=
+        #     cpol(g_star_inv_scaled) * h_star - fpol(g_star_inv_scaled);
+        # Print(aua);
+    od;
+    return Ps;
+end;
+
+
+
+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);
+    ConvertToMatrixRep(A, F);
+    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;
+
+# 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;
+
+# find symplectic and symmetric matrices in Forms
+# returns bases [[symmetric forms], [symplectic forms]]
+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;
+
+    # computes if f = f^T or f = -f^T or none
+    # return 0 if f = -f^T
+    # return 1 if f = f^T
+    # return 2 if f <> f^T and f <> -f^T
+    TransposedEqual := function(f, F, n) 
+        local i, j, symmetric_possible, symplectic_possible, mat;
+        
+        if Characteristic(F) = 2 then
+        # - = + now
+            for i in [1..n] do
+                for j in [1..n] do
+                    if f[i, j] <> f[j, i] then
+                        return 2;
+                    fi;
+                od;
+            od;
+            return 1;
+        fi;
+        symmetric_possible := true;
+        symplectic_possible := true;
+        for i in [1..n] do
+            for j in [i..n] do
+                if symmetric_possible and f[i, j] <> f[j, i] then
+                    symmetric_possible := false;
+                    if not symplectic_possible then
+                        return 2;
+                    fi;
+                fi;
+                if symplectic_possible and f[i, j] <> -f[j, i] then
+                    symplectic_possible := false;
+                    if not symmetric_possible then
+                        return 2;
+                    fi;
+                fi;
+            od;
+        od;
+        if symmetric_possible then
+            return 1;
+        fi;
+        if symplectic_possible then
+            return 0;
+        fi;
+        return 2;
+    end;
+
+    if Size(Forms) = 0 then
+        return [];
+    fi;
+    if Size(Forms) = 1 then
+        transposed_equal_result := TransposedEqual(Forms[1], F, n);
+        if transposed_equal_result = 2 then
+            return [];
+        fi;
+        if transposed_equal_result = 1 then
+            return [[Forms[1]], []];
+        fi;
+        if transposed_equal_result = 2 then
+            return [[], [Forms[1]]];
+        fi;
+
+    fi;
+
+    if Characteristic(F) <> 2 then
+        # 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, [], 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));
+
+        # for mat in symmetric_base_vecs do
+        #     ConvertToMatrixRep(mat, F);
+        # od;
+
+        # for mat in symplectic_base_vecs do
+        #     ConvertToMatrixRep(mat, F);
+        # od;
+
+        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: ignore the diagonal entries for the symmetric matrices, for symplectic these need to be zero.
+    # 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 the nice thing about the current algorithm for char F <> 2 is that we do not have to solve many systems of equations, rather just check wether we have the zero matrix
+
+    ## TODO: this code is broken right now, FIX!! the issue seems to lie in the function FORMS_MatrixReorganize
+    char_2_eqs := [];
+    for form in Forms do
+        Add(char_2_eqs, FORMS_MatrixReorganize(form - TransposedMat(form), n, F, n));
+    od;
+    sol := NullspaceMatDestructive(char_2_eqs);
+    out := [];
+    for form in sol do
+        Add(out, FORMS_VectorReorganize(form, n, F, n));
+    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} TODO: THIS NEEDS SOME FURTHER INVESTIGATION
+# there must be a better way to compute these matrices..
+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;
+    if Size(Forms) = 0 then
+        return [];
+    fi;
+    p := Characteristic(F);
+    q := p^(DegreeOverPrimeField(F)/2);
+    
+    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.
+        # Print(ff);
+        l := FORMS_ScalarFormIdentifyCheck(Forms[1], F, n, hom, p, q);
+        if l = fail then
+            return [];
+        fi;
+        
+        return [Forms[1] * l];
+    fi;
+    # this is where it gets interesting
+    # Print("ahhh this needs work!\n");
+    # kind of okay case?
+       
+    ## TODO proof if gAg^* = cA for some c it must hold that c \in GF(q) (not in GF(q^2))
+
+    ## we use (A + A*) is a hermitian form if gAg* = A the problem here is that for A, B such that gA*g = A and gB*g = B we may loose information namely it may be the case that 1/2 (A + A*) = c/2 (B + B*) this is annoying.... one thing one could do is check whether these matrices <1/2 (A + A*), 1/2 (B + B*), ...> are lineraly independent. if that is the case, we know that all forms must have been found (but do we???). but what if not? then there might be another form... this is annoying. Then we may add matrices A, B and so on such that <C1,C2, ., D1, D2..> is a basis of F where C1 and so on are hermitian and D1, D2 and so on are not. We may write D1 = A + B where A is hermitian and B is not. we can then try to write B in terms of the other matrices??? does this help... idk :(
+    ## for char = 2 this can be a bad idea as it can make the diagonal disaapear..   oh well
+    if p <> 2 then
+        Base := MutableBasis(GF(q), [NullMat(n, n, GF(q))]);
+        # 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
+            # l := FORMS_ScalarFormIdentifyCheck(FF, F, n, hom, p, q);
+            # if l = fail then
+            tr_form := TransposedMat(FF^hom);
+            CloseMutableBasis(Base, FF + tr_form);
+            CloseMutableBasis(Base, gf_base * FF + hgf_base*tr_form);
+            # else
+                # CloseMutableBasis(Base, l * FF);
+            # fi;
+        od;
+
+        O := [];
+        baseVecs := BasisVectors(ImmutableBasis(Base));
+        # for FF in baseVecs do
+        #     Add(O, HermitianFormByMatrix(FF, F));
+        # od;
+
+        # if Size(O) = Size(Forms) then
+        #     return O;
+        # fi;
+        # if Size(baseVecs) = Size(Forms) then
+        #     return baseVecs;
+        # fi;
+        # for mat in baseVecs do
+        #     ConvertToMatrixRep(mat, F);
+        # od;
+        return baseVecs;
+    fi;
+    Print("char 2 case is missing!!");
+    # # TODO this function definetly needs work!!!
+    # Print("Could not find a basis of hermitian Forms. Returned hermitian Forms and a Bigger space of matrices that contains all possible hermitian forms. \n");
+    # return [baseVecs, Forms];
+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#function(p, g, frob_base, frob_base_inv, F, n)
+        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 := [];
+    #(Image, spin_elem, frob_base_blocks, block_index, n, F)
+    
+    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));
+        od;
+    od;
+    return OutForms;
+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, needs_checking, A, failed_check, nspace, vec;
+    first := true;
+    needs_checking := false; ## maybe remove?
+
+    for i in [1..d] do
+        if needs_checking then
+            break;
+        fi;
+        h := Gens[i];
+        h_star := FORMS_CalculateAdjoint(h, unitary, hom, n, F);
+        for j in [1..n] do
+            # vec := RandomVector(F, n);
+            vec := g_res[4][j];
+            # Display(vec);
+            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);
+            
+            ConvertToMatrixRep(Conds, F);
+            if not first then
+                nspace := NullspaceMat(W * Conds);
+                ConvertToMatrixRep(nspace, F);
+                if Size(nspace) = 0 then 
+                    # Print("empty");
+                    return [];
+                fi;
+                W := nspace * W;
+            fi;
+            if first then
+                nspace := NullspaceMat(Conds);
+                ConvertToMatrixRep(nspace, F);
+                if Size(nspace) = 0 then 
+                    # Print("empty");
+                    return [];
+                fi;
+                W := nspace;
+                first := false;
+                # technically we still need checking here..
+            fi;
+            if Size(nspace) = 1 then
+                needs_checking := true;
+                break;
+            fi;
+        od;
+    od;
+    O := [];
+    # W := frob_base_inv_star * W;
+    # Print(WWW);
+    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);
+        ## This whole checking should be removed, it is stupid to always check!!!
+        if needs_checking then
+            failed_check := false;
+            for i in [1..d] do
+                if not failed_check and Gens[i] * A * FORMS_CalculateAdjoint(Gens[i], unitary, hom, n, F) <> Lambdas[i] * A then
+                    failed_check := true;
+                fi;
+            od;
+            if not failed_check then
+                Add(O, A);
+                # if not unitary then
+                #     Add(O, BilinearFormByMatrix(A, F));
+                # else
+                #     Add(O, A);
+                # fi;
+            fi;
+        fi;
+        if not needs_checking then
+            # if not unitary then
+            #     # throws errors since the forms package only accepts symmetric bilinear forms... i think this should be changed..
+            #     # or should i solve for symmetric matrices as a system of linear equations..?
+            #     Add(O, BilinearFormByMatrix(A, F));
+            # else
+            #     Add(O, A);
+            # fi;
+            Add(O, A);
+        fi;
+    od;
+    # if unitary then
+    #     return FORMS_FilterUnitaryForms(O, F, n, hom);
+    # fi;
+    return O;
+end;
+
+
+#! @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;
+        # Prüfen ob es sich um einen endlichen körper handelt??
+        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);
+        g_inv_frob := Inverse(g_res[4]);
+        #CalculateAdjoint := function(mat, mode, hom, n)
+        g_star_inv_unscaled := TransposedMat(Inverse(g_res[1]));
+        if unitary then
+            g_star_inv_unscaled := g_star_inv_unscaled^hom;
+        fi;
+        #* g_res[2]; 
+        # Todo the computation of this can probably be sped up by using the known information about g!!!
+        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);
+        # Prüfen ob es sich um einen endlichen körper handelt??
+        Gens := GeneratorsOfGroup(G);
+        # F := DefaultFieldOfMatrix(Gens[1]);
+        n := NrRows(Gens[1]);
+        d := Size(Gens);
+        hom := fail;
+
+        if d = 1 then
+            # Todo.... !!
+            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);
+        ConvertToMatrixRep(g_res[1], F);
+        ConvertToMatrixRep(g_res[4], F);
+        g_inv_frob := Inverse(g_res[4]);
+        #CalculateAdjoint := function(mat, mode, hom, n)
+        g_star_inv_unscaled := TransposedMat(Inverse(g_res[1]));
+        #* g_res[2]; 
+        # Todo the computation of this can probably be sped up by using the known information about g!!!
+        g_star_inv_scaled_frob := FrobeniusNormalForm(g_star_inv_unscaled); # hier ist eventuell noch ein bug mit den skalaren
+        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.
+#! The reason the the field must be given as an argument, is so that the Field automorphims of order two, which is used to compute the adjoined, is specified.
+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);
diff --git a/lib/recognition_new.gi b/lib/recognition_new.gi
index 5fb035e..f209137 100644
--- a/lib/recognition_new.gi
+++ b/lib/recognition_new.gi
@@ -355,7 +355,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/formspace/custom_test_functions.g b/tst/formspace/custom_test_functions.g
new file mode 100644
index 0000000..e15acf0
--- /dev/null
+++ b/tst/formspace/custom_test_functions.g
@@ -0,0 +1,148 @@
+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;
+
+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/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/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..8af73e1
--- /dev/null
+++ b/tst/interesting_groups.g
@@ -0,0 +1,373 @@
+## 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
+
+# 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,