Substring Search via Suffix Array Binary Search, in depth: lower and upper bound, prefix-aware comparison, min(l, r) skipping and O(m + log n) with LCP

Why occurrences form one range

Let the text be s of length n and the pattern p of length m. A position i is an occurrence when s[i:i+m] == p, which is the same as saying p is a prefix of the suffix starting at i. Now look at the suffixes in sorted order. All strings that start with p sort together: anything smaller than p on its first m characters comes before them, anything larger comes after. So the occurrences form one contiguous run of ranks, and two numbers describe the answer completely:

• lower bound: the first rank whose suffix, cut to m characters, is not less than p;
• upper bound: the first rank whose suffix, cut to m characters, is greater than p.

The count of occurrences is upper - lower, and the positions are sa[lower:upper]. If the two bounds are equal, p does not occur. Both are monotone predicates over ranks, which is exactly what binary search needs.

Comparing a suffix with a pattern

Every search step compares one suffix with p. Three outcomes matter: the suffix's first m characters are smaller than p, p is a prefix of the suffix, or they are larger. One more case trips people up: the suffix can end before p does. A suffix that is a proper prefix of p, such as an against ana, is smaller than p, because a string sorts before any longer string it is a prefix of.

The function below also takes a starting offset k: the caller promises that the first k characters already match, and the function returns how many characters match in total. Returning that count is what the faster searches later build on.

def compare(s, i, p, k=0):
    """Compare suffix s[i:] with p, given that the first k characters already match.
    Returns (matched, sign): sign < 0 if the suffix sorts before p (on its first m chars),
    0 if p is a prefix of the suffix, > 0 if it sorts after."""
    n, m = len(s), len(p)
    while k < m and i + k < n and s[i + k] == p[k]:
        k += 1
    if k == m:
        return k, 0           # p is a prefix of this suffix
    if i + k == n:
        return k, -1          # suffix ended first: it is a proper prefix of p, so smaller
    return k, (-1 if s[i + k] < p[k] else 1)

Avoid the tempting one-liner s[sa[mid]:sa[mid] + m] < p. It is correct, but in Python it allocates a new string of m characters at every step, and it hides the match length you need for the optimisations below.

Lower and upper bound

Both bounds use the same loop. Keep two ranks L and R with the invariant that suffix(L) is on the left side of the predicate and suffix(R) on the right, starting from the imaginary ranks -1 and n. The only difference between the bounds is which side the "p is a prefix" case falls on.

def search(s, sa, p, upper=False):
    """Lower bound (upper=False) or upper bound (upper=True) of p in suffix array sa."""
    L, R = -1, len(sa)       # imaginary sentinels: suffix(-1) < p, suffix(n) > p
    l = r = 0                # characters of p matched against suffix(L) and suffix(R)
    while R - L > 1:
        M = (L + R) // 2
        k, sign = compare(s, sa[M], p, min(l, r))
        if sign < 0 or (upper and sign == 0):
            L, l = M, k      # suffix(M) is on the left side
        else:
            R, r = M, k      # suffix(M) is on the right side
    return R

def occurrences(s, sa, p):
    lo, hi = search(s, sa, p), search(s, sa, p, upper=True)
    return hi - lo, sorted(sa[lo:hi])   # count, positions in text order

The half-open, sentinel-based form avoids the off-by-one edge cases of the lo <= hi style: the loop ends when L and R are adjacent, and R is the answer. An empty pattern matches everywhere, and the code returns lower bound 0 and upper bound n, as it should.

Worked example: ana in banana

Search for ana in banana, with sa = [5, 3, 1, 0, 4, 2].

The answer is ranks [1, 3), count 2, positions sorted([3, 1]) = [1, 3]. Look at the third lower-bound step: it starts comparing at offset 1, not 0, because both bounds already agreed with p on the first character. That saving is the next section.

Skipping characters already matched

Plain binary search costs O(m log n): log n steps, each comparing up to m characters. On repetitive data such as DNA, log files or source code, the suffixes near the answer share long prefixes with p, so the full m characters really are compared at nearly every step.

The fix rests on one fact about sorted strings. If suffix(L) <= suffix(M) <= suffix(R) and p lies between suffix(L) and suffix(R), then p agrees with suffix(M) on at least min(l, r) characters, where l and r are p's match lengths against the two bounds. Everything between two sorted strings shares their common prefix, and p shares at least min(l, r) characters with both. So each comparison may start at that offset, which is exactly what search above does by passing min(l, r).

This simple acceleration, due to Manber and Myers, makes most real searches close to O(m + log n). It does not change the worst case: if one bound matches a long prefix of p and the other matches nothing, min(l, r) stays at zero and the search re-reads the same characters step after step.

O(m + log n) with LCP information

To guarantee O(m + log n), use precomputed longest-common-prefix information between suffixes. Let lcp(a, b) be the common prefix length of the suffixes at ranks a and b. From the LCP array built by Kasai's algorithm, where LCP[i] compares ranks i - 1 and i, lcp(a, b) for a < b is the minimum of LCP[a+1..b], a range-minimum query. Suppose l >= r, and let q = lcp(L, M). Suffix(L) agrees with p on l characters and then falls below it:

• q > l: suffix(M) agrees with suffix(L) beyond position l, so it falls below p at the same place. Move L to M; l is unchanged. No characters read.
• q < l: suffix(M) leaves suffix(L) at position q, where suffix(L) still equals p, and it must go upward to stay sorted. So suffix(M) is above p. Move R to M with r = q. No characters read.
• q = l: compare from offset l. Every character that matches raises max(l, r) for good.

The case r > l is symmetric, using lcp(M, R). Because max(l, r) never decreases and each step reads at most one mismatching character, the total work is O(m + log n). The same routine serves both bounds.

class RangeMinLCP:
    """lcp(a, b) between ranks a < b in O(1) via a sparse table over the LCP array."""
    def __init__(self, lcp):
        self.n, self.t = len(lcp), [list(lcp)]
        j = 1
        while (1 << j) <= self.n:
            prev, h = self.t[-1], 1 << (j - 1)
            self.t.append([min(prev[i], prev[i + h]) for i in range(self.n - (1 << j) + 1)])
            j += 1

    def between(self, a, b):
        if a < 0 or b >= self.n:          # a sentinel shares nothing
            return 0
        lo, hi = a + 1, b
        j = (hi - lo + 1).bit_length() - 1
        return min(self.t[j][lo], self.t[j][hi - (1 << j) + 1])

def search_lcp(s, sa, rmq, p, upper=False):
    L, R, l, r = -1, len(sa), 0, 0
    while R - L > 1:
        M = (L + R) // 2
        if l >= r:
            q = rmq.between(L, M)
            if q > l:
                L = M; continue
            if q < l:
                R, r = M, q; continue
            k, sign = compare(s, sa[M], p, l)
        else:
            q = rmq.between(M, R)
            if q > r:
                R = M; continue
            if q < r:
                L, l = M, q; continue
            k, sign = compare(s, sa[M], p, r)
        if sign < 0 or (upper and sign == 0):
            L, l = M, k
        else:
            R, r = M, k
    return R

The sparse table costs n log n integers, which is too much for big texts. The original Manber-Myers paper instead precomputes lcp only for the (L, M) and (M, R) pairs a binary search can actually visit, about 2n values, usually called LCP-LR. The search logic is identical; only the lookup changes.

Answering real queries

• Counting is upper - lower: two searches, no matter how many occurrences.
• Reporting costs O(occ) more; sort the range only if callers need text order.
• Leftmost occurrence without sorting: a range-minimum query over the suffix array values in [lower, upper).
• Prefix buckets: a table of the rank range for every possible first two bytes (65,536 entries) replaces the first several search steps with one lookup, and those early steps are the ones that miss cache.
• Many patterns: sort them and search in order; consecutive patterns share prefixes and their ranges narrow each other's starting interval.

Large texts and testing

On a 1 GB text, each search makes about 30 steps per bound, and each step touches a random entry of the suffix array and a random place in the text: roughly 60 cache misses, or page faults if both live in memory-mapped files. Memory is the other constraint: unsigned 32-bit entries cover texts up to 4 GiB, signed ones (common in C libraries, and needed for the -1 sentinel above) only 2 GiB; beyond that you need 40- or 64-bit entries, so the array alone is 4n to 8n bytes before any LCP data.

Compare in code units, not locale order. Searching bytes is fastest, and UTF-8 byte order agrees with code point order, so a UTF-8 suffix array answers Unicode substring queries correctly, as long as pattern and text are normalised the same way.

Test every search against brute force: random strings over a two- or three-letter alphabet (small alphabets produce long shared prefixes and duplicates), random patterns including the empty one and ones longer than the text, and assertions that search and search_lcp return identical bounds.

Failure modes

• Comparing whole suffixes instead of m-character prefixes gives a valid lower bound but a wrong upper bound, because every longer suffix starting with p compares greater than p.
• Treating an ended suffix as larger. A suffix shorter than p that matches it so far is smaller.
• Skipping without the invariant. Starting at max(l, r) instead of min(l, r) skips characters that were never verified and silently returns wrong ranges on repetitive text.
• Sentinel confusion. If the array was built with an appended $, search the same text with the sentinel, and make sure no pattern contains that character.
• Integer overflow in (L + R) // 2 in fixed-width languages; use L + (R - L) / 2.

Trade-offs

What to do next

1. Implement compare and search from this page and fuzz them against str.startswith on random two-letter strings.
2. Measure characters compared per search on your own data, with and without min(l, r) skipping.
3. If patterns are long and the text repetitive, add LCP information and confirm identical bounds.
4. For multi-gigabyte texts, plan entry width, memory mapping and a prefix bucket table before writing code.
5. Keep learning: suffix arrays in depth, the Z algorithm and Rabin-Karp hashing.