SA

· · 算法·理论

引入及定义

令字符串为 s,长度为 n,下标从 1 开始。s[i\dots n] 这一段后缀称为「后缀 i」或者「第 i 个后缀」。

后缀数组(Suffix Array, SA)是解决一类与字符串后缀有关问题的算法。

两个基本的数组为 sa_i,rk_i。其中 sa_i 表示排名为 i 的后缀是哪个,rk_i 表示后缀 i 的排名。这里的排名指,将 s 的所有的 n 个后缀取出,按字典序从小到大排序后,某个后缀的排名。

不难发现 sa_{rk_i}=rk_{sa_i}=i,即两数组互逆。

这个数组求解后有什么意义呢?注意到后缀 sa_i 和后缀 sa_{i+1} 的「相似度」是很高的,因为它们的 LCP(最长公共前缀)一般是很大的。通过这样的排序可以方便解决一些与后缀 LCP 有关的问题。

下面介绍后缀数组的求解方法,以及相关应用。

求解

求解后缀数组的算法中,最常见的是倍增+基数排序的 \mathcal O(n \log n) 解法。尽管存在 \mathcal O(n) 的科技,其冗长的代码和极大的常数往往不是人们的第一选择。我们介绍的是 \mathcal O(n \log n) 解法。

首先每个后缀的长度都是不同的,最小是 1,最大是 n。但在字典序比较时,如果一个串是另一个的真前缀,那么这个短的串是小于长的串的。因此在比较两个后缀时,我们在字符串的最后补上若干个空字符,且这个空字符是小于任意一个其它字符的。这样比较时可以默认两字符串长度相等,按位比较即可。

当然实现上并不需要体现这一点,因为数组默认是 0

考虑进行一个增的倍。从小到大枚举 k,并对 n 个长度为 2^k 的字串排序。如同上面所说,最后几个长度不满 2^k 的字串我们会在最后添加空字符。不难发现当 k \ge \log n 时, 排好序的长度 2^k 的字串,就是 s 的所有后缀。

而在某一轮中排序也很简单。我们知道 s[x,x+2^i-1] 可以看作 s[x,x+2^{i-1}-1]+s[x+2^{i-1},x+2^i-1]。所以比较 s[x,x+2^i-1]s[y,y+2^i-1] 时,我们将它们分成两半,以前半段为第一关键字,后半段为第二关键字比较即可。

怎样分别比较前半段和后半段?它们的长度都是 2^{i-1},上一轮倍增中已经求出排名了,排名大小可以反映字典序大小。

这里很重要的一点是,如果两个不同位置的字串完全相同,那么它们的排名应该相同。

于是可以写出这样的代码:

for (int k = 1; k <= n; k <<= 1) {
    iota(sa + 1, sa + n + 1, 1);
    sort(sa + 1, sa + n + 1,
        [&](int x, int y) {
            return rk[x] == rk[y] ? rk[x + k] < rk[y + k] : rk[x] < rk[y];
        });
    memcpy(b, rk, sizeof b);
  int m = 0;
    for (int i = 1; i <= n; ++ i )
        if (i != 1 && b[sa[i]] == b[sa[i] - 1] && b[sa[i] + k] == b[sa[i - 1] + k]) {
            rk[sa[i]] = m;
        } else {
            rk[sa[i]] = ++ m;
        }
}

由于内部的排序以及外部的倍增,总复杂度是 \mathcal O(n \log ^2n) 的。

注意到排序是双关键字排序,且每个关键字的值域都不超过 n。这提醒我们使用基数排序。

在基数排序中,首先将数据按第二关键字排序,然后按第一关键字做稳定排序。这里稳定是指不改变第一关键字相同的元素排序前后的相对位置关系。这很类似二维数点中将数据离线后按某一维排序,将这一维转化为时间维进行下一步对另一维的操作。

因为值域不超过 n,所以按某一关键字排序时可以直接捅排序。这里我们并不要得到排序后的数组,而是排序后与排序前一一对应的位置,即 sark 数组。

我们将所有数据按照某一关键字拍到桶(cnt)中后,对桶做一遍前缀和,此时 cnt_i 即表示 \le i 的数的个数。

那么所有等于 i 的数的排名即 cnt_{i-1}+1 \sim cnt_i

此时再遍历一遍原数组。对于某个数 x 而言,每次遇到它时,将当前位置的排名赋为 cnt_x,并接着将 cnt_i 减一。

因为要做稳定排序,但是上面遇到相同数时,越靠前的排名越大,所以最后赋排名时倒过来枚举。

即:

// 按 a[i] 排序,p[i] 表示排名为 i 的是第几个
memset(cnt, 0, sizeof cnt);
for (int i = 1; i <= n; ++ i ) cnt[a[i]] ++ ;
for (int i = 1; i <= m; ++ i ) cnt[i] += cnt[i - 1];        // m 是值域
for (int i = n; i; -- i ) p[cnt[a[i]] -- ] = i;

把两边这个东西替换掉上面的 std::sort 即可。

然后是一些常数优化。首先如果某一轮循环结束后 m=n,即 rk 数组的值域等于 n,这代表我们已经将「以每个点开头的一段前缀」完全区分开了,那么直接退出即可。

然后我们考虑按 rk_{i+k} 为关键字做的第一轮排序。首先当 n - k + 1 \le i \le n 时,s_{i+k} 是空的,即这些 i 应该是排序后最靠前的。

for (int i = n - k + 1; i <= n; ++ i ) id[ ++ nums] = i;

对于剩下的 i,注意到上一轮排名为 i 的后缀是 sa_i,那么这一轮是 sa_i - k。因为我们要比较 rk_{sa_i-k+k} 的大小。即:

for (int i = 1; i <= n; ++ i ) if (sa[i] > k) id[ ++ nums] = sa[i];

因此总代码为:

m = 130;

for (int i = 1; i <= n; ++ i ) cnt[rk[i] = s[i]] ++ ;
for (int i = 2; i <= m; ++ i ) cnt[i] += cnt[i - 1];
for (int i = n; i; -- i ) sa[cnt[rk[i]] -- ] = i;

for (int k = 1; k <= n; k <<= 1) {
    nums = 0;
    for (int i = n - k + 1; i <= n; ++ i ) id[ ++ nums] = i;
    for (int i = 1; i <= n; ++ i )
        if (sa[i] > k) id[ ++ nums] = sa[i] - k;

    memset(cnt, 0, sizeof cnt);
    for (int i = 1; i <= n; ++ i ) cnt[rk[i]] ++ ;
    for (int i = 2; i <= m; ++ i ) cnt[i] += cnt[i - 1];
    for (int i = n; i; -- i ) sa[cnt[rk[id[i]]] -- ] = id[i];

    memcpy(b, rk, sizeof b);
    m = 0;
    for (int i = 1; i <= n; ++ i )
        if (b[sa[i]] == b[sa[i - 1]] && b[sa[i] + k] == b[sa[i - 1] + k]) {
            rk[sa[i]] = m;
        } else {
            rk[sa[i]] = ++ m;
        }
    if (m == n) break;
}

高度数组

我们定义高度数组 h_i = \operatorname{LCP}(i-1,i)。其中 \operatorname{LCP}(i,j) 表示后缀 sa_i 和后缀 sa_j 的最长公共前缀(注意不是后缀 i 和后缀 j 的最长公共前缀,它是 \operatorname{LCP}(rk_i,rk_j))。接下来我们将介绍高度数组的性质。

引理 1:\operatorname{LCP}(i,j) = \min \operatorname{LCP}(i,k), \operatorname{LCP}(k,j)\forall i \le k \le j

即将所有后缀排序后,任意两个后缀的最长公共前缀,是它们之间的某个后缀与自己的最长公共前缀的最小值。

证明之前明确一个事实。从 x,y 往后数 k 个字符都相同,等价于 \operatorname{LCP}(x,y) \ge k

考虑分别证明不大于和不小于。

  1. 若 $i \le k \le j$,那么从 $s$ 的 $sa_i,sa_j,sa_k$ 位置往后数 $\operatorname{LCP}(i,j)$ 个字符形成的子串都是相同的。这个想一下就很清楚而且严谨证明也不难(夹逼定理(?))。 所以 $s[sa_i,sa_i + \operatorname{LCP}(i,j)-1] = s[sa_k,sa_k + \operatorname{LCP}(i,j)-1]$,那么一定有 $\operatorname{LCP}(i,k) \ge \operatorname{LCP}(i,j)$。$\operatorname{LCP}(k,j)$ 同理。
  2. 根据定义容易得到,从 $i,k$ 往后数 $\operatorname{LCP}(i,k)$ 个字符是相同的,从 $k,j$ 往后数 $\operatorname{LCP}(k,j)$ 个字符是相同的,所以从 $i,j$ 往后数 $\min \operatorname{LCP}(i,k),\operatorname{LCP}(k,j)$ 个字符是相同的。所以 $\operatorname{LCP}(i,j) \ge \min \operatorname{LCP}(i,k),\operatorname{LCP}(k,j)$。

由此可得:

引理 1 推论:\operatorname{LCP}(i,j) = \min \operatorname{LCP}(i,i+1),\operatorname{LCP}(i+1,i+2),\dots,\operatorname{LCP}(j-1,j)

进一步的 \operatorname{LCP}(i,j) \le \operatorname{LCP}(j-1,j)

这个是显然的,归纳一下即可。

引理 2:h_{rk_i} \ge h_{rk_{i-1}}-1

这里会有点绕。h_{rk_i} 其实就是后缀 i 和比他小的最大后缀的 \operatorname{LCP}。同理 h_{rk_{i-1}} 是后缀 i-1 和比他小的最大后缀的 \operatorname{LCP}

设「比后缀 i-1 小的最大后缀」是后缀 j。即 s[j\dots n] < s[i-1\dots n]。根据定义 h_{rk_{i-1}} 就是后缀 i-1 和后缀 j\operatorname{LCP}。稍作讨论:

根据引理 2 可以得到一个简单的线性求解高度数组的方法:i 从小到大枚举,依次求解 h_{rk_i}。初始将 h_{rk_i} 设为 h_{rk_{i-1}}-1,然后暴力加,直到不合法。

for (int i = 1, k = 0; i <= n; ++ i ) {
  if (rk[i] == 1) h[rk[i]] = 0;
  else {
    if (k) -- k;
    int j = sa[rk[i] - 1];
    while (i + k <= n && j + k <= n && s[i + k] == s[j + k]) k ++ ;
    h[rk[i]] = k;
  }
}