【模板】后缀数组 SA

posted on 2022-08-07 17:04:40 | under 模板 | source

目录

\(LCP(S, T)\) 是两个字符串 \(S, T\) 的最长公共前缀（The Longest Common Prefix）。

后缀数组

简单来说，后缀数组需要一个长为 \(n\) 的字符串 \(S\)，它求出了：

  • \(sa_i\)，将 \(S\) 的所有真后缀按字典序排序，排在第 \(i\) 位的后缀（的起始点，显然 \(1\leq sa_i\leq n\)，以下默认如此）。
  • \(rnk_i\)，第 \(i\) 个后缀的排名，显然有 \(rnk_{sa_i}=i\)。
  • \(hei_i\)，为 \(lcp(i-1,i)\)，即排序后第 \(i\) 小后缀和第 \(i-1\) 小后缀的最长公共前缀。此处 \(lcp(i,j)=LCP(S[sa_i,n],S[sa_j,n])\)，即这两个后缀的 LCP。

实现时使用倍增，令当前的排序长度为 \(2^j\)，每次将上一轮的 \(rnk\) 合并在一起就可以比较两个（第 \(i\) 个后缀可以被 \((rnk_i,rnk_{i+2^{j-1}})\) 表示，从而可以 \(O(1)\) 比较）。使用桶排序，总复杂度 \(O(n\log n)\)。

对于 \(hei_i\) 数组，我们求一个 \(h_i=hei_{rnk_i}\)（\(h_i\) 表示后缀 \(S[i, n]\) 与其在 \(sa\) 中上一位的后缀的 LCP），由题解得 \(h_{i-1}-1\leq h_i\)，于是考虑像 manacher 一样维护算出 \(h_i\)，进而得到 \(hei_i\)。显然是线性复杂度 \(O(n)\)。

code

点击查看代码 
template <int N, int M = 128>
struct SA {
  int sa[N + 10], rnk[N + 10], h[N + 10], hei[N + 10];
  void bucketsort(int n) {
    memset(hei, 0, sizeof hei);
    for (int i = 1; i <= n; i++) hei[rnk[h[i]]]++;
    for (int i = 1; i <= max(n, M); i++) hei[i] += hei[i - 1];
    for (int i = n; i >= 1; i--) sa[hei[rnk[h[i]]]--] = h[i];
  }
  void build(char *s) {
    int n = strlen(s + 1);
    rnk[n + 1] = 0;
    for (int i = 1; i <= n; i++) h[i] = i, rnk[i] = s[i];
    bucketsort(n);
    for (int j = 1; j <= n; j <<= 1) {
      int tot = 0;
      for (int i = n - j + 1; i <= n; i++) h[++tot] = i;
      for (int i = 1; i <= n; i++)
        if (sa[i] - j > 0) h[++tot] = sa[i] - j;
      bucketsort(n);
      tot = h[sa[1]] = 1;
      for (int i = 2; i <= n; i++)
        h[sa[i]] =
            rnk[sa[i - 1]] == rnk[sa[i]] && rnk[sa[i - 1] + j] == rnk[sa[i] + j]
                ? tot
                : ++tot;
      memcpy(rnk, h, sizeof h);
    }
    memset(hei, 0, sizeof hei);
    for (int i = 1, k = 0; i <= n; i++) {
      if (rnk[i] == 1) continue;
      int j = sa[rnk[i] - 1];
      k = max(k - 1, 0);
      while (max(i, j) + k <= n && s[i + k] == s[j + k]) k++;
      hei[rnk[i]] = k;
    }
  }
};
点击查看代码 
int n, sa[50010], rnk[50010], hei[50010];
char a[50010];
namespace suffixsa {
void getsa() {
  static int tmp[50010];
  auto bucketsort = [&](int v) {
    auto& buc = hei;
    for (int i = 0; i <= v; i++) buc[i] = 0;
    for (int i = 1; i <= n; i++) buc[rnk[tmp[i]]] += 1;
    for (int i = 1; i <= v; i++) buc[i] += buc[i - 1];
    for (int i = n; i >= 1; i--) sa[buc[rnk[tmp[i]]]--] = tmp[i];
  };
  for (int i = 1; i <= n; i++) tmp[i] = i, rnk[i] = a[i];
  bucketsort(128);
  rnk[n + 1] = 0;
  for (int j = 1; j <= n; j <<= 1) {
    int cnt = 0;
    for (int i = n - j + 1; i <= n; i++) tmp[++cnt] = i;
    for (int i = 1; i <= n; i++)
      if (sa[i] > j) tmp[++cnt] = sa[i] - j;
    bucketsort(n);
    cnt = tmp[sa[1]] = 1;
    for (int i = 2; i <= n; i++) {
      if (rnk[sa[i - 1]] != rnk[sa[i]] || rnk[sa[i - 1] + j] != rnk[sa[i] + j])
        cnt += 1;
      tmp[sa[i]] = cnt;
    }
    for (int i = 1; i <= n; i++) rnk[i] = tmp[i];
  }
}
void getlcp() {
  int h = 0;
  hei[1] = 0;
  for (int i = 1; i <= n; i++) {
    if (rnk[i] == 1) continue;
    int j = sa[rnk[i] - 1];
    if (h) h -= 1;
    while (max(i, j) + h <= n && a[i + h] == a[j + h]) h += 1;
    hei[rnk[i]] = h;
  }
}
};  // namespace suffixsa
点击查看代码（完全 vector） 
vector<int> get_sa(const string& a) {
  int n = a.size();
  vector<int> sa(n), rk(a.begin(), a.end()), h(n);
  auto bucketsort = [&](const auto& key) {
    vector<int> buc(max(n, 128));
    reverse(h.begin(), h.end());
    for (int i : h) buc[rk[i]] += 1;
    partial_sum(buc.begin(), buc.end(), buc.begin());
    for (int i : h) sa[--buc[rk[i]]] = i;
    vector<int> tmp(n);
    tmp[sa[0]] = 0;
    for (int i = 1; i < n; i++) tmp[sa[i]] = tmp[sa[i - 1]] + (key(sa[i - 1]) != key(sa[i]));
    rk = move(tmp);
  };
  for (int i = 0; i < n; i++) h[i] = i;
  bucketsort([&](int x) { return rk[x]; });
  for (int j = 1; j < n; j <<= 1) {
    h.clear();
    for (int i = n - j; i < n; i++) h.push_back(i);
    for (int i = 0; i < n; i++) if (sa[i] >= j) h.push_back(sa[i] - j);
    bucketsort([&](int x) { return make_pair(rk[x], x + j < n ? rk[x + j] : -1); });
  }
  return sa;
}
vector<int> get_lcp(const string& a, const vector<int>& sa) {
  int n = a.size();
  vector<int> rk(n), lcp(n - 1);
  for (int i = 0; i < n; i++) rk[sa[i]] = i;
  for (int i = 0, h = 0; i < n; i++) {
    if (!rk[i]) continue;
    if (h) h -= 1;
    int j = sa[rk[i] - 1];
    while (max(i, j) + h < n && a[i + h] == a[j + h]) h += 1;
    lcp[rk[i] - 1] = h;
  }
  return lcp;
}

证明 \(rnk_{sa_i+j}\) 只会越界到 \(n+1\)

您写出了这个代码，发现您的 \(rnk\) 开了一倍且写出了 rnk[sa[i-1]+j]==rnk[sa[i]+j] 这样的东西，这不会越界吗？

证明

定义一个空字符为 \(0\)，且不与其它任何字符相同。

假设我们枚举到 \(n=8,j=4\)，分类讨论：

  • 有三个后缀（绿色）超出了原串范围。但是，这些后缀的 \(rnk\) 互不相同，你找不到其它任何一个后缀和它们相等，因为空字符和其它字符都不相同，所以 rnk[sa[i-1]]==rnk[sa[i]] 必然不成立，&& 立即短路，不会越界。
  • 有一个刚好到达串尾的后缀（橙色）看起来很危险，它有可能和前面一个后缀相等，但是它的起始位置是 \(n-j+1\)，右移 \(j\) 位是 \(n+1\)，只会越界到 \(n+1\)，因此我们只需要初始化 \(rnk_{n+1}=0\) 即可。
  • 剩余的四个后缀没有超出原串范围（红色），没有越界风险。

归纳证明，易得其它情况亦成立，故 \(rnk_{sa_i+j}\) 只会越界到 \(n+1\)，原命题得证。

解决方法

由上文得，我们只需要令 \(rnk_{n+1}=0\)。保险起见，你的字符集最好不要出现 \(0\)。

构造数据

只需要一个长为 \(2^j\) 的后缀，它在原串的出现次数大于一，就会出现以上情况。一个形如 aaaaaaa......aaa 或 ababababab......ababab 可以尽量多地卡到上限。

求两个后缀的 LCP

性质：有相同前缀的后缀，在 \(sa\) 中是一段连续区间。

LCP Lemma：\(lcp(i,j)=\min\limits_{i\leq k\leq j}\{lcp(i,k),lcp(k,j)\}\)。
LCP Theorem：\(lcp(i,j)=\min\limits_{i<k\leq j}lcp(k-1,k)=\min\limits_{i<k\leq j}hei_k\)。

于是问题变为在 \(hei_i\) 上做 RMQ（\(hei[rnk_l+1,rnk_r]\) 的最小值）。使用 ST 表可以做到 \(O(n\log n)-O(1)\)。

求 \(S[l,r]\) 在原串的出现次数

性质：有相同前缀的后缀，在 \(sa\) 中是一段连续区间。
性质：\(S[l,r]\) 等价于 \(S[l,n]\) 这个后缀的一段长为 \(r-l+1\) 的前缀。

考虑在 \(sa\) 中找到 \(S[l,n]\) 这段后缀（\(rnk_l\)），在左/右两边分别二分出最小/大的 \(x,y\)，使 \(lcp(x,rnk_l)\geq r-l+1\)，\(lcp(rnk_l,y)\geq r-l+1\)（这样 \([x,y]\) 区间都有了这一段前缀），与 0x01 的 LCP 求法类似，复杂度 \(O(n\log n)-O(\log n)\)。

本质不同的子串个数（位置不同算相同）

Recall：\(hei_i\)，为 \(lcp(i-1,i)\)，即排序后第 \(i\) 小后缀和第 \(i-1\) 小后缀的最长公共前缀。

考虑又一次把子串看成后缀的前缀，在 \(sa\) 意义下，后缀 \(sa_i\) 的长度为 \(n-sa_i+1\)，这些都是可以统计的子串。但是，有一些子串是重复算的，有 \(hei_i\) 个前缀一定在上一个 \(sa_{i-1}\) 计算过了，于是我们要求的其实是：

\[\sum\limits_{1\leq i\leq n}(n-sa_i+1-hei_i). \]

复杂度 \(O(n\log n)\)。

循环同构的最小表示法

最小表示法：字符串 \(S\) 的所有循环同构中，字典序最小的。

考虑倍长原串，这样 \(S\) 的循环同构就能表示成 \(S[i,i+n-1] (1\leq i\leq n)\)，对它做后缀排序，取 \(sa\) 中最小的在 \([1,n]\) 范围内的后缀（多那么一点不影响答案），这就是 \(S\) 的最小表示法。复杂度 \(O(n\log n)\)。

\(S,T\) 的最长公共子串

考虑造一个串 \(R=S+\text{'\#'}+T\)，求出它的后缀数组，于是有一个 naive 的想法：求出任意两个后缀的 \(LCP\)。但很显然是假的。

众所周知，前缀 \(\min\) 有单调性，于是我们可以贪心的选两个（在 \(sa\) 中）相邻的后缀，它们的 \(LCP\) 一定最优。但是！如果你取两个都来自 \(S\) 的前缀，那就寄了，想象一个把 \(S\) 的后缀染成白色，\(T\) 的后缀染成黑色，我们应该统计的是，在 \(sa\) 中，两个相邻的异色的后缀，的 \(LCP\)（即 \(hei\)）。复杂度 \(O(n\log n)\)。

判定后缀数组

有字符串 \(s\) 和后缀数组排列 \(p\)，令 \(p_0=n+1\)，\(s_{n+1}\) 小于 \(s\) 中其它字符。则 \(p\) 为 \(s\) 的后缀数组 \(\iff\) 对于所有 \(i\)，在 \(p\) 上：

  • 若 \(p_i+1\) 在 \(p_{i+1}+1\) 之前，则要求 \(s_i\leq s_{i+1}\)。
  • 若 \(p_i+1\) 在 \(p_{i+1}+1\) 之后，则要求 \(s_i< s_{i+1}\)。
posted @ 2022-11-06 19:17  caijianhong  阅读(48)  评论(0)  编辑  收藏  举报