【学习笔记】后缀数组SA

38 篇文章 0 订阅
13 篇文章 0 订阅

前言

建议先学习 基数排序 。本文讲解倍增法。想学奇技淫巧 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,,sn1) 。称 ( a 0 , a 1 , a 2 , … , a m − 1 ) (a_0,a_1,a_2,\dots,a_{m-1}) (a0,a1,a2,,am1) s s s 的后缀,当且仅当 m ≤ n m\le n mn ∀ 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+nm

更简单一点:后缀是包含最后一个字符的子串。

在后面的叙述中,将以第 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,,sn1),称为:后缀 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+w1) 的大小,就是 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+2wn )的呢?认为其第二关键字是 − ∞ -\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+w1>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(x1)] 。其中 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+1minxheight(i)height(x)=h(x)

如果是这样的话,假设后缀 ( x − 1 ) (x-1) (x1) 与后缀 ( y − 1 ) (y-1) (y1) 有一个 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(x1)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;
	}
}

例题

洛谷的板题

传送门 to luogu

为了让你有地方找题解——题解中也讲得很好。

#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 xyω(也就是错开了至少 ω \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(x1) 的前缀,一定是已经被计算过的。而且,这也是最长的一个 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=0n1height(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 算法基础教程》——于瑞国主编,北京大学出版社

后文

听说常数更小的线性做法,我看不懂,但我大受震撼。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值