好多年一直听到这个东西,不久前还研究过一次,但始终没有真正理解,短短的二十行代码,每次都是无功而返,就算打了模板最后还是忘得一干二净,今天难得大花了一比时间先把get_sa给稍微弄清楚了,这里仅对于代码谈谈自己的理解(至于大部分的思想可以参照网上其他的博客和罗穗骞的论文)。
后缀数组,我们要做的第一步也就是给后缀排名,但为了优化复杂度,代码还是大有讲究的,利用倍增的思想以及基数排序,妙哉!
首先讲一下每个数组的用法:
x数组:也就是rank数组,x[i]表示以第i个字符为开头的后缀的排名,实时更新,记录的是名次
sa数组:可以视为rank的反数组,sa[i]表示排名第i的后缀是哪一个,记录的是位置,实时更新,也是最终有用的结果
y数组:在倍增的前半部分作为按第二关键字排序的临时sa数组,在后半部分作临时的x数组,(最迷的一个数组)
c数组:可以视为一个桶,基数排序专用(不会基数排序的可以参考一下网上的博客,其实只要知道思想就行了),而且是按x数组记录的排名来排序的
这里还要申明一点,x数组记录的排名是可以相同的,sa数组却强行让所有的字符串排名两两不同
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#define For(i,a,b) for(i=(a);i<=(b);++i)
#define ll long long
#define rep(i,a,b) for(i=(a);i>=(b);--i)
#define mm(a,b) memset(a,b,sizeof(a))
#define inf 999999999
using namespace std;
const int maxn = 1000100;
char s[maxn];
int c[maxn] ,sa[maxn] ,wx[maxn] ,wy[maxn] ,M = 130;
void get_sa(){
int i,j;int len = strlen(s) - 1;
int *x = wx, *y = wy ,k ,p;//写成指针形式方便下面swap
For(i, 0, len) ++ c[x[i] = s[i]];//从这里开始三行以及下面的长相相似的那几行,都是基数排序,因为只有一个字母,就首先把每一个数值按字典序丢进桶里
//这里的x数组其实就是rank数组,这里刚好可以根据字典序
For(i, 1, M)c[i] = c[i-1] + c[i];//记录前缀和,基数排序
rep(i, len, 0)sa[-- c[x[i]]] = i;//由于桶的性质,排名是递增的,因此可以前缀和之后定义sa数组(可以强行举例模拟一下)前缀和是多少就意味着有多少个数小于等于
for(k = 1;k <= len; k<<=1){//枚举倍增的长度,意味着此时需要更新的sa值是以i开头向后延伸k<<1个长度的字符串的排名
p = 0;
for(i = len;i >= len-k+1; --i)y[p++] = i;//y是临时的sa数组,这里按第二关键字排序,并且由于以i开头的后缀在这个范围内第二关键字为空,因此排名自然最小,因此目前排名最小
For(i, 0, len)if(sa[i] >= k)y[p++] = sa[i] - k;//这里枚举的是第二关键字的开头,因此比如字符串ab,目前按b排名,但实际上却是把整个字符串暂时按第二关键字排序,因此把排名记录到a的位置上。
//开始按第一关键字排序
For(i, 0, M)c[i] = 0;//清空桶
For(i, 0, len)c[x[i]] ++;//这里x记录的还是上一轮的rank,因为要把每一个后缀排名,所以可以无脑加进桶
For(i, 1, M)c[i] = c[i-1] + c[i];//同上
rep(i, len, 0)sa[--c[x[y[i]]]] = y[i];//这里有点迷,但发现结构同上,其实就是要把第二关键字的影响给带到第一关键字来,意味着如果第一关键字相同就可以利用桶的性质做到按第二关键字排序的结果排序,比如说一组数17,35,34,43排序,首先按照个位排序,可以得到3,4,5,7,再将数字的十位伴随过来就成了43 34 35 17,然后按十位排序,又由于是位比个位更影响结果,因此首要按十位的顺序拍好,只有十位相同的才按个位排序17,34,35,43就可以得到正确的排名
swap(x, y);//注意已经交换了数组
x[sa[0]] = 0;p = 1;//p记录已经有几个不同的rank值
For(i, 1, len)x[sa[i]] = y[sa[i]] == y[sa[i-1]] && y[sa[i]+k] == y[sa[i-1]+k] ? p - 1:p ++;//更新x值,这里只是为了特判,如果两个字符串完全相同,就rank值相同,否则p++,至于更新的原理,就是由于与当前字符串sa[i]相似度最高的串一定只会是sa[i-1],sa[i+1],因此只需要考虑sa[i-1]这个串是否和sa[i]相同,另外若当前两个字符串相同,则这两个串的左半边和右半边的rank值也一定分别相同
if(p > len)break;//若每个字符串的排名两两不同,则已排好排名,直接break;
M = p;//由于桶的大小只取决于rank值得大小,因此可以剪枝
}
}
int main(){
#ifndef ONLINE_JUDGE
freopen("THUNDER.in","r",stdin);
freopen("THUNDER.out","w",stdout);
#endif
int i,j;
scanf("%s",s);
get_sa();
int tmp = strlen(s) - 1;
For(i, 0, tmp)printf("%d ",sa[i]+1);printf("\n");
return 0;
}