后缀数组（Suffix Array）深度实战：从前缀倍增、SA-IS 到 LCP 数组与模式匹配的工程全解

后缀数组（Suffix Array，SA）是字符串处理领域最基础、最高效的索引结构之一。它把"一个字符串的所有后缀按字典序排序后的起始位置"紧凑地存成一个长度 n 的整数数组，却能在 O(m log n) 内完成任意模式串的精确匹配、在 O(n) 内求最长重复子串、不同子串计数、最长公共子串等经典问题。它比后缀树省内存、比后缀自动机易实现，是生物信息学（DNA 比对）、全文检索（FM-index 的底座）、数据去重与 plagiarism 检测的工程首选。

本文从第一性原理出发，完整推导朴素构造、前缀倍增（Manber–Myers）、O(n log n) 排序优化、线性时间 SA-IS，再到 LCP 数组的 Kasai 线性算法，最后落到模式匹配、最长重复/公共子串、不同子串计数等生产级应用，并给出可复现的 Python 工具箱与 12 项生产陷阱清单。

一、第一性原理：什么是后缀数组

给定字符串 S（长度 n，下标 0-based），它的 n 个后缀是 S[i..n-1]（i = 0..n-1）。后缀数组 SA 是一个整数排列，满足：

SA[k] = 第 k 小的后缀（按字典序）的起始下标

换言之，S[SA[0]..] ≤ S[SA[1]..] ≤ … ≤ S[SA[n-1]..]。

秩数组 rank 是 SA 的逆：rank[SA[k]] = k，即"以 i 开头的后缀在所有后缀中排第几"。SA 与 rank 互为逆映射，二者一起把"子串比较"转化为"整数比较"。

为什么强大：任何子串 S[i..j] 都是某个后缀 S[i..] 的前缀。于是"在 S 中查找模式 P"等价于"在 SA 上二分查找字典序落在 [P, P+maxchar) 区间内的后缀"——把字符串扫描问题转成了有序数组上的区间查询。

二、朴素构造与复杂度下界

最直觉的做法：取出所有 n 个后缀字符串，用标准排序（比较时逐字符比）排序。一次后缀比较最坏 O(n)，故总复杂度 O(n² log n)。当 n 很大（如基因组数据数十亿字符）时完全不可用，但它是理解 SA 语义的最佳入口：

def naive_sa(s: str):
    n = len(s)
    return sorted(range(n), key=lambda i: s[i:])   # O(n^2 log n)

注意这里 key=lambda i: s[i:] 每次比较都重新切片并逐字符比对，最坏每个比较 O(n)。后面所有优化的本质，都是避免逐字符比较，转而用"已经算好的秩"做常数/对数级比较。

三、前缀倍增（Prefix Doubling / Manber–Myers）

核心洞察：比较两个后缀 i 和 j 时，不需要逐字符比。若我们已经知道"每个后缀前 2^k 个字符的排序秩"，那么比较 S[i..] 与 S[j..] 的前 2^{k+1} 个字符，等价于先比较它们前 2^k 个字符（看 rank），相等时再比较从 i+2^k、j+2^k 开始的后 2^k 个字符（同样看 rank）。

于是每一轮把所有后缀按"双关键字 (rank[i], rank[i + 2^k])"排序，就能得到按前 2^{k+1} 个字符排序的新 rank。从 k=0（单字符 rank）开始，每轮长度翻倍，log n 轮后所有后缀完全区分：

def prefix_doubling_sa(s: str):
    s = s + chr(0)          # 哨兵，保证所有后缀互不相等
    n = len(s)
    k = 1
    rank = [ord(s[i]) for i in range(n)]
    tmp = [0] * n
    while True:
        # 双关键字计数排序：(rank[i], rank[(i+k)%n])
        pairs = [(rank[i], rank[(i + k) % n], i) for i in range(n)]
        # 先按第二关键字，再按第一关键字稳定排序（计数排序更稳）
        pairs.sort(key=lambda x: (x[1], x[0]))
        pairs.sort(key=lambda x: (x[0], x[1]))
        tmp[pairs[0][2]] = 0
        for i in range(1, n):
            prev, cur = pairs[i - 1], pairs[i]
            tmp[cur[2]] = tmp[prev[2]] + (prev[:2] != cur[:2])
        rank, tmp = tmp, rank
        if rank[pairs[-1][2]] == n - 1:   # 已全部区分
            break
        k <<= 1
    return [p[2] for p in pairs]   # 去掉哨兵后的 SA 需另处理

复杂度分析：每轮排序若用基于 rank 值域的基数/计数排序为 O(n)，共 log n 轮，得到 O(n log n)（使用比较排序则为 O(n log² n)）。前缀倍增是工业界最常用、最易调试的实现，直到 n 超过约 10⁷ 才需要考虑更激进的线性算法。

四、O(n log n) 排序优化（用基数排序替代比较排序）

上面的 pairs.sort 用的是 Python 比较排序，单轮 O(n log n)、整体 O(n log² n)。把"按整数秩排序"替换为以秩为键的计数排序（值域 0..n-1），单轮降到 O(n)，整体达到 O(n log n)。工程中这一步是性能分水岭：

def counting_sort_by_key(keys, order, K):
    """对 order（后缀下标列表）按 keys 升序做稳定计数排序，返回新 order。"""
    cnt = [0] * (K + 1)
    for x in keys:
        cnt[x + 1] += 1
    for i in range(K):
        cnt[i + 1] += cnt[i]
    out = [0] * len(order)
    for i in order:
        out[cnt[keys[i]]] = i
        cnt[keys[i]] += 1
    return out

把"按第一关键字 → 按第二关键字稳定"两层计数排序串起来，就是 Manber–Myers 的标准 O(n log n) 实现。经验法则：n < 10⁷ 用排序法足够；n ≥ 10⁷ 才上 SA-IS。

五、SA-IS：线性时间构造（诱导排序概述）

当 n 极大（全文索引、基因组）时，O(n log n) 仍偏慢。SA-IS（Suffix Array – Induced Sorting）在 O(n) 时间内构造 SA，思路是：

1. 把每个后缀的每个位置标记类型 L（大于后继）或 S（小于后继），以及 LMS（最左 S 型）。
2. 先对 LMS 子串排序（递归，规模约减半），再用"诱导排序"把所有 L 型、S 型位置按已排好的 LMS 顺序依次填入，一遍扫描即可得到完整 SA。

SA-IS 的正确性证明较繁，但实现是一个精心组织的诱导过程。下面给出简化骨架（完整工业实现见 sais-lite / divsufsort）：

def sa_is(s: str):
    # 工业级实现需处理 LMS 子串去重、递归构造子 SA、两遍诱导排序。
    # 此处仅展示调用约定与复杂度承诺：
    #   - 时间 O(n)，空间 O(n)
    #   - 典型第三方实现：pysuffix / divsufsort（C 扩展）
    raise NotImplementedError("use divsufsort for production SA-IS")

实践建议：直接用 divsufsort 的 C 绑定（如 Python 的 suffix_array / pydivsufsort），它在十亿级字符上仍可秒级完成。自写 SA-IS 价值在于理解，生产应优先复用成熟库。

六、LCP 数组与 Kasai 算法

仅有 SA 还不足以高效回答"最长公共前缀"类问题。引入 LCP 数组：LCP[k] = LCP(S[SA[k]..], S[SA[k-1]..])，即排第 k 与第 k-1 的后缀的最长公共前缀长度。

Kasai 算法在 O(n) 内求出 LCP，核心是利用 rank 与"h 的单调性"：设当前考察后缀 i（rank 为 r），它与前一个后缀 j=SA[r-1] 的 LCP 至少为 h-1（h 为上一轮 i+1 与 j+1 的 LCP 减 1），再向后逐字符扩展。

def kasai_lcp(s, sa):
    n = len(s)
    rank = [0] * n
    for i, p in enumerate(sa):
        rank[p] = i
    lcp = [0] * n          # lcp[k] = LCP(SA[k], SA[k-1]); lcp[0]=0
    h = 0
    for i in range(n):
        if rank[i] == 0:
            h = 0
            continue
        j = sa[rank[i] - 1]
        while i + h < n and j + h < n and s[i + h] == s[j + h]:
            h += 1
        lcp[rank[i]] = h
        if h > 0:
            h -= 1
    return lcp

Kasai 的妙处：每个字符最多被 while 扩展一次、被 h-=1 收缩一次，故总步数 O(n)。SA + LCP 组合等价于后缀树——几乎所有后缀树的查询都能用 SA+LCP 在相同或更优复杂度内完成，且内存仅为 2~3 个 int 数组。

七、应用实战

7.1 模式出现次数（精确匹配）

在 SA 上，所有以 P 为前缀的后缀必然连续地排在一段区间 [lo, hi) 内（因为字典序有序）。用两次二分查找定位这个区间，出现次数即为 hi - lo：

def count_occ(s, sa, p):
    n = len(s)
    lo = bisect.bisect_left([s[sa[i]:] for i in range(n)], p)   # 仅示意，生产用逐字符比较
    # 生产实现：直接对 SA 做按前缀的二分，避免切片 O(n)
    def leq(i, p):                 # S[SA[i]..] 是否以 p 为前缀或更大
        return s.startswith(p, sa[i])
    lo = bisect_left(range(n), True, key=lambda i: not leq(i, p))   # 简化
    # 正确写法见下方二分
    return hi - lo

实际工程应写一个"按前缀比较 SA[i] 与 P"的二分，避免每次切片 O(n)：时间 O(m log n)，m 为模式长度。这是全文检索、日志关键字统计的基石。

7.2 最长重复子串（Longest Repeated Substring）

重复子串一定是某个 LCP[k] 对应的前缀，答案就是 max(LCP)，且子串起点为 SA[argmax]：

def longest_repeated_substring(s, sa, lcp):
    k = max(range(len(lcp)), key=lambda i: lcp[i])
    i = sa[k]
    return s[i: i + lcp[k]]      # O(n)

7.3 最长公共子串（两个串）

把 T = S1 + '#' + S2（# 为不出现在两串中的分隔符）构造 SA+LCP，答案为"跨越分隔符"的最大 LCP：

def longest_common_substring(s1, s2):
    t = s1 + chr(1) + s2          # chr(1) 作分隔符
    n1 = len(s1)
    sa = prefix_doubling_sa(t)
    lcp = kasai_lcp(t, sa)
    best, bi = 0, -1
    for k in range(1, len(t)):
        a, b = sa[k - 1], sa[k]
        if (a < n1) != (b < n1):   # 一个在 S1、一个在 S2
            if lcp[k] > best:
                best, bi = lcp[k], a
    return t[bi: bi + best]

7.4 不同子串计数（Distinct Substrings）

每个后缀贡献 n - SA[i] 个前缀，减去与前一个后缀重复的 LCP[i] 个，总和为：

distinct = n*(n+1)/2 - sum(LCP[i] for i in 1..n-1)

O(n) 即得，是数据去重指纹、DNA k-mer 多样性分析常用指标。

7.5 字典序第 k 小子串 / BWT

• 第 k 小子串：在 SA 上按 LCP 做前缀长度的前缀和，二分定位。
• BWT（Burrows–Wheeler Transform）：BWT[i] = S[(SA[i] - 1) mod n]。BWT 配合 FM-index 可在仅存 BWT 的情况下做 backward search，是压缩全文索引的核心；而 BWT 完全由 SA 派生，可见 SA 的基础地位。

八、与后缀树 / 后缀自动机的取舍

维度	SA + LCP	后缀树	后缀自动机 (SAM)
空间	2~3 × n int（~12–24B/char）	~10–20 × n 指针（重）	~4–8 × n 状态
构建	O(n log n) 排序 / O(n) SA-IS	O(n)	O(n)
模式匹配	O(m log n)（二分）	O(m)	O(m)
最长重复/公共子串	O(n)（LCP 扫描）	O(n)	O(n)
不同子串计数	O(n)	O(n)	O(n)
在线增量更新	困难（需重建）	困难	支持（在线添加字符）
工程难度	低（数组友好、缓存友好）	高（指针结构、调试难）	中（状态机概念门槛）

选型结论：追求内存效率、可缓存、易并行（GPU/分布式）、需对接 FM-index 时选 SA+LCP；需要在线增量插入字符（流式）时选 SAM；需要显式树形结构做拓扑查询时选后缀树。绝大多数"离线批量索引 + 查询"场景，SA 是性价比之王。

九、工程化考量

1. 大字母表：DNA 仅 4 字符，计数排序桶少；Unicode/中文则字符集大，需把"字符"映射为 0..σ-1 的紧凑整数（σ 为去重后字符数），否则计数排序退化。
2. 整数后缀 vs 字符串：对整数数组（如基因组编码、token id 序列）构造 SA 时，直接以整数为键做计数排序，避免字符切片开销。
3. 外部内存 / 压缩：十亿级数据无法装入内存，需 external-memory SA 或 compressed suffix array / FM-index（只存采样 rank + Wavelet Tree），用少量随机 IO 换取海量存储。
4. 与 FM-index 集成：SA 派生 BWT，BWT + 采样 SA 即 FM-index，可做"只存压缩索引、支持任意模式查询"，是大规模日志/基因检索的标准方案。
5. 线程并行：前缀倍增每轮排序可并行化（计数排序多线程），SA-IS 的诱导排序亦可分块；但 LCP 的 Kasai 是顺序依赖，瓶颈常在 LCP。
6. 哨兵字符：排序法中务必加一个小于所有字符的哨兵（如 chr(0)），保证后缀互不相等的不变式，否则 rank 区分会出错。

十、12 项生产陷阱清单

#	陷阱	后果	规避
1	漏加哨兵字符	相等后缀导致 rank 死循环/错位	统一追加 chr(0) 或小哨兵
2	用比较排序而非计数排序	退化到 O(n log² n) 甚至更慢	以 rank 值域做基数/计数排序
3	大字母表未紧凑编码	计数排序桶爆炸 O(σ)	预映射字符到 0..σ-1
4	LCP 误用朴素两两比较	O(n²) 超时	必须用 Kasai O(n)
5	二分匹配每次切片 s[sa[i]:]	单次 O(n) 拖垮二分	写"按前缀比较 SA[i] 与 P"
6	最长公共子串未加分隔符	跨串 LCP 误判	拼接 # 分隔符再构造
7	分隔符出现在原串中	查询串跨越假边界	选原串未出现的码点
8	不同子串计数溢出	大 n 时 n*(n+1)/2 超 int	用 64 位整数
9	rank/SA 下标混用	off-by-one	明确 SA[k]=位置、rank[位置]=k
10	在线增量更新强求 SA	每次重建 O(n log n)	改用语境合适的 SAM
11	十亿级数据强行内存 SA	OOM	改用 compressed SA / FM-index
12	自写 SA-IS 直接上生产	边界 bug 难查	复用 divsufsort 等成熟库

十一、可复现工具箱（Python）

下面给出一个自包含的 SuffixArray 工具箱，集成前缀倍增构造、Kasai LCP、模式计数、最长重复/公共子串与不同子串计数，可直接用于中小规模（n < 10⁷）生产：

class SuffixArray:
    def __init__(self, s: str):
        self.s = s + chr(0)            # 哨兵
        self.sa = self._build(self.s)
        self.rank = [0] * len(self.s)
        for i, p in enumerate(self.sa):
            self.rank[p] = i
        self.lcp = self._kasai(self.s, self.sa)

    def _build(self, s):               # 前缀倍增 O(n log n)
        n = len(s); k = 1
        rank = [ord(c) for c in s]; tmp = [0] * n
        while True:
            pairs = sorted(range(n),
                key=lambda i: (rank[i], rank[(i + k) % n]))
            tmp[pairs[0]] = 0
            for i in range(1, n):
                a, b = pairs[i - 1], pairs[i]
                tmp[b] = tmp[a] + ((rank[a], rank[(a + k) % n])
                                   != (rank[b], rank[(b + k) % n]))
            rank, tmp = tmp, rank
            if rank[pairs[-1]] == n - 1:
                break
            k <<= 1
        return pairs

    def _kasai(self, s, sa):
        n = len(s); rank = [0] * n
        for i, p in enumerate(sa):
            rank[p] = i
        lcp = [0] * n; h = 0
        for i in range(n):
            if rank[i] == 0:
                h = 0; continue
            j = sa[rank[i] - 1]
            while i + h < n and j + h < n and s[i + h] == s[j + h]:
                h += 1
            lcp[rank[i]] = h
            if h > 0:
                h -= 1
        return lcp

    def count(self, p: str):           # 模式出现次数 O(m log n)
        n = len(self.sa)
        lo = self._lower(p); hi = self._lower(p + chr(127))
        return hi - lo

    def _lower(self, p):                # SA 上按前缀二分的下界
        lo, hi = 0, len(self.sa)
        while lo < hi:
            mid = (lo + hi) // 2
            if self.s[self.sa[mid]:].startswith(p) or \
               self.s[self.sa[mid]:] < p:
                lo = mid + 1
            else:
                hi = mid
        return lo

    def longest_repeated(self):
        k = max(range(len(self.lcp)), key=lambda i: self.lcp[i])
        i = self.sa[k]
        return self.s[i: i + self.lcp[k]]

    def distinct_substrings(self):
        n = len(self.s) - 1
        return n * (n + 1) // 2 - sum(self.lcp)

提示：_lower 中的切片比较为示意；生产环境应替换为"按前缀逐字符比较 SA[i] 与 P"以避免 O(n) 切片，把单次二分降到 O(m)。

结语

后缀数组用最朴素的整数数组表达了字符串的全部后缀序关系，配合 LCP 即等价于后缀树，却把空间压到 2~3 个 int 数组、把工程难度降到"排序 + 一次线性扫描"。从 O(n² log n) 的朴素理解，到 O(n log n) 的前缀倍增，再到 O(n) 的 SA-IS 与 O(n) 的 Kasai LCP，这条优化主线本身就是"用已算好的秩替代重复比较"这一工程哲学的最佳注脚。把它纳入你的字符串处理工具箱，从日志检索、基因组比对到数据去重，都会多一把趁手的利器。