diff --git a/gap/nofoma.gd b/gap/nofoma.gd
index ef092ea..8429ef0 100644
--- a/gap/nofoma.gd
+++ b/gap/nofoma.gd
@@ -344,13 +344,14 @@ DeclareGlobalFunction("JordanChevalleyDecMatF");
 #! @Section The primary decomposition
 #! @Arguments A
 #! @Description
-#!  Returns a list containing two elements. The first element is
+#!  Returns a list containing three elements. The first element is
 #!  a base change matrix <M>B</M> such that <M>B</M><A>A</A><M>B^{-1}</M> is a
 #!  primary form of <A>A</A>, i.e., a block diagonal matrix where the minimal polynomials
 #!  of the the diagonal blocks are precisely the powers of irreducible factors
-#!  of the minimal polynomial of <A>A</A>. The second element is the size of each 
-#!  block. 
-#!  This function uses a modified version of Steel's algorithm.
+#!  of the minimal polynomial of <A>A</A>, in descending order. The second element is a list containing
+#!  the collected irreducible factors of the minimal polynomial of <A>A</A>, in the same order. The
+#!  last element is a list containing the the size of each block. 
+#!  The exact algorithm used in this function is described in <Cite Key ="Bon26"/>
 #! 
 #! @BeginExampleSession
 #! gap> A := [ [ Z(5)^2, 0*Z(5), Z(5)^2, Z(5)^3, Z(5) ], 
@@ -358,8 +359,8 @@ DeclareGlobalFunction("JordanChevalleyDecMatF");
 #! >    [ Z(5), Z(5)^0, 0*Z(5), Z(5)^0, 0*Z(5) ],
 #! >    [ Z(5)^0, Z(5)^0, Z(5)^0, 0*Z(5), Z(5)^3 ],
 #! >    [ Z(5), 0*Z(5), Z(5)^3, 0*Z(5), Z(5)^3 ] ];;
-#! gap> B := PrimaryDecomp(A);;
-#! gap> Display(B[2]);
+#! gap> B := PrimaryDecomposition(A);;
+#! gap> Display(B[3]);
 #! [ 1, 4 ]
 #! gap> Factors(MinimalPolynomial(A));
 #! [ x_1-Z(5)^0, x_1^4-x_1^3+Z(5)^3*x_1+Z(5)^3 ]
@@ -369,7 +370,7 @@ DeclareGlobalFunction("JordanChevalleyDecMatF");
 #! gap> MinimalPolynomial(PrimA{[2..5]}{[2..5]});
 #! x_1^4-x_1^3+Z(5)^3*x_1+Z(5)^3
 #! @EndExampleSession
-DeclareGlobalFunction("PrimaryDecomp");
+DeclareGlobalFunction("PrimaryDecomposition");
 
 #! @Chapter Auxiliary functions
 #! @Section Vectors and matrices and their associated polynomials
diff --git a/gap/nofoma.gi b/gap/nofoma.gi
index 341e937..32b1e7e 100644
--- a/gap/nofoma.gi
+++ b/gap/nofoma.gi
@@ -735,21 +735,23 @@ end);
 #Primary Decomposition for cyclic matrices
 #Jordan normal form will call this function if a cyclic matrix is detected
 BindGlobal("nfmPrimaryDecompositionforJNFCyclic", function(A, minpol, minpolfacs)
-    local vspan,n,w,mf,wspan,qi,k,COB,dims;
+    local vspan,n,w,mf,wspan,qi,k,COB,dims,elDivs;
     n := NrRows(A);
     vspan := nfmFindCyclicVectorNC(A);
     dims := [];
     COB := ZeroMutable(A);
     k := 0;
+    elDivs := [];
     for mf in minpolfacs do
         qi := (mf[1])^(mf[2]);
+        Add(elDivs, qi);
         w := nfmPolyEvalFromSpan(vspan,Quotient(minpol,qi));
         wspan := nfmSpinUntil(w,A,Degree(mf[1])*mf[2]);
         CopySubMatrix(wspan, COB, [1..NrRows(wspan)], [k+1..k+NrRows(wspan)], [1..n], [1..n]);
         Add(dims, NrRows(wspan));
         k := k + NrRows(wspan);
     od;
-    return [COB, dims];
+    return [COB, elDivs, dims];
 end);
 
 #TODO: recognise cyclic matrices
@@ -757,16 +759,17 @@ end);
 #Standalone version
 #Returns matrix B such that B*A*B^-1 is in primary decomposition form
 #along with dimensions of primary subspaces
-InstallGlobalFunction(PrimaryDecomp, function(A)
+InstallGlobalFunction(PrimaryDecomposition, function(A)
     local rank,F,n,m,f,w,p,j,i,wspan,gens,facs,L_i,qi,k,v,
-    COB,pot,gs,f2,dims,toAdd,dim,minpol;
+    COB,pot,gs,f2,dims,toAdd,dim,minpol,collected;
     rank := 0;
     n := NrRows(A);
     F := DefaultFieldOfMatrix(A);
     v := nfmGenerateNonZeroVector(n,A);
     minpol := MinimalPolynomial(F,A);
+    collected := Collected(Factors(minpol));
     if Degree(minpol) = n then
-      return nfmPrimaryDecompositionforJNFCyclic(A,minpol,Collected(Factors(minpol)));
+      return nfmPrimaryDecompositionforJNFCyclic(A,minpol,collected);
     fi;
     if IsZero(A) then
       return IdentityMat(n,F);
@@ -838,7 +841,7 @@ InstallGlobalFunction(PrimaryDecomp, function(A)
       CopySubMatrix(L_i, COB, [1..dim], [k+1..k+dim],[1..n], [1..n]);
       k := k + dim;
     od;
-    return [COB, dims];
+    return [COB, collected, dims];
 end);
 
 #Input: For matrix A with minimal polynomial p^m: m, vector v and p(A)
@@ -999,7 +1002,7 @@ end);
 #Input: Matrix A
 #Returns matrix B such that A^Inverse(B) is in Jordan normal form
 InstallGlobalFunction(JordanNormalForm, function(A)
-    local n,F,pol,minpol,primary,primarydims,crhr,cyclicdims,elDivs,
+    local n,F,pol,minpol,primary,primarydims,crhr,cyclicdims,elDivs,prim1,
     subsubCOB,COB,subA,cy,i,j,facOcc,subCOB,crcy,subsubA,prepreCOB,preCOB;
     F := DefaultFieldOfMatrix(A);
     n := NrRows(A);
@@ -1013,7 +1016,8 @@ InstallGlobalFunction(JordanNormalForm, function(A)
         return JordanNormalFormIrred(A,minpol);
     fi;
     if Degree(minpol) = n then
-        primary := nfmPrimaryDecompositionforJNFCyclic(A, minpol, facOcc);
+        prim1 := nfmPrimaryDecompositionforJNFCyclic(A, minpol, facOcc);
+        primary := [prim1[1],prim1[3]];
     else
         primary := nfmPrimaryDecompositionforJNF(A, minpol, facOcc);
     fi;
diff --git a/tst/primary.tst b/tst/primary.tst
index 201fd7a..a54a8dd 100644
--- a/tst/primary.tst
+++ b/tst/primary.tst
@@ -1,18 +1,18 @@
-gap> START_TEST("PrimaryDecomposition");
+gap> START_TEST("PrimaryDecompositionosition");
 gap> ReadPackage("nofoma", "tst/utils.g");
 true
 
 ##Test Primary decomposition
-gap> nfmCheckPrimaryDecomp(GF(5),50);
+gap> nfmCheckPrimaryDecomposition(GF(5),50);
 true
-gap> nfmCheckPrimaryDecomp(GF(25),25);
+gap> nfmCheckPrimaryDecomposition(GF(25),25);
 true
-gap> nfmCheckPrimaryDecomp(GF(7),150);
+gap> nfmCheckPrimaryDecomposition(GF(7),150);
 true
-gap> nfmCheckPrimaryDecompNonCyclic(GF(5),50);
+gap> nfmCheckPrimaryDecompositionNonCyclic(GF(5),50);
 true
-gap> nfmCheckPrimaryDecompNonCyclic(GF(5^5),10);
+gap> nfmCheckPrimaryDecompositionNonCyclic(GF(5^5),10);
 true
 
 ## Stop test
-gap> STOP_TEST("PrimaryDecomposition");
+gap> STOP_TEST("PrimaryDecompositionosition");
diff --git a/tst/representations.tst b/tst/representations.tst
index 6e1165d..21d8f95 100644
--- a/tst/representations.tst
+++ b/tst/representations.tst
@@ -58,19 +58,19 @@ gap> CheckJordanChev(bigm, JordanChevalleyDecMat(bigm, MinimalPolynomial(bigm)))
 [ Z(5)^0, x_1 ]
 
 ## Primary decomposition should work for all these matrix families
-gap> nfmCheckPrimaryDecompForMatrix(rat);
+gap> nfmCheckPrimaryDecompositionForMatrix(rat);
 true
-gap> nfmCheckPrimaryDecompForMatrix(ratm);
+gap> nfmCheckPrimaryDecompositionForMatrix(ratm);
 true
-gap> nfmCheckPrimaryDecompForMatrix(smallplain);
+gap> nfmCheckPrimaryDecompositionForMatrix(smallplain);
 true
-gap> nfmCheckPrimaryDecompForMatrix(gf2);
+gap> nfmCheckPrimaryDecompositionForMatrix(gf2);
 true
-gap> nfmCheckPrimaryDecompForMatrix(gf9);
+gap> nfmCheckPrimaryDecompositionForMatrix(gf9);
 true
-gap> nfmCheckPrimaryDecompForMatrix(bigplain);
+gap> nfmCheckPrimaryDecompositionForMatrix(bigplain);
 true
-gap> nfmCheckPrimaryDecompForMatrix(bigm);
+gap> nfmCheckPrimaryDecompositionForMatrix(bigm);
 true
 
 ## Jordan normal form should preserve the input matrix family
diff --git a/tst/utils.g b/tst/utils.g
index 944a957..39f72e1 100644
--- a/tst/utils.g
+++ b/tst/utils.g
@@ -57,14 +57,14 @@ end;
 #TODO: now that primdecomp has sorted blocks this can be tested more efficiently
 #check primary decomp function with field F and dimension n
 #checks if the submatrices have the correct minimal polynomials
-nfmCheckPrimaryDecomp := function(F, n)
+nfmCheckPrimaryDecomposition := function(F, n)
   local A,Prim,B,dim,k,minpolfacs,sub;
   A := Matrix(F,RandomInvertibleMat(n,F));
-  Prim := PrimaryDecomp(A);
+  Prim := PrimaryDecomposition(A);
   B := A^Inverse(Prim[1]);
   minpolfacs := Factors(MinimalPolynomial(F,A));
   k := 1;
-  for dim in Prim[2] do
+  for dim in Prim[3] do
     sub := ExtractSubMatrix(B,[k..k+dim-1],[k..k+dim-1]);
     k := k+dim;
     if not Factors(MinimalPolynomial(F,sub))[1] in minpolfacs then
@@ -74,14 +74,14 @@ nfmCheckPrimaryDecomp := function(F, n)
   return true;
 end;
 
-nfmCheckPrimaryDecompNonCyclic := function(F, n)
+nfmCheckPrimaryDecompositionNonCyclic := function(F, n)
   local A,Prim,B,dim,k,minpolfacs,sub;
   A := Matrix(F,nfmGenerateNonCyclicMatrix(F,n));
-  Prim := PrimaryDecomp(A);
+  Prim := PrimaryDecomposition(A);
   B := A^Inverse(Prim[1]);
   minpolfacs := Factors(MinimalPolynomial(F,A));
   k := 1;
-  for dim in Prim[2] do
+  for dim in Prim[3] do
     sub := ExtractSubMatrix(B,[k..k+dim-1],[k..k+dim-1]);
     k := k+dim;
     if not Factors(MinimalPolynomial(F,sub))[1] in minpolfacs then
@@ -108,14 +108,14 @@ nfmCheckInvariantFactorsForMatrix := function(A)
   return InvariantFactorsMat(A) = FrobeniusNormalForm(A)[1];
 end;
 
-nfmCheckPrimaryDecompForMatrix := function(A)
+nfmCheckPrimaryDecompositionForMatrix := function(A)
   local Prim,B,dim,k,minpolfacs,sub,AA;
-  Prim := PrimaryDecomp(A);
+  Prim := PrimaryDecomposition(A);
   AA := Matrix(A, Prim[1]);
   B := AA^Inverse(Prim[1]);
   minpolfacs := Factors(MinimalPolynomial(AA));
   k := 1;
-  for dim in Prim[2] do
+  for dim in Prim[3] do
     sub := ExtractSubMatrix(B,[k..k+dim-1],[k..k+dim-1]);
     k := k+dim;
     if not Factors(MinimalPolynomial(sub))[1] in minpolfacs then