/*
Suffix array O(n lg^2 n)
LCP table O(n)
*/
#include <cstdio>
#include <algorithm>
#include <cstring>
using namespace std;
#define REP(i, n) for (int i = 0; i < (int)(n); ++i)
namespace SuffixArray
{
const int MAXN = 1 << 21;
char * S;
int N, gap;
int sa[MAXN], pos[MAXN], tmp[MAXN], lcp[MAXN];
bool sufCmp(int i, int j)
{
if (pos[i] != pos[j])
return pos[i] < pos[j];
i += gap;
j += gap;
return (i < N && j < N) ? pos[i] < pos[j] : i > j;
}
void buildSA()
{
N = strlen(S);
REP(i, N) sa[i] = i, pos[i] = S[i];
for (gap = 1;; gap *= 2)
{
sort(sa, sa + N, sufCmp);
REP(i, N - 1) tmp[i + 1] = tmp[i] + sufCmp(sa[i], sa[i + 1]);
REP(i, N) pos[sa[i]] = tmp[i];
if (tmp[N - 1] == N - 1) break;
}
}
void buildLCP()
{
for (int i = 0, k = 0; i < N; ++i) if (pos[i] != N - 1)
{
for (int j = sa[pos[i] + 1]; S[i + k] == S[j + k];)
++k;
lcp[pos[i]] = k;
if (k)--k;
}
}
} // end namespace SuffixArray

概述

这是一个用于后缀数组构造的O(n log n)算法(或者更确切地说，如果使用的不是::sort，而是2遍桶排序)。

它的工作原理是首先排序2-g^(*)，然后是4-g，然后是8-g，依此类推，原始字符串S，所以在第i次迭代中，我们排序2ⁱ-g。显然，这样的迭代不可能超过log₂(n)，诀窍是在第i步中对2ⁱ-g进行排序，通过确保两个2ⁱ-g的每次比较都在O(1)时间内完成(而不是O(2ⁱ)时间)来简化。

它是怎么做到的?在第一次迭代中它对2g(又名双元)进行排序，然后执行所谓的字典重命名。这意味着它创建了一个新的数组(长度为n)，用于存储每个双元数据的排名。

字典重命名示例:假设我们有一个排序好的一些字母的列表。然后，我们通过从左到右分配排名(即字典名称)，从第0位开始，每当遇到新的双字符更改时，就增加排名。因此，我们分配的排名如下:

ab : 0
ab : 0   [no change to previous]
ca : 1   [increment because different from previous]
cd : 2   [increment because different from previous]
cd : 2   [no change to previous]
ea : 3   [increment because different from previous]

这些等级被称为字典名称。

现在，在下一次迭代我们对4克进行排序。这涉及到很多不同4克之间的比较。我们如何比较两个4克?我们可以逐个字符地比较它们。每次比较最多有4个操作。但是，我们使用在前面步骤中生成的排位表，通过查找包含在它们中的两个双元的排名来比较它们。该秩表示前一个2-g排序的字典顺序，因此，如果对于任何给定的4-g，它的第一个2-g的秩高于另一个4-g的第一个2-g的秩，那么在字典顺序上，它必须在前两个字符的某个位置大于。因此，如果对于两个4-g，第一个2-g的秩是相同的，它们在前两个字符中必须是相同的。换句话说，在排名表中两次查找足以比较两个4-g的所有4个字符。

排序后，我们再次创建新的字典名称，这次是4-g。

在第三次迭代中，我们需要按8g排序。同样，在前一步的字典排序表中进行两次查找足以比较两个给定8-g的所有8个字符。

等等。每次迭代i有两个步骤:

1. 按2ⁱ-grams排序，使用前一次迭代中的字典名称，以便在每个

   中分2步(即O(1)时间)进行比较

2. 创建新的词典名称

重复此操作，直到所有2ⁱ-g都不相同。如果发生这种情况，我们就完了。我们怎么知道它们是否都是不同的呢?嗯，字典名称是一个递增的整数序列，从0开始。因此，如果在迭代中生成的最高字典名称与n-1相同，则每个2ⁱ-gram必须具有自己的独特字典名称。

---

实现现在让我们看看代码来确认这一切。使用的变量如下:sa[]是我们正在构建的后缀数组。pos[]是排名查询表(即它包含字典名称)，具体来说，pos[k]包含前一步k- m-gram的字典名称。tmp[]是一个辅助数组，用来帮助创建pos[]。

我将在代码行之间给出进一步的解释:

void buildSA()
{
N = strlen(S);
/* This is a loop that initializes sa[] and pos[].
For sa[] we assume the order the suffixes have
in the given string. For pos[] we set the lexicographic
rank of each 1-gram using the characters themselves.
That makes sense, right? */
REP(i, N) sa[i] = i, pos[i] = S[i];
/* Gap is the length of the m-gram in each step, divided by 2.
We start with 2-grams, so gap is 1 initially. It then increases
to 2, 4, 8 and so on. */
for (gap = 1;; gap *= 2)
{
/* We sort by (gap*2)-grams: */
sort(sa, sa + N, sufCmp);
/* We compute the lexicographic rank of each m-gram
that we have sorted above. Notice how the rank is computed
by comparing each n-gram at position i with its
neighbor at i+1. If they are identical, the comparison
yields 0, so the rank does not increase. Otherwise the
comparison yields 1, so the rank increases by 1. */
REP(i, N - 1) tmp[i + 1] = tmp[i] + sufCmp(sa[i], sa[i + 1]);
/* tmp contains the rank by position. Now we map this
into pos, so that in the next step we can look it
up per m-gram, rather than by position. */
REP(i, N) pos[sa[i]] = tmp[i];
/* If the largest lexicographic name generated is
n-1, we are finished, because this means all
m-grams must have been different. */
if (tmp[N - 1] == N - 1) break;
}
}

关于比较函数函数sufCmp用于按字典顺序比较两个(2*gap)克。在第一次迭代中，它比较双元，在第二次迭代中是4克，然后是8克，以此类推。这是由全局变量gap控制的。

sufCmp的简单实现是这样的:

bool sufCmp(int i, int j)
{
int pos_i = sa[i];
int pos_j = sa[j];
int end_i = pos_i + 2*gap;
int end_j = pos_j + 2*gap;
if (end_i > N)
end_i = N;
if (end_j > N)
end_j = N;
while (i < end_i && j < end_j)
{
if (S[pos_i] != S[pos_j])
return S[pos_i] < S[pos_j];
pos_i += 1;
pos_j += 1;
}
return (pos_i < N && pos_j < N) ? S[pos_i] < S[pos_j] : pos_i > pos_j;
}

这将比较第i个后缀pos_i:=sa[i]开头的(2*gap)-gram与第j个后缀pos_j:=sa[j]开头的(2*gap)-gram。然后逐个字符进行比较，即比较S[pos_i]和S[pos_j]，然后比较S[pos_i+1]和S[pos_j+1]，以此类推。只要字符相同，它就会继续。一旦它们不同，如果第i个后缀中的字符小于第j个后缀中的字符，则返回1，否则返回0。(注意，返回int的函数中的return a<b意味着如果条件为真返回1，如果条件为假返回0)

return语句中看起来很复杂的条件是处理(2*gap)-g之一位于字符串末尾的情况。在这种情况下，pos_i或pos_j将在比较所有(2*gap)字符之前到达N，即使在此之前的所有字符都是相同的。如果第i个后缀在末尾，它将返回1，如果第j个后缀在末尾，它将返回0。这是正确的，因为如果所有字符都相同，较短的在字典上较小。如果pos_i已经到达末尾，则第i个后缀必须小于第j个后缀。

显然，这种朴素的实现是0(间隙)，即它的复杂性在(2*间隙)-g的长度上是线性的。然而，代码中使用的函数使用字典名称将其减少到O(1)(具体来说，减少到最多两次比较):

bool sufCmp(int i, int j)
{
if (pos[i] != pos[j])
return pos[i] < pos[j];
i += gap;
j += gap;
return (i < N && j < N) ? pos[i] < pos[j] : i > j;
}

可以看到，我们不是查找单个字符S[i]和S[j]，而是检查第i和第j个后缀的字典顺序。在之前的迭代中，为gap-grams计算了字典排序。因此，如果pos[i] < pos[j]，那么第i个后缀sa[i]必须以字典顺序小于sa[j]开头的gap-gram开头。换句话说，简单地通过查找pos[i]和pos[j]并比较它们，我们已经比较了两个后缀的第一个间隙字符。

如果排名相同，我们继续比较pos[i+gap]和pos[j+gap]。这与比较(2*gap)-g的下一个gap字符相同，即后半部分。如果排名再次相同，则两个(2*gap)-g相同，因此返回0。否则，如果第i个后缀小于第j个后缀，则返回1，否则返回0。

---

例子下面的示例说明了该算法是如何操作的，并特别说明了字典名称在排序算法中的作用。

要排序的字符串是abcxabcd。它需要三次迭代才能为此生成后缀数组。在每次迭代中，我将显示S(字符串)、sa(后缀数组的当前状态)以及{'ab','ab','ca','cd','cd','ea'}0和pos，它们表示字典名称。

首先，初始化:

S   abcxabcd
sa  01234567
pos abcxabcd

请注意，最初表示单字组的字典顺序的字典名称是如何与字符(即单字组)本身完全相同的。

第一次迭代:

排序sa，使用双元图作为排序标准:

sa  04156273

前两个后缀是0和4，因为它们是双字母'ab'的位置。然后是1和5(双字母'bc'的位置)，然后是6(双字母'cd')，然后是2(双字母'cx')。然后是7(不完整的大写字母"d")，然后是3(大写字母"xa")。显然，位置与顺序相对应，仅基于字符双引号。

生成字典名称:

tmp 00112345

如前所述，字典名称被分配为递增的整数。前两个后缀(都以双字母'ab'开头)为0，后两个后缀(都以双字母'bc'开头)为1，然后是2、3、4、5(每个都是不同的双字母)。

最后，我们根据sa中的位置进行映射，得到pos:

sa  04156273
tmp 00112345
pos 01350124

生成pos的方法是这样的:从左到右遍历sa，并使用pos中的条目来定义索引。使用tmp中相应的条目来定义该索引的值。所以pos[0]:=0,pos[4]:=0,pos[1]:=1,pos[5]:=1,pos[6]:=2等等。索引来自sa，值来自tmp)

第二个迭代:

我们再次对sa排序，并再次查看pos中的双元组(每个双元组代表原始字符串的两个双元组的序列)。

sa  04516273

请注意，与sa之前的版本相比，15的位置发生了怎样的变化。以前是15，现在是51。这是因为在之前的迭代中，pos[1]和pos[5]的重图是相同的(都是bc)，但现在pos[5]的重图是12，而pos[1]的重图是13。因此位置5位于位置1之前。这是因为字典名称现在每个都表示原始字符串的双元组:pos[5]表示bc,pos[6]表示'cd'。因此，它们一起代表bcd，而pos[1]代表bc,pos[2]代表cx，所以它们一起代表bcx，从字典顺序上看，CC_76确实比bcd大。

同样，我们通过从左到右筛选sa的当前版本并比较pos中相应的双元图来生成字典名称:

tmp 00123456

前两个条目仍然相同(都是0)，因为pos中对应的双元图都是01。其余部分是一个严格递增的整数序列，因为pos中的所有其他双元都是唯一的。

我们像以前一样执行到新的pos的映射(从sa取索引，从tmp取值):

sa  04516273
tmp 00123456
pos 02460135

第三个迭代:

我们再次对sa排序，取pos的双元(和往常一样)，现在每个双元表示原始字符串的4个双元的序列。

sa  40516273

你会注意到，现在前两个条目已经交换了位置:04变成了40。这是因为pos[0]处的双字母是02，而pos[4]处的双字母是01，后者显然在字典顺序上较小。深层原因是这两个分别代表abcx和abcd。

生成字典名称产生:

tmp 01234567

它们都不一样，最高的是7，也就是n-1。我们做完了，因为排序现在是基于m个不同的图。即使我们继续，排序顺序也不会改变。

---

改进建议用于在每次迭代中对2ⁱ-gram进行排序的算法似乎是内置的sort(或std::sort)。这意味着它是一个比较排序，在最坏的情况下需要O(n log n)时间，在每次迭代中。因为在最坏的情况下有log n次迭代，这使得它是一个O(n (log n)²)时间的算法。然而，排序可以通过使用两次桶排序来执行，因为我们用于排序比较的键(即前一步的字典顺序名称)形成一个递增的整数序列。因此，这可以改进为一个实际的O(n log n)时间的后缀排序算法。

---

评论我相信这是Manber和Myers在1992年的论文中提出的后缀数组构造的原始算法(链接在Google Scholar;它应该是第一个点击，它可能有一个PDF链接)。这(同时，但独立于Gonnet和Baeza-Yates的一篇论文)是将后缀数组(当时也称为部分数组)作为一种有趣的数据结构引入进一步研究的原因。

后缀数组构造的现代算法是O(n)，因此以上不再是可用的最佳算法(至少在理论，最坏情况复杂度方面不是)。

---

脚注

<一口>(*)这里的2g是指原始字符串中两个连续的字符的序列。例如，当S=abcde是字符串时，则ab、bc、cd、de是S的2g。同理，abcd和bcde是4克。通常，m-gram(对于正整数m)是由m个连续字符组成的序列。1克也叫一克，2克叫两克，3克叫三克。有些人继续用四边形、五角形等等。

注意S开头和位置为i的后缀是S的一个(n-i) gram。此外，每个m-gram(对于任何m)都是S的一个后缀的前缀。因此，对m-gram进行排序(对于m尽可能大)可以是对后缀进行排序的第一步。