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.