后缀数组

后缀数组转述

概念

要了解后缀数组,首先要知道字符串的后缀。

后缀

​ 对于字符串str,从str的某个位置pos到str的末尾的所有字符依次排列组成的新字符串。

即:

str.suffix(pos) = str.substr(pos,str.length-pos)

/其中substr是C++的库函数(头文件),用来求字符串的子串。用法为:string::substr(起始位置,子串长度)/

我们可以看出,对于一个特定的字符串,它的所有后缀都可以用起始位置pos唯一确定,这在下面会用到。

后缀数组(SuffixArray)

​ 由于生产生活以及提升气质的需要,我们需要用到后缀数组。

顾名思义: 1.它是一个数组 2.它与后缀有关

其实,它是存放字典序排名为i的后缀的起始位置(编号),

即“排第i位是是从哪开始的后缀”。

在实际运用中,我们还会用到它的逆运算,称为rank数组,str.rank[i]用来存放特定字符串的某个后缀按字典序在所有后缀子串中的排名,

即“从i开始的后缀排第几位”。

str.rank[i] = str.suffix(i).rankIn( All str.suffix)

互相转化的具体方法为:

For(i=0 To length-1)

suffixArray[rank[i]]=i

即:后缀i的排位所在的后缀就是i,且不会有人跟它抢……

例子

讲了这么多,我们先来举个例子。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
str=mycat

str.suffix={mycat,ycat,cat,at,t}

    id:        0       1     2    3    4

str.suffix(0)=mycat,str.suffix(1)=ycat,str.suffix(2)=cat”……

DictionaryOrderSort(str.suffix);

str.suffix={at ,cat,mycat,t,ycat}

    id:       3     2      0     4    1

str.suffixArray[]={3,2,0,4,1}

str.rank[]={2,4,1,0,3}

构造

​ 在知道了后缀数组与rank数组的概念之后,我们现在所需要的是知道如何构造。最简单的方法是将所有后缀取出,作为普通的字符串直接排序(暴力大法好)。但是我们注意到,这些字符串不是普通的字符串,这样排序不符合我们OIer的精神,所以我们要利用它们之间的关系来优化。

目前常见的方法有倍增算法、DC3等,我们此处介绍相对便于实现的倍增算法。

倍增算法(Doubling Algorithm)

​ 首先介绍一个显然但重要的性质:

对于两个字符串stra、strb,我们取出它们长度相等且为k的前缀

1
strak=stra.prefix(k-1),strbk=strb.prefix(k-1)//假设下标从0开始

则有:

strak<strbk→stra<strb

strak>strbk→stra>strb

strak=strbk→nothing!

​ 所以对于str的后缀,我们可以每次取出长度为k的前缀,先做出初步的排序,这时有可能会有排名相同的后缀,我们就扩大k,再一次排序。

1级优化

就是这里我们用到了倍增思想,让k→2k

1
2
3
For(k=1;k-前缀相等的字符串;k=2k)

Sort(All str.suffix.prefix(k-1))

**!!!**但是这里还要注意一个问题,所有后缀的长度并不是一样的,那么在比较到后面的时候一定会有后缀的长度小于k。我们的解决方法是:在后面补0。

2级优化

我们发现,每次组成的2k-新前缀都是前一次的两个k-前缀的组合。为了更大限度地利用前一次排序的结果,我们可以选用基数排序(应该没有人不会吧 具体方法在此不作介绍),将前一次的两个k-前缀的排名分别作为第一和第二关键字来排序。

上一张借来的图(以字符串aabaaaab为例):

3级优化

在基数排序的过程中,我们要对第一和第二关键字分别进行排序。优化的做法是,

  1. 先将第二关键字排序

  2. 再用第一关键字排序,第一关键字相等的按刚刚排好的第二关键字排

具体做法:

  1. 所有后缀的第一关键字在上一层已经出现,但是起始位置为len-k~len-1的后缀的第二关键字会超出字符串而需要补零(见1级优化),这时候这些第二关键字一定会排在最前面,我们只要一次放入即可。

  2. 另外一部分后缀的第二关键字一定也在上一层出现了,并且一定是从k开始(0~k-1不可能作为本层的第二关键字)且连续的。我们可以利用上一层的SuffixArray依次放入这些第二关键字(排名)即可。

  3. 接下来对第一关键字排序。由于处理的大多是字符,这里我们选用桶排序来提高速度。

于是这样,原来的基数排序的概念已经支离破碎了,每一层的时间复杂度降到O(len+maxRank)。

4级小优化

在刚开始k比较小的时候,因为当前排名相等的后缀数量比较少,所以最大排名(排名数量)p也比较小,桶排序用到的“桶”较少。我们令桶的最大编号为p,可以节省初始化作为桶的数组的时间。

实现

上面介绍的基本是算法层面的优化,但是还远远不止于此。接下来我们以一份流传甚广的后缀数组模板代码介绍代码中的优化。

下面是一份本人改写的后缀数组代码,优化的点都已经在注释中指出,欢迎指出错误或者提供可以hack的数据。

本人代码

  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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
#include<iostream> 
#include<cstring>
#include<string>
#include<cstdio>
#include<cstdlib>
using namespace std;
const int maxn=1000;
char ss[2*maxn+1];//原字符串 
int s[2*maxn+1];//转成int的新字符串 
int lib1[maxn],lib2[maxn];
int* Rankkey1=lib1;
int* key2=lib2;
//key2[i]表示按第二关键字排名为i的后缀起始位置,
//!!!但是Rankkey1[i]表示起始位置为i的第一关键字(排名),最终就是Rank数组 
int cnt[maxn];
template<typename T>
void Swap(T& a,T& b){T t=a;a=b,b=t;}

int height[maxn];
int h[maxn];
int lenss,maxalpha=0;
bool Equal(int* Rank1,int pos1,int pos2,int len)//用来判断两个后缀排名是否相同 
{
	return Rank1[pos1]==Rank1[pos2]
		 &&Rank1[pos1+len]==Rank1[pos2+len];
}

int sa[maxn];//suffixArray
int* Rank;
void initSuffixAndRankArray(int* str,int len,int maxs)
{                  
/* 参数含义 
str:转成int数组的字符串
len:串长
maxs:最大需要“桶”的数量(即数组最大值)
*/ 
	//对长度为1前缀进行桶排序 (计数排序)(!不是基数排序) 
	memset(cnt,0,(maxs+1)*sizeof(int));
	for(int i=0;i<=len-1;++i)//装桶 
		++cnt[str[i]];
	for(int i=1;i<=maxs;++i)//求出某个桶的排名 
		cnt[i]+=cnt[i-1]; 
	for(int i=len-1;i>=0;--i)
		--cnt[str[i]],//这样可以使所有排名都不一样,cnt[str[i]]就是当前str[i]的排名 
		sa[cnt[str[i]]]=i;//第i个后缀的排名对应的排位就是i的宝座(逃) 
	
	//初始化第一关键字为:1-前缀
	memcpy(Rankkey1,str,(len+1)*(sizeof(int)));
	int p=0;
	for(int k=1;p<=len-1-1;k=k*2,maxs=p)//从k=1,k=2,k=4倍增 将maxs=p是减少桶的数量的优化 
	{
		//cout<<p<<endl; 
		//将第二关键字排序
		p=0;//p用来存放排名,当p=len-1时所有后缀排名不同 
		for(int i=len-k;i<=len-1;++i)
			key2[p++]=i;
		/*从pos=(len-k)开始的后缀,由于 pos+k>=len>len-1,
		其第二关键字为0 (我们自动补 0)所以直接排在最前面*/
		
		for(int i=0;i<=len-1;++i)
			if(sa[i]>=k)
				key2[p++]=sa[i]-k;
		/*只有从k开始的k-前缀可以作为第二关键字,
		直接利用上一层sa排好的顺序放入即可*/
		//注意,我们在排第二关键字的全过程并没有直接用到该层第二关键字的真实值! 
		
		//将第一关键字排序
		#define tmp(i) Rankkey1[key2[i]]//第二关键字排在第i位的后缀的第一关键字 
		//接下来与开头基本相同 
		memset(cnt,0,(maxs+1)*sizeof(int));
		for(int i=0;i<=len-1;++i)//装桶 
			++cnt[tmp(i)];
		for(int i=1;i<=maxs;++i)
			cnt[i]+=cnt[i-1];
		for(int i=len-1;i>=0;--i)//此处i为第二关键字 
			--cnt[tmp(i)],
			sa[cnt[tmp(i)]]=key2[i];
			/*在第一关键字相同的时候,第二关键字大的先算到,
			cnt[tmp[i]]更大,排名更靠后。
			右边的key2[i]用来对应获取后缀起始位置*/ 
		
		//更新第一关键字排名
		/*判断程序能否结束的标准就是所有后缀的排名是不是都不一样
		由于p是所有后缀的种数(本算法中从0开始),那么当p=len-1时就可以结束了*/ 
		int* oldRankkey1=Rankkey1;
		Rankkey1=key2;//借利用已经没有用的key2的内存空间 
		Rankkey1[sa[0]]=0;
		p=0;
		for(int i=1;i<=len-1;++i)
		{
			if(Equal(oldRankkey1,sa[i-1],sa[i],k))
			{
				Rankkey1[sa[i]]=p;//如果相等,与之前的排名相同 
			}
			else Rankkey1[sa[i]]=++p;
		}
		key2=oldRankkey1;//还给key2原来Rankkey1的内存空间,实际相当于交换 
	}
}
void initHeight(int* str,int len)
{
	int i,j,k;
	for(i=k=0;i<=len-1;height[Rank[i++]]=k)
    for(k?k--:0,j=sa[Rank[i]-1];str[i+k]==str[j+k];k++);
}
void outsa()
{

	cout<<"SuffixArray:"<<endl;
	for(int i=1;i<=lenss-1;++i)
	{
		printf("%4d",i-1);
	}
	cout<<endl;
	for(int i=1;i<=lenss-1;++i)
	{
		printf("%4d",sa[i]);
	}
	cout<<endl;
}
void outRank()
{
	cout<<"RankArray:"<<endl;
	for(int i=0;i<=lenss-1-1;++i)
	{
		printf("%4d",i);
	}
	cout<<endl;
	for(int i=0;i<=lenss-1-1;++i)
	{
		printf("%4d",Rank[i]-1);
	}
	cout<<endl;
}
void outsuffix()
{
	string e;
	for(int i=1;i<=lenss-1;++i)
	{
		cout<<string(ss).substr(sa[i],lenss-sa[i])<<endl;
	}
}
void outheight(int len)
{
	for(int i=0;i<=len-1;++i)
	{
		cout<<height[i]<<' ';
	}
	cout<<endl;
}
int main()
{
	freopen("suffixArray.in","r",stdin);
	//freopen("suffixArray.out","w",stdout);
	scanf("%s",ss);
	int i;
	lenss=strlen(ss);
	//转字符串 
	for(i=0;i<=lenss-1;i++)
	{
		s[i]=ss[i]-'a'+1;
		if(s[i]>maxalpha)maxalpha=s[i];
	}
	s[lenss]=0;
	++lenss; 
	/*!!!之所以需要加一,是因为之前的Equal()中使用Rankkey1来判断排名
	但是这样会出现一个bug,一些后缀的第二关键字排名为0
	但是还有一些后缀长度不到k,后面的Rank1[pos+len]是初始化的0,
	这样就造成了错误判相等
	加长度后相当于"0"这个后缀长期占据第0名的位置,就可以避免这个问题
	但是输出suffixArray的时候需要从1开始 
	*/ 
	initSuffixAndRankArray(s,lenss,maxalpha);
		Rank=Rankkey1;
	//接下来是输出sa、Rank、按sa排名的后缀 
	initHeight(s,lenss-1);
	outheight(lenss-1);
	return 0;
}

模板代码

 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
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
using namespace std;
const int maxn=20010;
char re[maxn];
int wa[maxn],wb[maxn],wv[maxn],wss[maxn];
int rs[maxn],erank[maxn];
int* rank;
int cmp(int *r,int a,int b,int l)
{return r[a]==r[b]&&r[a+l]==r[b+l];}

void da(int *r,int *sa,int n,int m)
{
	int i,j,p,*x=wa,*y=wb,*t;
	for(i=0;i<m;i++) wss[i]=0;
	for(i=0;i<n;i++) wss[x[i]=r[i]]++;
	for(i=1;i<m;i++) wss[i]+=wss[i-1];
	for(i=n-1;i>=0;i--)sa[--wss[x[i]]]=i;
	
	for(j=1,p=1;p<n;j*=2,m=p)
	{
		for(p=0,i=n-j;i<n;i++) y[p++]=i;
		for(i=0;i<n;i++) if(sa[i]>=j) y[p++]=sa[i]-j;
		for(i=0;i<n;i++) wv[i]=x[y[i]];
		for(i=0;i<m;i++) wss[i]=0;
		for(i=0;i<n;i++) wss[wv[i]]++;
		for(i=1;i<m;i++) wss[i]+=wss[i-1];
		for(i=n-1;i>=0;i--) sa[--wss[wv[i]]]=y[i];
		for(t=x,x=y,y=t,p=1,x[sa[0]]=0,i=1;i<n;i++)
			x[sa[i]]=cmp(y,sa[i-1],sa[i],j)?p-1:p++;
	}
	rank=x;
	return;
}
int height[maxn];
void calheight(int *r,int n) // 这里的n是原串的本来长度,即不包括新增的0 
{
    int i,j,k=0;
    for(i=1;i<=n;i++) rank[erank[i]]=i; // 有后缀数组得到名次数组,排名第0的后缀一定是添加的0 
 
    for(i=0;i<n;height[rank[i++]]=k) // 以 i 开始的后缀总能够从以 i-1 开始的后缀中继承 k-1 匹配项出来 
    for(k?k--:0,j=erank[rank[i]-1];r[i+k]==r[j+k];k++); // 进行一个暴力的匹配,但是整个算法的时间复杂度还是O(n)的 
    return;
}
void outheight(int len)
{
	for(int i=0;i<=len-1;++i)
	{
		cout<<height[i]<<' ';
	}
	cout<<endl;
}
int main()
{
	freopen("suffixArray.in","r",stdin);
	//freopen("suffixArray.out","w",stdout);
	scanf("%s",re);
	int i,len=strlen(re);
	for(i=0;i<len;i++)
		rs[i]=re[i]-'a'+1;
	da(rs,erank,len+1,88);
	calheight(rs,len);
	outheight(len);
	return 0;
}

参考资料

0%