[字符串相关]后缀数组 - 倍增算法

我们约定,下文所提到的字符串下标都从 \(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

posted @ 2021-07-06 21:50  Dfkuaid  阅读(127)  评论(0编辑  收藏  举报