后缀数组
好吧 由于我是个蒟蒻 c++不常用 特别是字符串渣成翔……pascal还好说 c++么……表示第一次用C++搞这个哎 ToT 顺便贴一发
后缀数组的实现
本节主要介绍后缀数组的两种实现方法: 倍增算法和 3 DC3 算法,并对两种算
法进行了比较。可能有的读者会认为这两种算法难以理解,即使理解了也难以用
程序实现。本节针对这个问题,在介绍这两种算法的基础上,还给出了简洁高效
的代码。其中倍增算法只有 25 行,DC3 算法只有 40 行。
1 1.1 基本定义
子串:字符串 S 的子串 r[i..j],i≤j,表示 r 串中从 i 到 j 这一段,
也就是顺次排列 r[i],r[i+1],...,r[j]形成的字符串。
后缀:后缀是指从某个位置 i 开始到整个串末尾结束的一个特殊子串。字
符 串 r 的 从 第 i 个 字 符 开 始 的 后 缀 表 示 为 Suffix(i) , 也 就 是
Suffix(i)=r[i..len(r)]。
大小比较:关于字符串的大小比较,是指通常所说的“字典顺序 ” 比较,也
就是对于两个字符串 u、v,令 i 从 1 开始顺次比较 u[i]和 v[i],如果
u[i]=v[i]则令 i 加 1,否则若 u[i]<v[i]则认为 u<v,u[i]>v[i]则认为 u>v
(也就是 v<u ) ,比较结束。如果 i>len(u)或者 i>len(v)仍比较不出结果,那
么 若 len(u)<len(v) 则 认 为 u<v , 若 len(u)=len(v) 则 认 为 u=v , 若
len(u)>len(v)则 u>v。
从字符串的大小比较的定义来看, S 的两个开头位置不同的后缀 u 和 v 进
行比较的结果不可能是相等,因为 u=v 的必要条件 len(u)=len(v)在这里不可
能满足。
后缀数组: 后缀数组 SA 是一个一维数组, 它保存 1..n 的某个排列 SA[1] ,
SA[2],……,SA[n],并且保证 Suffix(SA[i]) < Suffix(SA[i+1]),1≤i<n 。
也就是将 S 的 n 个后缀从小到大进行排序之后把排好序的后缀的开头位置顺
次放入 SA 中。
名次数组:名次数组 Rank[i]保存的是 Suffix(i)在所有后缀中从小到大排
列的“名次 ” 。
简单的说,后缀数组是“排第几的是谁? ” ,名次数组是“你排第几? ” 。 容
易看出,后缀数组和名次数组为互逆运算。如图 1 所示。
设字符串的长度为 n。 为了方便比较大小, 可以在字符串后面添加一个字符,
这个字符没有在前面的字符中出现过,而且比前面的字符都要小。在求出名次数
组后,可以仅用 O(1)的时间比较任意两个后缀的大小。在求出后缀数组或名次
数组中的其中一个以后,便可以用 O(n)的时间求出另外一个。任意两个后缀如
果直接比较大小,最多需要比较字符 n 次,也就是说最迟在比较第 n 个字符时一
定能分出“胜负 ” 。
1.2倍增算法
倍增算法的主要思路是:用倍增的方法对每个字符开始的长度为 2 k 的子字
符串进行排序,求出排名,即 rank 值。k 从 0 开始,每次加 1,当 2 k 大于 n 以
后,每个字符开始的长度为 2 k 的子字符串便相当于所有的后缀。并且这些子字
符串都一定已经比较出大小,即 rank 值中没有相同的值,那么此时的 rank 值就
是最后的结果。每一次排序都利用上次长度为 2 k-1 的字符串的 rank 值,那么长
度为 2 k 的字符串就可以用两个长度为 2 k-1 的字符串的排名作为关键字表示,然
后进行基数排序, 便得出了长度为 2 k 的字符串的 rank 值。以字符串 “aabaaaab ”
为例, 整个过程如图 2 所示。其中 x 、 y 是表示长度为 2 k 的字符串的两个关键 字 。
具体实现:
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];}
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++) 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;
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;
}
待排序的字符串放在 r 数组中,从 r[0]到 r[n-1],长度为 n,且最大值小
于 m。为了函数操作的方便,约定除 r[n-1]外所有的 r[i]都大于 0, r[n-1]=0 。
函数结束后,结果放在 sa 数组中,从 sa[0]到 sa[n-1]。
函数的第一步,要对长度为 1 的字符串进行排序。一般来说,在字符串的题
目中,r 的最大值不会很大,所以这里使用了基数排序。如果 r 的最大值很大,
那么把这段代码改成快速排序。代码:
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;
这里 x 数组保存的值相当于是 rank 值。下面的操作只是用 x 数组来比较字
符的大小,所以没有必要求出当前真实的 rank 值。
接下来进行若干次基数排序,在实现的时候,这里有一个小优化。基数排序
要分两次,第一次是对第二关键字排序,第二次是对第一关键字排序。对第二关
键字排序的结果实际上可以利用上一次求得的 sa 直接算出, 没有必要再算一次 。
代码:
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;
其中变量j是当前字符串的长度, 数组y保存的是对第二关键字排序的结 果 。
然后要对第一关键字进行排序,代码:
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];
这样便求出了新的 sa 值。在求出 sa 后,下一步是计算 rank 值。这里要注
意的是,可能有多个字符串的 rank 值是相同的,所以必须比较两个字符串是否
完全相同, y 数组的值已经没有必要保存, 为了节省空间, 这里用 y 数组保存 rank
值。这里又有一个小优化,将 x 和 y 定义为指针类型,复制整个数组的操作可以
用交换指针的值代替,不必将数组中值一个一个的复制。代码:
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++;
其中 cmp 函数的代码是:
int cmp(int *r,int a,int b,int l)
{return r[a]==r[b]&&r[a+l]==r[b+l];}
这里可以看到规定 r[n-1]=0 的好处,如果 r[a]=r[b],说明以 r[a]或 r[b]
开头的长度为l的字符串肯定不包括字符r[n-1] , 所以调用变量 r[a+l]和r[b+l]
不会导致数组下标越界,这样就不需要做特殊判断。执行完上面的代码后,rank
值保存在 x 数组中,而变量 p 的结果实际上就是不同的字符串的个数。这里可以
加一个小优化, 如果 p 等于 n, 那么函数可以结束。因为在当前长度的字符串 中 ,
已经没有相同的字符串,接下来的排序不会改变 rank 值。例如图 1 中的第四次
排序,实际上是没有必要的。对上面的两段代码,循环的初始赋值和终止条件可
以这样写:
for(j=1,p=1;p<n;j*=2,m=p) {…………}
在第一次排序以后,rank 数组中的最大值小于 p,所以让 m=p。
整个倍增算法基本写好,代码大约 25 行。
算法分析:
倍增算法的时间复杂度比较容易分析。每次基数排序的时间复杂度为 O(n) ,
排序的次数决定于最长公共子串的长度,最坏情况下,排序次数为 logn 次,所
以总的时间复杂度为 O(nlogn)。
1.3 3 DC3 算法
DC3 算法分 3 步:
(1) 、先将后缀分成两部分,然后对第一部分的后缀排序。
将后缀分成两部分,第一部分是后缀 k (k 模 3 不等于 0 ) ,第二部分是后 缀
k(k 模 3 等于 0) 。先对所有起始位置模 3 不等于 0 的后缀进行排序,即对
suffix(1), suffix(2), suffix(4), suffix(5), suffix(7)……进行排序。做
法是将 suffix(1)和 suffix(2)连接,如果这两个后缀的长度不是 3 的倍数,那
先各自在末尾添 0 使得长度都变成 3 的倍数。然后每 3 个字符为一组,进行基数
排序,将每组字符“合并 ” 成一个新的字符。然后用递归的方法求这个新的字符
串的后缀数组。如图 3 所示。在得到新的字符串的 sa 后,便可以计算出原字符
串所有起始位置模 3 不等于 0 的后缀的 sa。要注意的是,原字符串必须以一个
最小的且前面没有出现过的字符结尾,这样才能保证结果正确(请读者思考为什
么) 。
(2) 、利用( 1 1 )的结果,对第二部分的后缀排序。
剩下的后缀是起始位置模 3 等于 0 的后缀, 而这些后缀都可以看成是一个字
符加上一个在(1)中已经求出 rank 的后缀,所以只要一次基数排序便可以求出
剩下的后缀的 sa。
(3) 、将( 1 1 )和(2 2 )的结果合并,即完成对所有后缀排序。
这个合并操作跟合并排序中的合并操作一样。每次需要比较两个后缀的大
小。分两种情况考虑,第一种情况是 suffix(3*i)和 suffix(3*j+1)的比较,可
以把 suffix(3*i)和 suffix(3*j+1)表示成:
suffix(3*i) = r[3*i] + suffix(3*i+1)
suffix(3*j+1) = r[3*j+1] + suffix(3*j+2)
其中 suffix(3*i+1)和 suffix(3*j+2) 的比较可以利用 ( 2 ) 的结果快速得 到 。
第二种情况是 suffix(3*i)和 suffix(3*j+2)的比较,可以把 suffix(3*i)和
suffix(3*j+2)表示成:
suffix(3*i) = r[3*i] + r[3*i+1] + suffix(3*i+2)
suffix(3*j+2) = r[3*j+2] + r[3*j+3] + suffix(3*(j+1)+1)
同样的道理,suffix(3*i+2)和 suffix(3*(j+1)+1) 的比较可以利用(2 )
的结果快速得到。所以每次的比较都可以高效的完成,这也是之前要每 3 个字符
合并,而不是每 2 个字符合并的原因。
具体实现:
#define F(x) ((x)/3+((x)%3==1?0:tb))
#define G(x) ((x)<tb?(x)*3+1:((x)-tb)*3+2)
int wa[maxn],wb[maxn],wv[maxn],ws[maxn];
int c0(int *r,int a,int b)
{return r[a]==r[b]&&r[a+1]==r[b+1]&&r[a+2]==r[b+2];}
int c12(int k,int *r,int a,int b)
{if(k==2) return r[a]<r[b]||r[a]==r[b]&&c12(1,r,a+1,b+1);
else return r[a]<r[b]||r[a]==r[b]&&wv[a+1]<wv[b+1];}
void sort(int *r,int *a,int *b,int n,int m)
{
int i;
for(i=0;i<n;i++) wv[i]=r[a[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--) b[--ws[wv[i]]]=a[i];
return;
}
void dc3(int *r,int *sa,int n,int m)
{
int i,j,*rn=r+n,*san=sa+n,ta=0,tb=(n+1)/3,tbc=0,p;
r[n]=r[n+1]=0;
for(i=0;i<n;i++) if(i%3!=0) wa[tbc++]=i;
sort(r+2,wa,wb,tbc,m);
sort(r+1,wb,wa,tbc,m);
sort(r,wa,wb,tbc,m);
for(p=1,rn[F(wb[0])]=0,i=1;i<tbc;i++)
rn[F(wb[i])]=c0(r,wb[i-1],wb[i])?p-1:p++;
if(p<tbc) dc3(rn,san,tbc,p);
else for(i=0;i<tbc;i++) san[rn[i]]=i;
for(i=0;i<tbc;i++) if(san[i]<tb) wb[ta++]=san[i]*3;
if(n%3==1) wb[ta++]=n-1;
sort(r,wb,wa,ta,m);
for(i=0;i<tbc;i++) wv[wb[i]=G(san[i])]=i;
for(i=0,j=0,p=0;i<ta && j<tbc;p++)
sa[p]=c12(wb[j]%3,r,wa[i],wb[j])?wa[i++]:wb[j++];
for(;i<ta;p++) sa[p]=wa[i++];
for(;j<tbc;p++) sa[p]=wb[j++];
return;
}
各个参数的作用和前面的倍增算法一样,不同的地方是 r 数组和 sa 数组的
大小都要是 3*n,这为了方便下面的递归处理,不用每次都申请新的内存空间。
函数中用到的变量:
int i,j,*rn=r+n,*san=sa+n,ta=0,tb=(n+1)/3,tbc=0,p;
rn 数组保存的是 (1) 中要递归处理的新字符串, san 数组是新字符串的 sa 。
变量 ta 表示起始位置模 3 为 0 的后缀个数, 变量 tb 表示起始位置模 3 为 1 的后
缀个数,已经直接算出。变量 tbc 表示起始位置模 3 为 1 或 2 的后缀个数。先按
(1)中所说的用基数排序把 3 个字符“合并 ” 成一个新的字符。为了方便操作 ,
先将 r[n]和 r[n+1]赋值为 0。
代码:
r[n]=r[n+1]=0;
for(i=0;i<n;i++) if(i%3!=0) wa[tbc++]=i;
sort(r+2,wa,wb,tbc,m);
sort(r+1,wb,wa,tbc,m);
sort(r,wa,wb,tbc,m);
其中 sort 函数的作用是进行基数排序。代码:
void sort(int *r,int *a,int *b,int n,int m)
{
int i;
for(i=0;i<n;i++) wv[i]=r[a[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--) b[--ws[wv[i]]]=a[i];
return;
}
基数排序结束后,新的字符的排名保存在 wb 数组中。
跟倍增算法一样,在基数排序以后,求新的字符串时要判断两个字符组是否
完全相同。代码:
for(p=1,rn[F(wb[0])]=0,i=1;i<tbc;i++)
rn[F(wb[i])]=c0(r,wb[i-1],wb[i])?p-1:p++;
其中 F(x)是计算出原字符串的 suffix(x)在新的字符串中的起始位置,c0
函数是比较是否完全相同,在开头加一段代码:
#define F(x) ((x)/3+((x)%3==1?0:tb))
inline int c0(int *r,int a,int b)
{return r[a]==r[b]&&r[a+1]==r[b+1]&&r[a+2]==r[b+2];}
接下来是递归处理新的字符串,这里和倍增算法一样,可以加一个小优化,
如果 p 等于 tbc,那么说明在新的字符串中没有相同的字符,这样可以直接求出
san 数组,并不用递归处理。代码:
if(p<tbc) dc3(rn,san,tbc,p);
else for(i=0;i<tbc;i++) san[rn[i]]=i;
然后是第(2)步,将所有起始位置模 3 等于 0 的后缀进行排序。其中对第
二关键字的排序结果可以由新字符串的 sa 直接计算得到,没有必要再排一次。
代码:
for(i=0;i<tbc;i++) if(san[i]<tb) wb[ta++]=san[i]*3;
if(n%3==1) wb[ta++]=n-1;
sort(r,wb,wa,ta,m);
for(i=0;i<tbc;i++) wv[wb[i]=G(san[i])]=i;
要注意的是 , 如果 n%3==1 , 要特殊处理 suffix(n-1),因为在 san 数组里并
没有suffix(n) 。 G(x)是计算新字符串的 suffix(x)在原字符串中的位置 , 和F(x)
为互逆运算。在开头加一段:
#define G(x) ((x)<tb?(x)*3+1:((x)-tb)*3+2)。
最后是第(3)步,合并所有后缀的排序结果,保存在 sa 数组中。代码:
for(i=0,j=0,p=0;i<ta && j<tbc;p++)
sa[p]=c12(wb[j]%3,r,wa[i],wb[j])?wa[i++]:wb[j++];
for(;i<ta;p++) sa[p]=wa[i++];
for(;j<tbc;p++) sa[p]=wb[j++];
其中 c12 函数是按(3)中所说的比较后缀大小的函数,k=1 是第一种情况 ,
k=2 是第二种情况。代码:
int c12(int k,int *r,int a,int b)
{if(k==2) return r[a]<r[b]||r[a]==r[b]&&c12(1,r,a+1,b+1);
else return r[a]<r[b]||r[a]==r[b]&&wv[a+1]<wv[b+1];}
整个 DC3 算法基本写好,代码大约 40 行。
算法分析:
假设这个算法的时间复杂度为 f(n)。容易看出第(1)步排序的时间为 O(n)
(一般来说,m 比较小,这里忽略不计),新的字符串的长度不超过 2n/3,求新
字符串的 sa 的时间为 f(2n/3),第(2)和第(3)步的时间都是 O(n)。所以
f(n) = O(n) + f(2n/3)
f(n) ≤ c × n + f(2n/3)
f(n) ≤ c × n + c × (2n/3) + c × (4n/9) + c × (8n/27) + …… ≤ 3c × n
所以 f(n) = O(n)
由此看出,DC3 算法是一个优秀的线性算法。
1 1 1 1 .4 4 4 4 倍增算法与 3 DC3 算法的比较
从时间复杂度、空间复杂度、编程复杂度和实际效率等方面对倍增算法与
DC3 算法进行比较。
时间复杂度:
倍增算法的时间复杂度为 O(nlogn),DC3 算法的时间复杂度为 O(n)。从常
数上看,DC3 算法的常数要比倍增算法大。
空间复杂度:
倍增算法和 DC3 算法的空间复杂度都是 O(n)。按前面所讲的实现方法,倍
增算法所需数组总大小为 6n,DC3 算法所需数组总大小为 10n。
编程复杂度:
倍增算法的源程序长度为 25 行,DC3 算法的源程序长度为 40 行。
实际效率:
测试环境:NOI-linux Pentium(R) 4 CPU 2.80GHz
(不包括读入和输出的时间,单位:ms)
从表中可以看出,DC3 算法在实际效率上还是有一定优势的。倍增算法容易
实现,DC3 算法效率比较高,但是实现起来比倍增算法复杂一些。对于不同的题
目,应当根据数据规模的大小决定使用哪个算法。
这个很明显不是我写的,源自罗穗骞的论文,话说以前听他讲过课,感觉很萌的样子……
另外关于倍增算法,下面有助于理解:
倍增算法
一、主要思路:倍增,s[i..i + 2k − 1]的排名通过s[i..i + 2k − 1 − 1]和s[i + 2k − 1..i + 2k − 1]的排名得到。
二、简要过程:已知每个长度为2k − 1的字符串的排名,则可作为每个长度为2k的字符串求排名的关键字xy,s[i..i + 2k − 1]第一关键字x为s[i..i + 2k − 1 − 1]的排名,第二关键字y为s[i + 2k − 1..i + 2k − 1]的排名。以字符串aabaaaab为例:
- k=0,对每个字符开始的长度为20 = 1的子串进行排序,得到rank[1..8]={1,1,2,1,1,1,1,2}
- k=1,对每个字符开始的长度为21 = 2的子串进行排序:由k=0的rank得关键字xy[1..8]={11,12,21,11,11,11,12,20},得到rank[1..8]={1,2,4,1,1,1,2,3}
- k=2,对每个字符开始的长度为22 = 4的子串进行排序:由k=1的rank得关键字xy[1..8]={14,21,41,11,12,13,20,30},得到rank[1..8]={4,6,8,1,2,3,5,7}
- k=3,对每个字符开始的长度为23 = 8的子串进行排序:由k=2的rank得关键字xy[1..8]={42,63,85,17,20,30,50,70},得到rank[1..8]={4,6,8,1,2,3,5,7}
下面是一道裸题 也就是bzoj1031
1031: [JSOI2007]字符加密Cipher
Time Limit: 10 Sec Memory Limit: 162 MBSubmit: 666 Solved: 259
[Submit][Status][Discuss]
Description
Input
Output
Sample Input
Sample Output
数据规模
对于40%的数据字符串的长度不超过10000。
对于100%的数据字符串的长度不超过100000。
HINT
Source
做法2*原串 然后直接跑后缀数组就行了 可以做模版
#include<cstdio>
#include<algorithm>
#include<iostream>
#define rep(i,n) for(int i=0;i<n;i++)
using namespace std;
const int maxn=100000*2+1;
int n,m,sa[maxn],ta[maxn],tb[maxn],*x=ta,*y=tb;
//ta、tb可有可无 直接换成x[maxn],y[maxn]也可以
char C[maxn];
bool cmp(int i,int j,int l)
{
return y[i]==y[j]&&y[i+l]==y[j+l];
}
//比较函数
void Sort()
{
static int w[maxn];
rep(i,m)w[i]=0;rep(i,n)w[x[y[i]]]++;
rep(i,m-1)w[i+1]+=w[i];
for(int i=n-1;i>=0;i--)sa[--w[x[y[i]]]]=y[i];
}
//基数排序
void Da()
{
int i,j,p;
rep(i,n)x[i]=C[i],y[i]=i;Sort();
rep(i,n)printf("%d ",y[i]);
printf("\n");
for(p=1,j=1;p<n;m=p,j*=2)
{
cout <<n-j<<" "<<j<<endl;
for(p=0,i=n-j;i<n;i++)y[p++]=i;
rep(i,n)if(sa[i]>=j)y[p++]=sa[i]-j;//第二关键字排序(根据上一轮推出)
Sort();
for(swap(x,y),i=1,p=1,x[sa[0]]=0;i<n;i++)
x[sa[i]]=cmp(sa[i-1],sa[i],j)?p-1:p++;//第二关键字排序
}
}
int main()
{
//freopen("in","r",stdin);
gets(C);n=strlen(C);
rep(i,n)C[n+i]=C[i];n<<=1;C[n++]=0;m=256;//原串翻两倍 m为基数排序上界 因为有标点等字符 故为256
Da();int s=n/2;
rep(i,n)if(sa[i]<s)printf("%c",C[s+sa[i]-1]);//为什么是s+sa[i]-1而不是sa[i]-1? 避免sa[i]==0
printf("\n");
}
贴一下WJMZBMR的代码 很诚实地说,我的后缀数组就是从这个代码学来的……写得很简洁,我大概分了部分……嗯就是这样。顺便Orz
写完随笔 感觉自己萌萌哒