Schreier–Sims algorithm

Prerequisites: Basic concepts of abstract algebra, Group theory, and Permutations

Introduction

The Schreier–Sims algorithm is an algorithm in computational group theory, named after the mathematicians Otto Schreier and Charles Sims. It can solve many problems in polynomial time, such as finding the order of a finite permutation group and checking whether a given permutation belongs to a given group. The Schreier–Sims algorithm was first introduced by Sims in 1970 based on Schreier's lemma. In 1981¹, Donald Knuth further improved its running time. A faster randomized version was later developed. Computer algebra systems such as GAP and Magma usually use highly optimized Monte Carlo versions² of this algorithm.

Overview

The Schreier–Sims algorithm mainly attempts to solve the following problem:

• Given a set of size $n$, a set $X$, and a set $S$ of permutations, how can we efficiently store the permutation group generated by $S$, namely $G=\langle S\rangle$, in a computer and perform a series of queries on that group?

Clearly, such a group $G$ may be very large, far larger than the sets $X$ and $S$. For example, the order of the $n$-th symmetric group $S_n=\langle(123\cdots n),(12)\rangle$ is $n!$, but it can be generated by only two permutations. Storing every element of the group is impractical.

Similar to constructing a linear basis of a vector space using Gaussian elimination, the Schreier–Sims algorithm seeks a “basis” for the finite permutation group $G$:

1. The input is a generating set for $G$, namely $S$, containing several permutations.
2. If $G$ is nontrivial, the action of $G$ always moves some point $\beta$, namely $|\beta^G|>1$.
3. Find the orbit of point $\beta$, $\Delta=\beta^{G}$, and for every point $\delta\in\Delta$, find in $G$ a permutation $t_{\delta}$ that moves $\beta$ to $\delta$.
4. Find a generating set for the stabilizer of point $\beta$, $G_{\beta}$, denoted $S'$.
5. Recursively call the algorithm on $G'=\langle S'\rangle$ until the trivial group $\{e\}$ is reached.

The reason this works is that the stabilizer of point $\beta$, $G_{\beta}$, is a subgroup of $G$, and its right cosets partition $G$. These cosets correspond one-to-one with the orbit of point $\beta$, $\beta^G$. The permutations $t_{\delta}$ found in step 3 are representatives of these cosets. The set $T$ of coset representatives is called a transversal of the stabilizer $G_{\beta}$. In other words, every element of $G$, denoted $g$, corresponds to a unique pair $(h,t)\in G_\beta\times T$, with $g=ht$. Therefore, storing the subgroup $G_{\beta}$ together with its transversal $T$ is enough to store all of $G$. The problem of storing the subgroup $G_{\beta}$ has already been solved: simply make the recursive call.

The implementation still has many details to organize, which is the main subject of this article. Before that, we first examine the result after calling the algorithm and see what structure it uses to store the group and what queries it can answer. Some concepts are needed first.

Stabilizer chain

Suppose the algorithm is called on $G=\{S\}$ a total of $k$ times. At the $i$-th call, the input is a set $S^{(i-1)}$ of permutations, the point found is $\beta_i$, the transversal obtained is $T_i$, and the generating set of the stabilizer is $S^{(i)}$. Let $G^{(i)}=\langle S^{(i)}\rangle$. The algorithm obtains the subgroup chain

Moreover, every subgroup in the chain is a stabilizer:

Thus, the Schreier–Sims algorithm can be viewed as computing a stabilizer chain.

Base and strong generating set

If a subset of $X$, denoted $B$, satisfies $G_B=\{e\}$, then $B$ is called a base of permutation group $G$. Clearly, the algorithm obtains a base $B=\{\beta_1,\cdots,\beta_k\}$. This means that the results of the permutation $g\in G$ acting on these points uniquely determine this permutation in $G$. The base output by the algorithm also satisfies

That is, every point in the base contains useful information about the elements of the group. Such a base is called nonredundant. The algorithm always outputs a nonredundant base, whose corresponding stabilizer chain is strictly descending.

The union of the generating sets obtained by the algorithm,

is also a generating set for $G$, and $\langle\bar S\cap G^{(i)}\rangle = G^{(i)}$ holds. A generating set satisfying this condition is called a strong generating set of $G$ relative to the base $B$. Thus, the Schreier–Sims algorithm can also be viewed as computing a base and strong generating set (BSGS) for group $G$. These two descriptions are equivalent, and this article does not distinguish between them.

The algorithm also obtains a sequence of orbits $\Delta_i=\beta_{i}^{G^{(i-1)}}$ and their corresponding transversals $T_i$. These orbits are called the fundamental orbits of $G$. Whenever this article mentions the stabilizer chain of group $G$ or a base and strong generating set, the corresponding fundamental orbits and transversals are assumed to have been obtained as well.

Data structure

This article provides a series of pseudocode listings. In the pseudocode, the stabilizer chain (or BSGS) is stored in a data structure $C$:

The members of this structure are the generating set $S$ of the current group, the fundamental orbit $\Delta$, the corresponding transversal $T$, and the stabilizer $C'$ stored as a nested substructure. The stabilizer $C'$ is of course another structure of the same kind. The whole group is stored as a layered structure, with each layer describing a group in the stabilizer chain. The innermost layer is an empty structure representing $G^{(k)}=\{e\}$.

In the pseudocode, the members can be accessed as $C.generators$, $C.orbit$, $C.transversal$, and $C.next$, respectively. The first element $C.orbit[0]$ of the orbit is by default the base point $\beta$, and the corresponding coset representative $C.transversal[\beta]$ is by default the identity transformation $e$. Although array indices are used here to access elements of the orbit and transversal, they need not actually be stored as arrays; they should be understood as providing a way to access the first orbit element and query the corresponding coset representative by an orbit point. Specific implementation details are discussed below.

Applications

Once a base and strong generating set has been obtained, a series of queries about the group can be answered. The most basic ones, and the ones most often encountered in algorithm contests, are finding the order of the group and deciding whether a permutation belongs to a given group.

Order of a group

If a base and strong generating set for group $G$ is known, Lagrange's theorem and the orbit-stabilizer theorem show that the order of group $G$ can be computed as

Therefore, multiplying the sizes of all transversals $T_i$ (or equivalently, the lengths of all fundamental orbits $\Delta_i$) gives the order of $G$.

Membership testing

Given a base and strong generating set for $G$, we can also decide whether a permutation $h$ belongs to $G$. This is the membership testing problem.

The problem can be solved recursively. To decide whether $h\in G^{(i-1)}$, first find, in $G^{(i)}$, the representative of the coset containing $h$, namely $t\in T_i$. If it exists, set $h'=ht^{-1}$. Then $h=h't$, and $h\in G^{(i-1)}$ is equivalent to $h'\in G^{(i)}$; the problem is reduced to testing $h'\in G^{(i)}$. If no such $t$ exists, or if the recursion has reached $G^{(k)}=\{e\}$ with $h\neq e$, then $h\notin G$. In fact, this process not only decides whether $h\in G$, but when $h\in G$ it also expresses $h$ as a product $t_k\cdots t_2t_1$ of transversal representatives, where $t_i\in T_i$. Such a representation exists and is unique for elements of $G$. This again proves the formula for the order above.

The process can be written as the following pseudocode:

As we will see below, membership testing is also an important component of the implementation of the Schreier–Sims algorithm discussed in this article.

Computing orbits, transversals, and stabilizers

To implement Schreier–Sims, we first need to solve the following subproblem:³

• Given a group $G$ and its generating set $S$, and a point $\beta$, how do we find the orbit $\beta^G$, the corresponding transversal $T$, and a generating set for the stabilizer $G_\beta$?

This is the problem addressed in this section.

Storing orbits and transversals

The orbit and transversal can be obtained by direct search. The pseudocode is:

In an implementation, either breadth-first search or depth-first search works. The number of states visited is $|S||T|$, so the time complexity is acceptable as long as the orbit and transversal are stored reasonably.

In the context of permutation groups, an orbit is simply a set of at most $n$ points. A Boolean array or an unordered set can be used for efficient lookup and insertion. Both operations take $O(1)$ time and the total space is $O(n)$. Depending on the implementation of the transversal, an additional mark for the first element may be needed.

The question is what data structure should store the transversal $T$.

Direct storage

The simplest method is to store every element in the transversal $T$, denoted $t$, directly. A single permutation is stored in one-line notation, requiring exactly $n$ space, so the space complexity of storing the transversal is $O(|T|n)$. The advantage is $O(1)$ access to an individual representative, at the cost of $O(|T|n)$ time to compute all representatives initially.

Schreier tree

Another common method is to use a tree structure to store the transversal. It is called a Schreier tree or Schreier vector⁴. The root is $\beta$, and the vertices are the elements of the orbit $\Delta$. Whenever a new vertex reached from $\delta$, namely $\gamma=\delta^s$, is found during the search, an edge is added from the old vertex $\delta$ to the new vertex $\gamma$, and the index (or pointer) of $S$'s permutation $s$ is stored on the edge. Since the generating set is already stored, the additional space for the transversal is $O(|T|)$. For large $n$, this can save considerable space and the initial computation takes $O(|T|)$ time. The drawback is that obtaining a representative requires traversing the path from the vertex to the root and recomputing the representative, so the time complexity depends heavily on the depth of the Schreier tree. In general, the depth may reach $O(n)$.

In practice, the time-space tradeoff should be chosen according to the situation. In algorithm contests, $n$ is usually not large and space is sufficient, while time is often the bottleneck. As we will see, the most time-consuming step of Schreier–Sims needs to access coset representatives repeatedly, so direct storage is often better. In some applications, however, $n$ may be large and memory may be tight, making Schreier trees useful. To avoid excessive depth, researchers have proposed many methods for rebuilding a shallow Schreier tree when the depth becomes too large. Interested readers can consult the references at the end.

In the pseudocode, we do not distinguish between implementations of the transversal. We only assume that the data structure storing the transversal $T$ supports accessing and modifying, for each orbit element $\delta\in\Delta$, the corresponding representative $T[\delta]$.

Schreier's lemma

Once the orbit $\beta^G$ and transversal $T$ have been obtained, Schreier's lemma provides a way to obtain a generating set for the stabilizer $G_\beta$.

Since the subgroup corresponding to $T$ is the stabilizer $G_\beta$, after finding $T$, combining the generating set of $G$, namely $S$, yields a generating set for $G_\beta$.

Algorithm

A small modification of the pseudocode above computes the orbit, transversal, and generating set of the stabilizer at the same time:

For each pair $(\delta,s)$ in the pseudocode, $t_\delta st_{\delta^s}^{-1}$ is a new Schreier generator only when it does not produce a new orbit element; otherwise it is the identity transformation $e$. Therefore, the algorithm generates at most

Schreier generators, including the initial identity transformation. For general cases this bound is tight.⁶ For actual finite groups, however, the bound is quite loose: most newly obtained Schreier generators are duplicates of earlier generators or can be composed from them.

Since the basic Schreier–Sims procedure can be implemented by recursively calling the orbit-and-stabilizer algorithm above, we already have a naive implementation. Without filtering, however, the number of Schreier generators grows exponentially: repeatedly applying $O(|S_{i}|)=O(|S_{i-1}||T_i|)$ shows that the innermost stabilizer generating set can reach $O(|S||G|)$. This is absurdly inefficient, since the innermost stabilizer is $\{e\}$.

Sims provided a method to restrict the growth of Schreier generators, guaranteeing that the final strong generating set $\bar S$ has size $O(n^2)$. This makes it possible to compute a base and strong generating set in polynomial time.

The Schreier–Sims algorithm

To solve the problem above, this section discusses an incremental implementation of Schreier–Sims, whose strong generating set has size $O(n^2)$.

Sifting

The core optimization over the naive algorithm is simple: before adding any Schreier generator to the stabilizer's generating set, the algorithm first performs sifting. Sifting means checking whether a new Schreier generator already belongs to the subgroup generated by the existing generators, and adding only generators that do not yet belong to it. The membership algorithm $\textrm{MembershipTest}(C,h)$ introduced above is enough.

However, this requires that a stabilizer chain (or BSGS) for the group generated by the Schreier generators obtained so far has already been built. Thus, whenever a new Schreier generator is added to a stabilizer's generating set, the inner stabilizer chain must be maintained dynamically for later sifting. Yet inserting one generator at the current level may produce several Schreier generators, causing several updates to the inner structure; each update of the inner structure may in turn trigger many updates deeper inside.

It may seem that the exponential explosion remains. It does not. Because sifting is performed first, a level is updated only when the generator being added truly enlarges the group stored at that level. Therefore, the number of updates at one level equals the number of strict growths of its stored group. A group of size $|G|$ has a subgroup chain of length at most $\log|G|$, because Lagrange's theorem ensures that the group size at least doubles with every increase in chain length. Thus, each level is updated at most $\log|G|$ times, and the final strong generating set $\bar S$ has size $|B|\log|G|$.

This estimate can be improved. The group $G$ here is known to be a subgroup of the $n$-th symmetric group $S_n$, so the subgroup chain length of $G$ is no greater than that of $S_n$. It can be proved⁷ that the length of a strictly increasing subgroup chain in $S_n$ is at most $3n/2$. Hence, each level is actually updated $O(n)$ times. The base size is also at most $n$, so the final strong generating set $\bar S$ has size $O(n^2)$.

The sifting method described here was proposed by Sims and is also called the Sims filter. A more complicated method proposed by Jerrum, called the Jerrum filter, guarantees a strong generating set of size $O(n)$. Interested readers can study it independently.

One small optimization is to have $\textrm{MembershipTest}(C,h)$ return the final “siftee”⁸ $h$ rather than a Boolean value (that is, use $\textbf{return }h$ instead of lines $10$ and $14$ of the pseudocode). If the siftee $h\neq e$, membership testing has failed, and we can add the siftee $h$ itself rather than the original $h$ to the current level. This siftee $h$ has already had several transversal factors removed, so fewer elements need to be moved and the local computation is reduced. This optimization has no effect on the overall complexity.

Procedure

We can now describe the Schreier–Sims procedure. First initialize an empty structure $C$ to store the stabilizer chain. Then add the generators to structure $C$ from $S$ one by one. The final structure $C$ is the stabilizer chain of $\langle S\rangle$. The pseudocode is:

The core is adding a new generator to the current structure $C$, namely $s$, through the subroutine $\textrm{Extend}(C,s)$. Before and after adding $s$, $C$ must remain a stabilizer chain. Thus, before adding $s$, we first sift it. If $s$ is not already in the group, compute the orbit, transversal, and Schreier generators incrementally. “Incrementally” means that computations already performed are not repeated, which is essential for the intended complexity.

Consider how to turn $\textrm{OrbitTransversalStabilizer}(S,\beta)$ into an incremental version. The search state space is $\Delta\times S$. After adding a new generator $s$, the state space becomes $\Delta'\times\left(S\cup\{s\}\right)$. Their difference is

Therefore, when adding the new generator $s$, first compute combinations of the new generator $s$ with the old orbit and corresponding transversal. If new orbit elements are found, then consider those elements with all generators, old and new. Repeat until the orbit no longer grows.

The pseudocode for adding a permutation to structure $C$, namely $g$, is:

This gives the complete Schreier–Sims algorithm.

Another implementation

The implementation above is correct, but lines $12\sim 19$ and $23\sim 30$ are somewhat repetitive. Knuth proposed a more concise recursive implementation in his paper. He treats the repeated part as updating the transversal (and orbit). Each update combines the transversal with all generators; depending on whether a new representative is produced, it recursively calls itself or the generator-adding procedure. The pseudocode is:

Comparing this pseudocode with the previous section shows that it is correct, and the two implementations have the same complexity.

Complexity

We need some notation to analyze the complexity. Let the length of a permutation be $n$ and the size $|S|$ of the generating set be $m$. Let the length of the resulting nonredundant base be $|B|$. In the $i$-th stabilizer $G_i$ from the outside in, let the number of generators be $|S_{i-1}|$, and let the size of the transversal (or the orbit length) be $|T_i|$. We now analyze the time complexity of the incremental implementation above. The algorithm has two main parts: sifting, and computing orbits, transversals, and Schreier generators.

The input generators and the Schreier generators obtained during the algorithm must all be sifted, so the total number of sifting operations is $O(\sum_{i=1}^{|B|}|S_{i-1}||T_i|+|S|)$. A single sifting operation computes products with $O(|B|)$ transversal representatives. Let $\tau$ be the time to multiply by one representative. With direct storage, $\tau\in O(n)$; with a Schreier tree, $\tau\in O(n^2)$. Thus, the total time for sifting is $O(\tau|B|\sum_{i=1}^{|B|}|S_{i-1}||T_i|+\tau|B||S|)$.

For orbit-related information, the incremental implementation computes each pair $(\delta,s)\in\Delta_i\times S_{i-1}$ in the state space at most once. Depending on the storage method, computing a transversal representative takes either $O(n)$ or $O(1)$ time, but in either case this does not exceed the cost of computing a Schreier generator. With direct storage this is $O(n)$; with a Schreier tree it is $O(n^2)$. Comparing this with the total complexity of sifting shows that computing this information does not exceed the time spent sifting. Therefore, the time complexity of Schreier–Sims is the $O(\tau|B|\sum_{i=1}^{|B|}|S_{i-1}||T_i|+\tau|B||S|)$ obtained above.

For space complexity, the final data structure stores $O(\sum_{i=1}^{|B|}|S_{i-1}|)$ generators and $O(\sum_{i=1}^{|B|}|T_i|)$ transversal representatives. With direct storage, both generators and representatives require $O(n)$ space; with a Schreier tree, generators require $O(n)$ space and representatives require only $O(1)$ space.

As stated earlier, the length of a strictly increasing subgroup chain in the symmetric group on $n$ points is $O(n)$, and this applies to both $|B|$ and $|S_i|$. Therefore, an implementation using direct storage has time complexity $O(n^5+mn^2)$ and space complexity $O(n^3)$. When $\log|G|\in O(n)$, a better estimate is time complexity⁹ $O(n^2\log^3|G|+mn\log|G|)$ and space complexity $O(n^2\log|G|)$. For random generating sets, experiments show that the actual complexity is significantly below $\Theta(n^5)$, approximately $\Theta(n^4)$.

Although a Schreier tree introduces an additional factor of $n$ in the time complexity compared with direct storage, when $n$ is large and the group itself is much smaller than the symmetric group on $n$ points, its space complexity is $O(n\log^2|G|)$, far smaller than the $O(n^2\log|G|)$ of direct storage. Such cases where Schreier-tree storage is better are rare in algorithm contests.

Reference implementation

Here is a reference implementation of Schreier–Sims. Since $n$ is small, the implementation directly fixes the base as $B=\{n,n-1,\cdots,1\}$ instead of choosing it algorithmically. This makes the permutations in the $k$-th layer from the inside out (excluding the empty structure) change only the first $k$ elements, simplifying later computation. Another optimization is that the stored transversal representatives are actually their inverse permutations, which simplifies permutation operations.

#include <iostream>
#include <vector>

// A permutation.
class Permutation {
  std::vector<int> perm;

  void shrink() {
    int m = perm.size();
    for (; m && perm[m - 1] == m - 1; --m);
    perm.resize(m);
  }

 public:
  Permutation() {}

  Permutation(const std::vector<int>& vec) : perm(vec) { shrink(); }

  int operator[](int i) const { return i < perm.size() ? perm[i] : i; }

  bool empty() const { return perm.empty(); }

  // First LHS then RHS.
  Permutation operator*(const Permutation& rhs) const {
    Permutation res;
    res.perm.resize(std::max(perm.size(), rhs.perm.size()));
    for (int i = 0; i < res.perm.size(); ++i) {
      res.perm[i] = rhs[(*this)[i]];
    }
    res.shrink();
    return res;
  }

  // First LHS^{-1} then RHS.
  Permutation operator/(const Permutation& rhs) const {
    Permutation res;
    res.perm.resize(std::max(perm.size(), rhs.perm.size()));
    for (int i = 0; i < res.perm.size(); ++i) {
      res.perm[(*this)[i]] = rhs[i];
    }
    res.shrink();
    return res;
  }

  // Inverse.
  Permutation inv() const {
    Permutation res;
    res.perm.resize(perm.size());
    for (int i = 0; i < res.perm.size(); ++i) {
      res.perm[(*this)[i]] = i;
    }
    return res;
  }
};

// A stabilizer chain (a.k.a., BSGS) for a group.
class PermutationGroup {
  size_t n, k;
  std::vector<bool> orbit;               // Orbit of the n-th point.
  std::vector<Permutation> generators;   // Generators.
  std::vector<Permutation> transversal;  // Inverse of coset representatives.
  PermutationGroup* next;                // Stabilizer.

  // Sift a permutation.
  void sift(Permutation& h) const {
    if (!n) return;
    int i = h[n - 1];
    if (!orbit[i]) return;
    h = h * transversal[i];
    next->sift(h);
  }

  // Add one more element into the transversal.
  void extend_transversal(Permutation t) {
    int i = t[n - 1];
    if (!orbit[i]) {
      ++k;
      orbit[i] = true;
      transversal[i] = t.inv();
      for (const auto& s : generators) {
        extend_transversal(t * s);
      }
    } else {
      next->extend(t * transversal[i]);
    }
  }

 public:
  PermutationGroup(int n)
      : n(n), k(1), orbit(n), transversal(n), next(nullptr) {
    if (!n) return;
    // Initialize the current layer.
    orbit[n - 1] = true;
    next = new PermutationGroup(n - 1);
  }

  // Destructor.
  ~PermutationGroup() {
    if (next) delete next;
  }

  // Add one more permutation into the group.
  void extend(Permutation g) {
    sift(g);
    if (g.empty()) return;
    generators.emplace_back(g);
    for (int i = 0; i < n; ++i) {
      if (orbit[i]) {
        extend_transversal(transversal[i] / g);
      }
    }
  }

  // Check whether a permutation belongs to the group.
  bool membership_test(Permutation h) const {
    sift(h);
    return h.empty();
  }

  // Return the size of the group.
  long long size() const { return n ? next->size() * k : 1LL; }
};

int main() {
  int n, m;
  std::cin >> n >> m;
  PermutationGroup group(n);
  // Read permutations and insert them to the group.
  std::vector<int> vec(n);
  for (; m; --m) {
    for (int& x : vec) {
      std::cin >> x;
      --x;  // Index starting at 0.
    }
    group.extend(Permutation(vec));
  }
  // Output the size of the group.
  std::cout << group.size();
  return 0;
}

Exercises

• LOJ 177. Order of a generated subgroup
• [WC2017] Board
• Permutations
• [Grand Prix of Yekaterinburg 2015] Problem H. Heimdall

References and notes

• Schreier–Sims algorithm - Wikipedia
• Sims, Charles C, Computational methods in the study of permutation groups, Computational Problems in Abstract Algebra, pp. 169–183, Pergamon, Oxford, 1970.
• Knuth, Donald E. Efficient representation of perm groups, Combinatorica 11 (1991), no. 1, 33–43.
• Ákos Seress, Permutation Group Algorithms, Cambridge University Press
• Alexander Hulpke's Notes on Computational Group Theory
• Derek Holt's Slides on The Schreier–Sims algorithm for finite permutation groups
• Martin Jaggi, Implementations of 3 Types of the Schreier–Sims Algorithm, MAS334 - Mathematics Computing Project, 2005
• Henrik Bäärnhielm. The Schreier–Sims algorithm for matrix groups