ENHANCE: Faster MaximalAbelianQuotient for subgroups of fp groups

lib/ghomfp.gi

@@ -800,7 +800,7 @@ local aug,w,p,pres,f,fam,opt;
   fi;
 
   TzOptions(pres).printLevel:=InfoLevel(InfoFpGroup); 
-  if ValueOption("quick")=true then
+  if ValueOption("cheap")=true then
     TzGo(pres);
   else
     TzEliminateRareOcurrences(pres,50);
@@ -1079,6 +1079,7 @@ local m,s,g,i,j,rel,gen,img,fin,hom,gens;
   od;
 
   if Length(m)>0 then  
+    m:=ReducedRelationMat(m);
     s:=NormalFormIntMat(m,25); # 9+16: SNF with transforms, destructive
     SetAbelianInvariants(f,AbelianInvariantsOfList(DiagonalOfMat(s.normal)));
 
@@ -1089,11 +1090,11 @@ local m,s,g,i,j,rel,gen,img,fin,hom,gens;
         Add(rel,g.(i)^s.normal[i][i]);
       od;
       Append(img,ListWithIdenticalEntries(Length(gen)-s.rank,0));
-      SetAbelianInvariants(f,img);
       g:=g/rel;
       fin:=false;
     else  
-      SetAbelianInvariants(f,DiagonalOfMat(s.normal));
+      # Not `AbelianInvariantsOfList' as the structure of the group is as
+      # given by the normal form
       g:=AbelianGroup(DiagonalOfMat(s.normal));
       fin:=true;
     fi;
@@ -1124,30 +1125,131 @@ local m,s,g,i,j,rel,gen,img,fin,hom,gens;
 end);
 
 InstallMethod(MaximalAbelianQuotient,
-        "for subgroups of finitely presented groups",
-        true, [IsSubgroupFpGroup], 0,
+        "for subgroups of finitely presented groups, fallback",
+        true, [IsSubgroupFpGroup], -1,
 function(U)
 local phi,m;
   # do cheaper Tietze (and thus do not store)
   phi:=AttributeValueNotSet(IsomorphismFpGroup,U:
     eliminationsLimit:=50,
     generatorsLimit:=Length(GeneratorsOfGroup(Parent(U)))*LogInt(IndexInWholeGroup(U),2),
-    quick); 
+    cheap); 
   m:=MaximalAbelianQuotient(Image(phi));
   SetAbelianInvariants(U,AbelianInvariants(Image(phi)));
   return phi*MaximalAbelianQuotient(Image(phi));
 end);
 
+InstallMethod(MaximalAbelianQuotient,
+        "subgroups of fp. abelian rewriting", true, [IsSubgroupFpGroup], 0,
+function(u)
+local aug,r,sec,t,expwrd,rels,ab,s,m,img,gen,i,j,t1,t2,tn;
+  if (HasIsWholeFamily(u) and IsWholeFamily(u))
+  # catch trivial case of rank 0 group
+   or Length(GeneratorsOfGroup(FamilyObj(u)!.wholeGroup))=0 then
+    TryNextMethod();
+  fi;
+
+  # get an augmented coset table from the group. Since we don't care about
+  # any particular generating set, we let the function chose.
+  aug:=AugmentedCosetTableInWholeGroup(u);
+
+  aug:=CopiedAugmentedCosetTable(aug);
+
+  r:=Length(aug.primaryGeneratorWords);
+  Info( InfoFpGroup, 1, "Abelian presentation with ",
+    Length(aug.subgroupGenerators), " generators");
+
+  # make vectors 
+  expwrd:=function(l)
+  local v,i;
+    v:=ListWithIdenticalEntries(r,0);
+    for i in l do
+      if i>0 then v:=v+sec[i];
+      else v:=v-sec[-i];fi;
+    od;
+    return v;
+  end;
+
+#  sec:=ShallowCopy(IdentityMat(r,1)); # initialize so next command works
+#
+#  sec:=List(GeneratorTranslationAugmentedCosetTable(aug),
+#    x->expwrd(LetterRepAssocWord(x)));
+#
+#  m:=sec;
+  # do GeneratorTranslation abelianized
+  sec:=ShallowCopy(IdentityMat(r,1)); # initialize so next command works
+
+  t1:=aug.tree[1];
+  t2:=aug.tree[2];
+  tn:=aug.treeNumbers;
+  if Length(tn)>0 then
+    for i in [Length(sec)+1..Maximum(tn)] do
+      sec[i]:=sec[AbsInt(t1[i])]*SignInt(t1[i])
+            +sec[AbsInt(t2[i])]*SignInt(t2[i]);
+    od;
+  fi;
+#  if sec<>m then Error("ZZZ");fi;
+
+  sec:=sec{aug.treeNumbers};
+
+  # now make relators abelian
+  rels:=[];
+  rels:=RewriteSubgroupRelators( aug, aug.groupRelators);
+  rels:=List(rels,expwrd);
+
+  if Length(rels)=0 then
+    Add(rels,ListWithIdenticalEntries(r,0));
+  fi;
+  rels:=ReducedRelationMat(rels);
+  s:=NormalFormIntMat(rels,25); # 9+16: SNF with transforms, destructive
+  SetAbelianInvariants(u,AbelianInvariantsOfList(DiagonalOfMat(s.normal)));
+  if r>s.rank then
+    # TODO: Reproduce creation of infinite abelian group
+    TryNextMethod();
+  else
+      # Not `AbelianInvariantsOfList' as the structure of the group is as
+      # given by the normal form
+    ab:=AbelianGroup(DiagonalOfMat(s.normal));
+  fi;
+  gen:=GeneratorsOfGroup(ab);
+  s:=s.coltrans;
+  img:=[];
+  for i in [1..Length(s)] do
+    m:=One(ab);
+    for j in [1..Length(gen)] do
+      m:=m*gen[j]^s[i][j];
+    od;
+    Add(img,m);
+  od;
+  aug.primaryImages:=img;
+  sec:=List(sec,x->LinearCombinationPcgs(img,x));
+  aug.secondaryImages:=sec;
+
+  m:=List(aug.primaryGeneratorWords,x->ElementOfFpGroup(FamilyObj(One(u)),x));
+  m:=GroupHomomorphismByImagesNC(u,ab,m,img:noassert);
+
+  # but give it `aug' as coset table, so we will use rewriting for images
+  SetCosetTableFpHom(m,aug);
+
+  SetIsSurjective(m,true);
+
+  return m;
+end);
+
 # u must be a subgroup of the image of home
 InstallGlobalFunction(
 LargerQuotientBySubgroupAbelianization,function(hom,u)
-local G,v,ma,mau,a,gens,imgs,q,k,co,aiu,aiv,primes,irrel;
+local G,v,ma,mau,a,gens,imgs,q,k,co,aiu,aiv,primes,irrel,pcgs,ind,i,tst;
   v:=PreImage(hom,u);
   aiu:=AbelianInvariants(u);
 
   G:= FamilyObj(v)!.wholeGroup;
-  aiv:=AbelianInvariantsSubgroupFpGroup( G, v );
-  if aiv=fail or aiu=aiv then
+  aiv:=AbelianInvariantsSubgroupFpGroup( G, v:cheap:=false );
+  if aiv=fail then
+    ma:=MaximalAbelianQuotient(v);
+    aiv:=AbelianInvariants(Image(ma,v));
+  fi;
+  if aiu=aiv then
     return fail;
   fi;
   # are there irrelevant primes?
@@ -1177,14 +1279,18 @@ local G,v,ma,mau,a,gens,imgs,q,k,co,aiu,aiv,primes,irrel;
   imgs:=List(gens,x->Image(mau,Image(hom,PreImagesRepresentative(ma,x))));
   q:=GroupHomomorphismByImages(a,Image(mau),gens,imgs);
   k:=KernelOfMultiplicativeGeneralMapping(q);
-  co:=ComplementClassesRepresentatives(a,k);
-  if Length(co)=0 then
-    co:=List(ConjugacyClassesSubgroups(a),Representative);
-    co:=Filtered(co,x->Size(Intersection(k,x))=1);
-    Sort(co,function(a,b) return Size(a)>Size(b);end);
-  fi;
-  Info(InfoFpGroup,2,"Degree larger ",Index(a,co[1]),"\n");
-  return PreImage(ma,co[1]);
+  # try to use indices
+  pcgs:=Pcgs(a);
+  ind:=List(InducedPcgs(pcgs,k),x->DepthOfPcElement(pcgs,x));
+  co:=TrivialSubgroup(a);
+  for i in Reversed(Difference([1..Length(pcgs)],ind)) do
+    tst:=ClosureSubgroup(co,pcgs[i]);
+    if Size(Intersection(tst,k))=1 then
+      co:=tst;
+    fi;
+  od;
+  Info(InfoFpGroup,2,"Degree larger ",Index(a,co));
+  return PreImage(ma,co);
 end);
 
 DeclareRepresentation("IsModuloPcgsFpGroupRep",

lib/matint.gd

@@ -559,3 +559,24 @@ DeclareGlobalFunction("DeterminantIntMat");
 ##  </ManSection>
 ##
 DeclareGlobalFunction("SNFofREF");
+
+#############################################################################
+##
+#O  ReducedRelationMat(<mat>)
+##
+##  <#GAPDoc Label="ReducedRelationMat">
+##  <ManSection>
+##  <Func Name="ReducedRelationMat" Arg='mat'/>
+##
+##  <Description>
+##  Let <A>mat</A> be a matrix that has been obtained as abelianized
+##  relations. Such matrices tend to have a particular form with some short
+##  vectors. This function runs a (quick) heuristic row reduction, 
+##  resulting in a matrix with the same Z-row space but fewer/shorter vectors,
+##  thus speeding up a subsequent SNF. It does not do a full HNF but should be
+##  much quicker.
+##  </Description>
+##  </ManSection>
+##  <#/GAPDoc>
+##
+DeclareGlobalFunction("ReducedRelationMat");

lib/matint.gi

@@ -1000,3 +1000,133 @@ InstallOtherMethod( AbelianInvariantsOfList,
     [ IsList and IsEmpty ],
     list -> [] );
 
+
+# Reduce a list of abelianized relations: Heuristic reduction without
+# making big vectors, iterate three times. Does not aim to do full HNF
+InstallGlobalFunction(ReducedRelationMat,function(mat)
+local n,zero,nv,new,pip,piv,i,v,p,w,g,nov,pin,now,rat,extra,clean,assign,try;
+
+  nv:=v->v*SignInt(v[PositionNonZero(v)]);
+  assign:=function(p,v)
+  local a,i,w,wn;
+    a:=v[p];
+    for i in [1..Length(pip)] do
+      if i<>p and IsInt(pip[i]) and mat[pip[i]][p]<>0 then
+        w:=mat[pip[i]]-QuoInt(mat[pip[i]][p],a)*v;
+        wn:=w*w;
+        if wn<=rat*pin[i] then
+          mat[pip[i]]:=nv(w);
+          pin[i]:=wn;
+        fi;
+      fi;
+    od;
+    mat[pip[p]]:=v;
+    # also try to reduce extra vectors
+    for i in [1..Length(extra)] do
+      w:=extra[i];
+      if not IsZero(extra[i]) then
+        wn:=w*w;
+        w:=w-QuoInt(w[p],a)*v;
+        if w*w<=rat*wn then
+          extra[i]:=w;
+        fi;
+      fi;
+    od;
+  end;
+
+  n:=NrCols(mat);
+  rat:=2; # growth ratio
+  zero:=ListWithIdenticalEntries(n,0);
+  mat:=Filtered(mat,x->not IsZero(x));
+  new:=Set(mat,nv); # kill duplicates
+  Info(InfoMatInt,1,"Reduce ",Length(mat)," to ",Length(new));
+  pip:=ListWithIdenticalEntries(n,fail);
+  piv:=[];
+  pin:=[];
+  mat:=[];
+  extra:=[];
+
+  # we once reduce and then go over the remainders again in case they were
+  # nice and short
+  for try in [1..3] do
+    SortBy(new, x -> - x*x); # reversed norm sort
+    i:=Length(new);
+    while i>0 do
+      v:=ShallowCopy(new[i]);
+      Info(InfoMatInt,3,"Process ",i);#" Norm:",v*v,"\n");
+      Unbind(new[i]); # take off stack
+      i:=i-1;
+      clean:=true;
+      p:=PositionNonZero(v);
+      while p<=n and pip[p]<>fail do
+        if v[p] mod piv[p]=0 then
+          # divides, reduce
+          #v:=v-QuoInt(v[p],piv[p])*mat[pip[p]];
+          AddRowVector(v,mat[pip[p]],-QuoInt(v[p],piv[p]));
+          p:=PositionNonZero(v,p);
+        elif clean and piv[p] mod v[p]=0 then
+          # swap and clean out
+          v:=nv(v);
+          Info(InfoMatInt,2,"Replace pivot ",piv[p],"@",p," to ",v[p]);
+          w:=mat[pip[p]];
+          #mat[pip[p]]:=v;
+          assign(p,v);
+          pin[p]:=v*v;
+          piv[p]:=v[p];
+          v:=w;
+          #v:=v-QuoInt(v[p],piv[p])*mat[pip[p]];
+          AddRowVector(v,mat[pip[p]],-QuoInt(v[p],piv[p]));
+          p:=PositionNonZero(v,p);
+        else
+          g:=Gcdex(v[p],piv[p]);
+          # form new pivot with gcd 
+          #w:=g.coeff2*mat[pip[p]]+g.coeff1*v; # automatically normed by Gcdex
+          w:=g.coeff2*mat[pip[p]];
+          AddRowVector(w,v,g.coeff1); # automatically normed by Gcdex
+          now:=w*w;
+          if (not clean) or now>rat*pin[p] then 
+            # only reduce a bit, not full gcd
+            #v:=v-QuoInt(v[p],piv[p])*mat[pip[p]];
+            AddRowVector(v,mat[pip[p]],-QuoInt(v[p],piv[p]));
+            p:=PositionNonZero(v,p);
+            clean:=false;
+          else
+            # replace with cgd pivot
+            Info(InfoMatInt,2,"Reduce pivot ",piv[p],"@",p," to ",g.gcd);
+            new[i+1]:=v; # keep old vectors to process
+            new[i+2]:=mat[pip[p]];
+            i:=i+2;
+            #mat[pip[p]]:=w;
+            assign(p,w);
+            piv[p]:=w[p];
+            pin[p]:=now;
+            p:=fail; # to bail out while loop
+          fi;
+
+        fi;
+      od;
+      if not clean then
+        # only reduced, did not do gcd
+        Add(extra,v);
+      elif p<=n then
+        # new pivot position
+        v:=nv(v); # norm (so we can compare with <)
+        pip[p]:=Length(mat)+1;
+        #Add(mat,v);
+        assign(p,v);
+        Info(InfoMatInt,1,"Added @",Length(mat));
+        piv[p]:=v[p];
+        pin[p]:=v*v;
+      fi;
+    od;
+    # now we've processed all. Clean out extra
+    new:=List(Filtered(Set(extra),x->not IsZero(x)),nv);
+    Info(InfoMatInt,1,"After ",try,": ",Length(extra)," to ",Length(new),
+      " new ones");
+    extra:=[];
+
+  od;
+
+  return Filtered(Concatenation(mat,new),x->not IsZero(x));
+
+end);