后缀数组(倍增)算法原理及应用
概述
后缀数组是指对某一字符串的所有后缀按照字典序排序后的结果,这一信息常用于处理各种字符串问题,维护该信息的算法之一:倍增法的复杂度为O(n \log n)。尽管有更优秀的DC3法,但是倍增法更加通用好写,因此这里只介绍倍增法求后缀数组。
原理与实现
最初的思路
我们知道,如果处理出所有后缀,再使用快速排序等排序方式,复杂度是O(n^2 \log n)的。
我们希望优化复杂度,可以从排序算法上下手。比快排更优秀的排序方法也是有的,比如计数排序。计数排序需要消耗值域规模的空间,如果用在此处空间应该是O(n)的,不成问题。并且,计数排序并非基于比较的排序算法,之前复杂度中计入的O(n)比较字符串字典序现在也不需要了。那么,我们单次排序的复杂度降为O(n)了。
可是,计数排序需要有一个数值来代表每个字符串,基于数值来比较大小。对于代表字符串的数值,我们首先可以想到Hash。不过字符串长了Hash要取模无法比较大小不说,光是计算一遍Hash都需要O(n^2)时间,白白浪费了精力对排序进行优化。
倍增?
接下来是一个很神奇的操作,即倍增法的精髓所在:倍增。
我们可以考虑倍增关键字的长度,采用双关键字排序的方法。比如,当长度为1的时候,两个关键字就是每个后缀的第一个和第二个字符。每次对长度倍增,最多当长度变为 \displaystyle \lfloor \frac{n}{2} \rfloor的时候,所有排序就完成了。这样来看,复杂度就是O(n \log n)了。
接下来,我们要考虑具体的实现问题,比如,如何完成双关键字的计数排序。这一段,不放我们来直接阅读分析实现代码。
代码分析
建议在阅读下面的内容的时候配合后缀数组算法总结 | 远航休息栈一文中的若干示例,动手模拟一下,加深理解。
首先我们给出整个程序的代码。这里用struct简单封装了一下。
struct SuffixArray {
char s[MAXN];
int n, m, sa[MAXN], t[MAXN], t1[MAXN], c[MAXN];
inline void init() {
memset(this, 0, sizeof(SuffixArray));
}
inline void calsa() {
int *x = t, *y = t1;
for(int i = 1; i <= m; i++) c[i] = 0;
for(int i = 1; i <= n; i++) c[x[i] = s[i]]++;
for(int i = 1; i <= m; i++) c[i] += c[i - 1];
for(int i = n; i >= 1; i--) sa[c[x[i]]--] = i;
for(int k = 1; k <= n; k <<= 1) {
int p = 0;
for(int i = n - k + 1; i <= n; i++) y[++p] = i;
for(int i = 1; i <= n; i++) if(sa[i] > k) y[++p] = sa[i] - k;
for(int i = 1; i <= m; i++) c[i] = 0;
for(int i = 1; i <= n; i++) c[x[y[i]]]++;
for(int i = 1; i <= m; i++) c[i] += c[i - 1];
for(int i = n; i >= 1; i--) sa[c[x[y[i]]]--] = y[i];
std::swap(x, y);
p = x[sa[1]] = 1;
for(int i = 2; i <= n; i++) {
x[sa[i]] = (y[sa[i - 1]] == y[sa[i]] && y[sa[i - 1] + k] == y[sa[i] + k])
? p : ++p;
}
if(p >= n) break;
m = p;
}
}
};
// Code by KSkun, 2018/5
变量声明
n
代表字符串长度,m
代表字符集大小。s
是字符串本体,内容从下标1开始。
sa
在排序中是一个排名数组,下标为排名,值为该排名对应的后缀首字母位置(后称“后缀编号”)。严格意义上说,这个数组实际上保存了每一次排序的结果,但如遇本次排序仍然存在两个后缀排序关键字相同,则长度大的(编号小的)后缀在前。SA计算完成后,该数组就是SA本体。
t
与t1
是临时数组,用于存放后面的x
和y
的信息。
c
是计数排序的计数数组,下标为计数的值,值为该值出现的次数。排序过程中,我们会对其求前缀和以便计算排名。
x
是一个中间量数组,意义为得到第一关键字的大小,对于一次排序,下标为代表后缀编号,值为象征对应后缀编号第一关键字大小的值(事实上可以视作排名)。
y
是一个中间量数组,意义为第二关键字的排名,对于一次排序,下标为排名,值为第二关键字对应的后缀的编号。
初始化
int *x = t, *y = t1;
for(int i = 1; i <= m; i++) c[i] = 0;
for(int i = 1; i <= n; i++) c[x[i] = s[i]]++;
for(int i = 1; i <= m; i++) c[i] += c[i - 1];
for(int i = n; i >= 1; i--) sa[c[x[i]]--] = i;
第1/8行分配t
与t1
两个内存空间给x
数组和y
数组。
第2/9行清空用于计数排序的计数数组c
。
第3/10行,统计该次第一关键字(首字母)的出现次数。
第4/11行,计算前缀和。
第5/12行,计算第一次排序的sa
。此次排序,仅区分首字母和后缀编号,首字母小的在前,首字母相同,后缀编号小(后缀长度大)的在前。
倍增:计算y
数组
int p = 0;
for(int i = n - k + 1; i <= n; i++) y[++p] = i;
for(int i = 1; i <= n; i++) if(sa[i] > k) y[++p] = sa[i] - k;
第2/15行,有一些后缀拆分后的第二关键字为空,这个循环将这些后缀的排名提到最前,以表示空字符的最高优先级。
第3/16行,对于满的第二关键字排序,利用某一后缀的后缀也是原串的后缀这一特点,由上一次的sa
计算而来。考虑那些排名靠前的后缀,如果往前补倍增长度位,那么这个后缀对应的长度前缀就变成了第二关键字。如我们有后缀aa
和abaa
,长度为2,计算abaa
的第二关键字的时候,实际上查到的是aa
的排名,并且往前补了2位。
倍增:计数排序计算x
数组
for(int i = 1; i <= m; i++) c[i] = 0;
for(int i = 1; i <= n; i++) c[x[y[i]]]++;
for(int i = 1; i <= m; i++) c[i] += c[i - 1];
for(int i = n; i >= 1; i--) sa[c[x[y[i]]]--] = y[i];
第2/18行,依然是在统计各个量的出现次数。y[i]
表示的是第二关键字排名为i的后缀编号,而x[i]
拿到的则是第一关键字排序时,后缀编号为i的后缀对应的象征第一关键字大小的值(在初始化过程中,我们将该值设置为首字母的值了,因此,对初始化还有一种理解,即y[i] = i
,并没对第二关键字排序)。因此x[y[i]]
的值反映的是代表第一关键字大小的值。
第4/20行,对i
逆序统计,其用意是:当第一关键字相同的时候,第二关键字排名大的,计算后sa
的值也会更大。
倍增:计算下一次的x
数组
std::swap(x, y);
p = x[sa[1]] = 1;
for(int i = 2; i <= n; i++) {
x[sa[i]] = (y[sa[i - 1]] == y[sa[i]] && y[sa[i - 1] + k] == y[sa[i] + k]) ? p : ++p;
}
第1/21行,把旧的x数组值换到y数组上。
第2/22行,规定当前排名第一的后缀x
值为1,以后的后缀x
值不会小于1。
第4/24行,考虑什么时候一个后缀和前面的后缀下一次的第一关键字(=这一次的第一+第二关键字,因为是倍增)相同,只有这一次的第一第二关键字都相等才可以判定下一次的第一关键字相等,否则应该比上一次的大。
倍增:边界条件
if(p >= n) break;
m = p;
第1/27行,当本次排序中已经产生了不小于n
种不同的第一关键字的值,排序便没有必要继续下去了,因为已经区分开了这些后缀。
第2/28行,本次排序更新了关键字的种类数,因此要更新m
值。
例题:【P3809】【模板】后缀排序 – 洛谷
把上面的板子打上去就可以了。
// Code by KSkun, 2018/5
#include <cstdio>
#include <cctype>
#include <cstring>
#include <algorithm>
const int MAXN = 1000005;
struct SuffixArray {
char s[MAXN];
int n, m, sa[MAXN], t[MAXN], t1[MAXN], c[MAXN];
inline void init() {
memset(this, 0, sizeof(SuffixArray));
}
inline void calsa() {
int *x = t, *y = t1;
for(int i = 1; i <= m; i++) c[i] = 0;
for(int i = 1; i <= n; i++) c[x[i] = s[i]]++;
for(int i = 1; i <= m; i++) c[i] += c[i - 1];
for(int i = n; i >= 1; i--) sa[c[x[i]]--] = i;
for(int k = 1; k <= n; k <<= 1) {
int p = 0;
for(int i = n - k + 1; i <= n; i++) y[++p] = i;
for(int i = 1; i <= n; i++) if(sa[i] > k) y[++p] = sa[i] - k;
for(int i = 1; i <= m; i++) c[i] = 0;
for(int i = 1; i <= n; i++) c[x[y[i]]]++;
for(int i = 1; i <= m; i++) c[i] += c[i - 1];
for(int i = n; i >= 1; i--) sa[c[x[y[i]]]--] = y[i];
std::swap(x, y);
p = x[sa[1]] = 1;
for(int i = 2; i <= n; i++) {
x[sa[i]] = (y[sa[i - 1]] == y[sa[i]] && y[sa[i - 1] + k] == y[sa[i] + k])
? p : ++p;
}
if(p >= n) break;
m = p;
}
}
} sa;
int main() {
scanf("%s", sa.s + 1);
sa.n = strlen(sa.s + 1); sa.m = 256;
sa.calsa();
for(int i = 1; i <= sa.n; i++) {
printf("%d ", sa.sa[i]);
}
return 0;
}
应用
求height
值
上文提到过,x
数组的值实际上是后缀编号对应的排名,我们将其内容复制一份进入rnk
数组待用。
下面给出height
数组的定义:代表该后缀与前一排名后缀的最长公共前缀(LCP)长度,下标为后缀编号。以下将该数组简写为hei
数组。
我们只需要对每个后缀找到紧接着后面排名的后缀,来统计相同的部分即可。需要注意的是,当我们向后枚举排名的时候,实际上可以看做比较的后缀各删去了第一位,可以用上次的k-1来作为基础,然后再比较前k-1位之后有没有更多的相同的部分,这样可以优化时间。
以下代码就是求hei
数组的代码,紧接着求sa
的代码。
memcpy(rnk, x, sizeof(rnk));
int k = 0;
for(int i = 1; i <= n; i++) {
int j = sa[rnk[i] + 1];
if(k) k--;
if(!j) continue;
while(s[i + k] == s[j + k]) k++;
hei[rnk[i]] = k;
}
求前缀信息
把字符串逆序一下扔进去求个SA就好了。
求前缀的LCS或后缀的LCP
要求两个排名不相邻的后缀的LCP怎么办?我们可以把上面的hei
数组建一个ST表,查询排名区间的最小值即可。hei
数组本身的含义就是LCP长度,这个性质有一定的传递性。在这里不太能说的明白,倒是可以采用观察一个sa
的例子验证一下的做法。
至于前缀的LCS,其实你发现逆序了以后就变成一个后缀LCP问题了。
这个问题的一个经典应用题目是:[NOI2016]优秀的拆分
参考资料
- 后缀数组算法总结 | 远航休息栈
- 《算法竞赛入门经典 训练指南》,刘汝佳,陈锋著,清华大学出版社,9787302291077
- 【后缀数组之height数组】 – LLGemini – 博客园
- NOI2016优秀的拆分 后缀数组 – CSDN博客