Suffix Arrays, in depth: prefix doubling in code, SA-IS in outline, pattern search, LCP-powered queries and the BWT connection

The structure in one picture

Definition and a worked example

Take s = banana, length 6. Its suffixes start at positions 0 to 5: banana, anana, nana, ana, na, a. Sorting them lexicographically gives a, ana, anana, banana, na, nana, so the suffix array is [5, 3, 1, 0, 4, 2]. Entry SA[r] is the start position of the suffix with rank r. The inverse array, rank[i], gives the sorted position of the suffix starting at i.

Why is this useful? Every substring of s is a prefix of some suffix. Sorting the suffixes therefore puts all occurrences of any substring next to each other: every suffix beginning with "ana" sits in one contiguous block of the array. Finding a pattern becomes finding a block, and finding a block in a sorted list is binary search. That is the central idea; the rest of the article is about building the array fast and using the block structure well.

Why naive construction is too slow

The obvious construction sorts the n suffixes with a comparison sort. That takes O(n log n) comparisons, but each comparison of two suffixes can cost O(n) characters, so the worst case is O(n² log n). On text with long repeats, such as DNA, logs or a file of identical lines, the worst case is the normal case. In Python the one-liner sorted(range(n), key=lambda i: s[i:]) also copies every suffix, which is O(n²) memory. It is perfect as a test oracle and useless beyond a few tens of thousands of characters.

Prefix doubling

Prefix doubling, due to Manber and Myers, sorts suffixes by their first k characters for k = 1, 2, 4, 8 and so on. The trick is that the first 2k characters of the suffix at i are the first k characters at i followed by the first k characters at i + k. If every suffix already has a rank by its first k characters, its rank by 2k characters is determined by the pair (rank[i], rank[i+k]), with a sentinel value when i + k runs off the end. Sorting by pairs of integers is cheap, and after at most log₂ n rounds all ranks are distinct, which means the suffixes are fully sorted.

def suffix_array(s):
    n = len(s)
    if n == 0:
        return []
    rank = [ord(ch) for ch in s]
    sa = list(range(n))
    k = 1
    while True:
        key = lambda i: (rank[i], rank[i + k] if i + k < n else -1)
        sa.sort(key=key)                     # O(n log n) per round
        new = [0] * n
        for j in range(1, n):
            new[sa[j]] = new[sa[j - 1]] + (key(sa[j]) != key(sa[j - 1]))
        rank = new
        if rank[sa[-1]] == n - 1:            # every rank distinct: fully sorted
            return sa
        k *= 2

Trace banana. The initial ranks are letter codes. The first round sorts by two-character prefixes, ba, an, na, an, na and a$ (where $ marks the end), leaving ties between the two copies of an and of na. The second round compares four-character prefixes, anan against ana$ and nana against na$$, and all six ranks become distinct, so the loop ends. With a comparison sort each round costs O(n log n), giving O(n log² n) overall. Replacing the sort with two stable counting sorts, first by the second key and then by the first, makes each round O(n) and the total O(n log n).

def counting_sort(order, key, m):
    count = [0] * (m + 1)
    for i in order:
        count[key[i] + 1] += 1
    for v in range(m):
        count[v + 1] += count[v]
    out = [0] * len(order)
    for i in order:                          # stable: preserves the incoming order
        out[count[key[i]]] = i
        count[key[i]] += 1
    return out

def suffix_array_radix(s):
    n = len(s)
    if n == 0:
        return []
    alphabet = sorted(set(s))
    rank = [alphabet.index(ch) + 1 for ch in s]   # 0 is reserved for "past the end"
    k = 1
    while True:
        m = max(rank) + 1
        second = [rank[i + k] if i + k < n else 0 for i in range(n)]
        sa = counting_sort(list(range(n)), second, m)
        sa = counting_sort(sa, rank, m)
        new = [0] * n
        new[sa[0]] = 1
        for j in range(1, n):
            a, b = sa[j - 1], sa[j]
            new[b] = new[a] + ((rank[a], second[a]) != (rank[b], second[b]))
        rank = new
        if rank[sa[-1]] == n:
            return sa
        k *= 2

Both functions were checked against the naive oracle on thousands of random strings over two- and four-letter alphabets, which is where repeats, and therefore bugs, are most common. Note the early exit: on text with few repeats the loop ends after a handful of rounds.

Linear time: SA-IS and DC3

Linear-time construction is possible, and production libraries use it. DC3 (also called the skew algorithm, by Kärkkäinen and Sanders) sorts the suffixes at positions not divisible by 3 recursively on a string of one-third the length, then merges in the rest. SA-IS (Nong, Zhang and Chan) classifies each suffix as S-type or L-type depending on whether it is smaller or larger than the next suffix. It sorts a small set of "leftmost S" substrings, recursing only if they are not unique, and then places all other suffixes by induced sorting: two linear scans that fill buckets from the already-placed suffixes. SA-IS is fast in practice and memory-frugal.

Unless you are writing a library, do not implement these. Use a maintained one, such as libdivsufsort or an SA-IS implementation in C or C++, or the suffix-array routines in a succinct-data-structure library, and keep a doubling implementation and the naive oracle for tests. Memory is the binding constraint at scale: 32-bit integers cost 4 bytes per character and are enough below about two billion characters; beyond that you need 64-bit entries, doubling the array, or external-memory construction.

Searching for patterns

To find a pattern p of length m, binary search for the first suffix whose first m characters are not less than p, and then for the first whose first m characters are greater than p. Everything between is an occurrence, so the count comes for free and the positions are the array entries in that range.

def find(s, sa, p):
    m = len(p)
    lo, hi = 0, len(sa)
    while lo < hi:                           # lower bound: first prefix >= p
        mid = (lo + hi) // 2
        if s[sa[mid]:sa[mid] + m] < p:
            lo = mid + 1
        else:
            hi = mid
    start, hi = lo, len(sa)
    while lo < hi:                           # upper bound: first prefix > p
        mid = (lo + hi) // 2
        if s[sa[mid]:sa[mid] + m] <= p:
            lo = mid + 1
        else:
            hi = mid
    return sorted(sa[start:lo])

# find("banana", [5, 3, 1, 0, 4, 2], "ana") -> [1, 3]

Each comparison costs up to m characters, so a search is O(m log n). Manber and Myers showed that by keeping the longest common prefix with the current bounds, and with precomputed LCP information, the search drops to O(m + log n). Simpler implementations get most of that benefit by skipping the prefix already matched against both bounds. Unlike KMP, which scans the whole text per pattern, a suffix array is built once and then answers many patterns in time that barely depends on the text length.

The LCP array and what it unlocks

The LCP array stores, for each rank r, the length of the longest common prefix between the suffixes at ranks r - 1 and r. Kasai's algorithm builds it in O(n) from the suffix array and its inverse. For banana it is [0, 1, 3, 0, 0, 2]: "a" and "ana" share 1 character, "ana" and "anana" share 3, "na" and "nana" share 2. Many classic problems become one pass over these two arrays.

• Distinct substrings. Each suffix contributes its length minus its LCP with the previous suffix, so the count is n(n+1)/2 minus the sum of the LCP array: 21 - 6 = 15 for banana.
• Longest repeated substring. The maximum LCP value and its position: 3, giving "ana".
• Longest common substring of two strings. Build the array for a + '#' + b with a separator in neither string, and take the maximum LCP between adjacent suffixes that start on opposite sides of the separator.
• Substrings repeated at least k times. A sliding window of k - 1 consecutive LCP values; the window minimum is a substring occurring k times.

Suffix arrays and their relatives

The connection to the Burrows-Wheeler transform is direct. With a unique terminator appended, the BWT character at rank r is s[SA[r] - 1], the character before each sorted suffix. Constructing a suffix array is therefore how most BWT and FM-index builders work, which is why the structure sits under read aligners in genomics and under block-sorting compressors.

Where suffix arrays are used

Beyond textbook problems, suffix arrays earn their keep in a few places. Genome read aligners build FM-indexes from a suffix array of a reference genome. Code-search and plagiarism tools find long shared substrings across files by concatenating them with separators. Machine-learning data pipelines use them for exact-substring deduplication: build a suffix array over a tokenised corpus and the LCP array reveals every span repeated beyond a threshold, which is how some language-model training sets have been cleaned of near-verbatim duplicates that encourage memorisation. Log analytics can use them to find recurring message templates.

A worked example: deduplicating 200 MB of support transcripts before fine-tuning. Concatenate the documents as bytes with a separator byte that never occurs in them, build the array with a C library (about 800 MB of 32-bit integers), compute LCP, and scan for runs of adjacent suffixes whose LCP exceeds a few hundred bytes. Each run identifies a repeated span and the documents containing it. Remove or mask every copy except one, then re-check counts. On a single machine this is minutes of work; a hash-of-shingles approach would miss repeats that do not align to shingle boundaries.

Failure modes

• Missing sentinel handling. Comparing past the end of the string, or using a sentinel that also appears in the text, silently misorders suffixes.
• Separator collisions. In multi-string arrays the separator must not occur in the inputs, and LCPs must not run across it.
• Bytes versus characters. Building on UTF-8 bytes and searching with code points, or the reverse, gives wrong offsets. Pick one representation.
• Integer overflow. 32-bit arrays fail past about two billion entries.
• Python at scale. Lists of Python ints cost far more than 4 bytes each; use NumPy arrays or a compiled library beyond a few million characters.
• Assuming dynamic updates. Suffix arrays are static; appending text means rebuilding, or using batching and merging.

Trade-offs

Choose a suffix array when the text is static, large and queried many times, and memory matters. Choose a suffix automaton or tree when text arrives online or queries need tree structure. Choose an FM-index when even n integers are too many. For a single pattern scan, classic string matching is simpler and needs no index. For prefix queries over a set of separate keys rather than substrings of one text, a trie fits better.

What to do next

1. Implement the naive oracle and the doubling algorithm, and test them against each other on random strings over a two-letter alphabet.
2. Add the two-bound binary search and check its results against a brute-force scan that tests s.startswith(p, i) at every position i.
3. Add Kasai's LCP, then solve distinct substrings and longest repeated substring.
4. Solve longest common substring of two strings with a separator.
5. Build the BWT from your array and invert it to check correctness.
6. Switch to a compiled library for real data, and measure memory against your 4n-byte estimate.