Suffix Array Construction, in depth: prefix doubling with radix passes, SA-IS induced sorting traced step by step, and building at scale

The construction contract

The output is a permutation sa of 0..n-1 such that the suffix starting at sa[0] is the smallest, and so on. Two conventions cause most off-by-one bugs.

The sentinel. Linear-time algorithms append a symbol $ that is smaller than every real symbol and occurs once, so no suffix is a prefix of another. Here the text is mapped with ord(ch) + 1 and the sentinel is 0. Its suffix always sorts first, so the public function drops slot 0.

The alphabet. Induced sorting allocates one bucket per symbol, so symbols must be small integers. For token vocabularies, remap the symbols that occur to a dense range first.

Keep a slow reference implementation from day one. It is obviously correct, which is its whole job:

def oracle(text):
    return sorted(range(len(text)), key=lambda i: text[i:])

Prefix doubling with radix passes

Prefix doubling, from Manber and Myers, rests on one observation. If every suffix has a rank ordering it by its first k characters, the pair (rank of i, rank of i + k) orders it by its first 2k. Start with k = 1, where the rank is the symbol, and double until all ranks are distinct: at most log2 n rounds. Two stable counting-sort passes per round, second key first, make each round linear, so the total is O(n log n):

def counting_sort(order, key, m):          # stable
    cnt = [0] * m
    for i in order:
        cnt[key[i]] += 1
    pos, acc = [0] * m, 0
    for v in range(m):
        pos[v], acc = acc, acc + cnt[v]
    out = [0] * len(order)
    for i in order:
        out[pos[key[i]]] = i
        pos[key[i]] += 1
    return out

def doubling_radix(text):
    n = len(text)
    if n == 0:
        return []
    letters = {ch: r + 1 for r, ch in enumerate(sorted(set(text)))}
    rank = [letters[ch] for ch in text]    # 0 means "past the end"
    k = 1
    while True:
        second = [rank[i + k] if i + k < n else 0 for i in range(n)]
        m = max(rank) + 1
        sa = counting_sort(counting_sort(range(n), second, m), rank, m)
        new, r = [0] * n, 1
        new[sa[0]] = 1
        for a, b in zip(sa, sa[1:]):
            r += rank[a] != rank[b] or second[a] != second[b]
            new[b] = r
        rank = new
        if r == n:                         # all ranks distinct
            return sa
        k *= 2

The early exit means text without long repeats finishes in a few rounds; a long run of one symbol needs all of them. Doubling is short enough to review in one sitting and makes an independent check for a faster builder. Its weakness is memory traffic: every round streams several n-sized arrays through random access, which makes it several times slower than induced sorting on large inputs.

SA-IS: L-type, S-type and LMS positions

SA-IS labels each position. Suffix i is S-type if it is smaller than suffix i + 1 and L-type if larger. One right-to-left pass computes it: t[i] < t[i+1] means S, t[i] > t[i+1] means L, and equal symbols inherit the label of i + 1. The sentinel is S. An LMS (leftmost-S) position is an S-type position whose left neighbour is L-type; the text between consecutive LMS positions, inclusive, is an LMS substring.

Two facts drive the algorithm. Within the bucket of suffixes starting with one symbol, L-type suffixes sort before S-type ones, so each bucket has an L head and an S tail. And the induced-sorting lemma: given the sorted LMS suffixes, one left-to-right scan places every L-type suffix and one right-to-left scan places every S-type suffix. For each suffix j passed, suffix j - 1 goes to the next free head (if L) or tail (if S) of its first symbol's bucket. Order propagates because suffix j - 1 is suffix j with one symbol prepended.

So the problem shrinks to sorting the LMS suffixes, of which there are at most n/2, since no two LMS positions are adjacent.

Worked example: mississippi step by step

Take mississippi with the sentinel appended, positions 0 to 11:

Buckets: $ at slot 0, i at 1 to 4, m at 5, p at 6 to 7, s at 8 to 11. Pass one seeds LMS positions 1, 4, 7 and 11 at bucket tails in text order, then runs both scans:

The LMS positions now appear in substring order: 11 ($), 7 (ippi$), then 1 and 4, both issi. Pass one sorts LMS substrings, not suffixes: 1 and 4 tie, and slots 3 to 4, 8 to 9 and 10 to 11 still hold their pairs in the wrong order.

Naming gives $ the name 0, ippi$ 1 and issi 2. In text order of positions 1, 4, 7, 11 the reduced string is [2, 2, 1, 0]. Names repeat, so we recurse. Its suffix array [3, 2, 1, 0] maps back to LMS order 11, 7, 4, 1, breaking the tie because issippi$ is smaller than ississippi$.

Dropping slot 0 gives [10, 7, 4, 1, 0, 9, 8, 6, 3, 5, 2], matching the oracle. The final S scan changed nothing only because every S-type suffix here is already a seeded LMS suffix.

The complete SA-IS implementation

The full algorithm favours clarity over memory; production code reuses the output array for nearly all working space.

def sais(t, k):
    """t: ints in [1, k) ending with a unique sentinel 0."""
    n = len(t)
    if n == 1:
        return [0]
    s = [False] * n                    # 1. types: True = S
    s[-1] = True
    for i in range(n - 2, -1, -1):
        s[i] = t[i] < t[i + 1] or (t[i] == t[i + 1] and s[i + 1])
    lms = lambda i: i > 0 and s[i] and not s[i - 1]
    cnt = [0] * k
    for ch in t:
        cnt[ch] += 1
    def bucket(tail):                  # start (or last slot) of each bucket
        b, acc = [0] * k, 0
        for ch in range(k):
            acc += cnt[ch]
            b[ch] = acc - 1 if tail else acc - cnt[ch]
        return b

    def induce(seeds):
        sa, b = [-1] * n, bucket(True)
        for i in reversed(seeds):      # 2a. seeds at bucket tails
            sa[b[t[i]]] = i; b[t[i]] -= 1
        b = bucket(False)
        for j in range(n):             # 2b. L-types, left to right
            i = sa[j] - 1
            if sa[j] > 0 and not s[i]:
                sa[b[t[i]]] = i; b[t[i]] += 1
        b = bucket(True)
        for j in range(n - 1, -1, -1): # 2c. S-types, right to left
            i = sa[j] - 1
            if sa[j] > 0 and s[i]:
                sa[b[t[i]]] = i; b[t[i]] -= 1
        return sa

    def same(a, b):                    # equal LMS substrings?
        if n - 1 in (a, b):
            return a == b
        j = 0
        while t[a + j] == t[b + j] and s[a + j] == s[b + j]:
            j += 1
            if lms(a + j) or lms(b + j):
                return lms(a + j) and lms(b + j)
        return False

    pos = [i for i in range(n) if lms(i)]
    names, name, prev = {}, -1, None   # 3. name sorted LMS substrings
    for i in induce(pos):
        if lms(i):
            if prev is None or not same(prev, i):
                name += 1
            names[i], prev = name, i
    reduced = [names[i] for i in pos]
    if name + 1 < len(pos):            # 4. recurse only if names repeat
        order = sais(reduced, name + 1)
    else:
        order = [0] * len(pos)         # names are unique: invert them
        for idx, nm in enumerate(reduced):
            order[nm] = idx
    return induce([pos[j] for j in order])   # 5. final pass

def suffix_array(text):
    t = [ord(ch) + 1 for ch in text] + [0]
    return sais(t, max(t) + 1)[1:]     # slot 0 holds the sentinel

LMS substrings are compared on symbol and type, since equal symbols with different types sort differently; the comparisons are linear overall because each position belongs to at most two LMS substrings. The reduced string ends with the sentinel's unique name 0, so the recursive call meets the same contract as the top-level one.

Building at scale: libraries, memory and cache

Python is right for understanding SA-IS and wrong for real data. In production, wrap a maintained C library and keep your own builder as a test double. Common choices are Yuta Mori's libdivsufsort, long the practical speed reference, and Ilya Grebnov's libsais, an induced-sorting implementation with optional OpenMP parallelism and 64-bit variants. Check each README for supported index widths.

Plan memory first. The output is n integers: with 32-bit entries, 4 bytes per symbol plus the text, and signed 32-bit indices stop near 2.1 billion symbols. Beyond that, 64-bit entries need 8n bytes, about 24 GB for a 3-billion-base genome. Speed is dominated by cache misses: induce scans read sequentially but write through bucket pointers that jump across the array, so libraries prefetch upcoming targets. When the array does not fit in RAM, external-memory variants such as eSAIS (Bingmann, Fischer and Osipov) adapt induced sorting to disk.

Testing a builder

1. Differential tests against the oracle on thousands of random strings of length 1 to 60 over tiny alphabets (a, ab, acgt), which force long repeats and deep recursion. Include length 1, a single repeated symbol, and periodic strings.
2. Cross-check two fast builders (SA-IS against doubling or a library) on inputs too large for the oracle.
3. Invariants on full-size output: the array is a permutation, and each adjacent pair is in order. Build the LCP array with Kasai's algorithm and compare only the first differing symbol, which keeps the check linear; the Kasai LCP article explains the construction.

Failure modes

• Missing or non-unique sentinel. Types are mislabelled and the output is a plausible but wrong permutation. Remap the input so 0 is free.
• Comparing LMS substrings by symbols only. Passes most random tests and fails where equal symbols carry different types.
• Signed 32-bit overflow past about 2.1 billion symbols, showing up as negative indices late in a long job. Choose the width from n before allocating.
• Assuming doubling's typical speed. Logs, versioned documents and genome collections are repetitive and need every round.
• Mismatched symbol definitions. UTF-8 bytes and code points sort differently; the builder and the search side must agree.

Trade-offs

For one pattern against one text, a Z-array or KMP table needs no index at all. When many queries hit the same text, build the array once.

What to do next

1. Write the oracle and a differential test harness before any builder.
2. Implement doubling with radix passes and pass 10,000 random small-alphabet strings.
3. Implement SA-IS and compare your mississippi trace with the tables above.
4. Time both on a 10-million-symbol file, repetitive and random.
5. Choose the index width from your largest n; budget text plus 4n or 8n bytes.
6. In production, wrap a library, keep your builder as the cross-check, and run the permutation and adjacent-pair invariants after every build.