How is GAP generating all subgroups?

GAP can return all subgroups of $S_4$ pretty much instantaneously:

gap> AllSubgroups(SymmetricGroup(4));
[ Group(()), Group([ (1,2)(3,4) ]), Group([ (1,3)(2,4) ]), Group([ (1,4)(2,3) ]), Group([ (3,4) ]), Group([ (2,3) ]),
  Group([ (2,4) ]), Group([ (1,2) ]), Group([ (1,3) ]), Group([ (1,4) ]), Group([ (2,4,3) ]), Group([ (1,3,2) ]),
  Group([ (1,4,2) ]), Group([ (1,4,3) ]), Group([ (1,4)(2,3), (1,3)(2,4) ]), Group([ (3,4), (1,2)(3,4) ]),
  Group([ (1,4), (1,4)(2,3) ]), Group([ (2,4), (1,3)(2,4) ]), Group([ (1,3,2,4), (1,2)(3,4) ]), Group([ (1,4,3,2), (1,
   3)(2,4) ]), Group([ (1,2,4,3), (1,4)(2,3) ]), Group([ (3,4), (2,4,3) ]), Group([ (1,4), (1,4,3) ]),
  Group([ (2,3), (1,3,2) ]), Group([ (1,2), (1,4,2) ]), Group([ (1,4)(2,3), (1,3)(2,4), (3,4) ]), Group([ (1,2)
  (3,4), (1,3)(2,4), (1,4) ]), Group([ (1,2)(3,4), (1,4)(2,3), (2,4) ]), Group([ (1,4)(2,3), (1,3)(2,4), (2,4,3) ]),
  Group([ (1,4)(2,3), (1,3)(2,4), (2,4,3), (3,4) ]) ]

$S_4$ has 24 elements ($4! = 24$).

A naive approach to generating all subgroups for a given group G would consider all possible subsets of G (i.e. power set of G). In this case, there are 16,777,216 possible subsets to consider ($2^{24} = 16777216$).

Of course, the order of each subgroup must evenly divide the order of G so we can narrow down the list of subsets to consider by using this fact (i.e. only consider subsets that have an order that evenly divides the order of G). For $S_4$ that narrows down the list to 3,587,174.

That's still a lot of sets to consider! I'd be surprised if GAP is computing the subgroups on the fly using this naive approach.

So I'm wondering how does GAP return the answer so quickly. Does GAP have the subgroups to $S_4$ (and other common groups) pre-computed?

GAP stores precomputed lists of subgroups of (automorphism groups) of a number of simple subgroups, the rest is computed on the fly. However it uses far more group theory than only considering subsets. (For reasons of memory use, GAP will calculate subgroup only up to conjugacy by default.)

The methods use (very roughly) follow the following strategy:

• If $G$ has an elementary abelian normal subgroup $N\lhd G$, every subgroup $U$ as an image $A=UN/N$ in $G/N$ and an intersection $B=U\cap N$. Recursively calculate subgroups of $G/N$, as well as of $N$ (this is a vector space) and for every pair $(A,B)$ with $A/N\le G/N$ and $B\le N$ such that $B\lhd A$ calculate complements to $N/B$ in $A/B$ using cohomology.

• If $G$ has no solvable normal subgroup, then the socle $S=Soc(G)$ is a direct product of simple groups. Calculate subgroups of $S$ as subdirect products of subgroups of the simple factors, and then for each subgroup consider subgroups of $N_G(S)/N$.

• The subgroups of the simple factors are stored.

Two articles describing the process (Caveat: author bias) are:

MR1806212 Cannon, John J.; Cox, Bruce C.; Holt, Derek F. Computing the subgroups of a permutation group. J. Symbolic Comput. 31 (2001), no. 1-2, 149–161.

MR3206359 Hulpke, Alexander Calculation of the subgroups of a trivial-fitting group. ISSAC 2013—Proceedings of the 38th International Symposium on Symbolic and Algebraic Computation, 205–210, ACM, New York, 2013.

A couple things:

1. I highly doubt they precompute S4. You're right the naive approach is bad, and they likely don't use that either. I won't try to develop a better algorithm here, but it seems like there might be a slightly less naive way to find all of the distinct generating sets via dynamic programming, which I'd imagine would be reasonably efficient in comparison.

2. GAP is open source (hosted on Github here), so it's actually not that difficult to look into how the literal computations themselves work. The "AllSubgroups" command is defined in this file. It appears to rely on in a key way the ConjugacyClassesSubgroups method from this file. That being said, I don't know the syntax of GAP, so someone else would have to analyze those particular files better.

3. Looking through the code, I saw a reference to a paper on better methods for computing normal subgroups, which is here. I think I saw a reference to them using the underlying method here, but this may no longer be the state of the art. Still, for the theoretical side of (part) of what GAP does, this is probably the best bet.