SA后缀数组的模板与相关题目

2018-04-08Chentiancai_nn?

开头的话

虽然这是SA学习文章,但我关于SA的某些教程应该是不会写的 也许以后有想法了再补。这里有个人的教程写的很不错,我的板子也是从那里抄的。
链接:五分钟搞懂后缀数组!后缀数组解析以及应用(附详解代码)
真是良心博主,真是大佬 .jpg
那么就先上模板代码

代码

这是倍增算法($O(nlogn)$),DC3($O(n)$)不会写,据说常数顶一只log,想了想倍增也够用了。毕竟好记好学。

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
//倍增构造SA,顺便用st表来计算RMQ,这样求lcp是O(1)
const int N = 100005;
int n,m; //n:字符串长度 m:倍增算法中用的
int sa[N*2],height[N*2],Rank[N*2],tax[N*2],tp[N*2]; //大小通常开两倍,不知道为什么,反正线性的内存不卡
int st[18][N];
int bin[18],Log[N];

inline void Rsort(){ //这是基数排序
for(int i = 1;i <= m;++i) tax[i] = 0;
for(int i = 1;i <= n;++i) ++tax[Rank[tp[i]]];
for(int i = 1;i <= m;++i) tax[i] += tax[i-1];
for(int i = n;i;--i) sa[tax[Rank[tp[i]]]--] = tp[i];
}
inline bool cmp(int *f,int x,int y,int w){return f[x] == f[y] && f[x+w] == f[y+w];}
inline void GetSA(){
for(int i = 1;i <= n;++i) Rank[i] = s[i],tp[i] = i;
m = 127; /*字符集大小*/Rsort();
for(int w = 1,p = 1,i;p < n;w += w,m = p){
for(p = 0,i = n-w+1;i <= n;++i) tp[++p] = i;
for(i = 1;i <= n;++i) if(sa[i] > w) tp[++p] = sa[i]-w;
Rsort(); swap(Rank,tp); Rank[sa[1]] = p = 1;
for(i = 2;i <= n;++i) Rank[sa[i]] = cmp(tp,sa[i],sa[i-1],w) ? p : ++p;
}
int j,k = 0;
for(int i = 1;i <= n;height[Rank[i++]] = k)
for(k = k ? k-1 : k,j = sa[Rank[i]-1];s[i+k] == s[j+k];++k);
}
inline void GetST(){
for(int i = 1;i <= n;++i) st[0][i] = height[i];
for(int j = 1;j <= Log[n];++j)
for(int i = 1;i+bin[j]-1 <= n;++i)
st[j][i] = Min(st[j-1][i],st[j-1][i+bin[j-1]]);
}
inline int rmq(int x,int y){
int t = Log[y-x+1];
return Min(st[t][x],st[t][y-bin[t]+1]);
}
inline int lcp(int x,int y){
x = Rank[x]; y = Rank[y];
if(x > y) swap(x,y);
return rmq(x+1,y);
}

练手题目

chapter 1

SA维护出了一个height数组,可以很方便的求出两个子串的lcp,用了st表维护之后查询就是$O(1)$。而求子串的lcp的题目有一些,当然此类题目hash算法也可以轻松做到而且常数更小,虽然查询复杂度是log(二分lcp的长度),至于卡…… 你告诉我Hash Killer III怎么过再说。(正常的双模数hash也相当难卡)。当然此类题目的难点不在SA上。
题目:
1、洛谷:P1117 [NOI2016]优秀的拆分
一上来就是比较难的题目,但冷静分析发现AABB其实可以直接看成两个AA然后,定义两个数组$f[i]$为以第$i$个字符为结尾的AA串个数,$g[i]$以第$i$个字符为开头的AA串的个数。答案就是$\Sigma f[i] \times g[i+1]$。
然后就是怎么求$f,g$这个问题。
我们发现如果有两个串完全相同且串长小于等于两个串的开头相差的距离。
那么设开头相差的距离为$L$串长为$M$则有$M-L+1$个长为$2\times L$的AA串这个随便举个例子自己手玩一下。
然后就是怎么计算出这个。我们可以先枚举AA串一半的长度即$L$
枚举$i4=k\times L,j=i+L$。
记$x=lcp(suf(i),suf(j))$, 记$y=lcs(pre(i-1),pre(j-1))$,
$lcp$为最长公共前缀,$lcs$为最长公共后缀。
$suf$为以第$i$个位置开头的后缀,$pre$为以第$i$个位置结尾的前缀。
如果$x+y \ge L$则得到了第一个串以$i-y$位置开头的两个串,从$i-y$到$i-y+(x+y-L+1)-1$为开头都有长为$L*2$的AA串,然后对于$f$和$g$的贡献是一段区间。这里直接差分即可。
$lcs$的求法就是串翻转求$lcp$。

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
#include<bits/stdc++.h>
using namespace std;
#define Min(x,y) ((x)<(y)?(x):(y))
typedef long long ll;
inline int read(){
int x = 0,f = 1;char ch = getchar();
while(ch < '0' || ch > '9'){if(ch == '-')f = -1;ch = getchar();}
while(ch >= '0' && ch <= '9'){x = x*10+ch-'0';ch = getchar();}
return x*f;
}

const int N = 60005;
char s[N];
int n,m,l,r;
int Log[N],bin[17],f[N],g[N];
struct SA{
char s[N];
int Rank[N],sa[N],height[N],tax[N],tp[N];
int st[17][N];

inline void Rsort(){
for(int i = 0;i <= m;++i) tax[i] = 0;
for(int i = 1;i <= n;++i) ++tax[Rank[tp[i]]];
for(int i = 1;i <= m;++i) tax[i] += tax[i-1];
for(int i = n;i;--i) sa[tax[Rank[tp[i]]]--] = tp[i];
}
inline bool cmp(int *f,int x,int y,int w){
return f[x] == f[y] && f[x+w] == f[y+w];
}
inline void GetSA(){
for(int i = 1;i <= n*2+5;++i) sa[i] = Rank[i] = height[i] = tp[i] = 0;
for(int i = 1;i <= m;++i) tax[i] = 0;
for(int i = 1;i <= n;++i) Rank[i] = s[i],tp[i] = i;
m = 127; Rsort();
for(int w = 1,p = 1,i;p < n;w += w,m = p){
for(p = 0,i = n-w+1;i <= n;++i) tp[++p] = i;
for(i = 1;i <= n;++i) if(sa[i] > w) tp[++p] = sa[i]-w;
Rsort(); swap(Rank,tp); Rank[sa[1]] = p = 1;
for(i = 2;i <= n;++i) Rank[sa[i]] = cmp(tp,sa[i],sa[i-1],w) ? p : ++p;
}
int j,k = 0;
for(int i = 1;i <= n;height[Rank[i++]] = k)
for(k = k ? k-1 : k,j = sa[Rank[i]-1];s[i+k] == s[j+k];++k);
}
inline void GetST(){
for(int i = 1;i <= n;++i) st[0][i] = height[i];
for(int j = 1;j <= Log[n];++j)
for(int i = 1;i+bin[j]-1 <= n;++i)
st[j][i] = Min(st[j-1][i],st[j-1][i+bin[j-1]]);
}
inline int rmq(int x,int y){
int t = Log[y-x+1];
return Min(st[t][x],st[t][y-bin[t]+1]);
}
inline int lcp(int x,int y){
x = Rank[x]; y = Rank[y];
if(x > y) swap(x,y);
return rmq(x+1,y);
}
}A,B;

inline int lcp(int x,int y){return A.lcp(x,y);}
inline int lcs(int x,int y){return B.lcp(n-x+1,n-y+1);}

signed main(){
bin[0] = 1; for(int i = 1;i < 16;++i) bin[i] = bin[i-1]<<1;
Log[0] = -1; for(int i = 1;i < N;++i) Log[i] = Log[i>>1]+1;
int T = read();
while(T--){
scanf("%s",s+1);
n = strlen(s+1);
for(int i = 1;i <= 2*n+5;++i) A.s[i] = B.s[i] = 0;
for(int i = 1;i <= n;++i) A.s[i] = B.s[n-i+1] = s[i];
A.GetSA(); A.GetST(); B.GetSA(); B.GetST();
for(int i = 1;i <= n;++i) f[i] = g[i] = 0;
for(int i = 1;i*2 <= n;++i)
for(int j = i*2;j <= n;j += i){
l = min(lcp(j,j-i),i); r = min(lcs(j-1,j-i-1),i-1);
int tmp = (l+r-i+1);
if(l+r >= i){
f[j+l-tmp]++; f[j+l]--;
g[j-i-r]++; g[j-i-r+tmp]--;
}
}
for(int i = 1;i <= n;++i) f[i] += f[i-1],g[i] += g[i-1];
ll ans = 0;
for(int i = 1;i < n;++i) ans += 1LL*f[i]*g[i+1];
printf("%lld\n", ans);
}
return 0;
}

SA

检测到你还在使用百度这个搜索引擎,
做为一个程序员,这是一种自暴自弃!

作环保的程序员,从不用百度开始!

Close