字符串后缀排序

Posted 2021-03-24 yifusuyi

后缀排序

Task

Description

给定一个字符串，要求按字典序升序输出它的所有后缀子串的第一个字符所在位置。

Requirements & Limitations

字符集大小为常数，要求时间复杂度 (O(n log n))，其中 (n) 为字符串长度

Algorithm

这就是大(ren)名(lei)鼎(zhi)鼎(hui)的后缀排序了。这里记录倍增法。

定义后缀 (i) 为首字符所在原字符串位置为 (i) 的后缀子串。

定义 (sa_i) 为排名为 (i) 的后缀的首字符所在位置，(rnk_i) 为后缀 (i) 的排名。也即对 (sa_i) 做排序相当于对原字符串做排序。

注意到 (sa) 和 (rnk) 是可以在 (O(n)) 时间内互相推出的，因为 排名为 (后缀 (i) 的排名) 的后缀为 (i)，形式化的说，(sa_{rnk_i} = i)，

考虑进行多次排序，每次排序固定一个长度，对所有后缀子串的该长度的前缀排序(长度不足补 (0))。我们倍增这个长度。

首先，我们对所有后缀子串的第一位进行排序，当字符相同时，我们要求首字符所在位置小的在前面，例如对于字符串 (ababa) 的所有后缀子串

1 ababa
2 baba
3 aba
4 ba
5 a

进行一次排序的结果为

1 ababa
3 aba
5 a
2 baba
4 ba

即我们 (sa) 数组当前的值为 ({1,~3,~5,~2,~4})，然后计算 (rnk) 数组，但是需要注意的是，两字符串第一位相同时，他们的 (rnk) 也应该相同。例如，(rnk_1 = rnk_3 = 1)。类似的得到 (rnk) 数组当前的值为 ({1,~2,~1,~2,~1})。

然后进行下一次排序。

由于第一次排序的长度为 (1)，所以本次排序的长度为 (1 imes 2 = 2)

对于上次的排序结果，我们它们的前两位进行排序。但是注意到一个性质：我们已经排好了第一位，现在对第二位进行大小比较，我们发现后缀 (i) 的第 (2) 位正是后缀 ((i + 1)) 的第 (1) 位，因此第二位的字符的相对大小关系我们也是知道的，具体的，比较 ((sa_i,~sa_j))，若 (rnk_{sa_i} eq rnk_{sa_j})，则他们的大小关系为 (rnk_{sa_i}) 和 (rnk_{sa_j}) 的大小关系。否则则说明它们的第一位是相同的，那么比较第二位，则它们的大小关系是 (rnk_{sa_i + 1}) 与 (rnk_{sa_j + 1}) 的大小关系。

排序以后的结果为

5 a
1 ababa
3 aba
4 ba
2 baba

即 (sa) 数组当前的值为 ({5,~1,~3,~4,~2})，而 (rnk) 当前的值为 ({2，~3，~2，~3，~1})

然后进行下一次排序。

由于上次的排序长度为 (2)，所以本次排序的长度为 (2 imes 2 = 4)

类似于上次的排序，我们已经知道了前两位的相对大小关系，现在要求前四位的相对大小关系，而后缀 (i) 的后两位正是后缀 ((i + 2)) 的前两位，在上次排序已经被排好了。于是他们后两位字符的相对大小关系也是知道的。

类似的，得到排序的结果为

5 a
3 aba
1 ababa
4 ba
2 baba

即 (sa = {5,~3,~1,~4,~2})，(rnk = {3,~5,~2,~4,~1})

至此， (rnk) 数组已经是一个 (1 sim n) 的排列，这意味着所有的后缀已经都排好了序，输出即可。

形式化的，我们倍增所排序的长度，设当前需要排序的长度为 (2len)，已经对 (len) 长度排好了序，且 (sa) 和 (rnk) 存储的是长度为 (len) 时的排序信息，那么对于两个后缀 ((i, j))，他们的前半部分大小关系为 ((rnk_i, rnk_j))，而注意到他们后半部分即为后缀 ((i + len)) （或后缀 ((j + len))）的前半部分，因此他们后半部分的大小关系为 ((rnk_{i + len}, rnk_{j + len}))。根据字符串比较规则，如果前半部分大小不同，则大小关系为前半部分大小关系，否则为后半部分大小关系。也就相当于对 ((rnk_i, rnk_{i + len})) 这个二元组进行排序。一直到所倍增的长度大于 (n) 时，则所有后缀的大小已被排好。

Algorithm (1)

Solution

我会 std::sort！

我们倍增长度，对二元组 ((rnk_i, rnk_{i + len})) 排序时，直接使用 std::sort。注意到我们倍增了 (O(log n)) 次，每次需要进行一次 (O(n log n)) 的排序，因此总时间复杂度 (O(n log^2 n))。但是依然能跑过你谷的模板题。

这个做法非常易于理解，因为这只是翻译了上面推导过程，代码也好记好写。

Code

#include <cstdio>
#include <algorithm>

const int maxn = 2000005;

int n;
char S[maxn];
int rnk[maxn], sa[maxn];
std::pair<int, int> MU[maxn];

int ReadStr(char *p);
bool cmp(const int &a, const int &b);

int main() {
  freopen("1.in", "r", stdin);
  n = ReadStr(S);
  for (int i = 1; i <= n; ++i) {
    rnk[i] = 1;
    sa[i] = i;
    MU[i].second = S[i];
  }
  for (int len = 1; len <= n; len <<= 1) {
    std::sort(sa + 1, sa + 1 + n, cmp);
    for (int i = 1; i <= n; ++i) {
      if (MU[sa[i]] == MU[sa[i - 1]]) {
        rnk[sa[i]] = rnk[sa[i - 1]];
      } else  {
        rnk[sa[i]] = i;
      }
    }
    for (int i = 1; i <= n; ++i) {
      MU[i].first = rnk[i];
      MU[i].second = rnk[i + len];
    }
  }
  for (int i = 1; i <= n; ++i) {
    qw(sa[i], ' ', true);
  }
  putchar('
');
  return 0;
}

int ReadStr(char *p) {
  auto beg = p;
  do *(++p) = IPT::GetChar(); while (((*p >= 'a') && (*p <= 'z')) || ((*p >= '0') && (*p <= '9')) || ((*p >= 'A') && (*p <= 'z')));
  *p = 0;
  return p - beg - 1;
}

inline bool cmp(const int &a, const int &b) {
  if (MU[a] != MU[b]) {
    return MU[a] < MU[b];
  } else {
    return a < b;
  }
}

Algorithm (2)

考虑是对双关键字进行排序，因此我们可以先对第二关键字做桶排序，然后对第一关键字做桶排序，在排序时，记录每个桶所压入的序号序列，取出是按照序列的顺序取出每个下标，作为排好序的序列。由于桶排序是稳定的排序，这样在第一关键字相同的时候先取出的是第二关键字较小的下标。这样我们得到了一个 (O(n)) 的基数排序，但是由于 std::vector 常数太大，这样的写法甚至比 (O(n log^2 n)) 的做法还要慢一倍，在你谷的板板题上会T三个点。

这个做法只是用这个大常数基数排序替换上面的 std::sort，虽然很慢但是他的复杂度已经正确了。

Code

其余部分相同，只是将 std::sort 换成 RadixSort，排序部分代码如下

std::vector<int>bk[maxn];

void RadixSort(int *const beg, int *const ed) {
  for (auto it = beg; it != ed; ++it) {
    bk[MU[*it].second].push_back(*it);
  }
  auto p = beg;
  for (int i = 0; i <= n; ++i) {
    for (auto u : bk[i]) {
      *(p++) = u;
    }
    bk[i].clear();
  }
  for (auto it = beg; it != ed; ++it) {
    bk[MU[*it].first].push_back(*it);
  }
  p = beg;
  for (int i = 0; i <= n; ++i) {
    for (auto u : bk[i]) {
      *(p++) = u;
    }
    bk[i].clear();
  }
}

Algorithm (3)

考虑使用人类智慧优化我们的基数排序部分。

我们再定义一个数组 (tp)，(tp_i) 代表 第二关键字 为 (i) 时字符串的位置。

（注意我们下方的假设为已经对 (len) 做好了排序，现在要排序 (2len)）

例如，(len = 2,~rnk_3 = 2)，则有 (tp_2 = 1)，因为后缀 (1) 的第二关键字为 (rnk_{1 + 2 = 3} = 2)，因此 (tp_2 = 1)。

形式化的，当 (i + len leq n) 时，(tp_{rnk_{i + len}} = i)。

而当 (i + len > n) 时，我们发现这一部分的 (i) 对应的第二关键字都是 (0)（因为后半部分是空串），比不是零的都要小，考虑到要求字典序相同时序号小的在前，我们只需要令 (tp_{i + len - n} = i) 即可。

当然，赋值时我们考虑枚举 (sa) 的值而不是 (rnk) 的（因为初始 (rnk) 的赋值会是字符集，需要一些额外处理），于是有

(forall i leq len,~tp_i = n - len + 1)，这一部分是对应第二关键字为 (0) 的部分。

(forall sa_i > len,~~tp_{++pos} = sa_i - len)。其中 (pos) 为计数器，初值为 (len)。注意我们是从小到大枚举 (i) 而不是枚举 (sa_i)。

现在用另一种方法解释第二个式子：

我们从小到大枚举了 (sa_i) ，即我们按照排名从小到大枚举了每个字符串，那么他们作为第二关键字的值当然是从小到大的（因为第二关键也是排名），所以我们按顺序记录这些位置，而这些位置作为第二关键字，对应第一关键字的位置是 (sa_i - len)，因为当前需要排序的长度为 (2len)，所以第一关键字与第二关键字的位置相差 (len)。

上面两句求 (tp) 的方法，是人类智慧基数排序的关键前提条件。

现在考虑进行基数排序。

基数排序的代码如下

void RadixSort() {
  for (int i = 0; i <= m; ++i) tax[i] = 0;
  for (int i = 1; i <= n; ++i) ++tax[rnk[i]];
  for (int i = 1; i <= m; ++i) tax[i] += tax[i - 1];
  for (int i = n; i; --i) sa[ tax[rnk[tp[i]]]-- ] = tp[i];
}

(tax) 代表桶，首先将桶清零，然后枚举每个字符串，将字符串的排名放入桶内。

然后对桶做一个前缀和。我们发现我们现在放入桶内的是第一关键字，这意味着从这个桶 (tax_i) 里出来的字符串所在的排名定大于 (sum_{x = 1}^{i - 1} tax_x)，因为这部分的第一关键字小于 (i)，他们的排名一定小于第一关键字为 (i) 的字符串。同理，他们的排名一定不大于 (sum_{x = 1}^i tax_x)。

然后进入第四行，我们倒序枚举第二关键字的值，用 (tp) 来确定这个第二关键字是哪个字符串的，然后用 (rnk) 找到这个字符串对应哪个桶。这个桶即为 (tax_{rnk_{tp_i}})。然后我们注意到从这个桶里出来的元素下标上界为当前的前缀和值，而我们是倒序枚举的第二关键字，所以它的排名应该是桶里剩下字符串中第二关键字最大也即排名最大的，因此它的排名就是 (tax_{rnk_{tp_i}})。完成后，将这个桶的值自减 (1)，代表将它的上界减一。这么做下去就可以完成基数排序。

好人类智慧啊QAQ

Code

#include <ctime>
#include <cstdio>
#include <cstring>
#include <vector>
#include <algorithm>

int n, m;
char S[maxn];
int rnk[maxn], sa[maxn], tp[maxn], tax[maxn];

void RadixSort();
int ReadStr(char *p);
bool cmp(const int &a, const int &b);

int main() {
  freopen("1.in", "r", stdin);
  n = ReadStr(S);
  for (int i = 1; i <= n; ++i) {
    rnk[i] = S[i];
    tp[i] = i;
  }
  m = 1000;
  RadixSort();
  for (int len = 1, p = 0; p < n; len <<= 1, m = p) {
    p = 0;
    for (int i = 1; i <= len; ++i) tp[++p] = n - len + i;
    for (int i = 1; i <= n; ++i) if (sa[i] > len) tp[++p] = sa[i] - len;
    RadixSort();
    std::swap(tp, rnk);
    rnk[sa[1]] = p = 1;
    for (int i = 2; i <= n; ++i) {
      rnk[sa[i]] = ((tp[sa[i - 1]] == tp[sa[i]]) && (tp[sa[i - 1] + len] == tp[sa[i] + len])) ? p : ++p;
    }
  }
  for (int i = 1; i <= n; ++i) {
    qw(sa[i], ' ', true);
  }
  putchar('
');
  return 0;
}

int ReadStr(char *p) {
  auto beg = p;
  do *(++p) = IPT::GetChar(); while (((*p >= 'a') && (*p <= 'z')) || ((*p >= '0') && (*p <= '9')) || ((*p >= 'A') && (*p <= 'z')));
  *p = 0;
  return p - beg - 1;
}

void RadixSort() {
  for (int i = 0; i <= m; ++i) tax[i] = 0;
  for (int i = 1; i <= n; ++i) ++tax[rnk[i]];
  for (int i = 1; i <= m; ++i) tax[i] += tax[i - 1];
  for (int i = n; i; --i) sa[ tax[rnk[tp[i]]]-- ] = tp[i];
}

Appreciation

(Algorithm~~3) 的代码参考于 @自为风月马前卒 的 blog。