诱导排序与SA-IS算法
riteme.site

诱导排序与 SA-IS 算法

SA-IS 是一种在实际运用中相当快速的、线性时间构建字符串的后缀数组的算法。本文将通过对原论文1进行一些翻译和总结来介绍该算法。最后会提供一个 SA-IS 算法的 C++ 实现。

1. 记号

现在先规定一些记号,以方便之后的讨论。

我们通常使用大写字母 $S, A, B, \dots$ 来表示一个字符串,小写字母 $a, b, c, \dots$ 表示单个的字符。字符通常被认为是一个整数。字符在一定条件下数量是有限的,我们将所有会用到的字符放入字符集 $\Sigma$ 中,此时 $|\Sigma|$ 表示字符集的大小。字符串之间可以任意连接,如 $AB$ 表示字符串 $A$ 与字符串 $B$ 依次连接而成的新字符串。字符串与单个字符之间也是如此。特别的,我们将空串记为 $\epsilon$

对于一个字符串我们可以访问其中的任意字符,用 $S[x]$ 表示 $S$ 中下标为 $x$ 的字符。下标从 $0$ 开始。我们用 $S[a, b]$ 表示 $S$ 的一个子串,即下标为 $a$ 的字符开始一直到下标为 $b$ 的字符依次连接而成的字符串,其中必须满足 $0 \le a \le b \lt |S|$。定义一个字符串 $S$ 的前缀为 $\text{prefix}(S, i) = S[0, i]$,其后缀为 $\text{suffix}(S, i) = S[i, |S| - 1]$

两个字符串 $A$$B$ 之间存在比较关系。当 $A$$B$ 的长度相同,且其中对应下标的字符也相同时,则为 $A = B$。当从左至右出现第一对不相同的字符时,则以这两个字符的大小关系作为 $A$$B$ 间的大小关系,这称为字典序。为了简化一些操作,定义 $\#$ 是字典序最小的字符,并且将其默认作为每个字符串的最后一个字符2,即作为 $S[|S|]$
另外,定义函数 $\text{lcp}(A, B)$ 表示 $A$$B$ 的最长公共前缀的长度,即两者所有相同的前缀中,最长的一个的长度。由此我们可以得知字典序的一个性质:

定理 1.1 $A < B$ 当且仅当 $A[\text{lcp}(A, B)] < B[\text{lcp}(A, B)]$

证明 这个结论显然成立,因为 $\text{prefix}(A, \text{lcp}(A, B) - 1) = \text{prefix}(B, \text{lcp}(A, B) - 1)$

2. SA-IS 算法

SA-IS 算法是基于诱导排序这种思想。基本想法就是将问题的规模缩小,通过解决更小的问题,获取足够信息,就可以快速的解决原始问题。从这里也可以看出,这一过程需要递归处理子问题。

在介绍 SA-IS 算法前,我们先来看一下该算法的基本框架:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
function SA-IS(S):
    t = bool[]
    S1 = int[]
    P = int[]
    bucket = int[]
    扫描倒序字符串确定每一个后缀的类型 -> t
    扫描t数组确定所有的 LMS 子串 -> P
    对所有的 LMS 子串进行诱导排序
    对每一个 LMS 子串重新命名,生成新的串 S1

    if S1 中的每一个字符都不一样:
        直接计算 SA1
    else
        SA1 = SA-IS(S1)  # 递归计算 SA1

    利用 SA1 来进行诱导排序,计算 SA
    return SA

现在你不会明白这里面都写了些什么东西,因为这只是一个简单的流程。接下来将会介绍其中的每一步并证明其正确性。

2.1 后缀类型

在第 6 行中,确定后缀的类型可能会令人无法理解,这是因为我们还没有定义什么是后缀的类型......

对于每一个后缀 $\text{suffix}(S, i)$,当 $\text{suffix}(S, i) < \text{suffix}(S, i + 1)$ 时,是 $S$ 型后缀。当 $\text{suffix}(S, i) > \text{suffix}(S, i + 1)$ 时,是 $L$ 型后缀。对于特殊的后缀 $\text{suffix}(S, |S|) = \#$,它默认为 $S$ 型。
例如,对于字符串 mmiissiissiippii,每一后缀的类型为:

1
2
3
      0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
S[i]: m  m  i  i  s  s  i  i  s  s  i  i  p  p  i  i  #
t[i]: L  L  S  S  L  L  S  S  L  L  S  S  L  L  L  L  S

既然我们需要得知每一个后缀类型,就需要一个快速的算法来计算。在此之前,我们发现它们有如下的性质:

引理 2.1 (后缀类型递推性质) 对于任意的 $i \in [0, |S| - 1]$:
如果 $t[i] = \text{S-type}$,当且仅当下面任意一项成立:

  1. $S[i] < S[i + 1]$
  2. $S[i] = S[i + 1]$$t[i + 1] = \text{S-type}$

如果 $t[i] = \text{L-type}$,当且仅当下面任意一项成立:

  1. $S[i] > S[i + 1]$
  2. $S[i] = S[i + 1]$$t[i + 1] = \text{L-type}$

证明 这里证明 $S$ 型的,$L$ 型是类似的。对于第一种情况,显然是成立的。对于第二种情况,我们设 $\text{suffix}(S, i) = aA, \text{suffix}(S, i + 1) = aB$。由于第一个字符是相同的,因此我们需要比较 $A$$B$ 的大小。因为它们是连续的后缀,所以 $A = \text{suffix}(S, i + 1), B = \text{suffix}(S, i + 2)$。由于我们是从右往左推出 $t$,所以 $A$$B$ 的关系实际上可以由 $t[i + 1]$ 给出。故 $t[i] = t[i + 1]$

因此,我们可以在 $\Theta(|S|)$ 的时间内,推出整个 $t$ 数组。

关于后缀类型,我们还可得出另外一个比较重要的性质:

引理 2.2 (后缀类型指导排序) 对于两个后缀 $A$$B$,如果 $A[0] = B[0]$$A$$S$ 型,$B$$L$ 型,则 $A > B$

证明 设 $A = abX, B = acY$,这里假设 $a \neq b$$a \neq c$。因为 $A$$S$ 型,所以可知 $a < b$。同理,$B$$L$ 型,可知 $a > c$。故 $ c < a < b$,所以 $A > B$。如果 $a = b$$a \neq c$,那么 $b = a > c$$A > B$,对于 $a = c$$a \neq b$ 同理。如果 $a = b = c$,则我们可以将 $A$$B$ 的第一个字符去掉,用新的后缀来进行比较。根据引理 2.1,去掉第一个字符后的后缀类型不变。因此我们可以通过这样的操作直到变为第一种情况。

2.2 LMS 子串

然而光有后缀类型,还不足以进行排序。因此我们在后缀类型的 $S$ 型中挑出特别的一类,记为 $*$ 型。$*$ 型是 $S$ 型的一种,它的特殊之处在于它要求它的左边的后缀必须是 $L$ 型的。依然以 mmiissiissiippii 为例:

1
2
3
4
      0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
S[i]: m  m  i  i  s  s  i  i  s  s  i  i  p  p  i  i  #
t[i]: L  L  S  S  L  L  S  S  L  L  S  S  L  L  L  L  S
            *           *           *                 *

可以将其理解为一连串的 $S$ 型中最靠左的一个。LMS (LeftMost S-type) 也正是这个意思。同时我们注意到,后缀 $\#$ 始终是 $*$ 型的。
对于每一个 $*$ 型所对应上的字符,我们称为 LMS 字符。上面的示例中,下标为 $2, 6, 10, 16$ 都是 LMS 字符。
位置相邻的两个 LMS 字符中间(包括这两个字符)所构成的子串,称为 LMS 子串。对于 mmiissiissiippii,其 LMS 子串依次为 iissiiissiiippii##

通过观察,我们发现 LMS 子串具有以下的性质:
引理 2.3 # 是最短的 LMS 子串。

引理 2.4 对于任意的非 # 的 LMS 子串,其长度大于 $2$

证明 因为两个 LMS 字符中间必定有一个 $L$ 型的后缀。

引理 2.5 (原串折半) 一个字符串中 LMS 子串的数量不超过 $\left\lceil{|S| / 2}\right\rceil$

证明 根据引理 2.4 可知。

引理 2.6 一个字符串的所有 LMS 子串的长度之和为 $O(|S|)$

此外,对于 LMS 子串间的大小比较,除了对每个字符的字典序进行比较外,还要对比每个字符所对应的后缀类型。$S$ 型的有更高的优先权,因为根据引理 2.2 可得,首字符相同的情况下,$S$ 型的后缀字典序更大。只有当每一个字符与其后缀类型都相同时,这两个 LMS 子串才被称为是相同的。
之所以这样定义,是因为在之后我们会利用 LMS 子串来进行对后缀的排序。如果缺少后缀类型的信息,就不能够胜任。
根据引理 2.6,我们可以利用基数排序$O(|S|)$ 的时间内对所有的 LMS 子串排序3。对 LMS 子串排序完后,我们按照字典序依次重新命名4,注意,如果两个 LMS 子串相同,则使用同样的名称5。这样给每个 LMS 子串命名后,按照其子串原有的顺序排出一个新串 $S_1$。继续以 mmiissiissiippii 为例:

1
2
3
4
5
6
      0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
S[i]: m  m  i  i  s  s  i  i  s  s  i  i  p  p  i  i  #
t[i]: L  L  S  S  L  L  S  S  L  L  S  S  L  L  L  L  S
            *           *           *                 *
新名称      2           2           1                 0
S1:   2 2 1 0

这样有什么用呢?我们发现,这实际上是将所有 $*$ 型的后缀进行了缩减,从而减小了问题的规模。对于这一点,我们有如下的引理:

引理 2.7 (问题缩减) $S_1$ 中两个后缀的字典序关系,就是 $S$ 中对应的 $*$ 型后缀的字典序关系。

证明 我们可以将 $S_1$ 视为是将 $*$ 后缀中不重合的部分进行切割并缩减。这样每一个 LMS 子串就可作为一个整体来进行比较。从而保持了这两者的一致性。

需要注意的是这里只是 $*$ 型后缀的字典序关系,与其它后缀无关。

2.3 从 SA1 诱导至 SA

从上面的引理 2.7 我们得知,只要获得了 $S_1$ 的后缀数组 $SA_1$,就可以得到所有 $*$ 型后缀的相对顺序。如果我们可以利用 $*$ 型后缀的相对顺序来对其它的 $L$ 型和 $S$ 型后缀6进行排序,就可以完成后缀数组的计算。

在这里我们先假定 $SA_1$ 已经计算出来,只需考虑如何计算 $SA$。在这之前,我们先观察一下后缀数组的形式。以 aabaaaab 为例,它的后缀数组是这样的:

1
2
3
4
5
6
7
8
9
#
aaaab#
aaab#
aab#
aabaaaab#
ab#
abaaaab#
b#
baaaab#

不难发现,首字母相同的后缀是连续排布的,这一点可以用反证法来证明。因此我们可以利用桶排序的思想,为每一个出现过的字符建立一个桶,用 $SA$ 数组来存储这些桶,每个桶之间按照字典序排列,这样就可以使后缀数组初步有序。
我们对每个后缀都赋予了一个后缀类型,那么在首字母一样的情况下,$S$ 型或 $L$ 型会连续分布吗?答案是肯定的。因为根据引理 2.2,首字母相同的后缀如果后缀类型不同,则相对顺序是确定的。因此易知不会出现 $S$ 型和 $L$ 型交替出现的情况。更进一步,由于 $L$ 型后缀更小,因此总是先排布 $L$ 型后缀,再排布 $S$ 型后缀。因此每一个字符的桶可以分为两部分,一个用于放置 $L$ 型后缀,另一个则用于 $S$ 型后缀。为了方便确定每一个桶的起始位置,$S$ 型后缀的桶的放置是倒序的。
但是如果首字母和后缀类型都一致,我们不能直接快速地判断大小关系。在这里就要利用到诱导排序了。

诱导排序的过程分为以下几步:

  1. $SA$ 数组初始化为每个元素都为 $-1$ 的数组。
  2. 确定每个桶 $S$ 型桶的起始位置。将 $SA_1$ 中的每一个 $*$ 型后缀按照 $SA_1$ 中的顺序放入相应的桶内。
  3. 确定每个桶 $L$ 型桶的起始位置。在 $SA$ 数组中从左往右扫一遍。如果 $SA[i] > 0$$t[SA[i] - 1] = \text{L-type}$,则将 $SA[i] - 1$ 所代表的后缀放入对应的桶中。
  4. 重新确定每个桶 $S$ 型桶的起始位置,因为所有的 $*$ 型后缀要重新被排序。由于 $S$ 型桶是逆序排放的,所以这次从右至左地扫描一遍 $SA$。如果 $SA[i] > 0$$t[SA[i] - 1] = \text{S-type}$,则将 $SA[i] - 1$ 所代表的后缀放入对应的桶中。

这样我们就可以完成从 $SA_1$ 诱导到 $SA$ 的排序工作。这里简单说明一下为什么这样做是正确的:首先对于所有的 $*$ 型后缀,都是有序排放的。从左至右扫描 $SA$ 数组实际上就是按照字典序扫描现有的已排序的后缀。对于两个相邻的 $L$ 型后缀 $A$$B$,这里假设 $|A| > |B|$,则必定有 $A > B$。由于 $B$ 会被先加入 $SA$ 中,所以我们保证了 $A$$B$ 之间的有序性。又因为 $L$ 型桶是从左往右的顺序加入的,所以所有的 $L$ 型后缀会逐步地按顺序加入到 $SA$ 中。最后所有的 $L$ 型后缀将会有序。
对于$S$型后缀,除了要注意是相反的顺序和需要重新对$*$型后缀排序外,其余的原理与$L$型的排序类似。

之前的讨论都是基于我们已知$SA_1$的情况下进行的。现在我们来考虑如何计算$SA_1$。由于$S_1$也是一个字符串,计算其后缀数组时可以考虑两种情况:

  1. 如果$S_1$中每一个字符都不一样,则可以直接利用桶排序直接计算$SA_1$
  2. 否则,递归计算$SA_1$。就如之前的算法框架所展示的一样。

2.4 对LMS子串排序

到这里,SA-IS算法几乎已经结束了,只是还有一个问题需要解决,就是对LMS子串的排序。
之前我们所提及的,我们可以利用基数排序。虽然可以在$O(|S|)$的时间内完成,但是事实上,这个基数排序不但常熟大,而且十分复杂(请想象一下对字符串进行基数排序......)。这个排序直接成为了整个算法的性能瓶颈。因此我们急切的需要一种新的算法来胜任这一任务。
这个算法依然是诱导排序

与之前从$SA_1$诱导到$SA$的算法一样,只是我们这里将第二步改为:

确定每个桶$S$型桶的起始位置。将每一个LMS子串的首字母按照任意顺序放入对应的桶中。

待算法完成,我们会获得一个$SA$数组,其中LMS子串之间是排好了序的。

为什么这个算法是正确的,我们需要扯到一个新的概念:LMS前缀。
规定LMS前缀函数$\text{pre}(S, i)$表示 (1) 如果$\text{suffix}(S, i)$$*$型的,则$\text{pre}(S, i) = S[i]$,即LMS字符的LMS前缀就是自己。 (2) 否则是从$S[i]$开始,到下一个LMS字符之间(包括首尾)的子串。同样的,按照$\text{suffix}(S, i)$的后缀类型的不同,LMS前缀也同样分为$S$型和$L$型。
例如,在mmiissiissiippii中,$\text{pre}(S, 2) = S[2] =$ i,而$\text{pre}(S, 3) = S[3, 6] =$ issi
因此上面的算法实际上是在对每个LMS前缀进行排序。

接下来将证明2-4步都是正确的,即每一步完成时$SA$数组中,LMS前缀都是有序的。

  1. 对于第二步,由于放入的LMS前缀都只有一个字符,因为桶的排列是按照字典序的,所以保证放置后一定有序。
  2. 对于第三步,当放入第一个$L$型LMS前缀时,$SA$数组必定是有序的(根据引理2.2)。假设我们已经放置了$k$$L$型LMS前缀,且它们在$SA$数组中保持有序,现在考虑放入的第$k + 1$个LMS前缀是否会保证有序。我们设这个LMS前缀为$\text{pre}(S, i)$,因为首字母不同的LMS前缀一定是保持有序的,因此我们只需要考虑它与其首字母相同的LMS前缀之间的关系。因为我们是从左至右扫描来使$L$型LMS前缀从小至大放置,那么对于所有之前放置的且首字母与其相同的LMS前缀$\text{pre}(S, j)$,应该都有$\text{pre}(S, j) < \text{pre}(S, i)$。假设我们存在一个这样的LMS前缀,使得$\text{pre}(S, j) > \text{pre}(S, i)$,由于$\text{pre}(S, j)[0] = \text{pre}(S, i)[0]$,所以我们得知$\text{pre}(S, j + 1) > \text{pre}(S, i + 1)$。而$\text{pre}(S, j + 1)$$\text{pre}(S, i + 1)$都是之前所加入过的,因此它们之间应当保持有序。而$\text{pre}(S, j) > \text{pre}(S, i)$告诉我们之前的$SA$数组不是有序的,与假设相反,故不存在这样的$\text{pre}(S, j)$。因此,当$\text{pre}(S, i)$放置后,$SA$数组保持有序。直到所有的$L$型LMS前缀加入完毕。
  3. 对于第四步的正确性的证明,与第三步的证明是类似的。读者可以自行推理一下。

这样,我们就完成了对这个诱导排序的正确性的证明。但我们需要注意的是,这只是对LMS前缀的排序,我们期望是对LMS子串进行排序。我们注意到每一个LMS子串都可以视为是一个LMS字符与一个LMS前缀连接而成的,因此我们只需要判断两个LMS字符和其后一位对应的LMS前缀的大小关系,就可以判断LMS子串之间的大小关系。
然而,我们其实可以做得更简单。下面的引理说明了上面的算法结束后的$SA$数组就是排好序的LMS子串的数组。

引理 2.8 (LMS子串排序) $SA$数组中每一个LMS字符所对应的LMS前缀已经按照其对应的LMS子串的顺序排好。

证明 对于首字母不同的LMS子串,这个结论是显然的。
对于首字母相同的LMS子串,它必定是由一个$L$型LMS前缀从右至左的推导过来的。由于不同的LMS子串必定对应着不同的递推的顺序,并且在第二步中,所有的$L$型LMS前缀已经排好了序。所以结论成立。

根据这个引理,我们可以直接对$SA$数组扫描一遍,就可以得到LMS子串的字典序,同时对它们进行命名。

3. 时空复杂度分析

在之前的讨论中,我们已经成功的运用诱导排序使每一步都是$\Theta(|S|)$。但是由于有一个递归的过程,时间复杂度似乎并不一定是线性的。
我们注意到,每次递归都是计算$S_1$的后缀数组。如果我们能够知道$S_1$的规模,就能够计算SA-IS的事件复杂度。
根据引理2.5,我们得知$|S_1| \le \left\lceil|S| / 2\right\rceil$。因此每一层的递归的问题规模都会减半。因此我们可以用以下的递归式来表示时间复杂度:
$$ T(n) = T(\left\lceil n/2\right\rceil) + \Theta(n) \tag{3.1} $$

求解可得:
$$ T(n) = \Theta(n) \tag{3.2} $$

因此总时间复杂度是$\Theta(n)$的。
对于这个递归式,我们可以理解为是递归了$O(\log n)$层,其中每一层的问题规模从小到大排序是$2^0, 2^1, 2^2, \dots, 2^{\left\lfloor\log n\right\rfloor}$。因此总复杂度就是对它们进行求和:
$$ \sum^{\left\lfloor\log n\right\rfloor}_{k=0} 2^k = 2 \times 2^{\left\lfloor\log n\right\rfloor} = \Theta(2^{\left\lfloor\log n\right\rfloor}) = \Theta(n) \tag{3.3} $$

对于空间复杂度的分析,与时间复杂度是如出一辙的。

4. 运行示例

下面将用aabaaaab作为输入字符串,展示计算其后缀数组的每一步。希望能通过这个示例能够更清晰地展现SA-IS算法的运作流程。由于涉及到对桶的操作,这里用@表示正在被处理的元素,而用^表示每个桶的起始位置。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
           0  1  2  3  4  5  6  7  8
S:         a  a  b  a  a  a  a  b  #
扫描后缀类型
t:         S  S  L  S  S  S  S  L  S
LMS characters:     *              *
         |# |       a         |  b  |  # 桶的名称
SA:       -1|-1 -1 -1 -1 -1 -1|-1 -1
对LMS子串进行排序
1. 放入LMS子串
SA:       08|-1 -1 -1 -1 -1 03|-1 -1
2. 从*型LMS前缀诱导到L型LMS前缀
SA:       08|-1 -1 -1 -1 -1 03|-1 -1
          @^  ^                 ^
SA:       08|-1 -1 -1 -1 -1 03|07 -1
           ^  ^             @      ^
SA:       08|-1 -1 -1 -1 -1 03|07 02  # pre(S, 6)不是L型的
           ^  ^                @   ^
SA:       08|-1 -1 -1 -1 -1 03|07 02
           ^  ^                   @^
SA:       08|-1 -1 -1 -1 -1 03|07 02
3. 从L型LMS前缀诱导到S型前缀
SA:       08|-1 -1 -1 -1 -1 03|07 02
           ^                 ^    @^
SA:       08|-1 -1 -1 -1 -1 01|07 02
           ^              ^    @   ^
SA:       08|-1 -1 -1 -1 06 01|07 02
           ^           ^    @      ^
SA:       08|-1 -1 -1 00 06 01|07 02
           ^        ^    @         ^
SA:       08|-1 -1 05 00 06 01|07 02  # 不存在pre(S, -1)
           ^     ^    @            ^
SA:       08|-1 -1 05 00 06 01|07 02
           ^     ^ @               ^
SA:       08|-1 04 05 00 06 01|07 02
           ^  ^ @                  ^
SA:       08|03 04 05 00 06 01|07 02
           ^ @^                    ^
SA:       08|03 04 05 00 06 01|07 02  # pre(S, 7)不是S型的
          @^  ^                    ^
SA:       08|03 04 05 00 06 01|07 02
扫描并重命名LMS子串
name:      1  2
S1:       2 1 0
由于每一个字符都不一样,直接计算SA1
SA1:     2 1 0
从SA1诱导到SA
         |$ |       a         |  b  |
SA:       -1|-1 -1 -1 -1 -1 -1|-1 -1
1. 按照SA1的原顺序放入(忽略S1最后的字符)
SA:       08|-1 -1 -1 -1 -1 03|-1 -1
           ^  ^                 ^
2. 从*型后缀诱导到L型后缀
SA:       08|-1 -1 -1 -1 -1 03|-1 -1
          @^  ^                 ^
SA:       08|-1 -1 -1 -1 -1 03|07 -1
           ^  ^             @      ^
SA:       08|-1 -1 -1 -1 -1 03|07 02
           ^  ^                @   ^
SA:       08|-1 -1 -1 -1 -1 03|07 02
           ^  ^                   @^
3. 从L型后缀诱导到S型后缀
SA:       08|-1 -1 -1 -1 -1 03|07 02
           ^                 ^    @^
SA:       08|-1 -1 -1 -1 -1 01|07 02
           ^              ^    @   ^
SA:       08|-1 -1 -1 -1 06 01|07 02
           ^           ^    @      ^
SA:       08|-1 -1 -1 00 06 01|07 02
           ^        ^    @         ^
SA:       08|-1 -1 05 00 06 01|07 02
           ^     ^    @            ^
SA:       08|-1 -1 05 00 06 01|07 02
           ^     ^ @               ^
SA:       08|-1 04 05 00 06 01|07 02
           ^  ^ @                  ^
SA:       08|03 04 05 00 06 01|07 02  # pre(S, 2)是L型
           ^ @^                    ^
SA:       08|03 04 05 00 06 01|07 02
          @^  ^                    ^
后缀数组构造完毕
SA: 8 3 4 5 0 6 1 7 2

return SA

5. 具体实现及性能对比

在原论文中说道SA-IS算法可以在100行左右的代码完成。我试了一下,基本符合这个要求。下面将会给出一个C++实现。需要注意的是,为了方便,这份代码中并没有回收分配的内存。

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
// 后缀类型
#define L_TYPE 0
#define S_TYPE 1

// 判断一个字符是否为LMS字符
inline bool is_lms_char(int *type, int x) {
    return x > 0 && type[x] == S_TYPE && type[x - 1] == L_TYPE;
}

// 判断两个LMS子串是否相同
inline bool equal_substring(int *S, int x, int y, int *type) {
    do {
        if (S[x] != S[y])
            return false;
        x++, y++;
    } while (!is_lms_char(type, x) && !is_lms_char(type, y));

    return S[x] == S[y];
}

// 诱导排序(从*型诱导到L型、从L型诱导到S型)
// 调用之前应将*型按要求放入SA中
inline void induced_sort(int *S, int *SA, int *type, int *bucket, int *lbucket,
                         int *sbucket, int n, int SIGMA) {
    for (int i = 0; i <= n; i++)
        if (SA[i] > 0 && type[SA[i] - 1] == L_TYPE)
            SA[lbucket[S[SA[i] - 1]]++] = SA[i] - 1;
    for (int i = 1; i <= SIGMA; i++)  // Reset S-type bucket
        sbucket[i] = bucket[i] - 1;
    for (int i = n; i >= 0; i--)
        if (SA[i] > 0 && type[SA[i] - 1] == S_TYPE)
            SA[sbucket[S[SA[i] - 1]]--] = SA[i] - 1;
}

// SA-IS主体
// S是输入字符串,length是字符串的长度, SIGMA是字符集的大小
static int *SAIS(int *S, int length, int SIGMA) {
    int n = length - 1;
    int *type = new int[n + 1];  // 后缀类型
    int *position = new int[n + 1];  // 记录LMS子串的起始位置
    int *name = new int[n + 1];  // 记录每个LMS子串的新名称
    int *SA = new int[n + 1];  // SA数组
    int *bucket = new int[SIGMA];  // 每个字符的桶
    int *lbucket = new int[SIGMA];  // 每个字符的L型桶的起始位置
    int *sbucket = new int[SIGMA];  // 每个字符的S型桶的起始位置

    // 初始化每个桶
    memset(bucket, 0, sizeof(int) * (SIGMA + 1));
    for (int i = 0; i <= n; i++)
        bucket[S[i]]++;
    for (int i = 1; i <= SIGMA; i++) {
        bucket[i] += bucket[i - 1];
        lbucket[i] = bucket[i - 1];
        sbucket[i] = bucket[i] - 1;
    }

    // 确定后缀类型(利用引理2.1)
    type[n] = S_TYPE;
    for (int i = n - 1; i >= 0; i--) {
        if (S[i] < S[i + 1])
            type[i] = S_TYPE;
        else if (S[i] > S[i + 1])
            type[i] = L_TYPE;
        else
            type[i] = type[i + 1];
    }

    // 寻找每个LMS子串
    int cnt = 0;
    for (int i = 1; i <= n; i++)
        if (type[i] == S_TYPE && type[i - 1] == L_TYPE)
            position[cnt++] = i;

    // 对LMS子串进行排序
    fill(SA, SA + n + 1, -1);
    for (int i = 0; i < cnt; i++)
        SA[sbucket[S[position[i]]]--] = position[i];
    induced_sort(S, SA, type, bucket, lbucket, sbucket, n, SIGMA);

    // 为每个LMS子串命名
    fill(name, name + n + 1, -1);
    int lastx = -1, namecnt = 1;  // 上一次处理的LMS子串与名称的计数
    bool flag = false;  // 这里顺便记录是否有重复的字符
    for (int i = 1; i <= n; i++) {
        int x = SA[i];

        if (is_lms_char(type, x)) {
            if (lastx >= 0 && !equal_substring(S, x, lastx, type))
                namecnt++;
            // 因为只有相同的LMS子串才会有同样的名称
            if (lastx >= 0 && namecnt == name[lastx])
                flag = true;

            name[x] = namecnt;
            lastx = x;
        }
    }  // for
    name[n] = 0;

    // 生成S1
    int *S1 = new int[cnt];
    int pos = 0;
    for (int i = 0; i <= n; i++)
        if (name[i] >= 0)
            S1[pos++] = name[i];

    int *SA1;
    if (!flag) {
        // 直接计算SA1
        SA1 = new int[cnt + 1];

        for (int i = 0; i < cnt; i++)
            SA1[S1[i]] = i;
    } else
        SA1 = SAIS(S1, cnt, namecnt);  // 递归计算SA1

    // 从SA1诱导到SA
    lbucket[0] = sbucket[0] = 0;
    for (int i = 1; i <= SIGMA; i++) {
        lbucket[i] = bucket[i - 1];
        sbucket[i] = bucket[i] - 1;
    }
    fill(SA, SA + n + 1, -1);
    for (int i = cnt - 1; i >= 0; i--)  // 这里是逆序扫描SA1,因为SA中S型桶是倒序的
        SA[sbucket[S[position[SA1[i]]]]--] = position[SA1[i]];
    induced_sort(S, SA, type, bucket, lbucket, sbucket, n, SIGMA);

    // 后缀数组计算完毕
    return SA;
}

SA-IS的实现还可以参考一篇很好的博文: A walk through the SA-IS Suffix Array Construction Algorithm,这篇文章的作者使用Python一步一步地介绍了如何实现一个基础的SA-IS算法。

我们看到SA-IS作为一个后缀排序的算法,代码量也并不小,然而实际上它速度非常快,其中一个重要的原因就是SA-IS算法几乎所有的操作都是顺序访问的,这样就可以很好地提高缓存命中率7。原论文中实测SA-IS算法击败了KA算法,不愧为目前为止速度最快的后缀数组的构建算法。

对中学OI界,构造后缀数组通常是使用倍增法和DC3算法。倍增法速度比DC3算法慢是众所周知的了,下面我将用不同规模的随机字符串来比较DC3算法与SA-IS算法。

数据规模 DC3 SA-IS
$10^4$ $0.015\text{s}$ $0.011\text{s}$
$10^5$ $0.055\text{s}$ $0.036\text{s}$
$10^6$ $0.365\text{s}$ $0.270\text{s}$
$2\times 10^6$ $0.997\text{s}$ $0.593\text{s}$
$5\times 10^6$ $2.803\text{s}$ $2.085\text{s}$
$10^7$ $4.012\text{s}$ $3.326\text{s}$

编译命令: g++ SAIS.cpp/DC3.cpp -o sais/dc3 -O3
运行环境: Ubuntu 14.04 LTS x64 / CPU: 2.0GHz 奔腾某CPU / 4GB内存 / 未启动X window

从上面的结果可以看出,SA-IS算法速度上明显优于DC3算法,并且数据规模越大,两者的速度差距越明显。
总而言之,SA-IS算法是一个相当不错的后缀数组的构建算法。似乎SA-IS算法所利用的诱导排序的思想还可以解决其它的一些字符串的问题,如计算LCP数组8,我并没有对此做过多的了解了。


  1. Nong, Ge; Zhang, Sen; Chan, Wai Hong (2009): Linear Suffix Array Construction by Almost Pure Induced-Sorting 

  2. 这里的 $\#$ 并不是指“井号”,而是一个特殊记号。字符串中不会出现这个字符。 

  3. 后面我们将会有更好的算法来对 LMS 子串排序,时间复杂度一致,但在实际运用中速度快很多。 

  4. 注意之后会创建新字符串$ S_1$,其字符集 $\Sigma$ 是不同的。 

  5. 判断两个 LMS 子串是否相等可以暴力判断,基于引理 2.6,可以证明其总复杂度为 $O(|S|)$。 

  6. 注意 $*$ 型后缀也属于 $S$ 型后缀,因此会对它们进行重排,从而确保新加入的 $S$ 型后缀的顺序是对的。 

  7. 即降低cache-missing。 

  8. Johannes Fischer (2011): Inducing the LCP-Array