Stabilizer chain gives wrong group size, not a divisor of actual size

[fingolfin]

While I would expect that a Monte-Carlo orbit stabilizer algorithm might sometimes underestimate the size of a group, I think it should always produce a divisor of the actual group size.

But Stefan Kohl reported that he sometimes encounters cases where this is not the case. Using his example, I compute a seed of the random source to reliably reproduce the issue:

gap> LoadPackage("genss");;
gap> G := Group((5,13)(8,14)(16,21)(18,22)(20,23),  (3,7)(8,12),(4,15)(6,16)(7,17)(9,18)(11,19)(12,20));;
gap> state := [ 43, [
>   189898554, 234529222, 126352480, 64129502, 218313656, 152728658, 215163821,
>   52767390, 147276579, 239260582, 188207164, 72447725, 101779050, 184338092,
>   96094684, 125971651, 75056515, 245721574, 163283225, 151620100, 174263367,
>   103006280, 91299698, 20291082, 240961624, 98550507, 149524538, 126499941,
>   244229841, 213694414, 15170774, 181886954, 40766256, 25841967, 37434579,
>   46798339, 142217533, 214018634, 166382803, 142691424, 108825002, 159053691,
>   1192096, 89158305, 27895964, 8507667, 68651150, 163843714, 202319118,
>   130760480, 145867311, 246634910, 24086403, 126617031, 24317613 ] ];;
gap> Reset(GlobalRandomSource, state);;
gap> S:=StabilizerChain(G);
<stabchain size=3456 orblen=3 layer=1 SchreierDepth=2>
 <stabchain size=1152 orblen=2 layer=2 SchreierDepth=1>
  <stabchain size=576 orblen=2 layer=3 SchreierDepth=1>
   <stabchain size=288 orblen=2 layer=4 SchreierDepth=1>
    <stabchain size=144 orblen=3 layer=5 SchreierDepth=2>
     <stabchain size=48 orblen=4 layer=6 SchreierDepth=2>
      <stabchain size=12 orblen=4 layer=7 SchreierDepth=2>
       <stabchain size=3 orblen=3 layer=8 SchreierDepth=2>
gap> Size(S);
3456
gap> Size(G);
4320
gap>

Here is a second example, again due to Stefan Kohl:

gap> LoadPackage("genss");;
gap> G := Group((4,15)(8,17)(12,19), (5,14)(9,16)(13,18), (6,7)(12,13));;
gap> state := [ 42, [
>   79989323, 27411214, 266168953, 1049796, 23280651, 202232511, 11005996,
>   253411464, 164097667, 28552034, 96638348, 115452439, 168282399, 26083525,
>   142194947, 22649614, 7236209, 168097464, 133401501, 19933020, 197193335,
>   190764831, 75671277, 11036977, 65983267, 92432952, 253721782, 181770330,
>   185331281, 169073647, 52541843, 169715992, 185394925, 155912256, 41518560,
>   234177554, 96338704, 78566396, 179705249, 148634593, 241641300, 125635167,
>   21979545, 118753100, 71778345, 132389610, 84882183, 44955068, 194977331,
>   258250357, 43755991, 166004159, 63606826, 151164300, 159967769 ] ];;
gap> Reset(GlobalRandomSource, state);;
gap> S:=StabilizerChain(G);
<stabchain size=72 orblen=2 layer=1 SchreierDepth=1>
 <stabchain size=36 orblen=2 layer=2 SchreierDepth=1>
  <stabchain size=18 orblen=2 layer=3 SchreierDepth=1>
   <stabchain size=9 orblen=3 layer=4 SchreierDepth=1>
    <stabchain size=3 orblen=3 layer=5 SchreierDepth=2>
gap> Size(S);
72
gap> Size(G);
96

[fingolfin]

OK, looking closer, and thinking a bit more about it, my previous thought that the computer size must be a divisor of the actual size was wrong -- resp. it only holds for some variants of the randomized Schreier-Sims, but certainly not the one implemened in genss.

The problem becomes visible in the second example:
There, the last two steps of the stab chain are recognized as a group of order 9, when in reality it should be a group of order 12, namely A4. The two layers of order 3 then correspond to different subgroups of order 3. The one at the bottom is correct there, but in layer 4, we are missing one generator.
This leads to an orbit of size 3 (instead of 4), and hence the result is off by a factor of 4/3.

So, all in all, this is not a bug in the code, it's an artifact of randomiziation, and the fact that the default error bound is only 1/1024. So if we run the computation several ten- or hundred thousand times, we simply will get such an error.

To make it less likely, the user can use the ErrorBound or VerifyElements options.

[fingolfin]

Hmm, on the other hand, the documentation for StabilizerChain explicitly states this:

However, the only error possible is that some of the subgroups in the
stabiliser chain are proper subgroups of the actual point stabilisers. So
the resulting group order is always a divisor of the actual group order and
if the two are equal, then the stabiliser chain is proved to be correct.

But I am not sure whether this is bug in the code, or in the documentation. Hrm.

[fingolfin]

One remark regarding the error distribution: Consider this test program

LoadPackage("genss");
gens := [(4,15)(8,17)(12,19), (5,14)(9,16)(13,18), (6,7)(12,13)];
N := Size(Group(gens));
i := 0;
repeat
  i := i + 1;
  G := Group(gens);
  S := StabilizerChain(G);
  n := Size(S);
  if i mod 1000 = 0 then Print("i = ",i,"\n"); fi;
until n <> N;

Running it eventually leads to a wrong answer (which by itself is not that surprising, we are looking at a Monte Carlo algorithm after all). I did this multiple times, and the number of iterations until an error occurred were as follows, starting in a fresh GAP session: 3259 (ouch), 36058, 39812, 96221, 49468.

Given that the default setting for ErrorBound is something like 1/1024, that seems not too bad...

[fingolfin]

Another remark: If I modify the above example to use S := StabilizerChain(G, rec(ImmediateVerificationElements:=10));, then even after over 9 million invocations there are no errors. Of course there will be errors eventualy even in this case... But hopefull not soon ;-).