Suffix array

Deep dive

How it works

In: a text of n bytes (build), then a pattern of m bytes · Out: all start offsets of the pattern (a contiguous SA range) · State: SA (4n bytes with u32), optional LCP (4n) and the text · Cost: build O(n) with SA-IS, query O(m log n)

Prefix doubling (Manber-Myers):

1. rank[i] = s[i]. k = 1.
2. Sort suffixes by the pair (rank[i], rank[i+k]), with "past the end" smallest. Now they are sorted by their first 2k bytes.
3. Re-rank: equal pairs get equal ranks. If all n ranks are distinct, stop; else k = 2k, repeat. At most ⌈log2 n⌉ rounds.

Kasai LCP: visit suffixes in text order. If suffix i shares h bytes with its SA predecessor, suffix i+1 shares at least h−1 with its own predecessor, so h decreases by at most 1 per step and the total work is O(n).

Search for p: lower = first suffix ≥ p, upper = first suffix after lower that does not start with p.

s = banana$
 r  sa  lcp  suffix
 0   5   0   a
 1   3   1   ana
 2   1   3   anana
 3   0   0   banana
 4   4   0   na
 5   2   2   nana
find "ana" -> r in 1..3 -> offsets 3, 1        longest repeat: max lcp = 3 -> "ana"

Code

/// Suffix array by prefix doubling: after round k, suffixes are sorted by their first 2k bytes.
pub fn suffix_array(s: &[u8]) -> Vec<u32> {
    let n = s.len();
    let mut rank: Vec<u32> = s.iter().map(|&c| c as u32).collect();
    let mut sa: Vec<u32> = (0..n as u32).collect();
    let mut keyed: Vec<(u64, u32)> = Vec::with_capacity(n);
    let mut k = 1;
    while n > 1 {
        // key(i) = (rank of first k bytes, rank of next k bytes + 1, or 0 past the end)
        keyed.clear();
        keyed.extend((0..n).map(|i| {
            let second = if i + k < n { rank[i + k] as u64 + 1 } else { 0 };
            (((rank[i] as u64) << 32) | second, i as u32)
        }));
        keyed.sort_unstable(); // radix sort here makes it O(n log n)
        let mut r = 0;
        for j in 0..n {
            if j > 0 && keyed[j].0 != keyed[j - 1].0 { r += 1; }
            sa[j] = keyed[j].1;
            rank[sa[j] as usize] = r; // safe to overwrite: keys are already materialized
        }
        if r as usize == n - 1 { break; } // all ranks distinct: fully sorted
        k *= 2;
    }
    sa
}

/// Kasai: lcp[r] = longest common prefix of suffixes sa[r-1] and sa[r] (lcp[0] = 0). O(n).
pub fn lcp_array(s: &[u8], sa: &[u32]) -> Vec<u32> {
    let n = s.len();
    let mut rank = vec![0u32; n];
    for (r, &i) in sa.iter().enumerate() { rank[i as usize] = r as u32; }
    let (mut lcp, mut h) = (vec![0u32; n], 0usize);
    for i in 0..n { // text order: h drops by at most 1 per step, so total work is O(n)
        let r = rank[i] as usize;
        if r == 0 { h = 0; continue; }
        let j = sa[r - 1] as usize;
        while i + h < n && j + h < n && s[i + h] == s[j + h] { h += 1; }
        lcp[r] = h as u32;
        h = h.saturating_sub(1);
    }
    lcp
}

/// Range of SA positions whose suffixes start with p; sa[range] are the match offsets.
pub fn find_all(s: &[u8], sa: &[u32], p: &[u8]) -> std::ops::Range<usize> {
    let suffix = |i: u32| &s[i as usize..];
    let lo = sa.partition_point(|&i| suffix(i) < p);
    let hi = lo + sa[lo..].partition_point(|&i| suffix(i).starts_with(p));
    lo..hi
}

/// Longest substring that occurs at least twice.
pub fn longest_repeat<'a>(s: &'a [u8], sa: &[u32], lcp: &[u32]) -> &'a [u8] {
    match (0..lcp.len()).max_by_key(|&r| lcp[r]) {
        Some(r) if lcp[r] > 0 => &s[sa[r] as usize..][..lcp[r] as usize],
        _ => &[],
    }
}

Complexity and constants

• Prefix doubling with comparison sort: O(n log² n) worst case; with radix sort per round, O(n log n). Rounds = log2(longest repeat), so typical text needs few rounds and periodic text needs about log2 n.
• Measured (code above, one Apple M-series core): 10 MB random 4-letter text builds in 1.0 s; 10 MB of abcabc… (24 rounds) in 4.2 s. LCP: 0.19 s and 0.01 s.
• Query on 10 MB random text, 12-byte patterns: 0.9-1.3 µs across runs (about 46 comparisons, each a random access into SA and text).
• Memory peak of the code: text + rank (4n) + SA (4n) + keyed pairs (16n) = 25n bytes. SA-IS needs about 5n-6n plus the text; libsais and divsufsort fit in 5n plus a small constant workspace for byte alphabets.

Efficiency

• Cost: query about 2·log2(n) suffix comparisons, 0.9-1.3 µs for 12-byte patterns on 10 MB random text (about 46 comparisons, each a random access); build by prefix doubling (the code above) 1.0 s for 10 MB random 4-letter text, 4.2 s for 10 MB of abcabc…; LCP 0.19 s (measured, one Apple M-series core). Memory: SA 4n bytes; the build code peaks at 25n, SA-IS at about 5n-6n plus the text. Bound: random memory access per comparison; construction is sorting-bound.
• Levers:
  • SA-IS construction (libsais, Go index/suffixarray) — linear time, about 5n-6n peak memory versus 25n here — any input beyond toy size
  • LCP-LR array (Manber-Myers) — O(m + log n) per query instead of O(m log n) — long patterns
  • cached top-level table of every 64th suffix's first bytes — fewer random accesses per query (order of) — hot, large indexes
  • sampled SA (every 32nd offset, FM-index style) — memory down, query time up — texts too big for 4n bytes
  • OpenMP-parallel build (libsais) — build scales with cores — multi-GB inputs
• Naive → standard → state of the art:
  • sort suffixes by full comparison · O(n² log n) on repetitive text
  • prefix doubling with radix sort · O(n log n); 1.0-4.2 s per 10 MB (measured above)
  • SA-IS (libsais) with LCP-LR or FM-index · O(n) build, O(m + log n) query

Implementation notes

• Use a library for big inputs: libsais (SA-IS, also builds BWT and LCP, OpenMP parallel) or libdivsufsort. Go's standard index/suffixarray uses SA-IS.
• Sentinel: algorithms often append a unique smallest symbol ($). Comparing slices with past-the-end as smallest, as above, gives the same order without a sentinel.
• Speed up queries: store an LCP-LR array (Manber-Myers) to skip common prefixes, getting O(m + log n) per query. Or keep a small top-level table of every 64th suffix's first bytes, which is cache-friendly like Eytzinger search.
• Large texts: offsets need u40 or u64 above 4 GB. Sampling the SA (keep every 32nd offset, the FM-index style) trades query time for memory.
• BWT from SA: bwt[r] = s[sa[r] − 1] (and the sentinel for sa[r] = 0).
• Disk-resident SA: binary search does about log2 n random reads. Cache the top levels in RAM, or use a B-tree over sampled suffixes.

Pitfalls

• u32 offsets overflow silently above 4 GiB of text.
• Prefix doubling on highly repetitive input (logs, genomes with repeats) takes log2 n rounds, each a full sort. Benchmark on real data, or use SA-IS.
• Comparing whole suffixes in a sort-based build (sort_by(|a, b| s[a..].cmp(&s[b..]))) is O(n² log n) on repetitive text.
• Kasai needs the inverse SA (another 4n bytes). Build LCP before you free it, or use the Φ-array variant.
• Unicode: the SA works on bytes. UTF-8 byte order matches code point order, so that is fine, but case-insensitive search needs normalized text.

Variants and alternatives

In production

Language libraries: Go index/suffixarray (SA-IS). Compression and differencing: zstd's dictionary builder (bundles divsufsort), bsdiff (qsufsort, Larsson-Sadakane). Genomics: BWA and Bowtie read aligners (FM-index built from suffix arrays); the Rust-Bio crate offers suffix arrays. Compressors: libsais (Ilya Grebnov, author of the bsc BWT compressor).

Implementations

Sources

• U. Manber, G. Myers, Suffix Arrays: A New Method for On-Line String Searches, SIAM Journal on Computing, 1993. paper
• T. Kasai, G. Lee, H. Arimura, S. Arikawa, K. Park, Linear-Time Longest-Common-Prefix Computation in Suffix Arrays and Its Applications, CPM, 2001. paper
• G. Nong, S. Zhang, W. H. Chan, Linear Suffix Array Construction by Almost Pure Induced-Sorting, DCC, 2009. paper · code
• N. J. Larsson, K. Sadakane, Faster Suffix Sorting, Theoretical Computer Science, 2007. paper
• P. Ferragina, G. Manzini, Opportunistic Data Structures with Applications, FOCS, 2000. paper
• The Go Authors, index/suffixarray package documentation. docs