再谈后缀数组：SA-IS

Arahc 11.0
  2022-03-23 08:00 2022-03-23 08:00 创建   2023-03-15 01:14:46 2023-03-15 01:14:46 更新    5.1k 字  22 分钟  

去年 8 月的时候教练亲自出马给我们讲了后缀数组，讲了用倍增法 O(nlog⁡n)\mathcal O(n\log n) 实现求 SA 的做法。这一次我们继续深入，用 SA-IS 算法做到 O(n)\mathcal O(n) 级别。

并实现常数和空间双重吊打后缀自动机，重振后缀数组荣光。

倍增法实现后缀数组博客连接：Link to my blog。

首先我们忽略以前那篇博客的一些人类智商天花板语句，比如：

O(n)\mathcal O(n) 版本的 SA 很复杂，码量大。

（事实上只是因为当时我只知道 DC3 可以 O(n)\mathcal O(n)。）

然后就可以正式开始了。

前排提示：你甚至可以不需要理解 SA 的倍增求法，只要你知道 sasa 数组是什么就能理解本博客，请放心阅读。

声明

  1. 若无特殊说明，字符串下标从 00 开始，大写字母表示字符串，小写字母表示字符，要处理的字符串为 SS。
  2. S[l,r]S[l,r] 表示 SS 在 [l,r][l,r] 区间内的子串，SiS_i 表示 SS 的下标为 ii 的项的字符。
  3. lcp(A,B)\mathrm{lcp}(A,B) 为 A,B 的最长公共前缀。
  4. preT(i)\mathrm{pre}_T(i) 表示字符串 TT 的以下标 ii 结尾的前缀，sufT(i)\mathrm{suf}_T(i) 表示字符串 TT 的以下标 ii 开头的后缀。为了方便，在不会引起歧义的地方可能会省略下标不写，记为 pre(i),suf(i)\mathrm{pre}(i),\mathrm{suf}(i)。
  5. 所有字符和字符串之间的大小比较均为字典序比较。

为了方便，我们在 SS 的末尾添加一个特殊字符（不妨记为 *\texttt{*}），并定义：*\texttt{*} 的字典序大于其它所有字符。

本篇博客大多参考自这个博客：Link to a blog。然而这个博客有少数地方有锅，甚至代码过不去 UOJ 的板子，但是思路非常清晰简明。

当然一定程度上还参考了一篇论文：Link to a thesis。虽然是英文的，但是初二英语水平就能看懂，不必担心。

后缀分类得 LMS

后缀类型

我们给 SS 的每一个后缀进行分类：

定义一个后缀 sufS(i)\mathrm{suf}_S(i) 是 S 类（S-type）的，当且仅当 suf(i)<suf(i+1)\mathrm{suf}(i)<\mathrm{suf}(i+1)。特别的，suf(∣S∣)=*\mathrm{suf}(\left|S\right|)=\texttt{*} 为 S 类。
定义一个后缀是 L 类（L-type）的，当且仅当其不是 S 类的。

根据定义，不难得到递推求后缀类型的方法：

后缀类型递推判定：

若 suf(i)\mathrm{suf}(i) 为 S 类（i<∣S∣i<\left|S\right|），当且仅当下面至少一项成立：

  1. Si<Si+1S_i<S_{i+1}。
  2. Si=Si+1S_i=S_{i+1} 且 suf(i+1)\mathrm{suf}(i+1) 为 S 类。

若 suf(i)\mathrm{suf}(i) 为 L 类（i<∣S∣i<\left|S\right|），当且仅当下面至少一项成立：

  1. Si>Si+1S_i>S_{i+1}。
  2. Si=Si+1S_i=S_{i+1} 且 suf(i+1)\mathrm{suf}(i+1) 为 L 类。

证明的话，若满足第一个条件显然成立；若满足第二个情况，相当于去掉这两个后缀的第一个字符，剩下的部分比较大小。

根据本引理，还可以得出一个推论：若一个后缀 TT 满足 T0=T1T_0=T_1，去掉其第一个字符，后缀类型不变。

于是可以用 Θ(n)\Theta(n) 的时间内，判断 SS 的每一个后缀的类型。

推出来有什么用？

后缀类型性质：

对于两个后缀 A,BA,B，若 A0=B0A_0=B_0 且 A,BA,B 后缀类型不同，则 S 类的那个后缀比 L 类的那个后缀大。

证明的话，不妨设 AA 为 S 类。记 A=abX,B=acYA=\texttt{abX},B=\texttt{acY}：

  1. a≠b,a≠ca\neq b,a\neq c，由两个字符串类型得：a<b,a>ca<b,a>c。因此 c<a<bc<a<b，因此 A>BA>B。
  2. a=b,a≠ca=b,a\neq c，同理可以得到 b=a>cb=a>c。a≠b,a=ca\neq b,a=c 同理。
  3. a=b=ca=b=c，同时去掉 A,BA,B 的第一个字符，其类型不变。重复进行次操作直到其不是第三种情况即可。

LMS

然而只有这样一个性质显然不足以排序。因此我们在 SS 的所有 S 类后缀中挑选一些特殊的类型，记为 P 类，P 类字符串集是 S 类字符串集的一个子集。确切地说，一个后缀 suf(i)\mathrm{suf}(i) 是 P 类，当且仅当其同时满足以下三个条件：

  1. 它是 S 类。
  2. i>0i>0。
  3. suf(i−1)\mathrm{suf}(i-1) 是 L 类。

据个例子，现在的字符串是 S=baaacbaabS=\texttt{baaacbaab}：

ii 0 1 2 3 4 5 6 7 8 9
SiS_i b\texttt{b} a\texttt{a} a\texttt{a} a\texttt{a} c\texttt{c} b\texttt{b} a\texttt{a} a\texttt{a} b\texttt{b} *\texttt{*}
Type L P S S L L P S L P

可以理解为，P 类是一串 S 类中最靠左的一个，即 Left-Most S-type\text{Left-Most S-type}。以下标 ii 开头的后缀是 P 类，则称这个下标 ii 上的字符是一个 LMS 字符。相邻两个 LMS 字符之间的部分（包括这两个字符）为 LMS 子串。特别的，*\texttt{*} 单独构成一个 LMS 子串。

例如上面这个例子，S[1,6],S[6,9],S[9,9]S[1,6],S[6,9],S[9,9] 是 LMS 子串。

显然有这样几个性质，证明非常简单，不必赘述：

  1. *\texttt{*} 是最短的 LMS 子串。
  2. 任意非 *\texttt{*} 的 LMS 子串长度 >2>2。
  3. 字符串 SS 的所有 LMS 子串的长度之和 ⩽2∣S∣\leqslant2\left|S\right|。
  4. 一个字符串的 LMS 子串数量 ⩽⌈∣S∣2⌉\leqslant \left\lceil\frac{\left|S\right|}{2}\right\rceil（长度折半性质）。

此外有一个重要引理：

不存在两个长度大于 11 的 LMS 子串 A,BA,B，满足 AA 是 BB 的真前缀，且 AA 的每个字符对应的后缀类型组成的字符串也是 BB 的每个字符对应的后缀类型组成的字符串的真前缀。
注：后文对于两个 LMS 子串 A,BA,B 是否是真前缀关系的描述，也是同时比较字符串本身和后缀类型串两个的。

据个例子，串 S=acbbccbbccbac*S=\texttt{acbbccbbccbac*}，其每个位置对应的后缀类型为 SLPSLLPSLLLPLP\texttt{SLPSLLPSLLLPLP}，此时 LSM 子串 S[2,6]=bbccbS[2,6]=\texttt{bbccb} 是 S[6,11]=bbccbaS[6,11]=\texttt{bbccba} 的真前缀，但是其后缀类型 PSLLP\texttt{PSLLP} 并不是 PSLLLP\texttt{PSLLLP} 的真前缀。

证明：若 AA 是 BB 的真前缀，则 BB 串的后缀类型中至少有三个 P，显然不符合 LMS 定义。

由简单性质第三点，我们可以用线性排序做法做到 O(∣S∣)\mathcal O(\left|S\right|) 做到给原串的所有 LMS 子串排序。至于这里的线性排序做法，考虑到基数排序常数较大，我们将会用到诱导排序（induced sort）。

无论如何，假定我们已经给这些 LMS 子串排好序了。将 LMS 子串按字典序大小编号，按原串内的顺序排列，得到排名数组 rkrk。这里还是以 S=baaacbaabS=\texttt{baaacbaab} 为例子：

LSM 子串（下标表示） S[1,6]S[1,6] S[6,9]S[6,9] S[9,9]S[9,9]
LSM 串 aaacba\texttt{aaacba} aab*\texttt{aab*} *\texttt{*}
字典序排名 1 2 0

于是 rk=⟨1,2,0⟩rk=\left\langle 1,2,0 \right\rangle。

得出 rkrk 有什么用呢？我们需要另外一个引理：

问题缩减引理：

rkrk 数组中两个后缀的大小关系，就是 SS 中 P 类后缀的大小关系。

证明的话，考虑暴力比较两个字符串字典序大小的过程，比对 rkrk 中的两个数，相当于比对 SS 中的两个 LMS 子串。因为上面所说的重要引理，被比较的 LMS 子串要么完全相同，要么存在一个位置上字符不同，要么存在一个位置上后缀类型不同。对于后缀类型不同的，S 类后缀字典序更大，而因为后缀类型性质，其字典序确实会更大。因此原命题成立。

诱导排序求 SA

如果可以用 P 类后缀的顺序对 S,L 类后缀排序，就可以得到后缀数组 sasa。

而我们发现，求出 SS 的 P 类后缀顺序，要求出 rkrk 的后缀数组，而求出它则需要求出 rkrk 的 P 类后缀顺序……这是一个递归的东西，因此后文为了方便，我们把递归 nn 层时要求的这个字符串叫做 SnS_n。例如这里的 rkrk 数组实际上为 S1S_1。而对于后缀数组 sasa，为了与下标区分，将 SnS_n 对应的后缀数组的第 ii 项记为 san(i)sa_n(i)。

还是对于这个串 S=baacbaabS=\texttt{baacbaab}，其后缀数组（包括添加的 *\texttt{*}）为：

排名 后缀
0 *\texttt{*}
1 aab*\texttt{aab*}
2 aacbaab*\texttt{aacbaab*}
3 ab*\texttt{ab*}
4 acbaab*\texttt{acbaab*}
5 b*\texttt{b*}
6 baab*\texttt{baab*}
7 baacbaab*\texttt{baacbaab*}
8 cbaab*\texttt{cbaab*}

显然，首字母相同的后缀连续排布。因此对于每一个首字母开一个桶。先考虑把每个桶内的后缀按照字典序排好。可以得出同一桶内 S 类和 L 类是连续分布的，因为首字母相同时，S 类总是比 L 类大。因此每个桶又可以分成两部分，先排 L 类，再排 S 类。如果后缀类型和首字母都一样？就需要用到诱导排序了。

诱导排序需要一个关于 P 类后缀的有序的后缀数组 sa1sa_1，考虑如何从 sa1sa_1 推出当前的 sasa（sa0sa_0）。诱导排序的过程为：

  1. 将 sasa 数组初始化为 −1-1。
  2. 确定每个桶的 S 类桶的第一名的排名（起始位置），将 sa1sa_1 的每一个 P 类后缀按照 sa1sa_1 的顺序入桶。为了方便，S 类后缀的桶的放置是倒序的。
  3. 确定每个桶的 L 类桶的第一名的排名（起始位置），将 sasa 数组从左到右扫一遍，若 sa(i)>0sa(i)>0 且 suf(sa(i)−1)\mathrm{suf}(sa(i)-1) 是 L 类后缀，则将这个后缀入桶。
  4. 重新确定每个桶的 S 类桶的起始位置，因为所有 P 类后缀需要重新排序。注意 S 类桶是倒序放的，所以扫的顺序是从右往左，如果 sa(i)>0sa(i)>0 且 suf(sa(i)−1)\mathrm{suf}(sa(i)-1) 是 S 类后缀，则将这个后缀入桶。

注意，这里说的桶是 sasa 数组这一“大桶”，因为首字母和后缀分类将其划分成“小桶”，类似于基数排序，提前算好每个小桶的大小，入桶时放入对应的小桶内。小桶间不会相互影响。后文的诱导排序的桶都是这样的。

为什么这样排序是正确的？

我们已经知道，后缀排序之后，首字母相同的在同一段，同一段内是先有一段 L 类，剩下的部分是一段 S 类。因此先排 L 再排 S 显然没有问题。

因此我们重点要证明对于两个首字母相同的 L 类后缀 A,BA,B（记 A=suf(i),B=suf(j)A=\mathrm{suf}(i),B=\mathrm{suf}(j)）满足 A<BA<B 时，AA 一定比 BB 先加入到 sasa 中。

考虑正常比较 A,BA,B 时，先比较首字母，如果首字母相同则比较下一个字符，也就是比较 suf(i+1),suf(j+1)\mathrm{suf}(i+1),\mathrm{suf}(j+1) 的大小关系，也就是递归比较。

这里因为两个后缀的首字母相同，因此两个后缀的大小关系就是两个后缀入桶的时间关系，而两个后缀入桶的时间关系，依赖于 A+1,B+1A+1,B+1 的关系，符合前面说的递归比较。而对于 S 类后缀，同理即可。

至此证明完毕。

既然可以完成 sa1sa_1 到 sasa 的诱导，如何求出 sa1sa_1？问题是类似的，因此可以递归：

  1. 若 S1S_1 中每个字符都不同，则用桶排序即可。
  2. 否则，递归到 S2S_2 并计算 sa2sa_2，诱导至 sa1sa_1。

最后的问题就是对 LMS 子串之间的排序。如果这一步用的是基数排序我们就可以得到 SA-KA 算法（基本上是 SA-IS 和 KA Algorithm 唯一的区别），但是众所周知基数排序常数很大，所以我们还是要灵活运用诱导排序。这里原博客绕了一点弯路，我们走直一点。

对 LMS 子串的排序用的诱导排序流程，如果沿用上面对 L 类和 S 类的那个诱导排序，则要求我们传入一个有序的 LMS 的 sasa（或许可以记作 lms1lms_1？）进去。如果我们传入一个乱序的 sa1sa_1 进去，诱导排序的结果显然不正确。

不过有一个结论：传入乱序的 sa1sa_1 时，其返回的后缀数组 sasa 中，关于 LMS 后缀的子序列的信息，是按照每个 LMS 后缀的开头的 LMS 子串为关键字进行排序的。

听起来挺离谱，尝试证明：

首先，如果一开始传入的 LMS 数组，如果只看其开头的第一个字母，别的字母不看，则在后缀数组中是有序的。

那么经过一轮诱导（L 类后缀解决）后，对于每个 L 类后缀，如果只看这个后缀开头到第一个 S 类字符的位置，别的地方不看，它们是有序的。

再经过一轮诱导（S 类后缀解决），对于每个 S 类后缀，如果只看这个后缀开头到第一个 LMS 字符所在的位置，别的地方不看，它们是有序的。

如果这个后缀本身就是一个 P 类后缀，那么其有序的部分，就恰好会是这个后缀开头的 LMS 子串的部分。

因此一次诱导排序之后可以给 LMS 子串排好序。

复杂度分析

根据长度折半性质，有复杂度式：

T(n)=T(⌈n2⌉)+Θ(n)T(n) = T\left(\left\lceil\frac{n}{2}\right\rceil\right) + \Theta(n)

不难得到，这个算法的复杂度是 O(n)\mathcal O(n) 的。应该是效率最高的 SA 算法（据说 DC1024 比 SA-IS 快，但是 DC1024 在竞赛中出现概率必然零）。虽然也是线性，但是相比其他 SA 算法常数减小明显，一是因为诱导排序本身比基数排序快，二是因为没什么寻址不连续的地方。

时间复杂度是 O(n)\mathcal O(n) 的，空间复杂度也是 O(n)\mathcal O(n) 的，不需要乘上字符集大小。

对比常数同样不大的 KA 算法（表格来自论文），运行任务为求出文本的后缀数组（模板 后缀排序）：

运行时间单位：秒；空间单位：MB。

测试数据 字符集大小 串长 运行时间（IS/KA） 空间（IS/KA） 递归深度（IS/KA）
《世界概况》 94 2473400 1.3 / 1.9 12.70 / 21.24 6 / 6
《圣经》 63 4047392 2.7 / 3.51 20.86 / 34.45 6 / 6
人类第 22 号染色体 4 34553758 24.7 / 33.4 178.09 / 289.97 6 / 8
大肠杆菌基因组 4 4638690 2.8 / 3.98 24.29 / 40.04 7 / 8
蛋白质序列数据库 66 109617186 94.6 / 132.8 9 554.58 / 930.06 7 / 9
古登堡计划 146 105277340 101.0 / 149.67 542.17 / 907.31 11 / 12

注意：因为这篇论文时代久远（2009），评测机配置也不怎么样（论文中写道评测机配置为 2.20GHz CPU；而截至 2022 年 3 月，CCF 系列赛评测机配置为 3.70GHz。且论文里写道评测机运行内存为 2GB，我们机房的电脑没加内存条是 4GB），因此虽然开了 O3 优化，还是不够贴切现在的感觉，具体时间请按照实际测试为准。

代码

下面这份代码可以通过 UOJ 的模板题。

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
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
#include<bits/stdc++.h>
using namespace std;
const int max_n=100005,max_m=53;
inline int read(){
    int x=0;bool w=0;char c=getchar();
    while(c<'0' || c>'9') w|=c=='-',c=getchar();
    while(c>='0' && c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return w?-x:x;
}
inline void write(int x){
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10^48);
}
struct SA{    
    inline bool isLMS(bool *tp,int x){
        return (x && tp[x] && !tp[x-1]);
    }
    inline void idsort(int *S,int *sa,bool *tp,int *buk,int *lbk,int *sbk,int n,int m){
        for(register int i=0;i<=n;++i)
            if(sa[i]>0 && !tp[sa[i]-1])
                sa[lbk[S[sa[i]-1]]++]=sa[i]-1;
        for(register int i=1;i<=m;++i)
            sbk[i]=buk[i]-1;
        for(register int i=n;i>=0;--i)
            if(sa[i]>0 && tp[sa[i]-1])
                sa[sbk[S[sa[i]-1]]--]=sa[i]-1;
    }
    inline bool equal(int *S,bool *tp,int x,int y){
        do{
            if(S[x]!=S[y]) return 0;
            ++x,++y;
        }while(!isLMS(tp,x) && !isLMS(tp,y));
        return S[x]==S[y];
    }
    inline int* sais(int *S,int n,int m){
        bool *tp=new bool[n+5];
        int *pos=new  int[n+5];
        int *rak=new  int[n+5];
        int *sa= new  int[n+5];
        int *buk=new  int[m+5];
        int *lbk=new  int[m+5];
        int *sbk=new  int[m+5];
        for(register int i=0;i<=m+2;++i)
            buk[i]=0;
        for(register int i=0;i<=n+2;++i)
            sa[i]=rak[i]=-1;
        lbk[0]=sbk[0]=0;
        for(register int i=0;i<=n;++i)
            ++buk[S[i]];
        for(register int i=1;i<=m;++i){
            buk[i]+=buk[i-1],
            lbk[i]=buk[i-1],
            sbk[i]=buk[i]-1;
        }
        tp[n]=1;
        for(register int i=n-1;i>=0;--i){
            if(S[i]<S[i+1])
                tp[i]=1;
            else if(S[i]>S[i+1])
                tp[i]=0;
            else
                tp[i]=tp[i+1];
        }
        int cnt=0;
        for(register int i=1;i<=n;++i)
            if(tp[i] && !tp[i-1])
                pos[cnt++]=i;
        for(register int i=0;i<cnt;++i)
            sa[sbk[S[pos[i]]]--]=pos[i];
        idsort(S,sa,tp,buk,lbk,sbk,n,m);
        int las=-1,tot=1;
        bool rep=0;
        for(register int i=1;i<=n;++i){
            int x=sa[i];
            if(!isLMS(tp,x)) continue;
            if(las>=0 && !equal(S,tp,x,las))
                ++tot;
            if(las>=0 && tot==rak[las])
                rep=1;
            rak[x]=tot,las=x;
        }
        rak[n]=0;
        int ps=0;
        int *sa1,*s1=new int[cnt+5];
        for(register int i=0;i<=n;++i)
            if(rak[i]>=0)
                s1[ps++]=rak[i];
        if(!rep){
            sa1=new int[cnt+5];
            for(register int i=0;i<cnt;++i)
                sa1[s1[i]]=i;
        }
        else
            sa1=sais(s1,cnt-1,tot);
        lbk[0]=sbk[0]=0;
        for(register int i=1;i<=m;++i){
            lbk[i]=buk[i-1],
            sbk[i]=buk[i]-1;
        }
        for(register int i=0;i<=n+2;++i)
            sa[i]=-1;
        for(register int i=cnt-1;i>=0;--i)
            sa[sbk[S[pos[sa1[i]]]]--]=pos[sa1[i]];
        idsort(S,sa,tp,buk,lbk,sbk,n,m);

        delete(tp),delete(buk),delete(sbk),delete(lbk),delete(s1),delete(sa1);
        return sa;
    }
    int sa[max_n],rk[max_n],ht[max_n];
    inline void build(int *a,int n,int m){
        int *p=sais(a,n,m);
        for(register int i=0;i<=n;++i)
            sa[i]=p[i];
        for(register int i=0;i<=n;++i)
            rk[p[i]]=i;
        for(register int i=0,j=0;i<=n;++i) if(rk[i]){
            if(j) --j;
            while(a[i+j]==a[j+p[rk[i]-1]])
                ++j;
            ht[rk[i]]=j;
        }
    }
    inline int& operator [] (const int &i){
        return sa[i];
    }
}sa;

int n,a[max_n];

signed main(){
    ios::sync_with_stdio(false),cin.tie(0);
    string s;
    cin>>s;n=s.size();
    for(register int i=0;i<n;++i)
        a[i]=s[i]-'a'+2;
    a[n]=1;
    sa.build(a,n,31);
    for(register int i=1;i<=n;++i)
        write(sa[i]+1),putchar(' ');
    putchar('\n');
    for(register int i=2;i<=n;++i)
        write(sa.ht[i]),putchar(' ');
    return 0;
}
  • 标题: 再谈后缀数组：SA-IS
  • 作者: Arahc
  • 创建于 : 2022-03-23 08:00:00
  • 更新于 : 2023-03-15 01:14:46
  • 链接: https://arahc.github.io/2022/03/23/【笔记】再谈后缀数组：SA-IS/
  • 版权声明: 本文章采用 CC BY-NC-SA 4.0 进行许可。
评论
目录
再谈后缀数组：SA-IS
  1. 声明
  2. 后缀分类得 LMS
    1. 后缀类型
    2. LMS
  3. 诱导排序求 SA
  4. 复杂度分析
  5. 代码