前言
建议先学习 基数排序
。本文讲解倍增法。想学奇技淫巧
D
C
3
DC3
DC3 的,建议自己看
2009
\rm 2009
2009 年集训队论文。
概述
一种非常普遍的应用于字符串问题的 思想(我姑且把它称为思想,而非算法或者数据结构)。——之后我们会知道,它可以很轻易地推广到树上。
思想
就是 把所有的后缀拿出来排序,然后就可以进行很多骚操作。比如,子串就是一个后缀的前缀。接下来,是正式的后缀数组。
什么是后缀数组
需要几个定义——
后缀
将字符串记作 ( s 0 , s 1 , s 2 , … , s n − 1 ) (s_0,s_1,s_2,\dots,s_{n-1}) (s0,s1,s2,…,sn−1) 。称 ( a 0 , a 1 , a 2 , … , a m − 1 ) (a_0,a_1,a_2,\dots,a_{m-1}) (a0,a1,a2,…,am−1) 是 s s s 的后缀,当且仅当 m ≤ n m\le n m≤n 且 ∀ x ∈ [ 0 , m ) , a x = s x + n − m \forall x\in[0,m),\;a_x=s_{x+n-m} ∀x∈[0,m),ax=sx+n−m 。
更简单一点:后缀是包含最后一个字符的子串。
在后面的叙述中,将以第 x x x 位开头的后缀,即 ( s x , s x + 1 , s x + 2 , … , s n − 1 ) (s_x,s_{x+1},s_{x+2},\dots,s_{n-1}) (sx,sx+1,sx+2,…,sn−1),称为:后缀 x x x 。
后缀数组
如果一个后缀在所有后缀中是第 x x x 大,就将其排名称为 x x x 。排名从零开始编号。
定义这些函数(也即数组):
- rank ( x ) \text{rank}(x) rank(x) 表示,后缀 x x x 的排名是多少。
- 与 rank ( x ) \text{rank}(x) rank(x) 正好相反, s a ( x ) sa(x) sa(x) 表示排名为 x x x 的后缀是后缀 s a ( x ) sa(x) sa(x) 。
对于 s a ( x ) sa(x) sa(x) 的理解一定要明晰,不然很容易学的头大。
我的记忆方法,就是牢记这一点: s a sa sa 中相邻的长得很像; s a sa sa 是排序后的后缀。
怎么求后缀数组
暴力排序
暴力排序复杂度是多少呢? O ( n log n ) \mathcal O(n\log n) O(nlogn) ?也挺快的啊?
但是这是 字符串比较,一次比较可能是 O ( n ) \mathcal O(n) O(n) 的,总复杂度就会趋近于 O ( n 2 log n ) \mathcal O(n^2\log n) O(n2logn) 。
基数排序
对于一个字符串,我们完全也可以进行基数排序嘛!
这样一来,我们只需要排
n
n
n 次,每次
O
(
n
)
\mathcal O(n)
O(n) ,总复杂度
O
(
n
2
)
\mathcal O(n^2)
O(n2) 。还是很慢。
倍增算法
对于一个 s s s 的子串 ( s l , s l + 1 , s l + 2 , … , s r ) (s_l,s_{l+1},s_{l+2},\dots,s_{r}) (sl,sl+1,sl+2,…,sr) ,它是由什么构成的?
设其长度不为一(即 l ≠ r l\ne r l=r ),令 m = ⌊ l + r 2 ⌋ m=\lfloor\frac{l+r}{2}\rfloor m=⌊2l+r⌋ ,那么 ( s l , s l + 1 , s l + 2 , … , s m ) (s_l,s_{l+1},s_{l+2},\dots,s_{m}) (sl,sl+1,sl+2,…,sm) 和 ( s m + 1 , s m + 2 , s m + 3 , … , s r ) (s_{m+1},s_{m+2},s_{m+3},\dots,s_{r}) (sm+1,sm+2,sm+3,…,sr) 拼接起来就是原子串 ( s l , s l + 1 , s l + 2 , … , s r ) (s_l,s_{l+1},s_{l+2},\dots,s_{r}) (sl,sl+1,sl+2,…,sr) 。多美妙啊!
然后,神奇的算法就出炉了:每次考虑将前 w w w 个字符拿出来进行排序。也就是说,只考虑每个后缀的前 w w w 个字符,将此作为依据进行大小关系比较。
假设上一次排序,也就是长度为 w 2 \frac{w}{2} 2w 时,结果是 s a sa sa (此时的 s a sa sa 类似于哈希值),那么 ( s x , s x + 1 , s x + 2 , … , s x + w − 1 ) (s_x,s_{x+1},s_{x+2},\dots,s_{x+w-1}) (sx,sx+1,sx+2,…,sx+w−1) 的大小,就是 s a ( x ) sa(x) sa(x) 作为第一关键字、 s a ( x + w 2 ) sa(x+\frac{w}{2}) sa(x+2w) 作为第二关键字进行排序。
没有后半截(满足 x + w 2 ≥ n x+\frac{w}{2}\ge n x+2w≥n )的呢?认为其第二关键字是 − ∞ -\infty −∞ 。毕竟空串是最小的。
代码如下:
int sa[MaxN], rnk[MaxN];
int tmp[MaxN], bucket[MaxN]; // 辅助数组
void getSA(char s[],int n){
int m = 255; // 字符集大小
for(int i=0; i<n; ++i)
rnk[i] = s[i]; // 第一次的Hash值为ASCII
/* 桶排序的过程:始 */
for(int i=0; i<=m; ++i)
bucket[i] = 0; // 清零好习惯
for(int i=0; i<n; ++i)
bucket[rnk[i]] ++;
for(int i=1; i<=m; ++i)
bucket[i] += bucket[i-1];
for(int i=n-1; ~i; --i)
sa[--bucket[rnk[i]]] = i;
/* 桶排序的过程:终 */
for(int w=1,p; w<n; w<<=1){
p = 0; // 桶的指针
/* 按照第二关键字排序 */
for(int i=n-w; i<n; ++i)
tmp[p ++] = i; // 第二关键字为-infty
for(int i=0; i<n; ++i)
if(sa[i] >= w) // 第二关键字为i
tmp[p ++] = sa[i]-w;
/* 此时tmp内的元素按照第二关键字有序 */
/* 桶排序的过程:始 */
for(int i=0; i<=m; ++i)
bucket[i] = 0; // 清零好习惯
for(int i=0; i<n; ++i)
bucket[rnk[i]] ++;
for(int i=1; i<=m; ++i)
bucket[i] += bucket[i-1];
for(int i=n-1; ~i; --i)
sa[--bucket[rnk[tmp[i]]]] = tmp[i];
/* 桶排序的过程:终 */
swap(tmp,rnk); // 目的是将rnk拷贝到tmp中
rnk[sa[0]] = p = 0; // 开始生成新的Hash值
for(int i=1; i<n; ++i){
# define HASH(x) ((x) < n ? tmp[(x)] : -1)
if(HASH(sa[i-1]+w) != HASH(sa[i]+w)) ++ p;
else if(HASH(sa[i-1]) != HASH(sa[i])) ++ p;
rnk[sa[i]] = p; // 与上一个不同,则哈希值加一
# undef HASH // 引用原有rnk,或者是-infty
}
if(p == n-1) break; else m = p;
}
}
再贴一份没有注释并且压行(也没有乱压行啦)的代码。注意这次是从 1 1 1 开始编号的。
# define rep(i,a,b) for(int i=(a); i<=(b); ++i)
# define drep(i,a,b) for(int i=(a); i>=(b); --i)
void getBuc(int n,int m){
memset(buc+1,0,m<<2);
rep(i,1,n) ++ buc[rnk[i]];
rep(i,2,m) buc[i] += buc[i-1];
}
void getSA(char s[],int n){
tmp[n+1] = 0; rep(i,1,n) rnk[i] = s[i];
int m = *max_element(rnk+1,rnk+n+1);
getBuc(n,m); if(m == n) ++ m;
rep(i,1,n) sa[buc[rnk[i]]--] = i;
for(int w=1,p=0; m!=n; w<<=1,m=p,p=0){
getBuc(n,m); rep(i,n-w+1,n) tmp[++ p] = i;
rep(i,1,n) if(sa[i] > w) tmp[++ p] = sa[i]-w;
drep(i,n,1) sa[buc[rnk[tmp[i]]]--] = tmp[i];
memcpy(tmp+1,rnk+1,n<<2); rnk[sa[1]] = p = 1;
rep(i,2,n) rnk[sa[i]] = (tmp[sa[i]] == tmp[sa[i-1]]
&& tmp[sa[i]+w] == tmp[sa[i-1]+w]) ? p : (++ p);
}
}
易错点:看上去,更新
r
a
n
k
\rm rank
rank 数组时,使用了
i
+
w
i+w
i+w 下标,需要将
t
m
p
\rm tmp
tmp 的长度设为
2
n
2n
2n 吗?事实上,当
i
+
w
−
1
>
n
i+w-1>n
i+w−1>n 时,以
i
i
i 开头的长度为
w
w
w 的子串,它不会与别的子串相同(因为它末尾是空字符,是独一无二的),根据 &&
短路性就不会进行判断。所以最多会用到
n
+
1
n+1
n+1 的位置。
后缀数组的兄弟
height \text{height} height 数组
定义 height ( x ) = lcp [ s a ( x ) , s a ( x − 1 ) ] \text{height}(x)=\text{lcp}[sa(x),\;sa(x{\rm-}1)] height(x)=lcp[sa(x),sa(x−1)] 。其中 lcp \text{lcp} lcp 表示最长公共前缀。
那么 lcp [ s a ( x ) , s a ( y ) ] \text{lcp}[sa(x),sa(y)] lcp[sa(x),sa(y)] 就有了快速算法——不妨设 x < y x<y x<y,那么 lcp [ s a ( x ) , s a ( y ) ] = min k = x + 1 y height ( k ) \text{lcp}[sa(x),sa(y)]=\min_{k=x+1}^{y}\text{height}(k) lcp[sa(x),sa(y)]=k=x+1minyheight(k)
为什么呢?请感性理解:设 lcp = ϑ \text{lcp}=\vartheta lcp=ϑ ,那么 ∀ k ∈ ( x , y ] \forall k\in(x,y] ∀k∈(x,y],毫无疑问 s a ( k ) sa(k) sa(k) 的前 ϑ \vartheta ϑ 位与 s a ( x ) sa(x) sa(x) 的前 ϑ \vartheta ϑ 位是一样的。
根本原因,是因为
s
a
sa
sa 中相邻的两个长得很像。说了是感性理解,你去厨房拿刀干嘛。
如何快速求解 height \text{height} height 数组呢?令 h ( x ) = height[rank ( x ) ] h(x)=\text{height[rank}(x)] h(x)=height[rank(x)],那么 h ( x ) h(x) h(x) 是否存在什么规律呢?
考虑 h ( x ) h(x) h(x) 的意义:所有字典序小于后缀 x x x 的后缀,与后缀 x x x 的 lcp \text{lcp} lcp 的最大值。
毕竟,令
x
′
=
rank
(
x
)
x'=\text{rank}(x)
x′=rank(x) ,那么对于另一个
y
(
y
<
x
′
)
y\;(y<x')
y(y<x′),就会有
lcp
[
s
a
(
y
)
,
s
a
(
x
′
)
]
=
min
i
=
y
+
1
x
′
height
(
i
)
≤
height
(
x
′
)
=
h
(
x
)
\text{lcp}[sa(y),sa(x')]=\min_{i=y+1}^{x'}\text{height}(i)\le\text{height}(x')=h(x)
lcp[sa(y),sa(x′)]=i=y+1minx′height(i)≤height(x′)=h(x)
如果是这样的话,假设后缀 ( x − 1 ) (x-1) (x−1) 与后缀 ( y − 1 ) (y-1) (y−1) 有一个 lcp = ϑ ( ϑ ≠ 0 ) \text{lcp}=\vartheta\;(\vartheta\ne0) lcp=ϑ(ϑ=0),那么,显然后缀 x x x 与后缀 y y y 有一个 lcp = ϑ − 1 \text{lcp}=\vartheta-1 lcp=ϑ−1 。并且后缀 x x x 字典序大于后缀 y y y 。如图。
也就是说,
h
(
x
)
⩾
h
(
x
−
1
)
−
1
h(x)\geqslant h(x{-}1)-1
h(x)⩾h(x−1)−1
有了这一点之后,就可以直接暴力判断了——因为每次 h ( x ) h(x) h(x) 只会减少一,总复杂度是 O ( n ) \mathcal O(n) O(n) 。
代码如下:
void GetHeight() {
for(int i=0,j,k=0; i<n; ++i){
if(k) -- k; // 每次减一
if(rnk[i] == 0){ // 没有意义,赋值为0
height[0] = k = 0;
continue;
}
j = sa[rnk[i]-1]; // 与后缀j的lcp
while(s[j+k] == s[i+k]) ++ k;
height[rnk[i]] = k; // 暴力匹配即可
}
}
老规矩,再给一份没有注释的。同样注意下标是从 1 1 1 开始的。
void getHeight(int n){
for(int i=1,j,k=0; i<=n; ++i){
j = sa[rnk[i]-1]; if(k) -- k;
if(rnk[i] == 1) continue;
while(s[i+k] == s[j+k]) ++ k;
heit[rnk[i]] = k;
}
}
例题
洛谷的板题
为了让你有地方找题解——题解中也讲得很好。
#include <cstdio>
#include <iostream>
#include <vector>
#include <cstring>
using namespace std;
inline int readint(){
int a = 0, f = 1; char c = getchar();
for(; c<'0' or c>'9'; c=getchar())
if(c == '-') f = -1;
for(; '0'<=c and c<='9'; c=getchar())
a = (a<<3)+(a<<1)+(c^48);
return a*f;
}
inline void writeint(int x){
if(x < 0) putchar('-'), x = -x;
if(x > 9) writeint(x/10);
putchar('0'+x%10);
}
const int MaxN = 1000050;
int sa[MaxN], rnk[MaxN<<1], n;
char s[MaxN+5];
int bucket[MaxN], tmp[MaxN<<1];
# define x rnk
# define y tmp
void DeBug(){
printf("sa:");
for(int i=0; i<n; ++i)
printf(" %d",sa[i]);
printf("\nrank:");
for(int i=0; i<n; ++i)
printf(" %d",rnk[i]);
printf("\ntmp:");
for(int i=0; i<n; ++i)
printf(" %d",y[i]);
printf("\n");
}
void getSA(int m){
for(int i=0; i<=m; ++i)
bucket[i] = 0;
for(int i=0; i<n; ++i)
bucket[x[i]=s[i]] ++;
for(int i=n; i<(n<<1); ++i)
x[i] = y[i] = -1;
for(int i=1; i<=m; ++i)
bucket[i] += bucket[i-1];
for(int i=n-1; ~i; --i)
sa[--bucket[x[i]]] = i;
for(int w=1; w<n; w<<=1){
int p = 0;
for(int i=n-w; i<n; ++i)
y[p++] = i;
for(int i=0; i<n; ++i)
if(sa[i] >= w)
y[p++] = sa[i]-w;
for(int i=0; i<=m; ++i)
bucket[i] = 0;
for(int i=0; i<n; ++i)
++ bucket[x[i]];
for(int i=1; i<=m; ++i)
bucket[i] += bucket[i-1];
for(int i=n-1; ~i; --i)
sa[--bucket[x[y[i]]]] = y[i];
swap(x,y);
x[sa[0]] = p = 0;
for(int i=1; i<n; ++i){
if(y[sa[i]] != y[sa[i-1]] or y[sa[i]+w] != y[sa[i-1]+w])
++ p;
x[sa[i]] = p;
}
if((m = p) == n-1) break;
}
}
# undef x
# undef y
int main(){
while(true){
s[n] = getchar();
if(s[n] == EOF or s[n] == '\n')
break;
++ n;
}
getSA(255);
for(int i=0; i<n; ++i){
writeint(sa[i]+1);
putchar(' ');
}
putchar('\n');
return 0;
}
可重叠最长重复子串
即,一个子串,这个子串在字符串中出现了至少两次(可以有字符重叠),求最长的。
答案是 height \text{height} height 数组中的最大值。因为子串就是后缀的前缀,出现两次就是 lcp \text{lcp} lcp 。
不可重叠最长重复子串
二分一个答案 ω \omega ω,考虑如何检查。
尝试着从 s a sa sa 里面拿出两个后缀,则一定不能跨过一个小于 ω \omega ω 的 height \text{height} height 。相当于将 height \text{height} height 数组切成了很多段。
对于每一段,由于不能重叠,就必须满足 ∣ x − y ∣ ≥ ω |x-y|\ge \omega ∣x−y∣≥ω(也就是错开了至少 ω \omega ω 位)。所以我们只需要求一个最大的 x x x 和一个最小的 y y y 。
注意一点, max \max max 和 min \min min 都是用 s a sa sa 更新。
代码源于此处(似乎有小错误),核心如下。
bool check(int len) {
int tiny, big; // height[0] = 0
for(int i=0; i<n; ++i){
if(height[i] >= len){
getMax(big,sa[i]);
getMin(tiny,sa[i]);
}
else tiny = big = sa[i];
if(big-tiny >= len) return true;
}
return false;
}
字符串匹配(最长公共子串)
只需要把两个字符串接在一起,中间插入一个特殊字符。
然后查 height \text{height} height 。因为 height \text{height} height 反映 lcp \text{lcp} lcp 。如果两个后缀分居特殊字符两侧,其 lcp \text{lcp} lcp 就是一个公共子串。
本质不同的子串的数量
子串一定是某个后缀的前缀。考虑在 s a sa sa 数组中依次考虑每个后缀带来的前缀。
对于 s a ( x ) sa(x) sa(x) 来说,其前缀有 height ( x ) \text{height}(x) height(x) 个是重复的——他们同样是 s a ( x − 1 ) sa(x-1) sa(x−1) 的前缀,一定是已经被计算过的。而且,这也是最长的一个 lcp \text{lcp} lcp ,所以不会有更多重复了。
最终答案就是总子串数量减去重复的。也就是,
a n s = n ( n + 1 ) 2 − ∑ i = 0 n − 1 height ( i ) ans=\frac{n(n+1)}{2}-\sum_{i=0}^{n-1}\text{height}(i) ans=2n(n+1)−i=0∑n−1height(i)
思路源于此博客。感谢。膜拜。
#include <cstdio>
#include <iostream>
#include <vector>
#include <cstring>
using namespace std;
inline int readint(){
int a = 0, f = 1; char c = getchar();
for(; c<'0' or c>'9'; c=getchar())
if(c == '-') f = -1;
for(; '0'<=c and c<='9'; c=getchar())
a = (a<<3)+(a<<1)+(c^48);
return a*f;
}
inline void writeint(long long x){
if(x < 0) putchar('-'), x = -x;
if(x > 9) writeint(x/10);
putchar('0'+(x%10));
}
const int MaxN = 1000050;
int sa[MaxN], rnk[MaxN<<1], n;
char s[MaxN+5];
int bucket[MaxN], tmp[MaxN<<1];
# define x rnk
# define y tmp
void DeBug(){
printf("sa:");
for(int i=0; i<n; ++i)
printf(" %d",sa[i]);
printf("\nrank:");
for(int i=0; i<n; ++i)
printf(" %d",rnk[i]);
printf("\ntmp:");
for(int i=0; i<n; ++i)
printf(" %d",y[i]);
printf("\n");
}
void getSA(int m){
for(int i=0; i<=m; ++i)
bucket[i] = 0;
for(int i=0; i<n; ++i)
bucket[x[i]=s[i]] ++;
for(int i=n; i<(n<<1); ++i)
x[i] = y[i] = -1;
for(int i=1; i<=m; ++i)
bucket[i] += bucket[i-1];
for(int i=n-1; ~i; --i)
sa[--bucket[x[i]]] = i;
for(int w=1,p; w<n; w<<=1){
p = 0;
for(int i=n-w; i<n; ++i)
y[p++] = i;
for(int i=0; i<n; ++i)
if(sa[i] >= w)
y[p++] = sa[i]-w;
for(int i=0; i<=m; ++i)
bucket[i] = 0;
for(int i=0; i<n; ++i)
++ bucket[x[i]];
for(int i=1; i<=m; ++i)
bucket[i] += bucket[i-1];
for(int i=n-1; ~i; --i)
sa[--bucket[x[y[i]]]] = y[i];
swap(x,y);
x[sa[0]] = p = 0;
for(int i=1; i<n; ++i){
if(y[sa[i]] != y[sa[i-1]] or y[sa[i]+w] != y[sa[i-1]+w])
++ p;
x[sa[i]] = p;
}
if((m = p) == n-1) break;
}
for(int i=0; i<n; ++i)
rnk[sa[i]] = i;
}
# undef x
# undef y
int height[MaxN];
void getHeight(){
int k = 0;
for(int i=0,j; i<n; ++i){
if(k) -- k;
if(not rnk[i]){
height[0] = k = 0;
continue;
}
j = sa[rnk[i]-1];
while(i+k < n and j+k < n and s[j+k] == s[i+k])
++ k;
height[rnk[i]] = k;
}
}
int main(){
n = readint();
scanf("%s",s);
getSA(255);
getHeight();
long long ans = 0;
for(int i=0; i<n; ++i)
ans += n-sa[i]-height[i];
writeint(ans); putchar('\n');
return 0;
}
字符串问题
一道省选题中使用了后缀数组,尽管这并不是难点。
参考
- 大佬的博客与学弟的博客。说来惭愧,比不过学弟了!
- 洛谷题解。——在
例题
中有题目链接。 - 《 A C M / I C P C \tt ACM/ICPC ACM/ICPC 算法基础教程》——于瑞国主编,北京大学出版社
后文
听说常数更小的线性做法,我看不懂,但我大受震撼。