[学习笔记]后缀数组
参考:
后缀数组 最详细讲解
上面一篇是转载这一篇的:
一、
后缀:suff(i),后缀要排序
sa[i],排名为i的后缀开始位置
rk[i],i开始位置的后缀的排名。
rk[sa[i]]=i,sa[rk[i]]=i
后缀要排序。
暴力排序:O(n^2logn)
可以利用倍增+基数排序优化到O(nlogn)
以下图片转载自:后缀数组 最详细讲解
为了避免常数过大,不能用vector当桶来用。
以及各种优化具体代码实现在下面。
二、
不能光排序。没有用。
设LCP(i,j)为,排名为i和排名为j的后缀的最长公共前缀长度。(注意i,j是排名)
即:LCP(suff[sa[i]],suff[sa[j]])(这个LCP是最长公共前缀的意思)
性质:
1.
LCP(i,j)=LCP(j,i)
LCP(i,i)=len(sa[i])=n-sa[i]+1
这个太平凡了。
2.LCP(i,k)=min(LCP(i,j),LCP(j,k)) 对于任意1<=i<=j<=k<=n,都满足
证明:
设p=min(LCP(i,j),LCP(j,k))
那么有:LCP(i,j)>=p,LCP(j,k)>=p
那么必然有LCP(i,k)>=p
假如LCP(i,k)=q>p
那么有:s[i+p+1]=s[k+p+1],并且s[i+p+1]!=s[j+p+1],s[j+p+1]!=s[k+p+1]
那么,i,k的字典序显然更接近。不满足:1<=i<=j<=k<=n
所以LCP(i,k)<=p
所以,LCP(i,k)=p
证毕
3.LCP(i,k)=min(LCP(j,j-1)) 对于所有的1<i<=j<=k<=n取一个min
有:LCP(i,k)=min(LCP(i,i+1),LCP(i+1,k)),然后不断拆下去,即可证明。
三、LCP求法
根据以上性质,
人们设计了:height[i]表示,LCP(suff[sa[i]],suff[sa[i-1]])
设h[i]表示,height[rk[i]]
关键性质:
h[i]>=h[i-1]-1;
证明:
设sa[rk[i-1]-1]=k
h[i-1]=p
那么有:suff[k],suff[i-1]的前p位相同。
1.s[k]!=s[i-1]
即p=0,显然成立。
2.s[k]=s[i-1]
suff[k],suff[i-1]的前p位相同。
所以,suff[k+1],suff[i]的前p-1位相同,
LCP(suff[i],suff[k+1])=p-1
由于k+1对于i位置来讲,suff[k+1]比较平凡,不如suff[sa[rk[i]-1]]和suff[i]相似度更大。
所以,一定有:h[i]>=p-1=h[i-1]-1
证毕
有了这个性质,就可以算出height来了。或者h
然后利用LCP(i,k)=min(LCP(j,j-1)) 对于所有的1<i<=j<=k<=n取一个min
这一条就是LCP(i,k)=min(height[j])对于所有的1<i<j<=k<=n取一个min
经常求出height之后,用RMQ等可以直接找到LCP(i,k)
代码:
#include<bits/stdc++.h> #define reg register int #define il inline #define numb (ch^'0') using namespace std; typedef long long ll; il void rd(int &x){ char ch;x=0;bool fl=false; while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true); for(x=numb;isdigit(ch=getchar());x=x*10+numb); (fl==true)&&(x=-x); } namespace Miracle{ const int N=1e6+6; char s[N]; int n,m; int x[2*N],y[2*N],sa[N];//不开两倍,下面求x的时候会RE(当然求x的时候,对sa[i]+k和n+1取min也可以,总之要是0 int c[N],num; int rk[N],hei[N]; void get_SA(){ m=122;//m表示剩下的排名种类个数 for(reg i=1;i<=n;++i) ++c[x[i]=s[i]];//c是桶,x是第一关键字 for(reg i=2;i<=m;++i) c[i]+=c[i-1];//c做前缀和,得到第一关键字为x的后缀排名最少是多少。 for(reg i=1;i<=n;++i) sa[c[x[i]]--]=i;//处理当前sa for(reg k=1;k<=n;k<<=1){ num=0; for(reg i=n-k+1;i<=n;++i) y[++num]=i;//y[i]表示第二关键字排在 第i位的后缀,起始位置是多少 //这里直接把第二关键字是0的放进去,因为它们肯定在最前面 for(reg i=1;i<=n;++i) if(sa[i]>k) y[++num]=sa[i]-k;//如果满足sa[i]>k,就有i会作为sa[i]-k的第二关键字。加进去 for(reg i=1;i<=m;++i) c[i]=0;//清空桶 for(reg i=1;i<=n;++i) ++c[x[i]];//累计第一关键字个数 for(reg i=2;i<=m;++i) c[i]+=c[i-1];//前缀和 for(reg i=n;i>=1;--i) sa[c[x[y[i]]]--]=y[i],y[i]=0;//处理当前sa,注意倒序这样第一关键字相同的,第二关键字一定升序排列 swap(x,y);//复制一下,便于后面处理x[i] num=1; x[sa[1]]=1; for(reg i=2;i<=n;++i) x[sa[i]]=(y[sa[i]]==y[sa[i-1]])&&(y[sa[i]+k]==y[sa[i-1]+k])?num:++num;//手动比较相邻排名是否真的不同,从而处理排名相同问题 if(num==n) break;//排好了 m=num;//重新设置最大关键字 //最后两步可以减小常数 } } void get_HEI(){ for(reg i=1;i<=n;++i) rk[sa[i]]=i;//刷rank int k=0; for(reg i=1;i<=n;++i){ if(rk[i]==1) continue;//h[sa[1]]=0 if(k) --k;//h[i]>=h[i-1]-1 int j=sa[rk[i]-1]; while(j+k<=n&&i+k<=n&&s[j+k]==s[i+k]) ++k;//暴力匹配,因为每次k最多-1,k最多是n,所以复杂度O(n) hei[rk[i]]=k; } } int main(){ scanf("%s",s+1); n=strlen(s+1); get_SA(); for(reg i=1;i<=n;++i){ printf("%d ",sa[i]); }return 0; } } int main(){ Miracle::main(); return 0; } /* Author: *Miracle* Date: 2018/11/14 21:24:52 */
四、例题
子串是所有后缀的所有前缀
对于出现两次的子串,必然起始位置的lcp大于等于长度。
某个子串每多出现一次,就会在某个hei[i]中出现。
所以,个数就是:n*(n-1)/2+n-∑hei[i]
每个出现多次的子串,都一定会被减掉多余的次数次。
更好的考虑方法是:
每个后缀贡献的前缀:就是n-sa[i]+1-height[i]
总共的前缀个数,减去和上一个排名的相同的前缀(即子串)个数。
剩下的一定都是之前没有出现过的。
所以所有之前出现出现过的都减去了。
当然,这个式子∑一下就是上面的那个。
BZOJ 4310: 跳蚤
求出本质不同的子串个数。
二分排名mid
从小到大枚举sa,求出第mid大的本质不同的子串出现位置。
具体来讲,每个位置的sa贡献的不同子串个数是:n-sa[i]+1-height[i]
而且新出现的一定是没有出现的最小的那几个。
所以,如果贡献的子串小于k个,就k-=贡献,继续找。否则就是剩下的第k个位置。
(当然,可以对贡献做一个前缀和,然后二分即可。O(n)->O(logn))
之后贪心砍
从后往前,由于是子串的最大子串,所以,每当得到一个子串比第mid个子串大,必须砍掉。
不能不砍,而如果在之前砍的话,留下的部分字典序一定更大。不优。
所以贪心正确。
比较两个串字典序大小,利用后缀数组直接RMQ求LCP即可。
注意,二分的mid是long long级别。re也要开long long
#include<bits/stdc++.h> #define reg register int #define il inline #define numb (ch^'0') using namespace std; typedef long long ll; il void rd(int &x){ char ch;x=0;bool fl=false; while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true); for(x=numb;isdigit(ch=getchar());x=x*10+numb); (fl==true)&&(x=-x); } namespace Miracle{ const int N=1e5+5; int n,m,p; char s[N]; int sa[N],rk[N],hei[N]; int x[2*N],y[2*N]; int c[N],num; void SA(){ m=230; for(reg i=1;i<=n;++i) ++c[x[i]=s[i]]; for(reg i=2;i<=m;++i) c[i]+=c[i-1]; for(reg i=n;i>=1;--i) sa[c[x[i]]--]=i; for(reg k=1;k<=n;k<<=1){ num=0; for(reg i=n-k+1;i<=n;++i) y[++num]=i; for(reg i=1;i<=n;++i) if(sa[i]>k) y[++num]=sa[i]-k; for(reg i=1;i<=m;++i) c[i]=0; for(reg i=1;i<=n;++i) ++c[x[i]]; for(reg i=2;i<=m;++i) c[i]+=c[i-1]; for(reg i=n;i>=1;--i) sa[c[x[y[i]]]--]=y[i],y[i]=0; swap(x,y); num=1;x[sa[1]]=1; for(reg i=2;i<=n;++i){ x[sa[i]]=(y[sa[i]]==y[sa[i-1]])&&(y[sa[i]+k]==y[sa[i-1]+k])?num:++num; } if(num==n) break; m=num; } } void HEI(){ for(reg i=1;i<=n;++i) rk[sa[i]]=i; int k=0; for(reg i=1;i<=n;++i){ if(rk[i]==1) continue; if(k) --k; int j=sa[rk[i]-1]; while(j+k<=n&&i+k<=n&&s[j+k]==s[i+k]) ++k; hei[rk[i]]=k; } } int f[N][20]; int lg[N]; int LCP(int i,int j){ if(i>j) swap(i,j); if(i==j) return n-sa[i]+1; int len=j-i; int ret=min(f[i+1][lg[len]],f[j-(1<<lg[len])+1][lg[len]]); return ret; } bool cmp(int s1,int d1,int s2,int d2){//s1<=s2 return 1 int lcp=LCP(rk[s1],rk[s2]); //cout<<" lcp "<<lcp<<endl; int l1=d1-s1+1,l2=d2-s2+1; if(l1<=l2){ if(lcp>=l1) return true; if(lcp<l1) return s[s1+lcp]<=s[s2+lcp]; } else{ if(lcp>=l2) return false; if(lcp<l2) return s[s1+lcp]<=s[s2+lcp]; } } ll l,r; int st,nd; bool che(ll mid){ //cout<<" mid ------------------------"<<mid<<endl; int tot=0; st=1,nd=0; ll re=mid; while(re>n-sa[st]+1-hei[st]){ re-=n-sa[st]+1-hei[st]; ++st; } nd=sa[st]+hei[st]-1+re; st=sa[st]; //cout<<" st "<<st<<" nd "<<nd<<endl; int rr=n; for(reg i=n;i>=1;--i){ if(tot>p) return false; while(i>=1&&cmp(i,rr,st,nd)) --i; ++tot; rr=i; //cout<<" tot "<<tot<<" : "<<i<<" "<<rr<<endl; ++i; } return tot<=p; } int main(){ scanf("%d",&p); scanf("%s",s+1); n=strlen(s+1); SA();HEI(); for(reg i=1;i<=n;++i) { lg[i]=(i>>(lg[i-1]+1))?lg[i-1]+1:lg[i-1]; f[i][0]=hei[i]; r+=n-sa[i]+1-hei[i]; } for(reg j=1;(1<<j)<=n;++j){ for(reg i=1;i+(1<<j)-1<=n;++i){ f[i][j]=min(f[i][j-1],f[i+(1<<(j-1))][j-1]); } } l=1; ll ans=0; // for(reg i=1;i<=n;++i){ // cout<<i<<" : "<<rk[i]<<" "<<sa[i]<<" "<<hei[i]<<endl; // } // cout<<" judge "<<cmp(4,5,1,4)<<endl; while(l<=r){ ll mid=(l+r)>>1; if(che(mid)){ ans=mid;r=mid-1; } else l=mid+1; } ans=che(ans); for(reg i=st;i<=nd;++i) putchar(s[i]); return 0; } } int main(){ Miracle::main(); return 0; } /* Author: *Miracle* Date: 2018/11/15 21:46:33 */
P2852 [USACO06DEC]牛奶模式Milk Patterns
对于出现k次的子串,考虑sa序列,枚举最长子串的出现位置中排名最靠后的位置的起始位置,为了找到最长的子串,必然找前面连续k个位置所代表的后缀,
不断求LCP,也就是求LCP(i,i-k+1),就是这个子串长度
单调队列优化即可。
#include<bits/stdc++.h> #define reg register int #define il inline #define numb (ch^'0') using namespace std; typedef long long ll; il void rd(int &x){ char ch;x=0;bool fl=false; while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true); for(x=numb;isdigit(ch=getchar());x=x*10+numb); (fl==true)&&(x=-x); } namespace Miracle{ const int N=1000000+5; int s[N]; int n,m,p; int sa[N],rk[N],hei[N],c[N]; int x[2*N],y[2*N]; int num; void SA(){ num=0; m=1000000; for(reg i=1;i<=n;++i) ++c[x[i]=s[i]]; for(reg i=2;i<=m;++i) c[i]+=c[i-1]; for(reg i=1;i<=n;++i) sa[c[x[i]]--]=i; for(reg k=1;k<=n;k<<=1){ num=0; for(reg i=n-k+1;i<=n;++i) y[++num]=i; for(reg i=1;i<=n;++i) if(sa[i]>k) y[++num]=sa[i]-k; for(reg i=1;i<=m;++i) c[i]=0; for(reg i=1;i<=n;++i) ++c[x[i]]; for(reg i=2;i<=m;++i) c[i]+=c[i-1]; for(reg i=n;i>=1;--i) sa[c[x[y[i]]]--]=y[i],y[i]=0; swap(x,y); num=1; x[sa[1]]=1; for(reg i=2;i<=n;++i){ x[sa[i]]=(y[sa[i]]==y[sa[i-1]])&&(y[sa[i]+k]==y[sa[i-1]+k])?num:++num; } if(num==n) break; m=num; } // for(reg i=1;i<=n;++i){ // cout<<i<<" : "<<sa[i]<<endl; // } } void HEI(){ for(reg i=1;i<=n;++i) rk[sa[i]]=i; int k=0; for(reg i=1;i<=n;++i){ if(rk[i]==1) continue; if(k) --k; int j=sa[rk[i]-1]; while(j+k<=n&&i+k<=n&&s[j+k]==s[i+k]) ++k; hei[rk[i]]=k; } } int q[N],l,r; int wrk(){ l=1; int ret=0; for(reg i=2;i<=p-1;++i){ while(l<=r&&hei[q[r]]>=hei[i]) --r; q[++r]=i; } //if(l<=r) ret=max(ret,hei[q[l]]); for(reg i=p;i<=n;++i){ while(l<=r&&hei[q[r]]>=hei[i]) --r; q[++r]=i; while(l<=r&&q[l]<i-p+2)++l; if(l<=r) ret=max(ret,hei[q[l]]); } return ret; } int main(){ rd(n);rd(p); for(reg i=1;i<=n;++i) rd(s[i]); SA();HEI(); int ans=wrk(); printf("%d\n",ans); return 0; } } int main(){ Miracle::main(); return 0; } /* Author: *Miracle* Date: 2018/11/15 8:52:42 */
另一篇:[HAOI2016]找相同字符
求出height,显然要先求出LCP=r的个数以及最大值。然后前缀处理一下。
求LCP=r的个数及最大值:
1.分治。类似于“找相同字符”那个题,注意,mid左边选择一个最值,mid右边选择一个最值。维护两个指针扫过的位置的最值。
最后合并即可。
#include<bits/stdc++.h>
#define reg register int
#define int long long
#define il inline
#define numb (ch^'0')
using namespace std;
typedef long long ll;
il void rd(int &x){
char ch;x=0;bool fl=false;
while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true);
for(x=numb;isdigit(ch=getchar());x=x*10+numb);
(fl==true)&&(x=-x);
}
namespace Miracle{
const int N=300000+5;
const int inf=0x3f3f3f3f3f3f3f3f;
int n,m;
char s[N];
ll a[N];
struct node{
ll ans;
ll mx;
void init(){
ans=0;mx=-inf;
}
node operator +(const node &b){
node c;c.init();
c.ans=ans+b.ans;
c.mx=max(mx,b.mx);
return c;
}
}b[N];
int sa[N],rk[N],hei[N];
int y[2*N],x[2*N];
int c[N],num;
void SA(){
m=122;
for(reg i=1;i<=n;++i) ++c[x[i]=s[i]];
for(reg i=2;i<=m;++i) c[i]+=c[i-1];
for(reg i=n;i>=1;--i) sa[c[x[i]]--]=i;
for(reg k=1;k<=n;k<<=1){
num=0;
for(reg i=n-k+1;i<=n;++i) y[++num]=i;
for(reg i=1;i<=n;++i) if(sa[i]>k) y[++num]=sa[i]-k;
for(reg i=1;i<=m;++i) c[i]=0;
for(reg i=1;i<=n;++i) ++c[x[i]];
for(reg i=2;i<=m;++i) c[i]+=c[i-1];
for(reg i=n;i>=1;--i) sa[c[x[y[i]]]--]=y[i],y[i]=0;
swap(x,y);
num=1;x[sa[1]]=1;
for(reg i=2;i<=n;++i){
x[sa[i]]=(y[sa[i]]==y[sa[i-1]])&&(y[sa[i]+k]==y[sa[i-1]+k])?num:++num;
}
if(num==n) break;
m=num;
}
}
void HEI(){
for(reg i=1;i<=n;++i) rk[sa[i]]=i;
int k=0;
for(reg i=1;i<=n;++i){
if(rk[i]==1) continue;
if(k) --k;
int j=sa[rk[i]-1];
while(j+k<=n&&i+k<=n&&s[j+k]==s[i+k]) ++k;
hei[rk[i]]=k;
}
}
void divi(int l,int r){
// cout<<" divi "<<l<<" and "<<r<<" ---------------------------"<<endl;
if(l==r) return;
int mid=(l+r)/2;
// cout<<" mid "<<mid<<endl;
ll mxl=-inf,mil=inf;
ll mxr=-inf,mir=inf;
int mi=0x3f3f3f3f;
int ptr=mid+1;
for(reg i=mid+1;i<=r;++i){
mi=min(mi,hei[i]);
while(ptr>=l+1&&hei[ptr]>=mi){
//tmp.upda(a[sa[ptr-1]]);
mil=min(mil,a[sa[ptr-1]]);
mxl=max(mxl,a[sa[ptr-1]]);
ptr--;
}
b[mi].ans+=(mid-ptr+1);
if(mil!=inf)b[mi].mx=max(b[mi].mx,mil*a[sa[i]]);
if(mxl!=-inf)b[mi].mx=max(b[mi].mx,mxl*a[sa[i]]);
// cout<<" ii "<<i<<" ptr "<<ptr<<endl;
// for(reg i=0;i<=2;++i){
// cout<<i<<" : "<<b[i].val[0]<<" "<<b[i].val[1]<<" || "<<b[i].val[2]<<" "<<b[i].val[3]<<endl;
// }
}
// cout<<" half "<<endl;
mi=inf;
ptr=mid;
for(reg i=mid;i>=l;--i){
mi=min(mi,hei[i+1]);
while(ptr+1<=r&&hei[ptr+1]>mi){
//tmp.upda(a[sa[ptr+1]]);
mir=min(mir,a[sa[ptr+1]]);
mxr=max(mxr,a[sa[ptr+1]]);
++ptr;
}
b[mi].ans+=(ptr-mid);
if(mir!=inf)b[mi].mx=max(b[mi].mx,mir*a[sa[i]]);
if(mxr!=-inf)b[mi].mx=max(b[mi].mx,mxr*a[sa[i]]);
}
// cout<<" end "<<endl;
// for(reg i=0;i<=2;++i){
// cout<<i<<" : "<<b[i].ans<<" "<<b[i].mx<<endl;
// }
divi(l,mid);
divi(mid+1,r);
}
int main(){
scanf("%lld",&n);
scanf("%s",s+1);
for(reg i=1;i<=n;++i)scanf("%lld",&a[i]);
SA();HEI();
// for(reg i=1;i<=n;++i){
// cout<<i<<" "<<sa[i]<<" "<<hei[i]<<endl;
// }
// cout<<"--------------------------------"<<endl;
//
// cout<<"hei : ";
// for(reg i=1;i<=n;++i){
// cout<<hei[i]<<" ";
// }cout<<endl;
// cout<<"a[i]: ";
// for(reg i=1;i<=n;++i){
// cout<<a[i]<<" ";
// }cout<<endl;
// cout<<" sa : ";
// for(reg i=1;i<=n;++i){
// cout<<sa[i]<<" ";
// }cout<<endl;
//
for(reg i=0;i<=n;++i)b[i].init();
divi(1,n);
for(reg i=n-1;i>=0;--i){
b[i]=b[i]+b[i+1];
}
for(reg i=0;i<=n-1;++i){
if(b[i].mx==-inf) b[i].mx=0;
printf("%lld %lld\n",b[i].ans,b[i].mx);
}
return 0;
}
}
signed main(){
Miracle::main();
return 0;
}
/*
Author: *Miracle*
Date: 2018/11/15 15:54:36
*/
2.并查集。按照height排序,每次合并i,i+1两个排名所代表的字符串集合。设lcp长度为k,两边的集合两两组合,lcp一定都为k。
把两棵树的size相乘,加入ans[k][0]。每个并查集维护集合内的a[sa[i]]的最大值和最小值,然后ans[k][1]=max(ans[k][1],mi1*mi2,mx1*mx2)
本质是考虑两个要求LCP的i,j,取得min的位置在哪里。在取得min的位置的那里,就会合并这两个集合并且统计答案。
3.SAM做法
利用Parent树DP,类似“差异”一题。对反串建SAM,在LCA处统计LCP个数以及最大值,最小值。然后按照len从后往前处理一下即可。
最后都要前缀处理一下。
五、总结
后缀数组可以对一个串求出height,
就可以快速支持查找i,j位置开始的lcp长度。即LCP(rk[i],rk[j]),然后用RMQ即可优化。
和区间问题比较像。但是不是直观地在原始字符串上处理区间,而是对sa序列进行区间处理。
由于是区间,又是找一个最小值,所以可以利用到的方法就比较多了。
分治,并查集,单调栈,单调队列都可以尝试。
并且相邻两个排名LCP长度是最长的,这个性质有时也会用到。