PrimaryDecomp recognises cyclic matrices (#62)

* PrimaryDecomp recognises cyclic matrices

Author: AliaBonnet

gap/nofoma.gi

@@ -627,32 +627,52 @@ BindGlobal("nfmPolyEvalFromSpan", function(span,pol)
     return resu;
 end);
 
-#TODO: recognise cyclic matrices 
-#Primary Decomposition using a modified version of Allan Steel's algorithm 
-#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) 
-    local rank,F,n,m,f,w,p,j,i,wspan,gens,facs,L_i,qi,k,v,
-    COB,pot,gs,f2,dims,toAdd,combined,dim;
+
+BindGlobal("factoriseByKnownFactors", function(kFacs,pol)
+  local fac,factup,i,resfacs,count,newpol,oldpol; 
+  resfacs := [];
+  oldpol := pol;
+  for factup in kFacs do 
+    fac := factup[1]; 
+    count := 0;
+    for i in [1..factup[2]] do 
+      newpol := Quotient(oldpol, fac);
+      if newpol = fail then
+        break;
+      fi;
+      count := count + 1;
+      oldpol := newpol;
+    od;
+    if not count = 0 then
+      Add(resfacs, [fac, count]);
+    fi;
+  od;
+  return resfacs;
+end);
+
+#TODO: if size of minpolfacs is one we can return identity 
+#Primary Decomposition using a modified version of Allan Steel's algorithm for use in the Jordan normal form function
+#Takes matrix along with its minimal polynomial and its factorised version as input
+#Returns matrix B such that B*A*B^-1 is in primary decomp. form, dimensions of primary subspaces
+#and factors of minimal polynomial in correct order
+BindGlobal("nfmPrimaryDecompositionforJNF", function(A, minpol, minpolfacs)
+    local rank,F,n,m,f,w,p,j,i,wspan,gens,facs,L_i,qi,k,v, 
+    COB,pot,gs,f2,toAdd,pos,dims,dim;
     rank := 0;
     n := NrRows(A);
     F := DefaultFieldOfMatrix(A);
     v := nfmGenerateRandomVector(F,n);
-    A := Matrix(F,A);
-    if IsZero(A) then 
-      return OneMutable(A);
-    fi;
-    gens := []; #Li_s as in Steel paper will go in here 
+    gens := []; #Li_s as in Steel paper will go in here (but without separate U)
     gs := []; #distinct factors of minimal polynomial 
+    dims := [];
     while not rank = n do 
         m := UnivariatePolynomial(F,SpinMatVector(A,v)[3]);
         p := One(PolynomialRing(F));
         for i in [1..Size(gens)] do 
             L_i := gens[i];
             if not IsOne(Gcd(m,gs[i])) then
                 f := m;
-                pot := 0; 
+                pot := 0;
                 for j in [1..n] do 
                     f2 := Quotient(f,gs[i]);
                     if f2 = fail then 
@@ -663,8 +683,8 @@ InstallGlobalFunction(PrimaryDecomp, function(A)
                 od; 
                 w := PolynomialToMatVec(A,f,v);
                 p := p * gs[i]^pot;
-                #minimal polynomial of w has degree smaller than or equal to n - degree(f)
-                wspan := nfmSpinUntil(w,A,n-Degree(f));
+                #minimal polynomial of w has degree smaller than or equal to minpol/f
+                wspan := nfmSpinUntil(w,A,Degree(minpol)-Degree(f));
                 toAdd := EcheloniseMat(Concatenation(wspan,gens[i]));
                 if not IsMatrix(toAdd) then 
                     toAdd := [toAdd];  # Convert vector to 1-row matrix
@@ -676,8 +696,9 @@ InstallGlobalFunction(PrimaryDecomp, function(A)
         v :=PolynomialToMatVec(A,p,v);
         m := Quotient(m,p);
         if not IsOne(m) then 
-            facs := Factors(m);
-            facs := Collected(facs);
+            #TODO: make this work
+            #facs := factoriseByKnownFactors(minpolfacs,m);
+            facs := Collected(Factors(m));
             for i in [1..Size(facs)] do 
                 qi := (facs[i][1])^(facs[i][2]);
                 w := PolynomialToMatVec(A,Quotient(m,qi),v);
@@ -686,7 +707,7 @@ InstallGlobalFunction(PrimaryDecomp, function(A)
                 Add(gs, facs[i][1]);
             od;
         fi;
-        #alles in eine matrix
+        # put everything into a single matrix
         COB := ZeroMatrix(F,n,n);
         k := 0;
         for i in [1..Size(gens)] do 
@@ -701,70 +722,70 @@ InstallGlobalFunction(PrimaryDecomp, function(A)
             v := nfmFindVectorNotInSubspaceNC(COB);
         fi;
     od;
-    #sorted blocks and count dims
     COB := ZeroMatrix(F,n,n);
-    combined := [];
-    for i in [1..Size(gs)] do 
-      Add(combined, [gens[i],gs[i]]); 
-    od;
-    Sort(combined, {v,w} -> v[2] < w[2]);
+    #sort primary subspaces back into order and collect dims
     k := 0;
-    dims := [];
-    for i in [1..Size(gs)] do 
-      L_i := combined[i][1];
-      dim := NrRows(L_i); 
+    for i in [1..Size(gens)] do 
+      pos := Position(gs, minpolfacs[i][1]);
+      L_i := gens[pos];
+      dim := NrRows(L_i); #?
       Add(dims, dim);
       CopySubMatrix(L_i, COB, [1..dim], [k+1..k+dim],[1..n], [1..n]);
       k := k + dim;
     od;
     return [COB, dims];
 end);
 
-BindGlobal("factoriseByKnownFactors", function(kFacs,pol)
-  local fac,factup,i,resfacs,count,newpol,oldpol; 
-  resfacs := [];
-  oldpol := pol;
-  for factup in kFacs do 
-    fac := factup[1]; 
-    count := 0;
-    for i in [1..factup[2]] do 
-      newpol := Quotient(oldpol, fac);
-      if newpol = fail then
-        break;
-      fi;
-      count := count + 1;
-      oldpol := newpol;
+#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;
+    n := NrRows(A);
+    vspan := nfmFindCyclicVectorNC(A);
+    dims := [];
+    COB := ZeroMutable(A);
+    k := 0;
+    for mf in minpolfacs do 
+        qi := (mf[1])^(mf[2]);
+        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;
-    if not count = 0 then
-      Add(resfacs, [fac, count]);
-    fi;
-  od;
-  return resfacs;
+    return [COB, dims];
 end);
 
-#TODO: if size of minpolfacs is one we can return identity 
-#Primary Decomposition using a modified version of Allan Steel's algorithm for use in the Jordan normal form function
-#Takes matrix along with its minimal polynomial and its factorised version as input
-#Returns matrix B such that B*A*B^-1 is in primary decomp. form, dimensions of primary subspaces
-#and factors of minimal polynomial in correct order
-BindGlobal("nfmPrimaryDecompositionforJNF", function(A, minpol, minpolfacs)
-    local rank,F,n,m,f,w,p,j,i,wspan,gens,facs,L_i,qi,k,v, 
-    COB,pot,gs,f2,toAdd,pos,dims,dim;
+#TODO: recognise cyclic matrices 
+#Primary Decomposition using a modified version of Allan Steel's algorithm 
+#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) 
+    local rank,F,n,m,f,w,p,j,i,wspan,gens,facs,L_i,qi,k,v,
+    COB,pot,gs,f2,dims,toAdd,combined,dim,minpol;
     rank := 0;
     n := NrRows(A);
     F := DefaultFieldOfMatrix(A);
     v := nfmGenerateRandomVector(F,n);
-    gens := []; #Li_s as in Steel paper will go in here (but without separate U)
+    A := Matrix(F,A);
+    minpol := MinimalPolynomial(F,A);
+    if Degree(minpol) = n then 
+      return nfmPrimaryDecompositionforJNFCyclic(A,minpol,Collected(Factors(minpol)));
+    fi;
+    if IsZero(A) then 
+      return IdentityMat(n,F);
+    fi;
+    gens := []; #Li_s as in Steel paper will go in here 
     gs := []; #distinct factors of minimal polynomial 
-    dims := [];
     while not rank = n do 
         m := UnivariatePolynomial(F,SpinMatVector(A,v)[3]);
         p := One(PolynomialRing(F));
         for i in [1..Size(gens)] do 
             L_i := gens[i];
             if not IsOne(Gcd(m,gs[i])) then
                 f := m;
-                pot := 0;
+                pot := 0; 
                 for j in [1..n] do 
                     f2 := Quotient(f,gs[i]);
                     if f2 = fail then 
@@ -775,8 +796,8 @@ BindGlobal("nfmPrimaryDecompositionforJNF", function(A, minpol, minpolfacs)
                 od; 
                 w := PolynomialToMatVec(A,f,v);
                 p := p * gs[i]^pot;
-                #minimal polynomial of w has degree smaller than or equal to minpol/f
-                wspan := nfmSpinUntil(w,A,Degree(minpol)-Degree(f));
+                #minimal polynomial of w has degree smaller than or equal to n - degree(f)
+                wspan := nfmSpinUntil(w,A,n-Degree(f));
                 toAdd := EcheloniseMat(Concatenation(wspan,gens[i]));
                 if not IsMatrix(toAdd) then 
                     toAdd := [toAdd];  # Convert vector to 1-row matrix
@@ -788,9 +809,8 @@ BindGlobal("nfmPrimaryDecompositionforJNF", function(A, minpol, minpolfacs)
         v :=PolynomialToMatVec(A,p,v);
         m := Quotient(m,p);
         if not IsOne(m) then 
-            #TODO: make this work
-            #facs := factoriseByKnownFactors(minpolfacs,m);
-            facs := Collected(Factors(m));
+            facs := Factors(m);
+            facs := Collected(facs);
             for i in [1..Size(facs)] do 
                 qi := (facs[i][1])^(facs[i][2]);
                 w := PolynomialToMatVec(A,Quotient(m,qi),v);
@@ -799,7 +819,7 @@ BindGlobal("nfmPrimaryDecompositionforJNF", function(A, minpol, minpolfacs)
                 Add(gs, facs[i][1]);
             od;
         fi;
-        # put everything into a single matrix
+        #alles in eine matrix
         COB := ZeroMatrix(F,n,n);
         k := 0;
         for i in [1..Size(gens)] do 
@@ -814,40 +834,25 @@ BindGlobal("nfmPrimaryDecompositionforJNF", function(A, minpol, minpolfacs)
             v := nfmFindVectorNotInSubspaceNC(COB);
         fi;
     od;
+    #sorted blocks and count dims
     COB := ZeroMatrix(F,n,n);
-    #sort primary subspaces back into order and collect dims
+    combined := [];
+    for i in [1..Size(gs)] do 
+      Add(combined, [gens[i],gs[i]]); 
+    od;
+    Sort(combined, function(v,w) return v[2] < w[2]; end);
     k := 0;
-    for i in [1..Size(gens)] do 
-      pos := Position(gs, minpolfacs[i][1]);
-      L_i := gens[pos];
-      dim := NrRows(L_i); #?
+    dims := [];
+    for i in [1..Size(gs)] do 
+      L_i := combined[i][1];
+      dim := NrRows(L_i); 
       Add(dims, dim);
       CopySubMatrix(L_i, COB, [1..dim], [k+1..k+dim],[1..n], [1..n]);
       k := k + dim;
     od;
     return [COB, dims];
 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;
-    n := NrRows(A);
-    vspan := nfmFindCyclicVectorNC(A);
-    dims := [];
-    COB := ZeroMutable(A);
-    k := 0;
-    for mf in minpolfacs do 
-        qi := (mf[1])^(mf[2]);
-        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];
-end);
-
 #Input: For matrix A with minimal polynomial p^m: m, vector v and p(A)
 #Returns: v, length r of v, vp^(r-1)(A)
 BindGlobal("GetMinPolPowerWithVec", function(m,v,Ainp)