后缀数组详解+模板
后缀数组
注
SA[] 第几名是谁
后缀数组:后缀数组 SA 是一个一维数组, 它保存 1..n 的某个排列 SA[1] ,SA[2],……,SA[n],并且保证 Suffix(SA[i]) < Suffix(SA[i+1]),1≤i<n 。也就是将 S 的 n 个后缀从小到大进行排序之后把排好序的后缀的开头位置顺次放入 SA 中。
Rank[] 谁是第几名名次数组:名次数组 Rank[i]保存的是 Suffix(i)在所有后缀中从小到大排列的“名次 ” 。
r[]:原始数据j当前字符串的长度,每次循环根据2个j长度的字符串的排名求得2j长度字符串的排名.
y[]:指示长度为2j的字符串的第二关键字的排序结果,通过存储2j长字符串的第一关键字的下标进行指示.
wv[]:2j长字符串的第一关键字的排名序号.
ws[]:计数数组,计数排序用到.
x[]:一开始是原始数据r的拷贝(其实也表示长度为1的字符串的排名),之后表示2j长度字符串的排名.
p:不同排名的个数.
片段
1.对长度为1的字符串进行排序(函数的第一步)
for(i=0;i<m;i++) ws[i]=0;
for(i=0;i<n;i++) ws[x[i]=r[i]]++;
for(i=1;i<m;i++) ws[i]+=ws[i-1];
for(i=n-1;i>=0;i--) sa[--ws[x[i]]]=i;
①用的是基数排序,也可以使用其它的排序
②r[]存储原本输入的字符串,x[]是对r[]的ASCII呈现(便于排序)
③m是一个估计数字,代表ASCII最大值,在循环中做边界
④n在这里是字符串的长度+1,后面的加加减减有所体现(貌似不介意直接用字符串的长度)
⑤最后一行比较难懂,但实践证明它确实是正确的,sa[i]=j表示第i名是j。
ws[i]是对第i及之前字符出现次数的累加,越往后ws[i]越大,而且对应的字符数值越大,举个例子,如果某一字符串为aaabaa,则a出现的次数为5,b出现的次数为1,按上述原理,可以看做ws[a]=5,ws[b]=6,固然a都在前5名,b在第六名。
对aabaaaab进行输出后为801345627,按照sa的定义对应起来
aabaaaab ~
23845679 1 非常正确
理解了这个,最后一行就能明白了
2.进行若干次基数排序
因为前面排序的名次可能有重复,所以要再进行若干次,直到所有的名次都不再相同
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++) Ws[i]=0;
for(i=0; i<n; i++) Ws[wv[i]]++;
for(i=1; i<m; i++) Ws[i]+=Ws[i-1];
for(i=n-1; i>=0; i--) sa[--Ws[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++;
}
相对于上面函数的第一步来说,这一坨代码更加复杂了
①从最外层循环可以看出,j是处于倍增状态的,代表正在比较的每一小段字符串的长度
②循环内的第一行,循环了j-1次,是对后面几个数的提前处理(其第二关键字都为0)如图
即所有加0的数
③第二行,再翻上去看一眼sa的作用。首先要明白这一行抛弃了一些东西,
由于是对第二关键字的排序,第一关键字先不看,所以有一条件if(sa[i]>=j)
这条语句后面y[p++]=sa[i]-j,要减去j也是因为这个
到这里,第二关键字的排序就完成了
④开始第一关键字的排序
假设需要排序的数为92 71 10 80 63 90
那么y[]=3 4 6 2 1 5 即对第二关键字排序后名次递增所对应的序号
x[]=10 80 90 71 92 63 即对第二关键字排序的结果
for(i=0; i<n; i++) wv[i]=x[y[i]];将x[]数组拷贝到wv[]中
⑤剩下的基数排序就与对长度为1的字符串进行排序一样了
完整的代码(参考理解)
#include <cstdio>
#include <iostream>
#include <cstring>
#define LL long long
#define ULL unsigned long long
using namespace std;
const int MAXN=100010;
//以下为倍增算法求后缀数组
int wa[MAXN],wb[MAXN],wv[MAXN],Ws[MAXN];
int cmp(int *r,int a,int b,int l)
{return r[a]==r[b]&&r[a+l]==r[b+l];}
/**< 传入参数:str,sa,len+1,ASCII_MAX+1 */
void da(const char r[],int sa[],int n,int m)
{
int i,j,p,*x=wa,*y=wb,*t;
for(i=0; i<m; i++) Ws[i]=0;
for(i=0; i<n; i++) Ws[x[i]=r[i]]++;//以字符的ascii码为下标
for(i=1; i<m; i++) Ws[i]+=Ws[i-1];
for(i=n-1; i>=0; i--) sa[--Ws[x[i]]]=i;
/*cout<<"SA"<<endl;;
for(int i=0;i<n+1;i++)cout<<sa[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++) Ws[i]=0;
for(i=0; i<n; i++) Ws[wv[i]]++;
for(i=1; i<m; i++) Ws[i]+=Ws[i-1];
for(i=n-1; i>=0; i--) sa[--Ws[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++;
}
return;
}
int sa[MAXN],Rank[MAXN],height[MAXN];
//求height数组
/**< str,sa,len */
void calheight(const char *r,int *sa,int n)
{
int i,j,k=0;
for(i=1; i<=n; i++) Rank[sa[i]]=i;
for(i=0; i<n; height[Rank[i++]]=k)
for(k?k--:0,j=sa[Rank[i]-1]; r[i+k]==r[j+k]; k++);
// Unified
for(int i=n;i>=1;--i) ++sa[i],Rank[i]=Rank[i-1];
}
char str[MAXN];
int main()
{
while(scanf("%s",str)!=EOF)
{
int len=strlen(str);
da(str,sa,len+1,130);
calheight(str,sa,len);
puts("--------------All Suffix--------------");
for(int i=1; i<=len; ++i)
{
printf("%d:\t",i);
for(int j=i-1; j<len; ++j)
printf("%c",str[j]);
puts("");
}
puts("");
puts("-------------After sort---------------");
for(int i=1; i<=len; ++i)
{
printf("sa[%2d ] = %2d\t",i,sa[i]);
for(int j=sa[i]-1; j<len; ++j)
printf("%c",str[j]);
puts("");
}
puts("");
puts("---------------Height-----------------");
for(int i=1; i<=len; ++i)
printf("height[%2d ]=%2d \n",i,height[i]);
puts("");
puts("----------------Rank------------------");
for(int i=1; i<=len; ++i)
printf("Rank[%2d ] = %2d\n",i,Rank[i]);
puts("------------------END-----------------");
}
return 0;
}