「字符串相关」后缀数组 - 倍增算法

我们约定，下文所提到的字符串下标都从 \(0\) 开始。下文中的“第 \(x\) 后缀”指的是从下标 \(x\) 开始的后缀。

对于后缀数组的学习，本人建议可以自己随便写一个字符串，然后按照相应的过程进行模拟，会很方便的理解每一步的作用以及原理。

#1.0 何为后缀数组

#1.1 后缀树

我们知道，像 \(\texttt{AC}\) 自动机 这样的多模式串匹配需要预处理的是模式串，但是在很多时候，我们会先知道文本串，而非模式串，对于每个模式串的询问要求在线回答，那怎么办？我们引入后缀树这个概念。

所谓后缀树，就是将一个字符串的所有后缀建成一棵树，设 \(S='abaca'\)，\(S\) 的后缀树为：

$ 符号是为了标记每一个后缀的结束位置，每个叶节点下面的数字表示该后缀从原字符串的第几位开始，我们已经把所有后缀按字典序从左向右摆放（默认 $ 符号的字典序最小），显然从上面寻找模式串是非常方便的，但我们发现整个后缀树的结构非常庞大，有没有合适的替代品呢？

#1.2 后缀数组

将上面的叶结点下的数字从左往右取出放在数组中，便得到了后缀数组（\(\texttt{Suffix Array}\)），

由此便不难得到后缀数组的定义：所有的后缀按字典序排序后的编号数组。

#2.0 后缀数组的构建

#2.1 朴素算法

一个很简单的算法便是直接采用快排 sort()，快排的时间复杂度为 \(O(n\log n)\)，但由于字符串比较需要按位比较，所以实际时间复杂度为 \(O(n^2\log n)\)，显然不能接受。

#2.2 倍增算法

这里采用倍增的双关键字排序。大致流程如下：

1. 给所有后缀的第一个字符排序；
2. 给所有后缀的前两个字符排序，实际是双关键字排序；
3. 枚举长度 \(l=2^k\)，给所有后缀的前 \(l\) 个字符排序，实际依旧是双关键字排序，因为可以借助上一次排序的结果（\(l=2^k=2^{k-1}+2^{k-1}\)），即错开 \(l-1\) 位进行双关键字排序；重复这一步直到所有后缀的排名不同。

显然，上面的过程最多只需要 \(\log n\) 次，如果这里的双关键字排序采用快排的话，时间复杂度为 \(O(n\log^2 n)\)

#2.3 基数排序优化

我们注意到最多只有 \(n\) 个后缀，所以可以采用基数排序（我们常用的桶排序便是单关键字的基数排序），这里需要采用双关键字的基数排序。

直接讲解并没有非常好讲，所以这里先贴出代码，我们来一行一行的分析。

char s[N];
int sa[N], t[N], t2[N], c[N], n;

inline void build_sa(int m) {
    //m 是编号的最大值
    int i, *x = t, *y = t2;
    for (i = 0; i < m; ++ i) c[i] = 0;
    for (i = 0; i < n; ++ i) c[x[i] = s[i]] ++;
    for (i = 1; i < m; ++ i) c[i] += c[i - 1];
    for (i = n - 1; i >= 0; i --) sa[-- c[x[i]]] = i;
    for (int k = 1; k <= n; k <<= 1) {
        int p = 0;
        for (i = n - k; i < n; i ++) y[p ++] = i;
        for (i = 0; i < n; i ++) if (sa[i] >= k)
          y[p ++] = sa[i] - k;
        for (i = 0; i < m; ++ i) c[i] = 0;
        for (i = 0; i < n; ++ i) c[x[y[i]]] ++;
        for (i = 1; i < m; ++ i) c[i] += c[i - 1];
        for (i = n - 1; i >= 0; i --)
          sa[-- c[x[y[i]]]] = y[i];
        swap(x, y); p = 1; x[sa[0]] = 0;
        for (i = 1; i < n; ++ i)
          x[sa[i]] = y[sa[i-1]] == y[sa[i]] && y[sa[i - 1] + k] == y[sa[i] + k] ? p - 1 : p ++;
        if (p >= n) break;
        m = p;
    }
}

UPD on 10.30.21: 上述代码可能存在问题，建议使用下述代码（思路一致）

char s[N];
int sa[N], rk[N], oldrk[N], px[N], id[N], cnt[N], n;

bool comp(int x, int y, int w) {
    return oldrk[x] == oldrk[y] && oldrk[x + w] == oldrk[y + w];
}

inline void build_sa(int m) {
    int i, p, w;
    for (i = 1; i <= n; ++ i) ++ cnt[rk[i] = s[i]];
    for (i = 1; i <= m; ++ i) cnt[i] += cnt[i - 1];
    for (i = n; i >= 1; -- i) sa[cnt[rk[i]] --] = i;
    for (w = 1;; w <<= 1, m = p) {
        for (p = 0, i = n; i > n - w; --i) id[++ p] = i;
        for (i = 1; i <= n; ++ i) if (sa[i] > w) id[++ p] = sa[i] - w;
        memset(cnt, 0, sizeof(cnt));
        for (i = 1; i <= n; ++ i) ++ cnt[px[i] = rk[id[i]]];
        for (i = 1; i <= m; ++ i) cnt[i] += cnt[i - 1];
        for (i = n; i >= 1; -- i) sa[cnt[px[i]] --] = id[i];
        memcpy(oldrk, rk, sizeof(rk));
        for (p = 0, i = 1; i <= n; ++ i)
          rk[sa[i]] = comp(sa[i], sa[i - 1], w) ? p : ++ p;
        if (p == n) {for (int i = 1; i <= n; ++ i) sa[rk[i]] = i; break;}
    }
}

---

首先是第一部分

for (i = 0; i < m; ++ i) c[i] = 0;
for (i = 0; i < n; ++ i) c[x[i] = s[i]] ++;
for (i = 1; i < m; ++ i) c[i] += c[i - 1];
for (i = n - 1; i >= 0; i --) sa[-- c[x[i]]] = i;

这里是先对每个后缀的首字母进行排序，x[i] 在这里表示的是第 i 个字符的编号（此时的第一关键字），sa[i] 表示的是按第一关键字排序后的第 i 个后缀在原字符串的起始位置。

---

然后是循环中的这一部分

int p = 0;
for (i = n - k; i < n; i ++) y[p ++] = i;
for (i = 0; i < n; i ++) if (sa[i] >= k)
  y[p ++] = sa[i] - k;

这一部分是按第二关键字排序，y[i] 储存的是按第二关键字排序后第 i 个后缀在原字符串上的开始位置。

其中第一个 for 是将没有第二关键字的后缀的第二关键字设置为无限小，所以要放在 y[] 数组的最前面。

考虑第二个 for 为什么可以直接借助 sa[] 数组进行赋值，我们知道，sa[i] 表示的是在长度为 \(\frac k 2\) 时，按当时的双关键字排序后的第 i 个后缀在原字符串的起始下标，这时对应的也就是在此时长度为 \(k\) 的情况下按第一关键字排序后的第 i 个后缀在原字符串的起始下标。

而 sa[i] - k 就相当于按第一关键字排序后的第 i 个后缀在原字符串的下标的第前 k 个位置，那么如果将 sa[i] 看作第二关键字，那么 sa[i] - k 便是所对应的第一关键字在原字符串中的下标位置，那么数组 y[] 中的数恰好便是按第二关键字排序后的第 i 个后缀在原字符串的起始下标。

---

循环中的第二大部分梦幻四重循环

for (i = 0; i < m; ++ i) c[i] = 0;
for (i = 0; i < n; ++ i) c[x[y[i]]] ++;
for (i = 1; i < m; ++ i) c[i] += c[i - 1];
for (i = n - 1; i >= 0; i --) sa[-- c[x[y[i]]]] = y[i];

这里 x[i] 表示的是第 \(i\) 后缀在按第一关键字排序后对应的 sa[] 数组上的排名（注意这样的排名的大小只是相对性的，可能并不连续）。

在整个算法的大部分时间里，x[] 都是从后缀映射到排名，sa[] 则是从排名映射到后缀。

不难发现这里和第一部分几乎一模一样。这里是结合两个关键字的分别排序，构建出当前的双关键字排序后的结果，也是下一次长度为 \(2k\) 时按第一关键字排序后的结果。

原因不难想到：y[] 本身具有顺序性，即第二关键字小的在前面，而本身这个排序则是将第一关键字小的排在前面，而不改变第一关键字相同的数之间的相对位置，综合起来，得到的答案便是按双关键字排序后的结果。

---

这是循环中的最后一部分

swap(x, y); p = 1; x[sa[0]] = 0;
for (i = 1; i < n; ++ i)
  x[sa[i]] = y[sa[i-1]] == y[sa[i]] && y[sa[i - 1] + k] == y[sa[i] + k] ? p - 1 : p ++;
if (p >= n) break;
m = p;

注意在开始前我们先交换了指针 x 和指针 y 所指向的地址，所以此时 x[i] 指的是按第二关键字排序后的第 i 个后缀在原字符串的起始下标，这个数据对于我们来说已经没有用处了；y[i] 表示的是第 \(i\) 后缀在按第一关键字排序后对应的 sa[] 数组上的排名（注意这样的排名的大小只是相对性的，可能并不连续）。

这里是统计按现在的结果，还会有多少个不同的排名，如果排名数达到了后缀的总数，那么就没有必要继续进行下去了，如果不够，那么同时为各个不同的后缀重新编号（记录按双关键字排序后对应的 sa[] 数组上的排名），这样编号最大为 \(p-1\)。

中间的三目运算符写成 if 如下

if (y[sa[i-1]] == y[sa[i]] && y[sa[i - 1] + k] == y[sa[i] + k])
  x[sa[i]] = p - 1;
else x[sa[i]] = p ++;

这里判等的条件分别是按长度为 \(2k\) 时的第一关键字相等（y[sa[i-1]] == y[sa[i]]）和第二关键字相等（y[sa[i - 1] + k] == y[sa[i] + k]），至于为什么 sa[] 的下标是连续的：两个双关键字都相等，而 sa[i] 中是当前双关键字排序后时排名为 \(i\) 的后缀的起始位置，双关键字排序的要求便是先满足第一关键字从小到大，再满足第二关键字从小到大，所以如果两个都相等，那么在 sa[] 中的排名必然是相邻的。

笔者中午在宿舍里思考这里时脑子抽了，是这样断的句：(((x[sa[i]] = y[sa[i-1]]) == y[sa[i]]) && (y[sa[i - 1] + k)] == y[sa[i] + k]) ? p - 1/*错误地以为是 p --*/ : p ++，当时愣是搞不明白为啥这样做是对的，吹起床哨了才发现我™又断错句了. QnQ

---

整体来说比较复杂，一遍不懂就需要多啃几遍。

#3.0 后缀数组的应用

#3.1 在线模式串匹配

采用二分查找的方式，在文本串中寻找是否存在一个位置开始的后缀的前缀是模式串。

我们可以二分所属的后缀的排名，并与模式串比较大小，如果小了就让排名增大，反之亦然。

设模式串长度为 \(m\)，文本串长度为 \(n\)，那么单次询问时间复杂度为 \(O(m\log n).\)

int m;
int cmp_suffix(char* pattern, int p) {
    return strncmp(pattern, s + sa[p], m);
}

int find(char* P) {
    m = strlen(P);
    if (cmp_suffix(P, 0) < 0) return -1;
    if (cmp_suffix(P, n - 1) > 0) return -1;
    int l = 0, r = n - 1;
    while (r >= l) {
        int mid = l + (r + l) / 2;
        int res = cmp_suffix(P, mid);
        if (!res) return mid;
        if (res < 0) r = mid - 1;
        else l = mid + 1;
    }
    return -1;
}

参考资料

[1] 刘汝佳, 陈锋. 算法竞赛入门经典：训练指南. 北京: 清华大学出版社, 2012.

[2] 后缀数组——倍增算法 - hyl天梦

[3] 后缀数组简介 - OI Wiki