Top

后缀数组DC3模板


后缀数组一些性质

  • suffix(j)和suffix(k)的最长公共前缀为height【rank【j】+1】到height【rank【k】】中的最小值;

  • 待更新

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
/*
POJ-2406
*/
#include<stdio.h>
#include<string.h>
#include<map>
#include<algorithm>
using namespace std;
#define REP(i,n) for(i=0;i<(n);++i)
#define UPTO(i,l,h) for(i=(l);i<=(h);++i)
#define DOWN(i,h,l) for(i=(h);i>=(l);--i)
const int maxn=1e6+10;
const int mod=1e9+7;
typedef long long ll;
template <typename T, int LEN>
struct suffixarray{
//rank[i] 第i个后缀的排名; SA[i] 排名为i的后缀位置; Height[i] 排名为i的后缀与排名为(i-1)的后缀的LCP
int str[LEN*3],sa[LEN*3];
int rank[LEN],height[LEN];
int id[LEN];
int best[LEN][20];
int len;
bool equal(int *str, int a, int b){
return str[a]==str[b]&&str[a+1]==str[b+1]&&str[a+2]==str[b+2];
}
bool cmp3(int *str, int *nstr, int a, int b){
if(str[a]!=str[b])return str[a]<str[b];
if(str[a+1]!=str[b+1])return str[a+1]<str[b+1];
return nstr[a+b%3]<nstr[b+b%3];
}
void radixsort(int *str, int *sa, int *res, int n, int m){
int i;
REP(i,m)id[i]=0;
REP(i,n)++id[str[sa[i]]];
REP(i,m)id[i+1]+=id[i];
DOWN(i,n-1,0)res[--id[str[sa[i]]]]=sa[i];
}
void dc3(int *str, int *sa, int n, int m){
#define F(x) ((x)/3+((x)%3==1?0:one))
#define G(x) ((x)<one?(x)*3+1:((x)-one)*3+2)
int *nstr=str+n, *nsa=sa+n, *tmpa=rank, *tmpb=height;
int i,j,k,len=0,num=0,zero=0,one=(n+1)/3;
REP(i,n)if(i%3)tmpa[len++]=i;
str[n]=str[n+1]=0;
radixsort(str+2, tmpa, tmpb, len, m);
radixsort(str+1, tmpb, tmpa, len, m);
radixsort(str+0, tmpa, tmpb, len, m);
nstr[F(tmpb[0])]=num++;
UPTO(i,1,len-1)
nstr[F(tmpb[i])]=equal(str,tmpb[i-1],tmpb[i])?num-1:num++;
if(num<len)dc3(nstr,nsa,len,num);
else REP(i,len)nsa[nstr[i]]=i;
if(n%3==1)tmpa[zero++]=n-1;
REP(i,len)if(nsa[i]<one)tmpa[zero++]=nsa[i]*3;
radixsort(str, tmpa, tmpb, zero, m);
REP(i,len)tmpa[nsa[i]=G(nsa[i])]=i;
i=j=0;
REP(k,n)
if(j>=len||(i<zero&&cmp3(str,tmpa,tmpb[i],nsa[j])))sa[k]=tmpb[i++];
else sa[k]=nsa[j++];
}
void initSA(T *s, int n,int m){
int i,j,k=0;
str[len=n]=0;//末尾增加一个0,这样就省去一些特殊情况的讨论,也就是最后一个mod 3刚好等于0
REP(i,n)str[i]=s[i];
dc3(str,sa,n+1,m); //可以切换成dc3
REP(i,n)sa[i]=sa[i+1];//第0小的默认为最后一个字符0,所以答案向前移动一位,da算法不用
//da(str,sa,n,m);
REP(i,n)rank[sa[i]]=i;
REP(i,n)//计算height数组
{
if(k)--k;
if(rank[i])for(j=sa[rank[i]-1];str[i+k]==str[j+k];++k);
else k=0;
height[rank[i]]=k;
}
}
void initRMQ(){
int i,j;
int m=(int)(log(len*1.0)/log(2.0));
REP(i,len)best[i][0]=height[i];
for(j=1;j<=m;++j)
for(i=0;i+(1<<j)-1<len;++i)
best[i][j]=min(best[i][j-1],best[i+(1<<(j-1))][j-1]);
}
int RMQ(int l, int r){
int k=int(log(r-l+1.0)/log(2.0));
return min(best[l][k],best[r-(1<<k)+1][k]);
}
int LCPSA(int a, int b){//查询区间RMQ(i,j)
a=rank[a],b=rank[b];
if(a>b)swap(a,b);
return RMQ(a+1,b);
}
};
suffixarray<char,maxn> msa;
map<int ,int > mymap; //计算m,m表示不同字符个数,如果是字母直接用256
char s[maxn];
int main(){
while(scanf("%s",s)){
if(s[0]=='.')break;
int len= strlen(s);
msa.initSA(s,len,256);
int flag = 0;
for (int i = 1; i <= len; i++) {
if (len % i == 0 && msa.rank[0] == msa.rank[i] + 1 && msa.height[msa.rank[0]] == len - i) {
printf("%d\n", len / i);
flag = 1;
break;
}
}
if (!flag)printf("1\n");
}

return 0;
}

一些模板利用

求重复次数最多的连续重复子串

例题:SPOJ - REPEATS-Repeats

1
2
3
4
5
6
7
8
9
10
11
12
13
msa.initRMQ();
for(int i=1;i<=len;i++){ //枚举长度
for(int j=0;j+i<len;j+=i){ //+i极大的降低了复杂度
//但也产生了需要向前比较的问题
ans=msa.LCPSA(j,j+i); //公共后缀的长度
k=j-(i-ans%i); //前推到k位置
ans=ans/i+1; //出现次数
if(k>=0&&msa.LCPSA(k,k+i)>=i)ans++;
//printf("L=%d,R=%d\n",i,ans);
Max=max(Max,ans);
}
}
printf("%d\n",Max);

求最小循环节

例题:POJ-2406 - Power Strings

1
2
3
4
5
6
7
8
9
int flag = 0;
for (int i = 1; i <= len; i++) {
if (len % i == 0 && msa.rank[0] == msa.rank[i] + 1 && msa.height[msa.rank[0]] == len - i) {
printf("%d\n", len / i);
flag = 1;
break;
}
}
if (!flag)printf("1\n");

求两串中长度大于k的公共子串有多少个。

例题:POJ - 3415 - Common Substrings

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
ll solve(int k){
int top=0;
ll ans=0,tot=0;
for(int i=1;i<n;i++){
if(msa.height[i]<k)tot=0,top=0;
else{
int cnt=0;
if(msa.sa[i-1]<lena){ //说明不在同一串中
cnt++;
tot+=msa.height[i]-k+1;
}
while(top&&msa.height[i]<=sta[top-1][0]){
top--;
tot+=(msa.height[i]-sta[top][0])*sta[top][1];
cnt+=sta[top][1];
}
sta[top][0]=msa.height[i],sta[top++][1]=cnt;
if(msa.sa[i]>lena)ans+=tot;
}
}

tot=top=0;
for(int i=1;i<n;i++){
if(msa.height[i]<k)tot=0,top=0;
else{
int cnt=0;
if(msa.sa[i-1]>lena){ //说明不在同一串中
cnt++;
tot+=msa.height[i]-k+1;
}
while(top&&msa.height[i]<=sta[top-1][0]){
top--;
tot+=(msa.height[i]-sta[top][0])*sta[top][1];
cnt+=sta[top][1];
}
sta[top][0]=msa.height[i],sta[top++][1]=cnt;
if(msa.sa[i]<lena)ans+=tot;
}
}
return ans;
}
文章版权为Anoyer博客所有,转载请以链接形式标明本文地址