String Algorithms

The simplest approach to finding pattern P in text T is to slide a window of length m across T and check character-by-character at each position. This is O(nm) in the worst case.

The worst case is not contrived. Consider T = aaaaaaaaab and P = aaaab. At every position, we match m-1 characters before a mismatch. On a genomic sequence with long repetitive regions, this is catastrophic. A human genome has ~3×10⁹ characters; reading it with a naive matcher and a 100-character pattern could require ~3×10¹¹ comparisons — minutes of runtime for a task we want done in seconds.

Where is the wasted work? When we shift the window by 1 and start comparing from scratch, we throw away all the comparison information we just gathered. If we already know that T[1..5] = abcab, and we're now aligning P at position 2, we already know what T[2..5] looks like — we compared it a moment ago. Every fast string algorithm exploits this observation in a different way.

The naive algorithm is still perfectly correct and often fast enough in practice (most modern text editors use Boyer-Moore variants, and for short patterns on ASCII text, cache effects make it competitive). Understanding its failure modes motivates the entire field.

def naive_search(text: str, pattern: str) -> list[int]:
    n, m = len(text), len(pattern)
    matches = []
    for i in range(n - m + 1):        # O(n) shifts
        if text[i:i+m] == pattern:     # O(m) comparison
            matches.append(i)
    return matches

# Worst case: text = 'a'*1000 + 'b', pattern = 'a'*500 + 'b'
# => ~500,000 character comparisons

KMP (Knuth, Morris, Pratt, 1977) achieves O(n + m) worst-case time by never moving the text pointer backward. The insight: after a partial match followed by a mismatch, the pattern itself tells us the minimum safe shift.

The failure function (also called the prefix function or pi array). For pattern P of length m, pi[i] is the length of the longest proper prefix of P[0..i] that is also a suffix. Example:

P =  a  b  a  b  c  a
i =  0  1  2  3  4  5
pi = 0  0  1  2  0  1

At i=3 (P[0..3] = abab), the longest proper prefix-suffix is ab (length 2). This means: if we're at position j in the pattern and hit a mismatch, we can resume comparison from pi[j-1] in the pattern (not from 0), because we know the first pi[j-1] characters of P already match the text just before our current position.

Intuition via analogy. Imagine you're reading a book aloud to a stenographer, and they ask you to repeat a phrase whenever you make a mistake. KMP tells you the minimum number of words to re-read — you don't need to go back to the beginning of the phrase, just to the last valid restart point.

Correctness argument. The invariant of the KMP matcher is: at each step, the text character T[i] is being compared against pattern character P[j]. If they match, advance both. If they mismatch and j > 0, set j = pi[j-1] and retry (do NOT advance i). If j = 0 and mismatch, advance i. The text pointer i never goes backward — it advances exactly n times. Each step either advances i or decreases j (which can only happen m times total since j is bounded by m). So the total work is O(n + m).

Building pi in O(m) time. The same logic applies recursively. When building pi[i], we're asking: is P[pi[i-1]] == P[i]? If yes, pi[i] = pi[i-1]+1. If no, fall back through the chain until we find a match or reach 0.

def build_pi(pattern: str) -> list[int]:
    m = len(pattern)
    pi = [0] * m
    k = 0  # length of previous longest prefix-suffix
    for i in range(1, m):
        while k > 0 and pattern[k] != pattern[i]:
            k = pi[k - 1]  # fall back
        if pattern[k] == pattern[i]:
            k += 1
        pi[i] = k
    return pi

def kmp_search(text: str, pattern: str) -> list[int]:
    n, m = len(text), len(pattern)
    pi = build_pi(pattern)
    matches = []
    j = 0  # number of chars matched so far
    for i in range(n):
        while j > 0 and text[i] != pattern[j]:
            j = pi[j - 1]  # fall back without moving i
        if text[i] == pattern[j]:
            j += 1
        if j == m:  # full match
            matches.append(i - m + 1)
            j = pi[j - 1]  # prepare for next match
    return matches

# Example:
# kmp_search('abababcab', 'abab') => [0, 2]

Rabin-Karp (Rabin and Karp, 1987) takes a completely different approach: instead of comparing characters, compare hash values. Compute a hash of P and of each length-m window of T. A mismatch in hashes proves no match; a hash match triggers a full verification.

The rolling hash trick. Naively, computing a new hash from scratch for each window takes O(m) time, which is no better than the naive algorithm. The breakthrough is a polynomial rolling hash (Rabin fingerprint) that can be updated in O(1):

Let the hash of a string s[0..m-1] be:

H = s[0]·b^(m-1) + s[1]·b^(m-2) + ... + s[m-1]·b^0   (mod q)

where b is a base (e.g. 256 for ASCII, 31 for lowercase letters) and q is a large prime.

When sliding the window one step right:

H_new = (H_old - s[0]·b^(m-1)) · b + s[m]   (mod q)

Removing the leftmost character costs O(1) (precompute b^(m-1) mod q once).

Spurious hits and correctness. If H(window) ≠ H(P), there is definitely no match. If H(window) = H(P), there might be a match — hash collisions are possible. Rabin-Karp always performs a full O(m) verification on a hash match. In the worst case (adversarial text with many hash collisions) this degrades to O(nm), but for random strings and a random prime q, the expected collision probability per window is 1/q. With q ≈ 10⁹, the expected extra work is negligible.

Where Rabin-Karp excels: multi-pattern search. If you're searching for k different patterns simultaneously (all of the same length), you can precompute a hash set of all k patterns and check each window against the set in O(1) expected time. This gives O(n + km) total, far better than running k KMP searches. This is the basis of plagiarism detection tools like MOSS.

Connection to randomized algorithms. Rabin-Karp is a beautiful example of Las Vegas randomization: it is always correct (it verifies on hash match), and its expected runtime is O(n + m). Choosing a random prime q from a large space makes adversarial inputs infeasible to construct.

def rabin_karp(text: str, pattern: str, base: int = 256, mod: int = 10**9+7) -> list[int]:
    n, m = len(text), len(pattern)
    if m > n:
        return []
    
    h_pat = 0  # hash of pattern
    h_win = 0  # hash of current window
    high_pow = pow(base, m - 1, mod)  # base^(m-1) mod q
    
    for i in range(m):
        h_pat = (h_pat * base + ord(pattern[i])) % mod
        h_win = (h_win * base + ord(text[i])) % mod
    
    matches = []
    for i in range(n - m + 1):
        if h_win == h_pat:
            if text[i:i+m] == pattern:  # verify to rule out spurious hits
                matches.append(i)
        if i < n - m:
            # Roll the hash: remove leftmost, add new rightmost
            h_win = (h_win - ord(text[i]) * high_pow) % mod
            h_win = (h_win * base + ord(text[i + m])) % mod
            h_win %= mod  # ensure non-negative
    return matches

Boyer-Moore (1977) is the algorithm that powers grep and most production string search tools. Unlike KMP and Rabin-Karp, it scans the pattern right-to-left and can skip large sections of text, achieving sublinear average-case behavior — it examines far fewer than n characters on typical inputs.

Right-to-left scanning. Align P against T[0..m-1]. Compare P[m-1] against T[m-1], then P[m-2] against T[m-2], etc. This seems counterintuitive, but it enables dramatic skipping.

Bad-character rule. Suppose we find a mismatch at position j in P (comparing with T[i]). The mismatched text character is T[i] = c. If c does not appear in P at all, we can shift P completely past position i — a shift of up to m characters! If c appears in P to the left of j, we align that rightmost occurrence of c in P with T[i] and continue. The bad-character table bc[c] = rightmost occurrence of character c in P.

Good-suffix rule. Suppose we've matched the suffix P[j+1..m-1] of P but hit a mismatch at j. We know that T[i+j+1..i+m-1] matches P[j+1..m-1]. The good-suffix rule asks: is there another occurrence of this matched suffix elsewhere in P? If so, shift P to align that occurrence with what we matched. If not, shift P entirely past the matched suffix.

Combining the rules. At each mismatch position, compute both shifts and take the maximum. In practice, Boyer-Moore skips an average of O(n/m) character comparisons on random text — examining only 1 in m characters! On large alphabets (like ASCII text), this is extremely fast.

Galil rule and worst-case O(n). The basic Boyer-Moore algorithm has O(nm) worst case on small alphabets (like binary strings). The Galil rule extends it to O(n) worst case by remembering which portion of the text has already been matched.

Why does industry use this? For interactive search (text editors, grep on files), average-case performance dominates. Boyer-Moore with the Galil rule is optimal in all cases and blazingly fast in practice. The GNU grep implementation uses a variant combining Boyer-Moore with Horspool's simplification.

def build_bad_char(pattern: str) -> dict:
    """rightmost occurrence of each character in pattern"""
    bc = {}
    for i, c in enumerate(pattern):
        bc[c] = i
    return bc

def boyer_moore_simple(text: str, pattern: str) -> list[int]:
    """Boyer-Moore with bad-character rule only (Horspool variant).
    O(nm) worst-case, O(n/m) average on large alphabets."""
    n, m = len(text), len(pattern)
    bc = build_bad_char(pattern)
    matches = []
    i = 0  # alignment position in text
    while i <= n - m:
        j = m - 1  # scan pattern right-to-left
        while j >= 0 and pattern[j] == text[i + j]:
            j -= 1
        if j < 0:  # full match
            matches.append(i)
            i += 1  # or use good-suffix rule
        else:
            # shift so rightmost occurrence of text[i+j] in P aligns with i+j
            bc_shift = j - bc.get(text[i + j], -1)
            i += max(1, bc_shift)
    return matches

KMP can be understood as implicitly building a Deterministic Finite Automaton (DFA) that reads text and transitions to accept states when the pattern is found. Making this explicit reveals deep connections to formal language theory.

The string-matching automaton. For a pattern P of length m over alphabet Σ, build a DFA with m+1 states {0, 1, ..., m}. State q represents 'we have matched the first q characters of P'. State m is accepting (pattern found). The transition function δ(q, c) answers: 'if we've matched q characters and see character c, how many characters have we now matched?'

Formal definition: δ(q, c) = the length of the longest prefix of P that is a suffix of P[0..q-1] + c.

Building the automaton. For each state q ∈ {0..m} and character c ∈ Σ:

• If c == P[q] (the next expected character), δ(q, c) = q + 1
• Otherwise, δ(q, c) = δ(pi(q), c) recursively — use the failure function!

Building the full table takes O(m|Σ|) time and space. Matching then runs in O(n) time with zero backtracking, just table lookups.

Connection to KMP. KMP is an implicit DFA — the failure function encodes the fall-back transitions without materializing the whole table. When |Σ| is large (e.g., Unicode), the DFA table is impractical; KMP's implicit representation saves space.

Multiple patterns: Aho-Corasick. When matching k patterns simultaneously, Aho-Corasick (1975) builds a single DFA-like automaton (a trie with failure links) for all patterns. Text is scanned once, with all k patterns matched simultaneously, in O(n + Σ|P_i| + output_size) time. This is the algorithm inside many virus scanners, network intrusion detection systems, and command-line tools like fgrep.

def build_automaton(pattern: str, alphabet: str) -> list[dict]:
    """δ[state][char] = next_state. O(m * |alphabet|) time."""
    m = len(pattern)
    pi = build_pi(pattern)  # reuse KMP failure function
    delta = [{} for _ in range(m + 1)]
    
    for q in range(m + 1):
        for c in alphabet:
            if q < m and c == pattern[q]:
                delta[q][c] = q + 1
            elif q == 0:
                delta[q][c] = 0
            else:
                # fall back: delta[q][c] = delta[pi[q-1]][c]
                delta[q][c] = delta[pi[q - 1]].get(c, 0)
    return delta

def automaton_search(text: str, pattern: str) -> list[int]:
    alpha = set(text) | set(pattern)
    delta = build_automaton(pattern, alpha)
    m = len(pattern)
    q = 0
    matches = []
    for i, c in enumerate(text):
        q = delta[q].get(c, 0)
        if q == m:
            matches.append(i - m + 1)
    return matches

All the algorithms above require O(m) time per query for a fixed text T. What if T is fixed (a database, a genome reference) and we want to answer many queries quickly? We want a text index — preprocess T once, then answer each query in O(m log n) or better.

The suffix array (Manber and Myers, 1990). The suffix array SA is a permutation of {0, 1, ..., n-1} such that SA[i] = j means the j-th suffix of T (T[j..n-1]) is the i-th in lexicographic order. For T = banana$ ($ is a lexicographically smallest sentinel):

Suffixes:          Sorted:
0: banana$         0: $          → SA[0]=6
1: anana$          1: a$         → SA[1]=5
2: nana$           2: ana$       → SA[2]=3
3: ana$            3: anana$     → SA[3]=1
4: na$             4: banana$    → SA[4]=0
5: a$              5: na$        → SA[5]=4
6: $               6: nana$      → SA[6]=2

Pattern matching with a suffix array. All occurrences of P in T correspond to a contiguous range in SA (since all suffixes starting with P appear consecutively in sorted order). Binary search for the leftmost and rightmost positions where SA[i] is a suffix starting with P. This takes O(m log n) comparisons (each comparison is O(m)), giving O(m log n) total.

LCP arrays accelerate this. The LCP (Longest Common Prefix) array LCP[i] stores the length of the longest common prefix of SA[i-1] and SA[i]. With LCP, the binary search can be improved to O(m + log n). And many other queries (longest repeated substring, number of distinct substrings, etc.) become O(n) or O(m) with the LCP array.

O(n log n) construction. The naive approach (sort all suffixes as strings) takes O(n² log n). The DC3/Skew algorithm achieves O(n). A practical O(n log n) approach uses prefix-doubling: sort by 2^k-length prefixes iteratively, using radix sort at each level. The SA-IS algorithm achieves O(n) in practice and is the basis of most production tools.

Applications. Suffix arrays power bwa mem and bowtie2 (genomic alignment), full-text search engines (combined with Burrows-Wheeler Transform), and bioinformatics tools that search for k-mers across large genomes.

def build_suffix_array_naive(text: str) -> list[int]:
    """O(n^2 log n) — correct but slow. Production uses SA-IS O(n)."""
    n = len(text)
    suffixes = sorted(range(n), key=lambda i: text[i:])
    return suffixes

def search_suffix_array(text: str, sa: list[int], pattern: str) -> tuple[int, int]:
    """Binary search: returns [lo, hi) range in SA of matching suffixes.
    O(m log n) time."""
    n, m = len(text), len(pattern)
    
    # Find leftmost position
    lo, hi = 0, n
    while lo < hi:
        mid = (lo + hi) // 2
        if text[sa[mid]:sa[mid]+m] < pattern:
            lo = mid + 1
        else:
            hi = mid
    left = lo
    
    # Find rightmost position
    hi = n
    while lo < hi:
        mid = (lo + hi) // 2
        if text[sa[mid]:sa[mid]+m] > pattern:
            hi = mid
        else:
            lo = mid + 1
    right = lo
    
    return left, right  # SA[left:right] contains all match positions

# Usage:
# T = 'banana$'
# sa = build_suffix_array_naive(T)
# l, r = search_suffix_array(T, sa, 'ana')
# match positions: [sa[l], sa[l+1], ..., sa[r-1]] = [1, 3]

The Burrows-Wheeler Transform (BWT) is a reversible string transformation that rearranges T into a form that is more compressible AND supports exact-match searching — the foundation of modern genomic aligners.

Computing BWT. Conceptually: form all n cyclic rotations of T (append a sentinel $ first), sort them lexicographically (this is exactly the suffix array), and take the last column. The result BWT[i] = T[(SA[i] - 1 + n) mod n].

For T = banana$:

Sorted rotations:   Last column (BWT):
$banana             a
a$banan             n
ana$ban             n
anana$b             b
banana$             $  (last char of original)
na$bana             a
nana$ba             a

BWT = annb$aa — notice the runs of identical characters (a's, n's). This makes it highly compressible with run-length encoding.

The FM-Index (Ferragina-Manzini, 2000). By storing BWT plus a small set of auxiliary structures (the F column from sorted suffixes, occurrence counts, and checkpoint arrays), you get a data structure that:

1. Stores T in near-compressed size (close to the entropy of T)
2. Supports exact pattern matching in O(m) time via backward search

Backward search intuition. To find P = ana in T: start with the range of all BWT rows (the full sorted suffix array). Narrow it character by character, processing P right-to-left. At each step, the range shrinks to only those suffixes whose corresponding text positions start with the characters we've processed so far. This is an LF-mapping (Last-to-First) step, computable in O(1) with the auxiliary structures.

Real-world impact. bwa mem, bowtie2, SOAP2, and samtools all use FM-Index internally. The human reference genome (hg38, ~3GB) is stored as a ~1.5GB FM-Index that supports querying a 150-character read in microseconds — one of the most elegant space-time tradeoffs in all of computer science.

def bwt(text: str) -> str:
    """Compute BWT from text (must include $ sentinel at end)."""
    assert text[-1] == '$'
    n = len(text)
    sa = sorted(range(n), key=lambda i: text[i:] + text[:i])
    return ''.join(text[(i - 1) % n] for i in sa)

def ibwt(bwt_str: str) -> str:
    """Invert BWT to recover original text."""
    n = len(bwt_str)
    # Build F (first column) and LF-mapping
    f = sorted(bwt_str)
    # LF mapping: LF[i] = position of bwt_str[i] in F
    from collections import defaultdict
    count = defaultdict(int)
    lf = [0] * n
    f_rank = {}  # cumulative counts in F
    for i, c in enumerate(sorted(set(bwt_str))):
        f_rank[c] = sum(1 for x in f[:f.index(c)])
    char_count = defaultdict(int)
    for i, c in enumerate(bwt_str):
        lf[i] = f_rank[c] + char_count[c]
        char_count[c] += 1
    # Walk back from $ to reconstruct
    result = []
    i = bwt_str.index('$')
    for _ in range(n - 1):
        i = lf[i]
        result.append(bwt_str[i])
    return ''.join(reversed(result))

# bwt('banana$') => 'annb$aa'
# ibwt('annb$aa') => 'banana'

The Z-array is a simpler alternative to KMP that is often easier to implement correctly and equally powerful. For string S of length n, Z[i] is the length of the longest substring starting at S[i] that is also a prefix of S.

Z[0] is conventionally defined as the full length n (or 0, depending on convention). For i > 0, Z[i] tells us how far S[i..] matches S[0..].

Example. For S = aabxaabaabx:

i:    0  1  2  3  4  5  6  7  8  9  10
S:    a  a  b  x  a  a  b  a  a  b  x
Z:   11  1  0  0  3  1  0  2  3  1  0

At i=4: S[4..6] = aab = S[0..2], so Z[4]=3.

O(n) construction via the Z-box. Maintain a window [L, R] — the rightmost segment [L, R] that matches a prefix of S. For each new i:

• If i > R: extend naively, update [L, R]
• If i ≤ R: Z[i] ≥ min(Z[i-L], R-i+1). Use this to avoid redundant comparisons.

Pattern matching with Z. Concatenate P + '$' + T (where $ is a character not in the alphabet). Compute the Z-array of this concatenated string. Any position i in the T-part with Z[i] ≥ m is a match. This is elegant and linear, often preferred in competitive programming for its simplicity.

Comparison with KMP. Both are O(n+m). KMP processes only T with pattern info encoded in pi. Z-algorithm processes the concatenation. KMP is more memory-efficient for online streaming (never needs to store T). Z-algorithm is simpler to implement correctly and reason about.

def z_function(s: str) -> list[int]:
    n = len(s)
    z = [0] * n
    z[0] = n
    l, r = 0, 0
    for i in range(1, n):
        if i < r:
            z[i] = min(r - i, z[i - l])
        while i + z[i] < n and s[z[i]] == s[i + z[i]]:
            z[i] += 1
        if i + z[i] > r:
            l, r = i, i + z[i]
    return z

def z_search(text: str, pattern: str) -> list[int]:
    combined = pattern + '$' + text
    z = z_function(combined)
    m = len(pattern)
    matches = []
    for i in range(m + 1, len(combined)):
        if z[i] >= m:
            matches.append(i - m - 1)  # position in original text
    return matches

# z_search('aababcaababc', 'aab') => [0, 5]

When you need to find occurrences of k patterns P₁, P₂, ..., Pₖ in text T simultaneously, running KMP k times takes O(k·n + Σ|Pᵢ|). Aho-Corasick achieves O(n + Σ|Pᵢ| + output_size) by building a single automaton for all patterns — text is scanned exactly once.

Construction (three steps).

1. Build a trie of all patterns. Each node represents a prefix of some pattern. Leaf nodes are tagged with the pattern they represent.
2. Add failure links (like KMP's failure function, but for the trie). The failure link of node v points to the node representing the longest proper suffix of v's string that is also a prefix in the trie.
3. Add output links (dictionary links): if a node v represents a prefix that is itself a suffix of another pattern in the trie, add a link to enable reporting all matching patterns.

Matching. Traverse the trie following characters of T. On a mismatch, follow failure links (never backtrack in T). At each node, report all patterns in the output closure. Total time: O(n + output_size) — linear in text length regardless of the number of patterns.

Real applications. Unix fgrep -f patterns.txt uses Aho-Corasick. Network intrusion detection systems (Snort, Suricata) match thousands of attack signatures against incoming packets using Aho-Corasick. Antivirus engines match millions of virus signatures using this algorithm. DNA motif scanners look for hundreds of transcription factor binding sites simultaneously.

Space. The trie has Σ|Pᵢ| nodes. The failure and output links add O(Σ|Pᵢ|) pointers. The DFA version (fully materializing all transitions) uses O(Σ|Pᵢ| · |Σ|) space, which can be large for Unicode.

from collections import deque, defaultdict

class AhoCorasick:
    def __init__(self):
        self.goto = [{}]  # goto[state][char] = next_state
        self.fail = [0]
        self.output = [set()]
    
    def add_pattern(self, pattern: str, pid: int):
        state = 0
        for c in pattern:
            if c not in self.goto[state]:
                self.goto.append({})
                self.fail.append(0)
                self.output.append(set())
                self.goto[state][c] = len(self.goto) - 1
            state = self.goto[state][c]
        self.output[state].add(pid)
    
    def build(self):
        q = deque()
        for c, s in self.goto[0].items():
            q.append(s)  # depth-1 nodes have fail -> 0
        while q:
            r = q.popleft()
            for c, s in self.goto[r].items():
                q.append(s)
                # Follow fail links to find longest matching suffix
                state = self.fail[r]
                while state != 0 and c not in self.goto[state]:
                    state = self.fail[state]
                self.fail[s] = self.goto[state].get(c, 0)
                if self.fail[s] == s:
                    self.fail[s] = 0
                self.output[s] |= self.output[self.fail[s]]
    
    def search(self, text: str) -> list[tuple]:
        state = 0
        results = []
        for i, c in enumerate(text):
            while state != 0 and c not in self.goto[state]:
                state = self.fail[state]
            state = self.goto[state].get(c, 0)
            for pid in self.output[state]:
                results.append((i, pid))
        return results

Bioinformatics is the domain where string algorithms have the most visible and consequential real-world impact. The central task: read mapping (also called alignment) — given millions of short DNA reads (150 characters each) produced by a sequencer, find where each read originated in the ~3 billion character human reference genome.

Scale. A single human whole-genome sequencing run produces ~300 million reads × 150 characters = 45 billion character comparisons if done naively. This must complete in hours on a laptop. No O(nm) algorithm survives.

Exact matching: FM-Index. Tools like bwa mem and bowtie2 build an FM-Index of the reference genome (~1.5 GB compressed). Each read is matched against the index in O(m) time using backward search. With ~300M reads of length 150, total time is O(300M × 150) = 45 billion elementary operations — feasible in hours on modern hardware.

Approximate matching and seeds. Real reads contain sequencing errors. The standard approach: seed-and-extend. Extract exact-match seeds (short substrings guaranteed to match if the read aligns within error tolerance), find seed locations with the FM-Index, then extend with banded Smith-Waterman dynamic programming (edit distance) in a small window.

k-mer indexing. For de novo assembly (no reference genome), tools use k-mer hash tables: index all length-k substrings. dBG (de Bruijn Graph) assemblers like SPAdes represent the genome as a path through a graph of (k-1)-mer overlaps, an O(n) construction that replaced O(n²) overlap-layout-consensus assemblers.

BLAST. The most-used bioinformatics tool, BLAST (Basic Local Alignment Search Tool) uses a word-based heuristic: find short exact matches (words), extend them with gapped alignment. The statistical framework (E-value = expected number of spurious matches) uses extreme-value distribution theory — an elegant union of string algorithms and statistics.

The Aho-Corasick connection. Scanning for multiple transcription factor binding site motifs (TFBS) across a genome is exactly multi-pattern matching. Tools like FIMO use algorithms equivalent to Aho-Corasick extended with PSSM (position-specific scoring matrices) for fuzzy matching.

# Simplified seed-and-extend in Python (illustrative, not production)

def build_kmer_index(text: str, k: int) -> dict:
    """Index all k-mers: kmer -> list of positions."""
    index = {}
    for i in range(len(text) - k + 1):
        kmer = text[i:i+k]
        index.setdefault(kmer, []).append(i)
    return index

def edit_distance(s: str, t: str) -> int:
    """Standard DP for edit distance."""
    m, n = len(s), len(t)
    dp = list(range(n + 1))
    for i in range(1, m + 1):
        prev, dp[0] = dp[0], i
        for j in range(1, n + 1):
            prev, dp[j] = dp[j], min(
                dp[j] + 1,           # deletion
                dp[j-1] + 1,         # insertion
                prev + (s[i-1] != t[j-1])  # sub/match
            )
    return dp[n]

def seed_and_extend(reference: str, read: str, k: int = 11, max_errors: int = 2) -> list[int]:
    """Find approximate matches of read in reference."""
    index = build_kmer_index(reference, k)
    candidates = set()
    for i in range(len(read) - k + 1):
        seed = read[i:i+k]
        for ref_pos in index.get(seed, []):
            start = ref_pos - i
            if 0 <= start <= len(reference) - len(read):
                candidates.add(start)
    matches = []
    for pos in candidates:
        window = reference[pos:pos+len(read)]
        if edit_distance(window, read) <= max_errors:
            matches.append(pos)
    return sorted(matches)