Fix (rare) wrong result for stabilizer computations

lib/ctblgrp.gi

@@ -2056,11 +2056,16 @@ end );
 ##
 #F  DixonSplit(<D>) . .  calculate matrix,split spaces and obtain characters
 ##
-InstallGlobalFunction( DixonSplit, function(D)
-local r,i,j,ch,ra,bsm,
+InstallGlobalFunction( DixonSplit, function(arg)
+local D,r,i,j,ch,ra,bsm,
       gens;
 
-  bsm:=BestSplittingMatrix(D);
+  D:=arg[1];
+  if Length(arg)>1 then
+    bsm:=arg[2];
+  else
+    bsm:=BestSplittingMatrix(D);
+  fi;
   if bsm<>fail then
     SplitStep(D,bsm);
   fi;

lib/ghomfp.gi

@@ -1225,7 +1225,11 @@ local aug,r,sec,t,expwrd,rels,ab,s,m,img,gen,i,j,t1,t2,tn;
     Add(img,m);
   od;
   aug.primaryImages:=img;
-  sec:=List(sec,x->LinearCombinationPcgs(img,x));
+  if ForAll(img,IsOne) then
+    sec:=List(sec,x->img[1]);
+  else
+    sec:=List(sec,x->LinearCombinationPcgs(img,x));
+  fi;
   aug.secondaryImages:=sec;
 
   m:=List(aug.primaryGeneratorWords,x->ElementOfFpGroup(FamilyObj(One(u)),x));

lib/grplatt.gi

@@ -3510,29 +3510,53 @@ local divs,limit,mode,l,process,done,bound,maxer,prime;
   if limit<20 then limit:=1;fi;
   maxer:=function(sub)
   local m,a,b,len,sz,i,j,k,r,tb;
-    #Print("call maxer for ",Size(sub)," |l|=",Length(l)," |process|=",Length(process),"\n");
-    if bound>1 and prime^(Length(AbelianInvariants(sub))-3)>bound then
-      i:=0;
-      repeat
-        m:=BoundedIndexAbelianized(G,sub,bound*prime^i);
-        i:=i+1;
-      until Size(m)<Size(sub);
-      if Size(m)^2<Size(sub) then
-        m:=MaximalSubgroupClassReps(sub:cheap:=false);
-      else
-        # add maxes not containing derived
-        if IsSolvableGroup(sub) then
-          i:=IsomorphismPcGroup(sub);
-          a:=DerivedSubgroup(Image(i));
-          a:=Filtered(MaximalSubgroupClassReps(Image(i)),
-            x->not IsSubset(x,a));
-          a:=List(a,x->PreImage(i,x));
+    Info(InfoLattice,1,"call maxer for ",Size(sub)," |l|=",Length(l),
+      " |process|=",Length(process));
+    if bound>1 then
+      # nonabelian indices
+      a:=CompositionSeries(sub);
+      m:=1;
+      for i in [2..Length(a)] do
+        if IndexNC(a[i-1],a[i])>m and 
+          not HasAbelianFactorGroup(a[i-1],a[i]) then
+          m:=Maximum(IndexNC(a[i-1],a[i]),m);
+        fi;
+      od;
+      # big (hope simple) bits
+      if m>=10^7 and (not IsPerfectGroup(sub)) 
+        # proper abelian factor
+        and IndexNC(sub,PerfectResiduum(sub)) in [2..bound] 
+        # not all abelian is direct factor
+        and IndexNC(sub,
+          ClosureGroup(PerfectResiduum(sub),RadicalGroup(sub)))>1 then
+        m:=MaximalSubgroupClassReps(
+          ClosureGroup(PerfectResiduum(sub),RadicalGroup(sub)):cheap:=false);
+      elif bound>1 and prime^(Length(AbelianInvariants(sub))-3)>bound then
+        i:=0;
+        repeat
+          m:=BoundedIndexAbelianized(G,sub,bound*prime^i);
+          i:=i+1;
+        until Size(m)<Size(sub);
+        if Size(m)^2<Size(sub) then
+          m:=MaximalSubgroupClassReps(sub:cheap:=false);
         else
-          a:=DerivedSubgroup(sub);
-          a:=Filtered(MaximalSubgroupClassReps(sub),
-            x->not IsSubset(x,a));
+          # add maxes not containing derived
+          if IsSolvableGroup(sub) then
+            i:=IsomorphismPcGroup(sub);
+            a:=DerivedSubgroup(Image(i));
+            a:=Filtered(MaximalSubgroupClassReps(Image(i)),
+              x->not IsSubset(x,a));
+            a:=List(a,x->PreImage(i,x));
+          else
+            a:=DerivedSubgroup(sub);
+            a:=Filtered(MaximalSubgroupClassReps(sub),
+              x->not IsSubset(x,a));
+          fi;
+          m:=Concatenation([m],a);
         fi;
-        m:=Concatenation([m],a);
+
+      else
+        m:=MaximalSubgroupClassReps(sub:cheap:=false);
       fi;
     else
       m:=MaximalSubgroupClassReps(sub:cheap:=false);

lib/oprt.gi

@@ -1355,19 +1355,21 @@ local   orb,  stb,  rep,  p,  q,  img,  sch,  i,d,act,r,
       fi;
       p:=stopat;
       while p<=Length(orb) do
-	img := act( orb[ p ], acts[ i ] );
-	MakeImmutable(img);
-	q:=LookupDictionary(dict,img);
-	if q<=p then
-	  sch := getrep( p ) * gens[ i ] / getrep( q );
-	  if not sch in stb then
-	    stb:=ClosureSubgroupNC(stb,sch);
-	    if Size(G)/Size(stb)=Length(orb) then
-	      p:=Length(orb); #done
-	    fi;
-	  fi;
-	fi;
-	p:=p+1;
+        i:=1;
+        while i<=Length(gens) do
+          img := act( orb[ p ], acts[ i ] );
+          MakeImmutable(img);
+          q:=LookupDictionary(dict,img);
+          sch := getrep( p ) * gens[ i ] / getrep( q );
+          if not sch in stb then
+            stb:=ClosureSubgroupNC(stb,sch);
+            if Size(G)/Size(stb)=Length(orb) then
+              p:=Length(orb);i:=Length(gens); #done
+            fi;
+          fi;
+          i:=i+1;
+        od;
+        p:=p+1;
       od;
       #Error("after");
     fi;

lib/twocohom.gi

@@ -1515,7 +1515,7 @@ InstallGlobalFunction(FpGroupCocycle,function(arg)
 local r,z,ogens,n,gens,str,dim,i,j,f,rels,new,quot,g,p,collect,m,e,fp,old,sim,
       it,hom,trysy,prime,mindeg,fps,ei,mgens,mwrd,nn,newfree,mfpi,mmats,sub,
       tab,tab0,evalprod,gensmrep,invsmrep,zerob,step,simi,simiq,wasbold,
-      mon,ord,mn,melmvec;
+      mon,ord,mn,melmvec,killgens,frew,fffam,ofgens,rws,formalinverse;
 
   # function to evaluate product (as integer list) in gens (and their
   # inverses invs) with corresponding action mats
@@ -1549,6 +1549,17 @@ local r,z,ogens,n,gens,str,dim,i,j,f,rels,new,quot,g,p,collect,m,e,fp,old,sim,
     return w;
   end;
 
+  # formal inverse of module element
+  formalinverse:=function(w)
+  local l,i;
+    l:=[];
+    for i in Reversed(LetterRepAssocWord(w)) do
+      if IsEvenInt(i) then Add(l,i-1);
+      else Add(l,i+1);fi;
+    od;
+    return AssocWordByLetterRep(FamilyObj(w),l);
+  end;
+
   r:=arg[1];
   z:=arg[2];
   ogens:=GeneratorsOfGroup(r.presentation.group);
@@ -1562,6 +1573,14 @@ local r,z,ogens,n,gens,str,dim,i,j,f,rels,new,quot,g,p,collect,m,e,fp,old,sim,
     # will allow for reduced multiplication. Do so before the group, so there is no
     # issue about overwriting variables that might be needed later
     mon:=Image(r.monhom);
+    ofgens:=FreeGeneratorsOfFpMonoid(mon);
+    fffam:=FamilyObj(One(FreeMonoidOfFpMonoid(mon)));
+    frew:=ReducedConfluentRewritingSystem(mon);
+    # generators that get killed. Do not use in RHS
+    killgens:=Filtered(GeneratorsOfMonoid(mon),x->UnderlyingElement(x)<>
+      ReducedForm(frew,UnderlyingElement(x)));
+    killgens:=Union(List(killgens,x->LetterRepAssocWord(UnderlyingElement(x))));
+
     str:=List(GeneratorsOfMonoid(mon),String);
     mn:=Length(str);
     for i in [1..dim] do
@@ -1577,17 +1596,8 @@ local r,z,ogens,n,gens,str,dim,i,j,f,rels,new,quot,g,p,collect,m,e,fp,old,sim,
       Add(rels,[gens[i]*gens[i+1],One(f)]);
       Add(rels,[gens[i+1]*gens[i],One(f)]);
     od;
-    # rules with tails
-    for i in [1..Length(r.presentation.monrulpos)] do
-      new:=RelationsOfFpMonoid(mon)[r.presentation.monrulpos[i]];
-      new:=List(new,x->MappedWord(x,FreeGeneratorsOfFpMonoid(mon),gens{[1..mn]}));
-      new[2]:=new[2]*melmvec(z{[(i-1)*dim+1..i*dim]});
-      Add(rels,new);
-    od;
 
-    if r.presentation.killrelators<>[] then Error("TODO killrelators");fi;
-
-    # elementary abelian
+    # module is elementary abelian
     for i in [mn+1,mn+3..Length(str)-1] do
       if prime=2 then
         Add(rels,[gens[i]^2,One(f)]);
@@ -1602,16 +1612,40 @@ local r,z,ogens,n,gens,str,dim,i,j,f,rels,new,quot,g,p,collect,m,e,fp,old,sim,
         Add(rels,[gens[j]*gens[i+1],gens[i+1]*gens[j]]);
       od;
     od;
+
     # module rules
     for i in [1..mn] do
-      for j in [mn+1..Length(str)] do
-        new:=r.module.generators[QuoInt(i+1,2)];
-        if IsEvenInt(i) then new:=new^-1; fi;
-        new:=new[QuoInt(j-mn+1,2)];
-        if IsEvenInt(j) then new:=-new;fi;
-        Add(rels,[gens[j]*gens[i],gens[i]*melmvec(new)]);
-      od;
+      if not i in killgens then
+        for j in [mn+1..Length(str)] do
+          new:=r.module.generators[QuoInt(i+1,2)];
+          if IsEvenInt(i) then new:=new^-1; fi;
+          new:=new[QuoInt(j-mn+1,2)];
+          if IsEvenInt(j) then new:=-new;fi;
+          Add(rels,[gens[j]*gens[i],gens[i]*melmvec(new)]);
+        od;
+      fi;
     od;
+
+    # rules with tails
+    for i in [1..Length(r.presentation.monrulpos)] do
+      new:=RelationsOfFpMonoid(mon)[r.presentation.monrulpos[i]];
+      new:=List(new,x->MappedWord(x,ofgens,gens{[1..mn]}));
+      new[2]:=new[2]*melmvec(z{[(i-1)*dim+1..i*dim]});
+      Add(rels,new);
+    od;
+
+    # Any killed generators -- use just the same word expression
+    if r.presentation.killrelators<>[] then 
+      for i in Filtered(killgens,IsOddInt) do
+        new:=AssocWordByLetterRep(fffam,[i]);
+        new:=[new,ReducedForm(frew,new)];
+        new:=List(new,x->MappedWord(x,
+          ofgens,gens{[1..mn]}));
+        Add(rels,new);
+        RemoveSet(killgens,i);
+      od;
+    fi;
+
     if HasReducedConfluentRewritingSystem(mon) then
       ord:=OrderingOfRewritingSystem(ReducedConfluentRewritingSystem(mon));
       if HasLevelsOfGenerators(ord) then
@@ -1642,10 +1676,26 @@ local r,z,ogens,n,gens,str,dim,i,j,f,rels,new,quot,g,p,collect,m,e,fp,old,sim,
       fi;
     od;
 
+    # temporary build to be able to reduce
+    mon:=FactorFreeMonoidByRelations(f,rels);
+    rws:=KnuthBendixRewritingSystem(FamilyObj(One(mon)),ord);
+
+    # handle inverses that get killed
+    for i in killgens do
+      new:=AssocWordByLetterRep(fffam,[i]);
+      new:=ReducedForm(frew,new); # what is it in the factor?
+      new:=MappedWord(new,ofgens,gens{[1..mn]});
+      # now left multiply with non-inverse generator
+      j:=ReducedForm(rws,gens[i-1]*new); #must be a tail.
+      j:=ReducedForm(rws,formalinverse(j)); # invert
+      Add(rels,[gens[i],new*j]);
+    od;
+
     mon:=FactorFreeMonoidByRelations(f,rels);
-    new:=KnuthBendixRewritingSystem(FamilyObj(One(mon)),ord);
-    MakeConfluent(new);  # will add rules to kill inverses, if needed
-    SetReducedConfluentRewritingSystem(mon,new);
+    rws:=KnuthBendixRewritingSystem(FamilyObj(One(mon)),ord:isconfluent);
+    ReduceRules(rws);
+    MakeConfluent(rws);  # will add rules to kill inverses, if needed
+    SetReducedConfluentRewritingSystem(mon,rws);
   fi;
 
   # now construct the group

tst/testbugfix/2021-09-17-Stabilizer.tst

@@ -0,0 +1,13 @@
+# Fix GitHub issue #4659, reported by Daniel Pook-Kolb
+gap> p_gens := [PermutationMat((1,2), 6), PermutationMat((1,2,3,4,5,6), 6)];;
+gap> p_group := Group(p_gens);;
+gap> s1 := DiagonalMat([-1, -1, 1, 1, 1, 1]);;
+gap> s_gens := Orbit(p_group, s1);;
+gap> ps_gens := Concatenation(p_gens, s_gens);;
+gap> ps_group := Group(ps_gens);;
+gap> v := [1,2,1,1,0,0];;
+gap> orbit := Orbit(ps_group, v);;
+gap> Size(ps_group);;
+gap> stab := Stabilizer(ps_group, v);;
+gap> Size(orbit) * Size(stab) = Size(ps_group);
+true