String Algorithms and Pattern Matching Questions

The idea in plain words

Edit distance (Levenshtein distance) is the minimum number of single-character insertions, deletions and substitutions that turn string a into string b. The textbook table D[i][j] (the distance between the first i characters of a and the first j of b) has n x m cells, which for two 100,000-character strings is 10 billion cells.

The insight: an alignment is a path through that table from the top-left to the bottom-right. Every insertion or deletion moves the path one diagonal away from the main diagonal (where i = j) and costs 1. So a path of total cost at most k never leaves the band where |i - j| <= k. If you expect the answer to be small, compute only that band, about (2k + 1) cells per row instead of m, for O(k * n) time.

The recurrence, with a tiny table

Each cell is the cheapest of three moves:

• Delete a[i-1]: D[i-1][j] + 1 (the first i-1 characters of a already match the first j of b; drop one more character of a).
• Insert b[j-1]: D[i][j-1] + 1 (the first i of a match the first j-1 of b; add the missing character of b).
• Match or substitute: D[i-1][j-1] plus 0 if a[i-1] == b[j-1], else plus 1.

The first row and column are D[0][j] = j and D[i][0] = i (turning a string into an empty one, or the reverse, costs its length). By hand for kitten to sitting: D[1][1] compares k with s, a mismatch, so it is min(D[0][1] + 1, D[1][0] + 1, D[0][0] + 1) = min(2, 2, 1) = 1. D[1][2] (k against si) is min(D[0][2] + 1, D[1][1] + 1, D[0][1] + 1) = min(3, 2, 2) = 2. D[2][2] (ki against si) has i == i, so it is D[1][1] + 0 = 1. The program prints the whole table; the band |i - j| <= 3 is the part that is not a dot.

The brute force and why it is not enough

The full table costs O(n * m) time. It can be done in O(min(n, m)) memory by keeping two rows, but the time is still quadratic. When the true distance d is small (two versions of the same file, two reads of the same genome region), almost all of that work fills in cells far from the diagonal that cannot be on an optimal path.

The catch: you do not know d in advance

The band width k must be at least the true distance, otherwise the band may cut off the optimal path. Two facts handle that:

• If abs(len(a) - len(b)) > k the answer is certainly more than k (each length difference needs at least one insertion or deletion), so such a band is rejected immediately.
• Inside a band of half-width k, a returned value that is at most k is exact, and a value above k means "the distance is more than k". Any value above k is clamped to k + 1, which also keeps the numbers small.

So double the band: try k = 1, 2, 4, 8, ..., stop at the first k whose result is at most k. The pass that succeeds has k < 2d (for d >= 1), because the previous pass with k/2 failed. For d = 12 the bands tried are 1, 2, 4, 8 (all too small) and 16 (which is below 2 * 12 = 24). Each pass costs O(k * n), so the costs are n, 2n, 4n, 8n, 16n: a geometric series, where each term is double the one before, so the whole sum is less than twice the last term. The total is dominated by the last pass: O(d * n) overall, within a small constant of one pass with the perfect k.

What to store

Only the band. For each row i keep the cells with i - k <= j <= i + k. To get the distance you need only the previous row, so memory is O(k). If you also need the alignment itself (which edits to make), keep all rows of the band for O(k * n) memory, or recompute with Hirschberg's divide-and-conquer (find which cell of the middle row the best path passes through using only two rows, then solve the two halves separately) to get O(k) memory at about twice the time.

Code

def full_edit_distance(a, b):
    prev = list(range(len(b) + 1))
    for i in range(1, len(a) + 1):
        cur = [i] + [0] * len(b)
        for j in range(1, len(b) + 1):
            cur[j] = min(prev[j] + 1, cur[j - 1] + 1, prev[j - 1] + (a[i - 1] != b[j - 1]))
        prev = cur
    return prev[-1]


def show_table(a, b, k):
    """Print the full table D[i][j]; cells outside the band |i - j| <= k print as '.'."""
    n, m = len(a), len(b)
    D = [[0] * (m + 1) for _ in range(n + 1)]
    for i in range(n + 1):
        for j in range(m + 1):
            if i == 0 or j == 0:
                D[i][j] = i + j
            else:
                D[i][j] = min(D[i - 1][j] + 1, D[i][j - 1] + 1, D[i - 1][j - 1] + (a[i - 1] != b[j - 1]))
    print("      " + "  ".join(["-"] + list(b)))
    for i in range(n + 1):
        row = [str(D[i][j]) if abs(i - j) <= k else "." for j in range(m + 1)]
        print(f"  {a[i - 1] if i else '-'}   " + "  ".join(row))
    edits, i, j = [], n, m
    while i or j:
        if i and j and D[i][j] == D[i - 1][j - 1] + (a[i - 1] != b[j - 1]):
            if a[i - 1] != b[j - 1]:
                edits.append(f"replace {a[i - 1]!r} with {b[j - 1]!r}")
            i, j = i - 1, j - 1
        elif i and D[i][j] == D[i - 1][j] + 1:
            edits.append(f"delete {a[i - 1]!r}")
            i -= 1
        else:
            edits.append(f"insert {b[j - 1]!r}")
            j -= 1
    print("  edits (last position first):", "; ".join(edits))


def banded_edit_distance(a, b, k):
    """Edit distance if it is <= k, else None. Only cells with |i - j| <= k are computed."""
    n, m = len(a), len(b)
    if abs(n - m) > k:
        return None, 0
    INF = k + 1                       # any value above k means "too far"
    prev = {j: j for j in range(0, min(m, k) + 1)}     # row 0
    cells = len(prev)
    for i in range(1, n + 1):
        lo, hi = max(0, i - k), min(m, i + k)
        cur = {}
        for j in range(lo, hi + 1):
            if j == 0:
                best = i
            else:
                best = min(
                    prev.get(j, INF) + 1,                                  # delete a[i-1]
                    cur.get(j - 1, INF) + 1,                               # insert b[j-1]
                    prev.get(j - 1, INF) + (a[i - 1] != b[j - 1]),         # match / substitute
                )
            cur[j] = min(best, INF)
        cells += len(cur)
        prev = cur
    d = prev.get(m, INF)
    return (d if d <= k else None), cells


def edit_distance_doubling(a, b, verbose=False):
    """Try k = 1, 2, 4, ... until the banded answer fits inside the band."""
    k, total_cells = 1, 0
    while True:
        d, cells = banded_edit_distance(a, b, k)
        total_cells += cells
        if verbose:
            print(f"  band k={k}: cells={cells}, result={d}")
        if d is not None:
            return d, k, total_cells
        k *= 2


if __name__ == "__main__":
    print("table for kitten -> sitting, cells with |i - j| <= 3 shown (band k = 3):")
    show_table("kitten", "sitting", 3)
    print(banded_edit_distance("kitten", "sitting", 3))
    print(banded_edit_distance("kitten", "sitting", 2))
    print(banded_edit_distance("ab", "abcdef", 2))      # length gap 4 > 2: rejected before any cell is computed
    print(edit_distance_doubling("", "abcde"))          # bands 1, 2, 4 are rejected for free; k = 8 answers
    import random
    random.seed(1)
    n = 2000
    a = "".join(random.choice("ACGT") for _ in range(n))
    b = list(a)
    for pos in random.sample(range(n), 12):
        b[pos] = random.choice("ACGT")
    del b[500]; b.insert(1200, "G")
    b = "".join(b)
    d, k, cells = edit_distance_doubling(a, b, verbose=True)
    print("distance", d, "final band", k, "cells computed", cells, "full table cells", (n + 1) * (len(b) + 1))
    print("matches full DP:", d == full_edit_distance(a, b))

Run in a python:3.12 container, it prints:

table for kitten -> sitting, cells with |i - j| <= 3 shown (band k = 3):
      -  s  i  t  t  i  n  g
  -   0  1  2  3  .  .  .  .
  k   1  1  2  3  4  .  .  .
  i   2  2  1  2  3  4  .  .
  t   3  3  2  1  2  3  4  .
  t   .  4  3  2  1  2  3  4
  e   .  .  4  3  2  2  3  4
  n   .  .  .  4  3  3  2  3
  edits (last position first): insert 'g'; replace 'e' with 'i'; replace 'k' with 's'
(3, 40)
(None, 31)
(None, 0)
(5, 8, 6)
  band k=1: cells=6001, result=None
  band k=2: cells=9999, result=None
  band k=4: cells=17989, result=None
  band k=8: cells=33945, result=None
  band k=16: cells=65761, result=12
distance 12 final band 16 cells computed 133695 full table cells 4004001
matches full DP: True

Worked trace

• "kitten" to "sitting" has distance 3, and the bottom-right cell of the printed table is 3. The edit list read from the table is: replace k with s, replace e with i, insert g. The band with k = 3 is every cell where the row and column numbers differ by at most 3: the printed table shows 4, 5, 6, 7, 7, 6 and 5 cells in its seven rows, 40 in total, against the 56 cells of the full 7 by 8 table. With band k = 3 the function returns 3 after computing those 40 cells. With k = 2 it returns None after 31 cells, correctly reporting "more than 2". The next two printed lines show the length-difference check: "ab" against "abcdef" with k = 2 returns (None, 0) with no cell computed, and "" against "abcde" is rejected for k = 1, 2 and 4 at no cost, then answers 5 with k = 8 after 6 cells.
• The 2,000-character DNA-like pair is built by applying up to 12 random substitutions (a random pick can repeat the same letter), one deletion and one insertion; the program reports the actual distance, 12. The doubling loop tries bands 1, 2, 4 and 8, each returning None, then band 16 returns 12. The cell counts per pass are 6,001; 9,999; 17,989; 33,945 and 65,761, adding up to 133,695.
• The full table would have 4,004,001 cells, so the banded search computed about 30 times fewer (4,004,001 / 133,695 = 29.9). The final line confirms the answer equals the full dynamic-programming result.

Complexity

• Time: each pass with half-width k computes at most 2k + 1 cells per row, so O(k * n). With doubling the passes cost O(1 * n), O(2 * n), O(4 * n) and so on; the sum is below twice the last pass, and the last pass has k < 2d, so the total is O(d * n).
• Space: the code keeps two rows of at most 2k + 1 entries, so O(d) at the end (a dictionary per row here; in a faster implementation two preallocated arrays of length 2k + 1).
• Going further: a typical interview needs the banded table with doubling; the diagonal-transition methods are deeper material. Landau and Vishkin's method (O(n + k^2) with a constant-time longest-common-extension query, that is, a precomputed answer to "how far do these two suffixes agree starting here?") and Myers' O(ND) difference algorithm (which counts insertions and deletions only, not substitutions) do not fill cells one by one. For each number of edits e and each diagonal they record only the furthest row reachable, sliding along matching characters for free, which beats the O(k * n) banded approach when d is small.

Edge cases and common mistakes

• Identical strings: d = 0 is found with k = 1 on the first pass.
• One or both strings empty: the distance is the other length; the length-difference check rejects small bands until k reaches it.
• Skipping the length-difference check: the answer stays correct, because the bottom-right cell then lies outside the band and reads as k + 1, which becomes None; but the code first fills the whole band for nothing (12 cells instead of 0 for "ab" against "abcdef" with k = 2, and a full pass of about n x (2k + 1) cells on long strings). The check is an early exit, not a correctness guard.
• Not capping values: anything above k is never needed, so cap at k + 1; that also makes the treatment of out-of-band neighbours consistent.
• Reading out-of-band cells: neighbours outside the band are treated as "greater than k" (the code uses k + 1 through prev.get(j, INF)), never as 0.
• Recomputing each pass from scratch: every pass recomputes its whole band; reusing the narrower band's rows is possible but not worth the complexity because the work is a geometric series.
• Thinking an answer of None is an error: it means the distance exceeds the band and the loop should double.

Verification

The banded function was compared with the full dynamic program on 3,000 random binary-alphabet string pairs (lengths 0 to 9) for every band from 0 to 5, and the doubling wrapper was compared with the full result on the same pairs. All agreed.

Running the code

Save the listing as banded_edit_distance.py (Python 3.12, standard library only), then: python banded_edit_distance.py

The idea in plain words

Every substring is a prefix of some suffix, so counting distinct substrings means counting distinct prefixes of suffixes. Two methods count them without listing them:

• Suffix array + LCP. Sort the suffixes. The suffix at sorted position r has len prefixes, but the first lcp[r] of them (where lcp[r] is the longest common prefix with the previous sorted suffix) already appear as prefixes of that earlier suffix. So it contributes len - lcp[r] new substrings. Sum over all suffixes.
• Suffix automaton (SAM): a smallest-possible automaton (a graph of states with letter-labelled edges) that accepts exactly the substrings of the string: start at the start state, follow one edge per letter, and the walk succeeds if and only if the letters form a substring. The end positions of a substring are the indexes where its occurrences end; in banana, na ends at positions {3, 5} and so does ana, while a ends at {1, 3, 5}. Each state stands for a group of substrings that end at the same set of positions, so na and ana share a state and a has its own. Every substring has exactly one end-position set, so it belongs to exactly one state: distinct states hold disjoint sets of substrings. Inside a state the substrings are all suffixes of the longest one, with every length from length[link[v]] + 1 up to length[v], where link[v] (suffix link) points to the state of the longest proper suffix that belongs to a different group. A state therefore holds length[v] - length[link[v]] distinct substrings. Sum over states except the start state.

The printed tables below show this on banana: state 9 holds na and ana (lengths 2 to 3, link to state 5, which is a), so it adds 3 - 1 = 2, and state 8 holds nana, anana and banana (lengths 4 to 6) and adds 3.

Brute force and why it is not enough

Put every substring s[i:j] into a set. There are about n^2/2 of them, each up to n characters long, so O(n^3) time and O(n^3) memory in the worst case. It is the reference the two methods are checked against.

The code

import random


def suffix_array(s):
    return sorted(range(len(s)), key=lambda i: s[i:])


def lcp_array(s, sa):
    n = len(s)
    rank = [0] * n
    for r, i in enumerate(sa):
        rank[i] = r
    lcp, h = [0] * n, 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
            h = max(h - 1, 0)
        else:
            h = 0
    return lcp


def distinct_by_suffix_array(s):
    """Every substring is a prefix of some suffix. Suffix at sa[r] adds
    (its length - lcp[r]) prefixes that no earlier sorted suffix already has."""
    sa = suffix_array(s)
    lcp = lcp_array(s, sa)
    return sum(len(s) - sa[r] - lcp[r] for r in range(len(s)))


class SuffixAutomaton:
    def __init__(self, s, trace=False):
        self.next = [{}]        # next[state][char] -> state
        self.link = [-1]        # suffix link
        self.length = [0]       # length of the longest substring in the state
        self.added = []         # added[i]: new distinct substrings created by appending s[i]
        last = 0
        for ch in s:
            cur = self._new_state(self.length[last] + 1)
            if trace:
                print(f"read {ch!r}: new state {cur}, longest length {self.length[cur]}")
            p = last
            while p != -1 and ch not in self.next[p]:
                self.next[p][ch] = cur
                if trace:
                    print(f"  state {p} has no {ch!r} edge: add {p} -{ch}-> {cur}")
                p = self.link[p]
            if p == -1:
                self.link[cur] = 0
                if trace:
                    print(f"  walked past the start state: link[{cur}] = 0")
            else:
                q = self.next[p][ch]
                if self.length[p] + 1 == self.length[q]:
                    self.link[cur] = q
                    if trace:
                        print(f"  p={p} already has {ch!r} -> {q}, and length[{p}]+1 == length[{q}]: link[{cur}] = {q}")
                else:
                    clone = self._new_state(self.length[p] + 1)
                    self.next[clone] = dict(self.next[q])
                    self.link[clone] = self.link[q]
                    if trace:
                        print(f"  p={p} already has {ch!r} -> {q}, but length[{p}]+1={self.length[p] + 1} != length[{q}]={self.length[q]}: "
                              f"clone {q} as state {clone} (length {self.length[clone]})")
                    while p != -1 and self.next[p].get(ch) == q:
                        self.next[p][ch] = clone
                        if trace:
                            print(f"  redirect edge {p} -{ch}-> {q} to {clone}")
                        p = self.link[p]
                    self.link[q] = self.link[cur] = clone
                    if trace:
                        print(f"  link[{q}] = link[{cur}] = {clone}")
            self.added.append(self.length[cur] - self.length[self.link[cur]])
            last = cur

    def _new_state(self, length):
        self.next.append({})
        self.link.append(0)
        self.length.append(length)
        return len(self.next) - 1


def distinct_by_automaton(s, trace=False):
    """State v holds length[v] - length[link[v]] distinct substrings."""
    sam = SuffixAutomaton(s)
    total = 0
    for v in range(1, len(sam.next)):
        add = sam.length[v] - sam.length[sam.link[v]]
        if trace:
            print(f"state {v}: longest len {sam.length[v]}, link {sam.link[v]} "
                  f"(len {sam.length[sam.link[v]]}) -> adds {add}")
        total += add
    return total


def substrings_by_state(s):
    """Group every distinct substring of s by the automaton state it leads to."""
    sam = SuffixAutomaton(s)
    groups = {}
    for i in range(len(s)):
        for j in range(i + 1, len(s) + 1):
            v = 0
            for ch in s[i:j]:
                v = sam.next[v][ch]
            groups.setdefault(v, set()).add(s[i:j])
    return sam, groups


def distinct_brute(s):
    return len({s[i:j] for i in range(len(s)) for j in range(i + 1, len(s) + 1)})


if __name__ == "__main__":
    s = "banana"
    sa = suffix_array(s)
    lcp = lcp_array(s, sa)
    print("suffix array method on banana")
    for r in range(len(s)):
        suffix = s[sa[r]:]
        print(f"  {suffix:7} len={len(suffix)} lcp={lcp[r]} new={len(suffix) - lcp[r]}")
    print("  total:", distinct_by_suffix_array(s))
    print("automaton method on banana")
    print("  total:", distinct_by_automaton(s, trace=True))
    print("brute force on banana:", distinct_brute(s))
    wrong = [len(s) - sa[r] - (lcp[r + 1] if r + 1 < len(s) else 0) for r in range(len(s))]
    print("mistake, next suffix's LCP:", wrong, "total", sum(wrong))
    sam_b = SuffixAutomaton(s)
    with_start = sum(sam_b.length[v] - sam_b.length[sam_b.link[v]] for v in range(len(sam_b.next)))
    print("mistake, start state included:", with_start)
    sam, groups = substrings_by_state(s)
    print("substrings held by each state of banana")
    for v in sorted(groups):
        held = sorted(groups[v], key=len)
        print(f"  state {v}: {held} (lengths {sam.length[sam.link[v]] + 1} to {sam.length[v]}, link {sam.link[v]})")
    print("building the automaton for 'abb' step by step")
    SuffixAutomaton("abb", trace=True)
    per_prefix = []
    running = 0
    for gain in SuffixAutomaton(s).added:
        running += gain
        per_prefix.append(running)
    print("distinct substrings of each prefix of banana:", per_prefix)
    assert per_prefix == [distinct_brute(s[:i]) for i in range(1, len(s) + 1)]

    rng = random.Random(9)
    for _ in range(500):
        t = "".join(rng.choice("abc") for _ in range(rng.randint(0, 14)))
        want = distinct_brute(t)
        assert distinct_by_suffix_array(t) == want, t
        assert distinct_by_automaton(t) == want, t
    print("500 random strings: both methods match brute force")

    for t in ("", "a", "aaaa", "abcd"):
        print(repr(t), distinct_by_suffix_array(t), distinct_by_automaton(t))

How the automaton is built

Each appended letter ch creates one new state cur for the whole new prefix (its longest string). The code then walks p back along suffix links from the previous end state, which visits the suffixes of the old prefix from longest to shortest:

• While p has no edge labelled ch, add one to cur: the new letter extends that suffix into a substring that did not exist before. Then move p to its suffix link.
• If p falls off the start (p == -1), the new state's link is the start state: no shorter suffix was seen before.
• Otherwise p already has an edge on ch to some state q. If length[p] + 1 == length[q], then q holds exactly the strings that end the same way as the new suffix, and cur simply links to q.
• If not, q also contains longer strings that do not end at the new position, so one state would mix two end-position sets. The fix is a clone: a copy of q limited to the shorter strings (length length[p] + 1), with the same edges and link. The edges that pointed to q from p and its suffix-link ancestors are redirected to the clone, and both q and cur link to the clone.

The abb run below is the smallest case that needs a clone. When the second b is read, state 0 already has b -> 2, but state 2 holds ab and b (length 2) while the new suffix b has length 1, so state 2 is cloned into state 4, which holds only b.

The abb automaton has 5 states including the start state, which is 2n - 1 for n = 3, the known upper bound. Each letter adds one new state and at most one clone, and the first two letters never need a clone, so the total including the start state is at most 1 + n + (n - 2) = 2n - 1 for n >= 2.

Worked trace

On banana the SA method prints each sorted suffix with its length, LCP and new substrings; the automaton method prints one line per state; the driver then prints which actual substrings each state holds, a step-by-step build of the automaton for abb, and the per-prefix counts the automaton gives for free. Complete output of the program:

suffix array method on banana
  a       len=1 lcp=0 new=1
  ana     len=3 lcp=1 new=2
  anana   len=5 lcp=3 new=2
  banana  len=6 lcp=0 new=6
  na      len=2 lcp=0 new=2
  nana    len=4 lcp=2 new=2
  total: 15
automaton method on banana
state 1: longest len 1, link 0 (len 0) -> adds 1
state 2: longest len 2, link 5 (len 1) -> adds 1
state 3: longest len 3, link 7 (len 2) -> adds 1
state 4: longest len 4, link 9 (len 3) -> adds 1
state 5: longest len 1, link 0 (len 0) -> adds 1
state 6: longest len 5, link 7 (len 2) -> adds 3
state 7: longest len 2, link 0 (len 0) -> adds 2
state 8: longest len 6, link 9 (len 3) -> adds 3
state 9: longest len 3, link 5 (len 1) -> adds 2
  total: 15
brute force on banana: 15
mistake, next suffix's LCP: [0, 0, 5, 6, 0, 4] total 15
mistake, start state included: 12
substrings held by each state of banana
  state 1: ['b'] (lengths 1 to 1, link 0)
  state 2: ['ba'] (lengths 2 to 2, link 5)
  state 3: ['ban'] (lengths 3 to 3, link 7)
  state 4: ['bana'] (lengths 4 to 4, link 9)
  state 5: ['a'] (lengths 1 to 1, link 0)
  state 6: ['nan', 'anan', 'banan'] (lengths 3 to 5, link 7)
  state 7: ['n', 'an'] (lengths 1 to 2, link 0)
  state 8: ['nana', 'anana', 'banana'] (lengths 4 to 6, link 9)
  state 9: ['na', 'ana'] (lengths 2 to 3, link 5)
building the automaton for 'abb' step by step
read 'a': new state 1, longest length 1
  state 0 has no 'a' edge: add 0 -a-> 1
  walked past the start state: link[1] = 0
read 'b': new state 2, longest length 2
  state 1 has no 'b' edge: add 1 -b-> 2
  state 0 has no 'b' edge: add 0 -b-> 2
  walked past the start state: link[2] = 0
read 'b': new state 3, longest length 3
  state 2 has no 'b' edge: add 2 -b-> 3
  p=0 already has 'b' -> 2, but length[0]+1=1 != length[2]=2: clone 2 as state 4 (length 1)
  redirect edge 0 -b-> 2 to 4
  link[2] = link[3] = 4
distinct substrings of each prefix of banana: [1, 3, 6, 9, 12, 15]
500 random strings: both methods match brute force
'' 0 0
'a' 1 1
'aaaa' 4 4
'abcd' 10 10

Check the arithmetic: the SA contributions are 1 + 2 + 2 + 6 + 2 + 2 = 15. The same number is n(n+1)/2 = 21 (the number of start and end position pairs) minus the sum of the LCP values, 0 + 1 + 3 + 0 + 0 + 2 = 6, so 21 - 6 = 15. The automaton contributions are 1 + 1 + 1 + 1 + 1 + 3 + 2 + 3 + 2 = 15. Both equal the brute-force count.

Complexity

The rows an interviewer expects are the prefix-doubling suffix array with a comparison sort (O(n log^2 n), short to write) or with radix sort (O(n log n)), and the suffix automaton (O(n)). SA-IS and DC3 are names of linear-time suffix-array constructions; know that they exist and that they are long to implement.

Kasai's LCP step is O(n) in every row. The SAM construction is linear in the number of transitions it adds and follows; with an array of size sigma per state (sigma = alphabet size) that is O(n sigma) memory, with a hash map or sorted edge list it is O(n) memory.

Handling large alphabets

• SA: replace each distinct symbol by its rank (sort the distinct symbols once, O(n log n)). After that the text is a list of integers in 0..n-1 and the algorithms do not care how large the original alphabet was.
• SAM: a table next[state][symbol] is only viable for tiny alphabets. For Unicode or word-level symbols store transitions in a hash map per state (as the code does) or in a small sorted list; construction becomes expected O(n) with hashing or O(n log sigma) with ordered maps.

When to prefer which

• SA + LCP: the string is fixed, memory per character matters (three integer arrays and no pointers, cache-friendly), and you also want other queries from the same index (pattern search, longest repeat, longest common substring).
• SAM: the string arrives as a stream or you need the answer after every appended character (the per-prefix counts above: 1, 3, 6, 9, 12, 15). Recomputing a suffix array for every prefix would cost far more. It is also the natural choice when you need per-substring occurrence counts.
• With a hard memory limit, SA-IS plus a compact LCP wins over a SAM stored as dictionaries.

Edge cases and common mistakes

• Empty string: 0. One character: 1. aaaa: 4 (substrings a, aa, aaa, aaaa).
• All distinct characters (abcd): n(n+1)/2 = 10, since every LCP is 0.
• Mistake: counting len - lcp with the next sorted suffix's LCP instead of the previous one. The grand total does not change, because the same LCP values are summed (the program prints total 15 for the wrong pairing too), but every suffix's own share is wrong: banana's shares become [0, 0, 5, 6, 0, 4] instead of [1, 2, 2, 6, 2, 2]. Anything that uses the shares one by one, such as finding the k-th distinct substring, breaks.
• Mistake: including the start state of the automaton (length 0) in the sum. Its suffix link is -1, which Python reads as the last state, so the total for banana comes out as 12 instead of 15 (the program prints it); in C++ or Java the same index is out of bounds.
• Mistake: forgetting to copy transitions when cloning a state (dict(self.next[q])), which silently breaks the automaton.
• Mistake: assuming the answer fits in 32 bits. For n = 100,000 it can be about 5 * 10^9, which is more than a signed 32-bit integer holds in C++ or Java (use a 64-bit type). Python integers do not overflow.

Running the code

Save the program as distinct.py and run: python distinct.py   (Python 3.12, standard library only)

The idea in plain words

Two strings can look identical on screen and still be different sequences of numbers. A code point is one Unicode number (U+00E9 is the single code point for "é"). A combining mark is a code point that attaches to the one before it (U+0301 is a combining acute accent), so "é" can also be stored as "e" followed by U+0301. A UTF-16 code unit is a 16-bit slot that Java and JavaScript strings are made of; a code point above U+FFFF does not fit in one slot and is stored as two (a surrogate pair), so such a character counts as 2 in String.length(). A grapheme cluster is what a reader sees as one character: it can be one code point, a base letter plus marks, or an emoji built from several code points joined by the invisible U+200D (ZWJ, zero-width joiner). Plain substring search compares numbers, so it fails in two ways: it misses a match when the same text is stored in a different form, and it reports a match that cuts a grapheme in half (searching "cafe" inside "cafe" + U+0301 "finds" a hit that strips the accent).

The fix has three parts: normalize both sides to one canonical form, search the normalized code points with an ordinary fast substring algorithm, and accept only hits that start and end on grapheme boundaries, keeping an offset map so results point back into the original text.

The brute force and why it is not enough

Comparing raw UTF-16 code units (Java String.indexOf) or raw code points (Python str.find) is the baseline. Run on "café" stored both ways, it gives these wrong answers: the composed pattern is not found in the decomposed text, and the pattern "cafe" is found at 0 even though the text's last visible character is "é". Counting length is also unreliable: the same visible "café" has length 4 or 5 depending on the form, and a family emoji has 5 code points, 8 UTF-16 units in Java, and 1 grapheme.

Terms used below

• Canonical equivalence: different code point sequences that must render and behave identically (U+00E9 vs "e" + U+0301). Normalization forms NFC (compose where possible) and NFD (decompose fully) map equivalent sequences to one sequence.
• Compatibility equivalence: looser equivalence for things like the ligature "ffi" vs "ffi" or full-width letters. NFKC / NFKD apply it; it loses information, so use it for search keys only, never for stored text.
• Case folding: Unicode's caseless-comparison mapping. It differs from lowercasing: "ß" folds to "ss".
• Extended grapheme cluster: the Unicode text-segmentation unit for "user-perceived character". The rules are written down in UAX #29 (Unicode Standard Annex #29, the Unicode Consortium's text-segmentation specification).
• Category Mn: the Unicode general category "Mark, nonspacing", which covers combining accents such as U+0301 and also marks that carry meaning in other scripts, such as the Devanagari vowel sign U+0941; unicodedata.category(ch) returns it.
• Collation-aware search: search that follows a language's sorting rules, so that for example a and ä can count as equal at a chosen strength level; it is a different thing from normalization.

Preprocessing steps

1. Decide the match mode up front: exact canonical match, case-insensitive, or accent-insensitive. This sets the key function.
2. Split the text into grapheme clusters, remembering each cluster's start and end offset in the original string.
3. Turn each cluster into a key: NFKC, then case fold, then NFC (and for accent-insensitive mode, decompose with NFD, drop combining marks, recompose). Do the pattern the same way.
4. Concatenate the keys into one folded string and record where each cluster's key starts and ends in it. Folding can expand a cluster ("ß" becomes "ss"), so folded offsets are not original offsets.
5. Run a normal substring search on the folded string (str.find, or KMP / Rabin-Karp for many patterns).
6. Keep a hit only if it starts at a cluster start and ends at a cluster end, then translate it to original offsets.

Representation: code points or graphemes

Commit to this: search on normalized code points, validate on grapheme boundaries. Searching a list of graphemes directly sounds cleaner, but a fold that expands one grapheme into several (ligatures, "ß") makes "office" fail to match "office" when keys are compared cluster by cluster. Searching the folded code points avoids that and lets you use any fast string search; the boundary check is what restores grapheme correctness. Use raw UTF-16 units only as an index into Java strings, never as the unit of meaning. If you need only canonical equivalence and not grapheme safety (for example, a database LIKE), normalizing both sides to NFC is enough.

Libraries

For everyday work, unicodedata (Python) or java.text.Normalizer (Java) for normalization plus the \X pattern for grapheme boundaries is enough, and that is what the code below uses. A typical interview expects normalization, a grapheme-aware search and the offset mapping; ICU4J, the Aho-Corasick and suffix-array structures that matter for many patterns or repeated queries, and inverted indexes are deeper material for large-scale or locale-specific requirements.

• Python: unicodedata.normalize, unicodedata.is_normalized (3.8+), unicodedata.category, str.casefold(). The standard library has no grapheme segmentation, so use the third-party regex module, whose \X matches one extended grapheme cluster (the standard re has no \X).
• Java: java.text.Normalizer (normalize, isNormalized), String.codePoints() and codePointCount, java.text.BreakIterator.getCharacterInstance, and \X in java.util.regex.Pattern (documented as the extended grapheme cluster matcher in the Java 21 API; the program below runs it on Java 21). For real case folding, collation-aware search, and newer Unicode data than the JDK ships, use ICU4J (Normalizer2, UCharacter.foldCase, StringSearch). Java's toLowerCase(Locale.ROOT) is lowercasing, not folding.

Code (Python)

import unicodedata
import regex   # pip install regex==2024.11.6  (third-party; the standard re module has no \X)


def normalize_key(s):
    """Matching key: NFKC, case fold, then NFC, so canonically equal text gets equal code points."""
    return unicodedata.normalize("NFC", unicodedata.normalize("NFKC", s).casefold())


def strip_marks(s):
    """Accent-insensitive key for Latin-script text: decompose, drop combining marks (category Mn), recompose."""
    d = unicodedata.normalize("NFD", s)
    return unicodedata.normalize("NFC", "".join(ch for ch in d if unicodedata.category(ch) != "Mn"))


def graphemes(text):
    """Split text into extended grapheme clusters with their offsets in the ORIGINAL string."""
    return [(m.group(), m.start(), m.end()) for m in regex.finditer(r"\X", text)]


def find_all(text, pattern, key=normalize_key, verbose=False):
    """Fold each grapheme of the text, search the folded string, keep only hits that start and
    end on grapheme boundaries, and report (start, end) offsets in the original text."""
    units = graphemes(text)
    folded, unit_start, unit_end = "", {}, {}
    for g, s, e in units:
        k = key(g)
        unit_start[len(folded)] = s            # folded offset where this grapheme's key begins
        folded += k
        unit_end[len(folded)] = e              # folded offset where it ends
    needle = "".join(key(g) for g, _, _ in graphemes(pattern))
    if verbose:
        print("  folded text:", [hex(ord(c)) for c in folded], "needle:", [hex(ord(c)) for c in needle])
        print("  grapheme starts in folded text -> original offset:", unit_start)
        print("  grapheme ends in folded text -> original offset:", unit_end)
    hits, i = [], folded.find(needle) if needle else -1
    while i != -1:
        j = i + len(needle)
        if i in unit_start and j in unit_end:
            hits.append((unit_start[i], unit_end[j]))
        i = folded.find(needle, i + 1)
    return hits


composed, decomposed = "caf\u00e9", "cafe\u0301"      # e-acute as 1 code point vs e + combining acute
print(composed == decomposed, len(composed), len(decomposed))
print("composed code points:  ", [hex(ord(c)) for c in composed])
print("decomposed code points: ", [hex(ord(c)) for c in decomposed])
print("naive find of 'caf\\u00e9' in the decomposed text:", decomposed.find(composed))
print("naive find of 'cafe' in the decomposed text:", decomposed.find("cafe"), "(false hit: it cut the accent off)")
print("grapheme find of the composed pattern:", find_all(decomposed, composed, verbose=True))
print("grapheme find of 'cafe':", find_all(decomposed, "cafe"))
print("accent-insensitive find of 'cafe':", find_all(decomposed, "cafe", key=lambda g: strip_marks(normalize_key(g))))

family = "\U0001F468\u200d\U0001F469\u200d\U0001F467"   # man ZWJ woman ZWJ girl
text = "x" + family + "y"
print("code points:", len(family), "graphemes:", len(graphemes(family)))
print("naive find of the man emoji:", text.find("\U0001F468"), "| grapheme find:", find_all(text, "\U0001F468"))
print("whole family found at:", find_all(text, family))
print("ligature office:", find_all("an o\ufb03ce", "office"), "| sharp s:", find_all("Stra\u00dfe", "STRASSE"))
print("partial ligature 'off' in 'office':", find_all("o\ufb03ce", "off"))

Run in a python:3.12 container with regex 2024.11.6, it prints:

False 4 5
composed code points:   ['0x63', '0x61', '0x66', '0xe9']
decomposed code points:  ['0x63', '0x61', '0x66', '0x65', '0x301']
naive find of 'caf\u00e9' in the decomposed text: -1
naive find of 'cafe' in the decomposed text: 0 (false hit: it cut the accent off)
  folded text: ['0x63', '0x61', '0x66', '0xe9'] needle: ['0x63', '0x61', '0x66', '0xe9']
  grapheme starts in folded text -> original offset: {0: 0, 1: 1, 2: 2, 3: 3}
  grapheme ends in folded text -> original offset: {1: 1, 2: 2, 3: 3, 4: 5}
grapheme find of the composed pattern: [(0, 5)]
grapheme find of 'cafe': []
accent-insensitive find of 'cafe': [(0, 5)]
code points: 5 graphemes: 1
naive find of the man emoji: 1 | grapheme find: []
whole family found at: [(1, 6)]
ligature office: [(3, 7)] | sharp s: [(0, 6)]
partial ligature 'off' in 'office': []

Reading the offset maps printed by the verbose call. The decomposed text cafe + U+0301 has 5 code points but 4 graphemes: c, a, f, and e plus its accent as one cluster. Its folded key is 4 code points: c a f and the single composed U+00E9. unit_start records where each cluster's key begins in the folded string, mapped to where the cluster begins in the original: {0: 0, 1: 1, 2: 2, 3: 3}. unit_end records where each key ends: {1: 1, 2: 2, 3: 3, 4: 5}. The last entry is the interesting one: folded offset 4 is original offset 5, because the final cluster spans two original code points. The needle (the composed pattern, 4 folded code points) is found at folded offset 0, so j = 0 + 4 = 4; both 0 and 4 are in the maps, and the hit translates to (unit_start[0], unit_end[4]) = (0, 5). For the needle cafe the folded string never contains cafe, because the key of the last cluster is é, not e, so there is no hit at all.

Reading the output: the naive search finds nothing for the composed pattern and a false hit for "cafe"; the grapheme-aware search finds the composed pattern at (0, 5), which are offsets in the original 5-code-point text; "cafe" matches only in accent-insensitive mode; the man emoji inside the family emoji is rejected because the hit ends inside a cluster, while the whole family matches at (1, 6); "office" matches the ligature and "STRASSE" matches "Straße" because the key function folds first; and "off" does not match inside the ligature because the hit would end in the middle of one cluster.

Code (Java)

import java.text.BreakIterator;
import java.text.Normalizer;
import java.util.Locale;
import java.util.regex.Matcher;
import java.util.regex.Pattern;

public class Uni {
    static String key(String s) {
        String folded = Normalizer.normalize(s, Normalizer.Form.NFKC).toLowerCase(Locale.ROOT);
        return Normalizer.normalize(folded, Normalizer.Form.NFC);
    }

    static int graphemeCount(String s) {
        Matcher m = Pattern.compile("\\X").matcher(s);   // \X: extended grapheme cluster (Java 9+)
        int n = 0;
        while (m.find()) n++;
        return n;
    }

    public static void main(String[] args) {
        String composed = "caf\u00e9", decomposed = "cafe\u0301";
        System.out.println(composed.equals(decomposed) + " " + composed.length() + " " + decomposed.length());
        System.out.println(key(composed).equals(key(decomposed)));
        System.out.println(Normalizer.isNormalized(decomposed, Normalizer.Form.NFC));
        System.out.println("indexOf cafe in decomposed: " + decomposed.indexOf("cafe"));
        System.out.println("lowercase of Straße equals strasse: " + "Stra\u00dfe".toLowerCase(Locale.ROOT).equals("strasse"));
        String family = "\uD83D\uDC68\u200D\uD83D\uDC69\u200D\uD83D\uDC67";
        System.out.println("UTF-16 units: " + family.length()
            + ", code points: " + family.codePointCount(0, family.length())
            + ", graphemes via \\X: " + graphemeCount(family));
        BreakIterator it = BreakIterator.getCharacterInstance(Locale.ROOT);
        it.setText(family);
        int bi = 0;
        while (it.next() != BreakIterator.DONE) bi++;
        System.out.println("graphemes via BreakIterator: " + bi);
    }
}

Run in an eclipse-temurin:21-jdk container, it prints:

false 4 5
true
false
indexOf cafe in decomposed: 0
lowercase of Straße equals strasse: false
UTF-16 units: 8, code points: 5, graphemes via \X: 1
graphemes via BreakIterator: 1

The Java build of the same pipeline is key, a \X or BreakIterator pass for boundaries, and String.indexOf on the key string. The Java key lowercases, which is weaker than the Python case fold: the line lowercase of Straße equals strasse: false shows that Java would not match "STRASSE" against "Straße" with this key. For true case folding in Java, use ICU4J (UCharacter.foldCase).

Performance trade-offs

• Normalization is linear in the text. Checking is_normalized / Normalizer.isNormalized first lets already-normalized text (most of the web) skip the rewrite.
• Grapheme segmentation is linear too but with a heavier constant than a raw scan, because it consults Unicode property tables for every code point. The shipped function segments the whole text. On a large corpus, normalize and fold once at ingest, store the folded key (and the offset map if you must highlight in the original), and keep the original for display.
• Searching the folded string costs what the chosen algorithm costs: str.find for one pattern, Aho-Corasick for many. The boundary check is O(1) per hit using the two offset maps, so cost is paid per candidate, not per character.
• Memory: the offset maps are O(number of graphemes). Folding can lengthen text ("ß" to "ss").
• Index first when queries repeat: build the suffix array or inverted index over folded keys, never over raw text.

Edge cases and common mistakes

• Searching raw text and calling it done: misses canonically equal text and returns hits that split a cluster.
• Lowercasing instead of case folding ("Straße" vs "STRASSE"), and ignoring locale rules such as Turkish dotted and dotless i (in Java, always pass a locale).
• Normalizing only the pattern, or normalizing to different forms on the two sides.
• Using the accent-insensitive key outside Latin text. It drops every category Mn mark, which is right for accents on Latin letters but changes words where the marks carry meaning: Devanagari कु (ka, U+0915, plus vowel sign U+0941) becomes क (ka), and Hebrew points, Arabic vowel marks and Thai tone marks are also Mn, so the key would treat different words as equal. Restrict the mode to the scripts where the marks are optional accents.
• Using NFKC for stored data: it rewrites "ffi", superscripts and full-width characters irreversibly.
• Reporting folded offsets to the UI. Normalization changes lengths, so map back to original offsets.
• Mixing index time and query time library versions: grapheme rules and emoji sequences change between Unicode versions, so an index built under one version can disagree with a query segmented under another.
• Slicing a Java string at length() / 2 or taking charAt of a string with supplementary characters: that splits a surrogate pair.

Running the code

Python: pip install regex==2024.11.6, then python uni.py (save the Python listing as uni.py)
Java:   save the Java listing as Uni.java, then javac Uni.java && java Uni (JDK 21)

The idea in plain words

Rabin-Karp slides a window of the pattern's length over the text and compares a number (a hash, a short fingerprint of the window) instead of comparing characters. The hash can be updated in constant time when the window moves one step (a rolling hash). Different windows can share a hash (a collision), so a fingerprint match is only a suspect. Two things make the program both correct and hard to attack:

1. Verify every hit. When the fingerprints match, compare the actual characters. A false positive then costs time but can never produce a wrong answer. This is what makes the result always correct.
2. Make the hash unpredictable. If the base and modulus are fixed in the source, anyone who has read the code can craft strings that collide. Pick the base at random on every run, and use two independent (modulus, base) pairs so a collision has to happen in both at once.

Why the textbook version is not enough

The common version uses one fixed base and lets the 64-bit unsigned integer overflow, which silently means arithmetic modulo 2^64. That hash has a known break: the Thue-Morse string (a two-letter sequence whose i-th letter is a or b according to whether the number of 1 bits in i is even or odd, starting abbabaab...) of length 2^11 and its letter-swapped complement have the same hash modulo 2^64 for every odd base, so changing the base constant does not help. The section "Why a power-of-two modulus fails" below shows the mechanism with a printed table. The program below shows it: the fixed-base version reports a match at position 0 that is not a match. It also skips verification, which is why it is wrong rather than merely slow.

Choosing the parameters

The hash of a window is the polynomial c_0 * B^(m-1) + c_1 * B^(m-2) + ... + c_(m-1) taken modulo a prime p, where m is the pattern length and B is the base.

• Modulus: a prime near 10^9. Two windows with different characters give two different polynomials. Their difference, seen as a function of the base B, is a polynomial that is not identically zero and has degree at most m-1. Modulo a prime p, such a polynomial has at most m-1 roots (values of B that make it zero), the same way a quadratic has at most 2 roots; this holds because arithmetic mod a prime lets you divide by any nonzero value, which a power-of-two modulus does not. The two windows collide only if B happens to be one of those roots. With B drawn uniformly at random (every value in range equally likely), the chance of that is at most (m-1)/(p-256), because the code draws the base from the p-256 values 256 to p-1. That is within one part in a million of (m-1)/p, so the figures below use (m-1)/p. Both values stay below 2^30, so a product of two of them stays below 2^60 and never overflows uint64_t.
• Base: random, fresh per call, not a constant. The adversary may know the moduli but not B, so no input can be built against it. Bases start at 256 so every byte value is a distinct digit.
• Two independent pairs. With p1 = 1,000,000,007 and p2 = 998,244,353 and independent random bases, the false-positive chance for one pair of windows is at most (m-1)^2 / (p1 * p2): the two bases are drawn independently, so the probability that both fall on a root is the product of the two single chances, (m-1)/p1 times (m-1)/p2. For m = 2048 that is (2047/1,000,000,007) * (2047/998,244,353) = 2.05e-6 * 2.05e-6 = 4.2e-12. The union bound says the chance that at least one of several bad events happens is at most the sum of their separate chances, so over n = 1,000,000 windows the chance that any window is a false positive is at most 1,000,000 * 4.2e-12 = 4.2e-6. All three figures are printed by the second program below.
• A single modulus 2^61 - 1 (a Mersenne prime, a prime of the form 2^k - 1) with a 128-bit multiply gives (m-1)/(2^61-1) = 8.9e-16 for m = 2048 with one pair, and is a good alternative where unsigned __int128 is available. The product of two 61-bit values needs up to 122 bits, so it is computed in unsigned __int128 and reduced with %: (u64)((u128)a * b % M61), shown in the second program. unsigned __int128 is a GCC and Clang extension, not standard C++.

Because every hit is verified, these probabilities affect speed only, never correctness.

The code

#include <cstdint>
#include <iostream>
#include <random>
#include <string>
#include <vector>
using namespace std;
using u64 = uint64_t;
using u128 = unsigned __int128;

constexpr u64 MOD1 = 1'000'000'007ULL;
constexpr u64 MOD2 = 998'244'353ULL;

struct Result {
    vector<size_t> positions;
    size_t fingerprint_hits = 0;   // windows whose double hash matched
    size_t false_positives = 0;    // hits rejected by the character check
};

Result rabin_karp(const string& text, const string& pat) {
    Result r;
    size_t n = text.size(), m = pat.size();
    if (m == 0 || m > n) return r;

    // Fresh random bases on every call: the input cannot be tuned to them.
    random_device rd;
    mt19937_64 gen(((u64)rd() << 32) ^ rd());
    u64 b1 = 256 + gen() % (MOD1 - 256);
    u64 b2 = 256 + gen() % (MOD2 - 256);

    u64 h1 = 0, h2 = 0, w1 = 1, w2 = 1;  // w = base^(m-1)
    for (size_t i = 0; i < m; i++) {
        h1 = (h1 * b1 + (unsigned char)pat[i]) % MOD1;
        h2 = (h2 * b2 + (unsigned char)pat[i]) % MOD2;
        if (i + 1 < m) { w1 = w1 * b1 % MOD1; w2 = w2 * b2 % MOD2; }
    }
    u64 t1 = 0, t2 = 0;
    for (size_t i = 0; i < m; i++) {
        t1 = (t1 * b1 + (unsigned char)text[i]) % MOD1;
        t2 = (t2 * b2 + (unsigned char)text[i]) % MOD2;
    }
    for (size_t i = 0;; i++) {
        if (t1 == h1 && t2 == h2) {
            r.fingerprint_hits++;
            if (text.compare(i, m, pat) == 0) r.positions.push_back(i);
            else r.false_positives++;
        }
        if (i + m >= n) break;
        u64 out = (unsigned char)text[i], in = (unsigned char)text[i + m];
        t1 = ((t1 + MOD1 - out * w1 % MOD1) % MOD1 * b1 + in) % MOD1;
        t2 = ((t2 + MOD2 - out * w2 % MOD2) % MOD2 * b2 + in) % MOD2;
    }
    return r;
}

// Textbook version: fixed base, arithmetic mod 2^64 (unsigned overflow), no verification.
vector<size_t> naive_fixed_hash(const string& text, const string& pat) {
    vector<size_t> out;
    size_t n = text.size(), m = pat.size();
    const u64 B = 131;
    u64 w = 1, h = 0, t = 0;
    for (size_t i = 0; i < m; i++) { h = h * B + (unsigned char)pat[i]; t = t * B + (unsigned char)text[i]; if (i + 1 < m) w *= B; }
    for (size_t i = 0;; i++) {
        if (t == h) out.push_back(i);
        if (i + m >= n) break;
        t = (t - (unsigned char)text[i] * w) * B + (unsigned char)text[i + m];
    }
    return out;
}

string thue_morse(int k, bool complement) {
    string s;
    for (size_t i = 0; i < (1u << k); i++) s += ((__builtin_popcountll(i) & 1) ^ complement) ? 'b' : 'a';
    return s;
}

int main() {
    string pat = thue_morse(11, false);
    string text = thue_morse(11, true) + "xyz" + pat;  // one real match, at index 2051
    auto naive = naive_fixed_hash(text, pat);
    cout << "fixed base, mod 2^64: reported";
    for (size_t p : naive) cout << ' ' << p;
    cout << '\n';
    Result r = rabin_karp(text, pat);
    cout << "random double hash: reported";
    for (size_t p : r.positions) cout << ' ' << p;
    cout << ", fingerprint hits " << r.fingerprint_hits << ", false positives " << r.false_positives << '\n';

    string t2 = "abracadabra abracadabra", p2 = "abra";
    Result q = rabin_karp(t2, p2);
    cout << "abra in text:";
    for (size_t p : q.positions) cout << ' ' << p;
    cout << '\n';
    Result e = rabin_karp("aaaa", "aaa");
    cout << "overlap aaa in aaaa:";
    for (size_t p : e.positions) cout << ' ' << p;
    cout << '\n';

    Result z = rabin_karp("abc", ""), y = rabin_karp("ab", "abc"), w = rabin_karp("abc", "abc");
    cout << "empty pattern: " << z.positions.size() << " positions, pattern longer than text: "
         << y.positions.size() << " positions, pattern equal to text:";
    for (size_t p : w.positions) cout << ' ' << p;
    cout << '\n';
}

Output

Run in a python:3.12 container (g++ 14.2.0), compiled with g++ -std=c++17 -O2, it prints:

fixed base, mod 2^64: reported 0 2051
random double hash: reported 2051, fingerprint hits 1, false positives 0
abra in text: 0 7 12 19
overlap aaa in aaaa: 0 1
empty pattern: 0 positions, pattern longer than text: 0 positions, pattern equal to text: 0

The first line is the break of the textbook hash: the pattern is the Thue-Morse string of length 2048, the text starts with its complement, then "xyz", then the pattern, so the only real match is at 2051. The fixed-base hash also reports 0, a false positive. The randomized double hash reports 2051 only, and its single fingerprint hit was the real match.

Why a power-of-two modulus fails

Letting uint64_t overflow means every result is reduced modulo 2^64, and 2^64 is not a prime: it has no division by even numbers. Write a for 0 and b for 1 (the real letters differ by exactly 1 as well, which only scales the difference). The hash of the Thue-Morse string of length 2^k minus the hash of its complement is, up to sign, the product (B - 1)(B^2 - 1)(B^4 - 1)...(B^(2^(k-1)) - 1). For k = 2 the strings are abba and baab, and the product is (B - 1)(B^2 - 1). For any odd B every factor is even, and each squaring adds another power of 2 to the factor: B^2 - 1 is divisible by 8, B^4 - 1 by 16, B^8 - 1 by 32, and so on. Multiplying k such factors piles up more and more factors of 2, and once the product has at least 64 factors of 2 it is 0 modulo 2^64. For k = 10 the guaranteed count is 1 + 3 + 4 + ... + 11 = 64, so every odd base fails from length 1024 up, and the base is irrelevant. The second program prints, for several bases, how many trailing zero bits (factors of 2) the hash difference has; "64 (equal)" means the two hashes are identical.

A prime modulus has no such trap: a product of nonzero values mod a prime is never zero, so only the at most m-1 roots from the bound above can cause a collision.

Small numbers you can follow by hand

The first part of the next program uses digits as characters, base B = 10 and the tiny prime modulus 13, so the whole rolling hash fits on paper. Pattern 23 hashes to 23 mod 13 = 10. The first window 62 is 62 mod 13 = 10: a collision, so the character check rejects it. Sliding to the next window removes the outgoing digit 6 (its weight is 10, so subtract 6 * 10), shifts, and adds the incoming 3: (62 - 60) * 10 + 3 = 23, hash 10, a true match. Sliding again: (23 - 20) * 10 + 5 = 35, hash 9, no match. The real code does this with two large primes instead of 13, so collisions become astronomically rare. The remaining parts print the Thue-Morse table and the probability figures quoted above.

#include <cstdint>
#include <cstdio>
#include <string>
using namespace std;
using u64 = uint64_t;
using u128 = unsigned __int128;

// Part 1: a hash you can follow with a pencil. Digits are the characters, base 10, modulus 13.
void toy() {
    const int B = 10, P = 13;
    string text = "6235", pat = "23";
    int m = pat.size(), w = B;  // w = B^(m-1) = 10
    int hp = 0, t = 0;
    for (char c : pat) hp = (hp * B + (c - '0')) % P;
    for (int i = 0; i < m; i++) t = (t * B + (text[i] - '0')) % P;
    printf("pattern %s hashes to %d\n", pat.c_str(), hp);
    for (int i = 0;; i++) {
        bool same = text.compare(i, m, pat) == 0;
        printf("window %s at %d: hash %d%s\n", text.substr(i, m).c_str(), i, t,
               t == hp ? (same ? "  hash match, characters match" : "  hash match, characters differ (false positive)") : "");
        if (i + m >= (int)text.size()) break;
        int out = text[i] - '0', in = text[i + m] - '0';
        t = ((t + P - out * w % P) % P * B + in) % P;
    }
}

// Part 2: why mod 2^64 breaks. Hash of the Thue-Morse string (a = 0, b = 1 digits) and of its complement.
string tm(int k, bool comp) {
    string s;
    for (size_t i = 0; i < (1u << k); i++) s += ((__builtin_popcountll(i) & 1) ^ comp) ? 'b' : 'a';
    return s;
}
u64 poly(const string& s, u64 B) { u64 h = 0; for (char c : s) h = h * B + (unsigned char)c; return h; }

void thue() {
    printf("k=2 strings: %s and %s\n", tm(2, false).c_str(), tm(2, true).c_str());
    const u64 bases[] = {3, 131, 65537};
    for (int k = 1; k <= 11; k++) {
        string x = tm(k, false), y = tm(k, true);
        printf("k=%2d length %4d:", k, 1 << k);
        for (u64 B : bases) {
            u64 d = poly(x, B) - poly(y, B);  // difference modulo 2^64
            printf("  base %5llu: %2d trailing zero bits%s", (unsigned long long)B,
                   d == 0 ? 64 : __builtin_ctzll(d), d == 0 ? " (equal)" : "");
        }
        printf("\n");
    }
}

// Part 3: modulus 2^61 - 1 with a 128-bit product.
const u64 M61 = (1ULL << 61) - 1;
u64 mulmod61(u64 a, u64 b) { return (u64)((u128)a * b % M61); }

int main() {
    toy();
    thue();
    printf("mulmod61(2^60, 2^60) = %llu\n", (unsigned long long)mulmod61(1ULL << 60, 1ULL << 60));
    double q = 2047.0;
    printf("per pair, one prime: %.2e and %.2e; both: %.2e; over 1e6 windows: %.2e; mod 2^61-1: %.2e\n",
           q / 1000000007.0, q / 998244353.0, q / 1000000007.0 * q / 998244353.0,
           1e6 * (q / 1000000007.0) * (q / 998244353.0), q / (double)M61);
}

Output (same g++ 14.2.0 and flags as above):

pattern 23 hashes to 10
window 62 at 0: hash 10  hash match, characters differ (false positive)
window 23 at 1: hash 10  hash match, characters match
window 35 at 2: hash 9
k=2 strings: abba and baab
k= 1 length    2:  base     3:  1 trailing zero bits  base   131:  1 trailing zero bits  base 65537: 16 trailing zero bits
k= 2 length    4:  base     3:  4 trailing zero bits  base   131:  4 trailing zero bits  base 65537: 33 trailing zero bits
k= 3 length    8:  base     3:  8 trailing zero bits  base   131:  8 trailing zero bits  base 65537: 51 trailing zero bits
k= 4 length   16:  base     3: 13 trailing zero bits  base   131: 13 trailing zero bits  base 65537: 64 trailing zero bits (equal)
k= 5 length   32:  base     3: 19 trailing zero bits  base   131: 19 trailing zero bits  base 65537: 64 trailing zero bits (equal)
k= 6 length   64:  base     3: 26 trailing zero bits  base   131: 26 trailing zero bits  base 65537: 64 trailing zero bits (equal)
k= 7 length  128:  base     3: 34 trailing zero bits  base   131: 34 trailing zero bits  base 65537: 64 trailing zero bits (equal)
k= 8 length  256:  base     3: 43 trailing zero bits  base   131: 43 trailing zero bits  base 65537: 64 trailing zero bits (equal)
k= 9 length  512:  base     3: 53 trailing zero bits  base   131: 53 trailing zero bits  base 65537: 64 trailing zero bits (equal)
k=10 length 1024:  base     3: 64 trailing zero bits (equal)  base   131: 64 trailing zero bits (equal)  base 65537: 64 trailing zero bits (equal)
k=11 length 2048:  base     3: 64 trailing zero bits (equal)  base   131: 64 trailing zero bits (equal)  base 65537: 64 trailing zero bits (equal)
mulmod61(2^60, 2^60) = 576460752303423488
per pair, one prime: 2.05e-06 and 2.05e-06; both: 4.20e-12; over 1e6 windows: 4.20e-06; mod 2^61-1: 8.88e-16

Read-out: the toy section shows one false positive (window 62) and one true match. In the table, bases 3 and 131 need k = 10 to reach 64 zero bits (their counts rise 13, 19, 26, 34, 43, 53 and then reach 64), matching the factor-of-2 count above, while base 65537 = 2^16 + 1 hits 64 already at k = 4 because B - 1 itself is divisible by 2^16. The last line gives the quoted probabilities: 2.05e-6 per prime, 4.2e-12 for both, 4.2e-6 over a million windows, 8.88e-16 for the 2^61 - 1 alternative. A coding round usually needs only the verified rolling hash with a random base modulo a prime; the factor-of-2 argument and the 128-bit modulus are deeper material that interviewers rarely ask for.

Trace on a small input

For text "aaaa" and pattern "aaa" the windows are at 0 and 1. Both have the same characters, so both fingerprints equal the pattern's, both are verified by text.compare, and both are reported: 0 1. Overlapping matches are found because the window advances one character at a time.

Expected cost

• Hashing: building the pattern hash and first window hash is O(m); each slide is O(1) with two modular steps, so O(n) for the slides. Total O(n + m).
• Verification: each hit costs O(m). A true match is always verified, so the cost is O(occurrences * m). A false hit happens with probability at most about 4.2e-12 per window for m = 2048, so the expected extra work from false hits is n * 4.2e-12 * m, which is negligible.
• Expected total when matches are rare: O(n + m). Worst case: O(n * m), reached when matches are plentiful, such as text of n copies of 'a' with a pattern of m copies of 'a', because every window really matches and is verified. If the question only needs "does a match exist", stop at the first verified hit.
• Space: O(1) beyond the output list.

Edge cases

• Empty pattern or pattern longer than the text returns no positions (the code returns early, and the program prints both cases, along with a pattern equal to the whole text).
• Overlapping matches ("aaa" in "aaaa") are reported.
• Subtracting the outgoing character adds the modulus first (t + MOD - out * w % MOD) so the value never goes negative on unsigned types.
• (unsigned char) casts keep bytes above 127 from becoming negative numbers.

Common mistakes

• Relying on uint64_t overflow as the modulus (breakable as shown) and a hard-coded base.
• Seeding the random generator with a constant or the clock second: the adversary only needs to guess it.
• Skipping verification and accepting a hash match as a match.
• Multiplying two 62-bit values in uint64_t for a large modulus, which overflows.

Running the code

g++ -std=c++17 -O2 -Wall -o rabin_karp rabin_karp.cpp && ./rabin_karp

Save the listing as rabin_karp.cpp. No other files or libraries are needed.

The idea in plain words

The suffix array of a string is the list of starting positions of all its suffixes, in the order the suffixes would appear in a dictionary. For banana the suffixes are banana, anana, nana, ana, na, a (starting at 0 to 5), and sorted they are a, ana, anana, banana, na, nana, so the suffix array is [5, 3, 1, 0, 4, 2]. With it you can find any substring by binary search, and it is the compact cousin of the suffix tree.

The doubling idea: do not compare whole suffixes; compare them by their first 1, then 2, then 4, then 8 characters. Suppose you already have a rank (an integer that says how many distinct smaller keys there are, so equal keys share a rank and ranks run 0, 1, 2, ... with no gaps) for every suffix based on its first k characters. The first 2k characters of the suffix starting at i are the first k characters at i followed by the first k characters at i + k. So the pair (rank[i], rank[i + k]) is enough to order suffixes by their first 2k characters, and one round doubles the prefix length. When all ranks are different, the order is final. That takes at most about log2(n) rounds, because k reaches n.

Brute force and why it is not enough

sorted(range(n), key=lambda i: s[i:]) is short, but every key is a copy of a suffix: n(n+1)/2 characters in total (5 billion for n = 100,000), and each comparison can read up to n characters, as in aaaa...a. Time is O(n^2 log n) in the worst case and memory O(n^2).

Approach

1. Sort positions by their first character and give each distinct character a rank (0, 1, 2, ...).
2. Repeat while some ranks are tied, with k = 1, 2, 4, ...:
   • Order the positions by their second key rank[i + k]. A suffix with k or fewer characters left has no second half, so it sorts first; those are the positions n-k .. n-1. The others come in the order of the existing array shifted back by k: for each j in sa with j >= k, take j - k.
   • Do a stable counting sort by the first key rank[i]. Stable means equal keys keep their earlier order, so the second-key order survives inside each tie. Counting sort on integers below n costs O(n), no comparisons: count how many items have each key, turn the counts into starting slots with a running sum, then drop each item into the next free slot of its key.
   • Recompute ranks: walk the array in order and increase the rank each time the pair (rank[i], rank[i + k]) differs from the previous one (use -1 when i + k is off the end).
3. Stop when the last suffix has rank n - 1, meaning all ranks are different.

The code

import random


def suffix_array(s: str, trace: bool = False, stats: dict | None = None) -> list[int]:
    """Prefix doubling with two counting-sort passes per round: O(n log n) time, O(n) space."""
    n = len(s)
    if n == 0:
        return []
    sa = sorted(range(n), key=lambda i: s[i])           # round 0: order by first character
    rank = [0] * n
    for idx in range(1, n):
        rank[sa[idx]] = rank[sa[idx - 1]] + (s[sa[idx]] != s[sa[idx - 1]])
    k = 1
    while rank[sa[-1]] < n - 1:                         # until every suffix has its own rank
        # Sort by the second key (rank of the suffix starting k later) first...
        by_second = list(range(n - k, n)) + [i - k for i in sa if i >= k]
        # ...then a stable counting sort on the first key (rank[i]) finishes the pair sort.
        count = [0] * n
        for i in by_second:
            count[rank[i]] += 1
        total = 0
        for r in range(n):
            count[r], total = total, total + count[r]
        for i in by_second:
            sa[count[rank[i]]] = i
            count[rank[i]] += 1

        def pair(i):                                     # (rank of first k, rank of next k or -1)
            return rank[i], (rank[i + k] if i + k < n else -1)

        new_rank = [0] * n
        for idx in range(1, n):
            new_rank[sa[idx]] = new_rank[sa[idx - 1]] + (pair(sa[idx]) != pair(sa[idx - 1]))
        rank = new_rank
        if trace:
            print(f"after sorting by first {2 * k} chars: sa={sa} rank={rank}")
        k *= 2
        if stats is not None:
            stats['rounds'] = stats.get('rounds', 0) + 1
    return sa


def naive_suffix_array(s: str) -> list[int]:
    return sorted(range(len(s)), key=lambda i: s[i:])


if __name__ == "__main__":
    word = "banana"
    sa = suffix_array(word, trace=True)
    print("suffix array of", word, "=", sa)
    for start in sa:
        print(f"  {start}: {word[start:]}")

    random.seed(4)
    for _ in range(2000):
        text = "".join(random.choice("ab") for _ in range(random.randint(0, 40)))
        assert suffix_array(text) == naive_suffix_array(text), text
    for text in ("a" * 300, "ab" * 150, "abc" * 100, "z", ""):
        assert suffix_array(text) == naive_suffix_array(text)
    print("random and repetitive checks passed")

    for name, text in [("a * 1024", "a" * 1024), ("a * 100000", "a" * 100000),
                       ("random, 100000 letters", "".join(random.choice("abcd") for _ in range(100000)))]:
        stats = {}
        suffix_array(text, stats=stats)
        print(f"{name}: {stats['rounds']} doubling rounds")

Output:

after sorting by first 2 chars: sa=[5, 1, 3, 0, 2, 4] rank=[2, 1, 3, 1, 3, 0]
after sorting by first 4 chars: sa=[5, 3, 1, 0, 4, 2] rank=[3, 2, 5, 1, 4, 0]
suffix array of banana = [5, 3, 1, 0, 4, 2]
  5: a
  3: ana
  1: anana
  0: banana
  4: na
  2: nana
random and repetitive checks passed
a * 1024: 10 doubling rounds
a * 100000: 17 doubling rounds
random, 100000 letters: 5 doubling rounds

Worked trace on banana

Round 0 ranks by first letter: a=0, b=1, n=2, so rank = [1, 0, 2, 0, 2, 0].

• k = 1 (sort by 2 characters). The pairs (rank[i], rank[i+1]) are: i=0 (1,0), i=1 (0,2), i=2 (2,0), i=3 (0,2), i=4 (2,0), i=5 (0,-1). Sorted: 5, then 1 and 3 (tied), then 0, then 2 and 4 (tied). The output line confirms sa=[5, 1, 3, 0, 2, 4] with rank=[2, 1, 3, 1, 3, 0]. Two ties remain, so another round.
• k = 2 (sort by 4 characters). The pairs now use the ranks from the k = 1 round, rank = [2, 1, 3, 1, 3, 0]: i=0 (2,3), i=1 (1,1), i=2 (3,3), i=3 (1,0), i=4 (3,-1), i=5 (0,-1). Tied suffixes 1 (anana) and 3 (ana) are separated by the rank two places ahead: suffix 3 reads rank[5] = 0 (the suffix a) and suffix 1 reads rank[3] = 1 (the suffix ana), so 3 sorts first. Suffix 4 (na) has only two characters left, so its second key is -1, and it sorts before suffix 2 (nana), whose second key is rank[4] = 3. The result is sa=[5, 3, 1, 0, 4, 2] with all ranks different (rank=[3, 2, 5, 1, 4, 0]), so the loop ends after 2 rounds. The printed table lists the suffixes in order.

The sort step, executed on banana

The clever line by_second = list(range(n - k, n)) + [i - k for i in sa if i >= k] builds the second-key order without sorting: the suffixes that have no second half come first, and every other suffix i is listed in the order of the suffix i + k that starts its second half, which is exactly the existing order sa shifted back by k. The next block runs those steps on banana, printing the intermediate arrays; its first round follows from the k = 1 pairs above and its second round from the k = 2 pairs.

s = "banana"
n = len(s)
sa = sorted(range(n), key=lambda i: s[i])
rank = [0] * n
for idx in range(1, n):
    rank[sa[idx]] = rank[sa[idx - 1]] + (s[sa[idx]] != s[sa[idx - 1]])
print("round 0: rank =", rank, "(a=0, b=1, n=2)")

k = 1
while rank[sa[-1]] < n - 1:
    print(f"\n--- k = {k}: sort by the first {2 * k} characters ---")
    tail = list(range(n - k, n))                        # suffixes with no second half
    shifted = [i - k for i in sa if i >= k]             # everything else, in the order of their second halves
    by_second = tail + shifted
    print("suffixes with no second half:", tail, "; the rest, from sa =", sa, "shifted back by", k, ":", shifted)
    print("by_second =", by_second, "(ordered by rank[i + k])")

    count = [0] * n
    for i in by_second:
        count[rank[i]] += 1
    print("count of each first-key rank:", count)
    total = 0
    for r in range(n):
        count[r], total = total, total + count[r]
    print("start slot of each rank (prefix sums):", count)
    for i in by_second:
        sa[count[rank[i]]] = i
        count[rank[i]] += 1
    print("sa after the stable counting sort =", sa)

    pair = lambda i: (rank[i], rank[i + k] if i + k < n else -1)
    print("pairs in sa order:", [(i, pair(i)) for i in sa])
    new_rank = [0] * n
    for idx in range(1, n):
        new_rank[sa[idx]] = new_rank[sa[idx - 1]] + (pair(sa[idx]) != pair(sa[idx - 1]))
    rank = new_rank
    print("new rank =", rank)
    k *= 2
print("\nsuffix array =", sa)

Output:

round 0: rank = [1, 0, 2, 0, 2, 0] (a=0, b=1, n=2)

--- k = 1: sort by the first 2 characters ---
suffixes with no second half: [5] ; the rest, from sa = [1, 3, 5, 0, 2, 4] shifted back by 1 : [0, 2, 4, 1, 3]
by_second = [5, 0, 2, 4, 1, 3] (ordered by rank[i + k])
count of each first-key rank: [3, 1, 2, 0, 0, 0]
start slot of each rank (prefix sums): [0, 3, 4, 6, 6, 6]
sa after the stable counting sort = [5, 1, 3, 0, 2, 4]
pairs in sa order: [(5, (0, -1)), (1, (0, 2)), (3, (0, 2)), (0, (1, 0)), (2, (2, 0)), (4, (2, 0))]
new rank = [2, 1, 3, 1, 3, 0]

--- k = 2: sort by the first 4 characters ---
suffixes with no second half: [4, 5] ; the rest, from sa = [5, 1, 3, 0, 2, 4] shifted back by 2 : [3, 1, 0, 2]
by_second = [4, 5, 3, 1, 0, 2] (ordered by rank[i + k])
count of each first-key rank: [1, 2, 1, 2, 0, 0]
start slot of each rank (prefix sums): [0, 1, 3, 4, 6, 6]
sa after the stable counting sort = [5, 3, 1, 0, 4, 2]
pairs in sa order: [(5, (0, -1)), (3, (1, 0)), (1, (1, 1)), (0, (2, 3)), (4, (3, -1)), (2, (3, 3))]
new rank = [3, 2, 5, 1, 4, 0]

suffix array = [5, 3, 1, 0, 4, 2]

Reading the output for k = 1: sa before the round is [1, 3, 5, 0, 2, 4] (sorted by first letter). Position 5 has no second half, so it goes first. The other positions are sa minus 1, giving [0, 2, 4, 1, 3]: suffix 0 comes first because the suffix after it (anana, starting at 1) leads the old order. count is [3, 1, 2, ...] because three suffixes start with a (rank 0), one with b and two with n; the running sum turns it into the starting slots [0, 3, 4, ...]; placing the items of by_second in order drops 5, 1, 3 into slots 0 to 2, then 0 into slot 3, then 2, 4 into slots 4 and 5. Because the placement walks by_second in order and is stable, 1 stays before 3 inside the tie, exactly as their second keys demand.

The last lines of the main listing's output show that the number of rounds follows the logarithm of the length: aaaa... of length 1024 takes 10 rounds, of length 100,000 takes 17, while a random 100,000-letter string over four letters finishes in 5 because ties disappear quickly. The listing asserts equality with a brute-force sort on thousands of random strings and on repetitive ones.

Complexity of the shipped code

• Time: the initial sorted on one character per position is O(n log n). Then at most ceil(log2 n) rounds, each with two linear passes plus the counting sort, so O(n) per round. Total O(n log n).
• Space: O(n): the arrays sa, rank, new_rank, by_second and count, each of n integers.
• If you replace the counting sort by sorted(..., key=pair) inside each round, the cost of each round is O(n log n) and the total becomes O(n log^2 n). That is the version most people write in interviews, and it is simpler to write.

Compared with the other methods

SA-IS in one pass of ideas

SA-IS is almost never implemented in an interview; the point of knowing it is that linear time is possible and roughly how. The labelling it starts from, computed for banana$:

s = "banana$"                      # $ is the sentinel, smaller than every letter
n = len(s)
kind = [""] * n
kind[n - 1] = "S"                  # by convention the sentinel suffix is S-type
for i in range(n - 2, -1, -1):
    if s[i] == s[i + 1]:
        kind[i] = kind[i + 1]      # equal letters: same type as the next suffix
    else:
        kind[i] = "S" if s[i] < s[i + 1] else "L"
lms = [i for i in range(1, n) if kind[i] == "S" and kind[i - 1] == "L"]
print("position:", " ".join(f"{i:2d}" for i in range(n)))
print("letter  :", " ".join(f"{c:>2s}" for c in s))
print("type    :", " ".join(f"{c:>2s}" for c in kind))
print("LMS positions (S-type with an L-type left neighbour):", lms)

Output:

position:  0  1  2  3  4  5  6
letter  :  b  a  n  a  n  a  $
type    :  L  S  L  S  L  L  S
LMS positions (S-type with an L-type left neighbour): [1, 3, 6]

The sentinel $ is the smallest character. Reading right to left: a at 5 is larger than $, so L; n at 4 is larger than the a after it, so L; a at 3 is smaller than n, so S; and so on. The LMS positions here are 1, 3 and 6.

SA-IS builds the suffix array in linear time by induced sorting. Add a sentinel smaller than every character at the end and label each suffix S-type if it is smaller than the suffix after it, and L-type if larger (the labelling listing above does this for banana$). An LMS position is an S-type position whose left neighbour is L-type (leftmost S of a run). The steps: (1) classify every suffix; (2) place the LMS suffixes in the ends of their buckets (a bucket holds all suffixes starting with one character), then scan left to right to induce the L-type suffixes from them and right to left to induce the S-type ones; this sorts the substrings between consecutive LMS positions (the LMS substrings); (3) give each LMS substring a name by its rank, and if names repeat, recurse on the string of names, which has at most n/2 symbols because LMS positions are never adjacent; (4) place the LMS suffixes in their true order and induce once more to fill in everything else. Each level does linear work on a string at most half as long, so the total is O(n). The price is complexity; use a library when speed matters.

Edge cases and common mistakes

• Empty string returns an empty list; one character returns [0] without entering the loop.
• Suffixes that run off the end must sort first within a tie. by_second lists them first and the stable counting sort keeps that order, so on banana at k = 2 na lands before nana; the trace above shows it. Listing them last in by_second is the wrong order: it returns [1, 3, 5, 0, 2, 4] for banana.
• The second key used when recomputing ranks is -1 for a suffix with no second half, so that it never equals a real rank (real ranks start at 0). The order itself comes from by_second, not from this value. Using 0 there makes such a suffix look equal to one whose second half has rank 0, and the random check in the listing then fails.
• Using the same array for old and new ranks while updating them; the code writes to a new rank list.
• Letting the counting-sort array be sized by the alphabet instead of by the number of distinct ranks (up to n).
• Assuming the loop runs a fixed log n rounds: it stops as soon as ranks are distinct, which is after 5 rounds for the random string above.
• Using s[i:] slices for comparisons inside the loop: it silently brings back the O(n^2) memory.

Running the code

# Python 3.12, standard library only
# save the main listing as suffix_array.py and the two small blocks as sort_step.py and sais_labels.py
python suffix_array.py
python sort_step.py
python sais_labels.py