Suffix arrays and trees

What this chapter covers

The last three chapters searched a text from scratch every single time. That's the right model when you scan a text once, but it's wasteful when you'll query the same text over and over — a genome you'll search for thousands of markers, a book you'll grep repeatedly, a log you'll probe all day. For that, you preprocess the text once into a searchable index, and then every query is fast. A suffix array is the classic such index, and its idea is elegant: sort all the suffixes of the text. Because every substring is the prefix of some suffix, the occurrences of any pattern form one contiguous block of the sorted suffixes — a block you find with binary search in O(m log n). This chapter builds the suffix array with prefix doubling, searches it, computes the companion LCP array that powers its richer uses, and shows the payoff: an index that costs 12 ms to build repays itself after just 50 queries.

A bit of history

Suffix arrays were introduced in 1990 by Udi Manber and Gene Myers as a space-efficient alternative to the suffix tree, a structure Peter Weiner had invented in 1973 (Donald Knuth called Weiner's linear-time construction "the algorithm of the year"). Suffix trees are powerful — they answer many string queries in optimal time — but they're memory-hungry, needing 15–20 bytes per character in practice, which made them impractical for the huge texts of genomics. Manber and Myers showed you could get most of a suffix tree's power from a plain sorted array of suffix start positions plus the LCP array, at a fraction of the memory (4 bytes per character). That space win is why suffix arrays, not suffix trees, became the backbone of DNA search. The build algorithms then improved for two decades: Manber-Myers was O(n log n), and in 2003 three groups independently found true linear-time O(n) constructions (the DC3/skew algorithm and others), with the practical SA-IS algorithm arriving in 2009. This chapter's prefix doubling is the Manber-Myers idea — the clearest to understand, and fast enough to see the point.

The intuition

Take the text and list all its suffixes — for banana, those are banana, anana, nana, ana, na, a. Now sort them alphabetically: a, ana, anana, banana, na, nana. The suffix array just records where each sorted suffix started in the original text — here [5, 3, 1, 0, 4, 2] — so it's n integers, not n strings. Here's why that sorted order is magical for search: every occurrence of a pattern in the text is the start of some suffix that begins with the pattern. And in sorted order, all the suffixes beginning with a given prefix sit together in one contiguous block. So to find every place ana occurs, you binary-search the sorted suffixes for the block that starts with ana — you find suffixes ana (position 3) and anana (position 1), a block of two, giving both occurrences. Two binary searches (for the block's start and end) locate all matches in O(m log n).

The clever part is building the sorted order without the O(n²) cost of comparing whole suffixes. Prefix doubling does it in rounds. Round one: sort the suffixes by their first character. Round two: sort by their first two characters — but you don't re-read them, you reuse round one's ranking, because the first two characters of suffix i are (rank of i by 1 char, rank of i+1 by 1 char). Round three: first four characters, from two overlapping two-character ranks. Each round doubles how many characters the order accounts for, so after log n rounds the suffixes are fully sorted, and each round is a single sort of integer pairs. That's O(n log² n) — the log n rounds times the sort — and it never compares raw strings after the first round.

Complexity: how it scales

Building the suffix array by prefix doubling is $O(n \log^2 n)$ (log n rounds, each an $O(n \log n)$ sort); the specialized SA-IS algorithm does it in true $O(n)$. Storage is $O(n)$ integers. Each query is $O(m \log n)$: two binary searches of $\log n$ steps, each comparing up to $m$ characters. The LCP array adds $O(n)$ via Kasai's algorithm. The whole value proposition is amortization — you pay the build once and every subsequent query is cheap — so the face-off measures cumulative time to answer $K$ queries, the suffix array (build once, then query) against rescanning the text for each query:

The suffix array starts behind — at a single query it's already spent 12 ms building the index while the rescan just scans once. But its per-query cost is tiny, so its line barely rises, while the rescan line climbs linearly: every query re-reads the whole 8000-character text. The two cross at about 50 queries, and past that the index pulls steadily ahead — by 2000 queries the rescan has re-read the text two thousand times while the index read it once. That crossover is the entire logic of indexing: preprocessing is a fixed cost you amortize, worth it exactly when you'll query enough times to repay the build. Search a text once, rescan; search it fifty times or a million, build the index.

What it's good at, what it isn't

Suffix arrays are the right tool when you query one fixed text repeatedly, or need whole-text structural queries. Genomics is the flagship: a reference genome is indexed once and searched for millions of reads (the FM-index behind BWA and Bowtie is a compressed suffix array). Full-text search engines, plagiarism detectors, and data-compression tools (the Burrows-Wheeler Transform, the heart of bzip2, is a sorted rotation closely related to the suffix array) all lean on them. And the suffix-array-plus-LCP combination answers longest-repeated-substring, distinct-substring-count, and longest-common-substring queries that single-pass matchers can't touch. When the text is fixed and the questions are many or structural, this is the tool.

Where it's the wrong choice is a one-shot search — if you'll scan a text only once, building an index is pure overhead; use Boyer-Moore or str.find and move on. It's also poorly suited to frequently changing text: the suffix array must be rebuilt when the text changes (there's no cheap incremental update), so it fits static or rarely-updated corpora, not live-edited documents. And the pure-Python $O(n \log^2 n)$ build here is far slower than a compiled SA-IS; for real genomic-scale texts you'd use pydivsufsort or a BWA-style index, not hand-rolled Python. The suffix array is an index — worth building only when you'll query enough to amortize it.

The data, or the inputs

The face-off fixes an 8000-character random text and answers a growing number of six-character queries two ways: building the suffix array once then binary-searching per query, versus rescanning the whole text on each query, timing both in Python so the comparison is fair. Correctness is checked hard — on hundreds of random strings the prefix-doubling suffix array must exactly equal the naive sort-the-suffixes reference, search_all must match Python's str.find for many patterns, and the LCP array must equal a brute-force longest-common-prefix of adjacent sorted suffixes. The animation binary-searches the six sorted suffixes of banana for the block starting with ana.

Build it, one function at a time

Building the suffix array by prefix doubling — rank on 1, then 2, then 4, … characters:

def build_suffix_array(s):
    """The suffix array: the start indices of all suffixes of `s`, sorted by the suffixes they
    begin. Built by PREFIX DOUBLING — sort suffixes by their first character, then by their first
    2, then 4, 8, ... reusing the previous round's ranks so each doubling is a single sort of
    (rank, next-rank) pairs. O(n log^2 n). After log n rounds every suffix has a unique rank and
    the array is fully sorted."""
    n = len(s)
    sa = list(range(n))
    rank = [ord(c) for c in s]                   # round 0: rank by first character
    tmp = [0] * n
    k = 1
    while True:
        def key(i):
            return (rank[i], rank[i + k] if i + k < n else -1)
        sa.sort(key=key)                          # sort suffixes by (rank, rank k ahead)
        tmp[sa[0]] = 0
        for i in range(1, n):
            tmp[sa[i]] = tmp[sa[i - 1]] + (key(sa[i]) != key(sa[i - 1]))
        rank = tmp[:]
        if rank[sa[-1]] == n - 1:                 # all ranks distinct → fully sorted
            break
        k *= 2
    return sa

Searching it — two binary searches bracket the contiguous block of matches:

def search_all(s, sa, pattern):
    """Every start position of `pattern` in `s`, via two binary searches on the suffix array.
    The occurrences form a contiguous range [lo, hi) of the sorted suffixes — those whose prefix
    is `pattern`. Each comparison looks at up to m characters, and there are log n of them, so
    this is O(m log n). Returns the match positions (sorted)."""
    n, m = len(s), len(pattern)
    if m == 0:
        return sorted(range(n + 1))

    def lower(bound_pattern):
        lo, hi = 0, n
        while lo < hi:
            mid = (lo + hi) // 2
            if s[sa[mid]:sa[mid] + m] < bound_pattern:
                lo = mid + 1
            else:
                hi = mid
        return lo

    start = lower(pattern)                        # first suffix whose prefix ≥ pattern
    # first suffix whose m-prefix is strictly greater than pattern
    lo, hi = start, n
    while lo < hi:
        mid = (lo + hi) // 2
        if s[sa[mid]:sa[mid] + m] <= pattern:
            lo = mid + 1
        else:
            hi = mid
    end = lo
    return sorted(sa[start:end])

The LCP array via Kasai's algorithm — the companion that unlocks the structural queries:

def kasai_lcp(s, sa):
    """The LCP array: lcp[i] is the length of the longest common prefix of the suffixes ranked i
    and i-1 in the suffix array (adjacent in sorted order). Kasai's algorithm computes it in O(n)
    by scanning suffixes in TEXT order and reusing the fact that dropping the first character can
    lower the LCP by at most one. LCP + suffix array together answer substring-count, longest-
    repeated-substring, and longest-common-substring queries."""
    n = len(s)
    rank = [0] * n
    for i, p in enumerate(sa):
        rank[p] = i
    lcp = [0] * n
    h = 0
    for i in range(n):
        if rank[i] > 0:
            j = sa[rank[i] - 1]
            while i + h < n and j + h < n and s[i + h] == s[j + h]:
                h += 1
            lcp[rank[i]] = h
            if h > 0:
                h -= 1                             # next suffix drops one leading char → LCP falls ≤ 1
        else:
            h = 0
    return lcp

Watch it work

Here are the six suffixes of banana, sorted: a, ana, anana, banana, na, nana. To find every occurrence of ana, we binary-search this sorted list for the block that starts with ana. Orange is the suffix being compared at the midpoint; the darker band is the range still in play. Watch it narrow: the midpoint suffix's first three characters are compared against ana, and the search moves left or right just like binary search on numbers — because sorted suffixes are sorted, alphabetically. It converges on a contiguous block of two suffixes, ana and anana, whose start positions (3 and 1) are exactly where ana occurs in banana. Every match, found by two binary searches over an index built once:

The complete code

The from-scratch tab is the suffix array — prefix-doubling build, binary-search query, and Kasai's LCP; the library tab is the naive suffix array (the correctness reference) and the rescan-each-query contrast, with the pydivsufsort call you'd use in production noted. Flip between them.

"""Suffix arrays — index the TEXT once, then find any pattern fast. The last three chapters
searched a text from scratch every time. But if you'll query the same text many times — a book,
a genome, a log file — it pays to preprocess it into a searchable index once and answer every
later query in a flash. A suffix array is that index, and it's beautifully simple: the sorted
order of all the text's suffixes.

Every substring of a text is a prefix of one of its suffixes, so if you sort all n suffixes, the
occurrences of any pattern land in a CONTIGUOUS block of that sorted order — the suffixes that
start with the pattern. Find that block with binary search and you've found every occurrence, in
O(m log n), no matter how many times you search. The suffix array stores just the n starting
positions in sorted-suffix order (n integers, not the suffixes themselves), so it's compact. This
chapter builds it with prefix doubling, searches it by binary search, and computes the companion
LCP array that unlocks its richer uses.
"""


# region: build
def build_suffix_array(s):
    """The suffix array: the start indices of all suffixes of `s`, sorted by the suffixes they
    begin. Built by PREFIX DOUBLING — sort suffixes by their first character, then by their first
    2, then 4, 8, ... reusing the previous round's ranks so each doubling is a single sort of
    (rank, next-rank) pairs. O(n log^2 n). After log n rounds every suffix has a unique rank and
    the array is fully sorted."""
    n = len(s)
    sa = list(range(n))
    rank = [ord(c) for c in s]                   # round 0: rank by first character
    tmp = [0] * n
    k = 1
    while True:
        def key(i):
            return (rank[i], rank[i + k] if i + k < n else -1)
        sa.sort(key=key)                          # sort suffixes by (rank, rank k ahead)
        tmp[sa[0]] = 0
        for i in range(1, n):
            tmp[sa[i]] = tmp[sa[i - 1]] + (key(sa[i]) != key(sa[i - 1]))
        rank = tmp[:]
        if rank[sa[-1]] == n - 1:                 # all ranks distinct → fully sorted
            break
        k *= 2
    return sa
# endregion


# region: search
def search_all(s, sa, pattern):
    """Every start position of `pattern` in `s`, via two binary searches on the suffix array.
    The occurrences form a contiguous range [lo, hi) of the sorted suffixes — those whose prefix
    is `pattern`. Each comparison looks at up to m characters, and there are log n of them, so
    this is O(m log n). Returns the match positions (sorted)."""
    n, m = len(s), len(pattern)
    if m == 0:
        return sorted(range(n + 1))

    def lower(bound_pattern):
        lo, hi = 0, n
        while lo < hi:
            mid = (lo + hi) // 2
            if s[sa[mid]:sa[mid] + m] < bound_pattern:
                lo = mid + 1
            else:
                hi = mid
        return lo

    start = lower(pattern)                        # first suffix whose prefix ≥ pattern
    # first suffix whose m-prefix is strictly greater than pattern
    lo, hi = start, n
    while lo < hi:
        mid = (lo + hi) // 2
        if s[sa[mid]:sa[mid] + m] <= pattern:
            lo = mid + 1
        else:
            hi = mid
    end = lo
    return sorted(sa[start:end])
# endregion


# region: lcp
def kasai_lcp(s, sa):
    """The LCP array: lcp[i] is the length of the longest common prefix of the suffixes ranked i
    and i-1 in the suffix array (adjacent in sorted order). Kasai's algorithm computes it in O(n)
    by scanning suffixes in TEXT order and reusing the fact that dropping the first character can
    lower the LCP by at most one. LCP + suffix array together answer substring-count, longest-
    repeated-substring, and longest-common-substring queries."""
    n = len(s)
    rank = [0] * n
    for i, p in enumerate(sa):
        rank[p] = i
    lcp = [0] * n
    h = 0
    for i in range(n):
        if rank[i] > 0:
            j = sa[rank[i] - 1]
            while i + h < n and j + h < n and s[i + h] == s[j + h]:
                h += 1
            lcp[rank[i]] = h
            if h > 0:
                h -= 1                             # next suffix drops one leading char → LCP falls ≤ 1
        else:
            h = 0
    return lcp
# endregion

"""The library counterpart, the reference, and the contrast. In production you'd build a suffix
array with a compiled library — `pydivsufsort` (SA-IS, true O(n)) or a suffix-automaton library —
never in pure Python for large texts:

    from pydivsufsort import divsufsort
    sa = divsufsort(text.encode())

`naive_suffix_array` below is the O(n^2 log n) reference — sort the actual suffix strings — used to
check the prefix-doubling build. And the real contrast for the *search* is Python's `str.find`:
without an index, every query rescans the whole text (O(n) per query), where the suffix array pays
O(n log n) once and then answers each query in O(m log n). The face-off times "index once, query
many" against "rescan every query."
"""


# region: naive
def naive_suffix_array(s):
    """Sort the suffixes directly by materializing them — simple and obviously correct, but
    O(n^2 log n) time and O(n^2) space. The trusted reference for the prefix-doubling build."""
    return sorted(range(len(s)), key=lambda i: s[i:])
# endregion


# region: rescan
def find_all_builtin(text, pattern):
    """All occurrences via str.find — no index, so it rescans the text each call (O(n) per query).
    Both the correctness reference and the 'no index' contrast in the face-off."""
    if not pattern:
        return sorted(range(len(text) + 1))
    out, i = [], text.find(pattern)
    while i != -1:
        out.append(i)
        i = text.find(pattern, i + 1)
    return out
# endregion

Scratch vs library

The lesson of this chapter is indexing as a strategy, and the face-off's crossover graph is its clearest statement: preprocessing trades a fixed upfront cost for cheap repeated queries, and it's worth it exactly when the number of queries exceeds the crossover. That's a decision pattern far bigger than string search — it's the same logic behind database indexes, sorting once to binary-search many times, and building a hash table before a batch of lookups. The suffix array is the string world's version: don't re-scan a text you'll interrogate a thousand times; sort its suffixes once and binary-search forever. The build details also reward study — prefix doubling's rank reuse is the same "reuse the previous round" thinking as dynamic programming, and Kasai's linear LCP is the same amortized argument as KMP. In production you'd use a compiled SA-IS (pydivsufsort) or an FM-index; building the array yourself is what turns "index the text" from an abstraction into the concrete image of sorted suffixes and a binary search converging on a block.

Where you'll actually meet it

Suffix arrays and their compressed descendants run the world's largest text searches. In genomics, read aligners like BWA and Bowtie index reference genomes with FM-indexes — compressed suffix arrays — to map billions of DNA reads, the workhorse of modern sequencing. Data compression uses the closely related Burrows-Wheeler Transform (sorted rotations) in bzip2 and in the compression inside those same genomic tools. Full-text search systems and code-search engines use suffix arrays or suffix automata for substring and regex queries over fixed corpora. Plagiarism and clone detectors find longest common substrings with them. Bioinformatics uses the LCP array for repeat-finding and sequence assembly. Anywhere a large, static text is searched or analyzed many times, a suffix array — usually compressed — is the index underneath.

Takeaways

A suffix array is the sorted order of a text's suffixes, stored as n start indices; because occurrences of any pattern form one contiguous block of that order, binary search finds all of them in O(m log n) after an O(n log n) build. Prefix doubling builds it by ranking on the first 1, 2, 4, … characters with rank reuse, and Kasai's algorithm adds the O(n) LCP array that unlocks longest-repeat, distinct-substring, and longest-common-substring queries — most of a suffix tree's power at a quarter of its memory. It's an index: worth building only when you'll query a fixed text enough to amortize the cost, which the face-off pegged at about 50 queries.

The strings tier ends with a return to matching, but for many patterns at once. Aho-Corasick builds a single automaton from a whole dictionary of patterns and finds all their occurrences in one linear pass over the text — it's KMP's failure function generalized from one pattern to a trie of thousands, and it's the engine behind virus scanners, network intrusion detection, and any tool that must spot a large fixed vocabulary in a stream.