SA 复习笔记

大家好，由于蒟蒻 tzc 最近被动态点分治这个学也学不会的毒瘤玩意儿虐得不轻，所以就准备换换脑筋来 Van 同样学也学不会的后缀数组了。

考虑一个非常经典的问题：【模板】后缀排序。

一些定义（very important，不要搞混）：

\(sa_i\) 表示字典序排在第 \(i\) 的字符串是谁，也就是一个排名 \(\to\) 下标的映射
\(rk_i\) 表示后缀 \(i\)（（注：在本文中，后缀 \(i\) 就是开头位置为 \(i\) 的后缀）的排名是啥，也就是一个下标 \(\to\) 排名的映射

字符串比较是 \(\mathcal O(n)\) 的，所以暴力排序复杂度是 \(n^2\log n\) 的，显然不是我们所满足的。

考虑模拟暴力排序的过程，我们以这 \(n\) 个后缀的第一位为第一关键字，再以它们的第二位为第二关键字，再以它们的第三位为第三关键字……如果后缀长度不足 \(n\) 就在它后面补上 \(c_0\)，其中 \(c_0\) 为某个小于 \(s\) 中所有字符的字符，比如说 \(!\)。

不难发现整个过程可以分“阶段”进行，我们记第 \(k\) 阶段为比较完前 \(k\) 个字符后的结果，设 \(rk_i(k)\) 表示将所有后缀的长度为 \(k\) 的前缀排序后，后缀 \(i\) 的排名。如果有相同那么它们的 \(rk\) 也应该相同。

很明显暴力排序的过程就是每次从第 \(k\) 阶段转移到第 \(k+1\) 阶段，具体方法是对每个 \(i\) 将 \((rk_i(k),s_{i+k})\) 当作一个 pair 进行排序，如果 \(i+k>n\) 那么 \(s_{i+k}=0\)，即上文中所说的 \(c_0\)。

那么我们考虑，每次能不能不从第 \(k\) 阶段转移到第 \(k+1\) 阶段，而是从第 \(k\) 阶段转移到第 \(2k\) 阶段呢？

答案是肯定的。

我们从第 \(k\) 阶段转移到第 \(k+1\) 阶段的思路是：一个长度为 \(k+1\) 个字符串可以断成一个长度为 \(k\) 的子串和一个字符，而这个长度为 \(k\) 的字符串的大小关系我们是知道的，单个字符的大小关系我们也是知道的，所以可以转移。

同样道理一个长度为 \(2k\) 个字符串可以断成两个长度为 \(k\) 的子串，而这两个长度为 \(k\) 的字符串的大小关系我们也是知道的，所以也可以转移。

具体来说，你将 \((rk_i(k),rk_{i+k}(k))\) 当作一个 pair 进行排序，若 \(i+k>n\) 则 \(rk_{i+k}(k)=0\)，就可以得到 \(rk_{i}(2k)\)。

算下时间复杂度：倍增 1 log，sort n log，总复杂度 \(n\log^2n\)。

for(int i=1;i<=n;i++) x[i]=mp(mp(s[i],0),i);
sort(x+1,x+n+1);int gr=0;//stage 1
for(int i=1;i<=n;i++){
	if(x[i].fi!=x[i-1].fi) gr++;//note that if s[i]=s[j],rk(1,i)=rk(1,j), so we have to do a process similar to discretion
	rk[x[i].se]=gr; 
}
for(int k=1;k<n;k<<=1){//from stage k to stage 2k
	for(int i=1;i<=n;i++){
		if(i+k<=n) x[i]=mp(mp(rk[i],rk[i+k]),i);//sort the pairs of (rk(i,k),rk(i+k,k))
		else x[i]=mp(mp(rk[i],0),i);//i+k>n
	}
	sort(x+1,x+n+1);gr=0;
	for(int i=1;i<=n;i++){
		if(x[i].fi!=x[i-1].fi) gr++;
		rk[x[i].se]=gr; 
	}
}

然而这个程序是 TLE 90 的，这说明有比 \(n\log^2n\) 更优的算法。那么怎么进一步优化呢？

不难发现在我们所有 pair 中的两个数最大不过 \(n\)，于是我们回忆起一个东西叫做 ~~gay♂~~ 基数排序。

基数排序是一种能够在 \(\mathcal O(n+m)\) 的时间内对某个数组进行排序的算法，其中 \(m\) 为数组的值域。

感觉 SA 鸽了这么久就是因为不会这玩意儿啊 qwq。

先考虑一维的基数排序，我们及一个桶 \(c_i\) 表示值为 \(i\) 的 \(a_i\) 有多少个。然后对 \(c_i\) 做一遍前缀和。

显然，如果 \(a_i\) 互不相同排名为 \(c_{a_i}\) 的数就是 \(a_i\)，因为 \(c_{a_i}\) 就是比 \(a_i\) 小的个数。

那如果有 \(a_i\) 相同怎么办呢？那显然排名在 \((c_{a_i-1},c_{a_i}]\) 之间的数都等于 \(a_i\)。

考虑从后到前扫一遍，每一次将 \(a_i\) 的排名设为 \(c_{a_i}\) 并令 \(c_{a_i}\) 减 \(1\)，这样你发现所有相同的 \(a_i\) 都被赋上了一个不同的排名，并且对于相同的 \(a_i\)，越靠后的排名越靠后，也就是我们对 \((a_i,i)\) 这个 pair 排了一遍序。

回过头来，如果是二维的，例如在后缀排序中遇到的这些 pair 怎么办呢？

考虑我们在刚才的基数排序中为什么要从后到前扫，因为这样一来如果有 \(a_i\) 相同的那么靠后的排名永远比考前的排名高，也就是按照第二维降序的顺序更新桶。那这样一来就好办了，先将所有 \(y_i\) 排个序，然后按照 \(y_i\) 从大到小的顺序更新桶。就可以实现 \(\mathcal O(n)\) 排序了。

注意到如果有相同的 pair 那么它们的 \(rk\) 也应该相同，所以排完序后还要做一遍类似于离散化的操作。

const int M=123;//the ASCII of 'z' is 122
pii x[MAXN+5];
int n,sa[MAXN+5],rk[MAXN+5],buc[MAXN+5],seq[MAXN+5];
void get_sa(){
	for(int i=1;i<=n;i++) buc[s[i]]++;
	for(int i=1;i<=M;i++) buc[i]+=buc[i-1];
	for(int i=n;i;i--) sa[buc[s[i]]--]=i;
	int gr=0;
	for(int i=1;i<=n;i++){
		if(s[sa[i]]!=s[sa[i-1]]) gr++;
		rk[sa[i]]=gr;
	}//calculate the ranks of stage 1
	for(int k=1;k<n;k<<=1){//from stage k to stage 2k
		for(int i=1;i<=n;i++){
			if(i+k<=n) x[i]=mp(rk[i],rk[i+k]);//sort the pairs of (rk(i,k),rk(i+k,k))
			else x[i]=mp(rk[i],0);
		}
		memset(buc,0,sizeof(buc));//clear the bucket
		for(int i=1;i<=n;i++) buc[x[i].se]++;//sort the values of the second dimension
		for(int i=1;i<=n;i++) buc[i]+=buc[i-1];
		for(int i=n;i;i--) seq[buc[x[i].se]--]=i;//seq[i] means the i-smallest index sorted by the second dimension
		memset(buc,0,sizeof(buc));
		for(int i=1;i<=n;i++) buc[x[i].fi]++;//sort the values of the first dimension
		for(int i=1;i<=n;i++) buc[i]+=buc[i-1];
		for(int i=n;i;i--) sa[buc[x[seq[i]].fi]--]=seq[i];//update the bucket in descending order of the second dimension
		gr=0;
		for(int i=1;i<=n;i++){
			if(x[sa[i]]!=x[sa[i-1]]) gr++;//note that if x[i]=x[j],rk(i)=rk(j), so we have to do a process similar to discretion
			rk[sa[i]]=gr;
		}
	}
}

这样过是过了，但是我们发现最慢的一个点跑了整整 1s！显然常数上天了，我们考虑进行常数优化。

首先我们发现其实不用用桶第二维进行排序。对于 \(i=n-k+1\sim n\) 显然它们的第二维都是 \(0\)，应该排在最前面；对于 \(i\leq n-k\) 我们可以用之前求得的 \(k\) 阶段的排名算出，所以 2 遍循环就行了。

现在最慢的一个点跑了 750ms，大约省去了 \(1/4\) 的常数

void get_sa(){
	for(int i=1;i<=n;i++) buc[s[i]]++;
	for(int i=1;i<=M;i++) buc[i]+=buc[i-1];
	for(int i=n;i;i--) sa[buc[s[i]]--]=i;
	int gr=0;
	for(int i=1;i<=n;i++){
		if(s[sa[i]]!=s[sa[i-1]]) gr++;
		rk[sa[i]]=gr;
	}
//	for(int i=1;i<=n;i++) printf("%d%c",rk[i],(i==n)?'\n':' ');
	for(int k=1;k<=n;k<<=1){
		for(int i=1;i<=n;i++){
			if(i+k<=n) x[i]=mp(rk[i],rk[i+k]);
			else x[i]=mp(rk[i],0);
		}
		int num=0;
		for(int i=n-k+1;i<=n;i++) seq[++num]=i;
		for(int i=1;i<=n;i++) if(sa[i]>k) seq[++num]=sa[i]-k;
		memset(buc,0,sizeof(buc));
		for(int i=1;i<=n;i++) buc[x[i].fi]++;
		for(int i=1;i<=n;i++) buc[i]+=buc[i-1];
		for(int i=n;i;i--) sa[buc[x[seq[i]].fi]--]=seq[i];
		gr=0;
		for(int i=1;i<=n;i++){
			if(x[sa[i]]!=x[sa[i-1]]) gr++;
			rk[sa[i]]=gr;
		}
	}
}

继续卡！我们发现，并不是每次的值域都能达到 \(n\) 那么大，所以可以记录一个 \(vmax\) 表示当前阶段的值域，更新桶的时候只用对 \([1,vmax]\) 做前缀和。

现在最慢的一个点跑了 600ms：

void get_sa(){
	int vmax=123;
	for(int i=1;i<=n;i++) buc[s[i]]++;
	for(int i=1;i<=vmax;i++) buc[i]+=buc[i-1];
	for(int i=n;i;i--) sa[buc[s[i]]--]=i;
	int gr=0;
	for(int i=1;i<=n;i++){
		if(s[sa[i]]!=s[sa[i-1]]) gr++;
		rk[sa[i]]=gr;
	} vmax=gr;
//	for(int i=1;i<=n;i++) printf("%d%c",rk[i],(i==n)?'\n':' ');
	for(int k=1;k<=n;k<<=1){
		for(int i=1;i<=n;i++){
			if(i+k<=n) x[i]=mp(rk[i],rk[i+k]);
			else x[i]=mp(rk[i],0);
		}
		int num=0;
		for(int i=n-k+1;i<=n;i++) seq[++num]=i;
		for(int i=1;i<=n;i++) if(sa[i]>k) seq[++num]=sa[i]-k;
		memset(buc,0,sizeof(buc));
		for(int i=1;i<=n;i++) buc[x[i].fi]++;
		for(int i=1;i<=vmax;i++) buc[i]+=buc[i-1];
		for(int i=n;i;i--) sa[buc[x[seq[i]].fi]--]=seq[i];
		gr=0;
		for(int i=1;i<=n;i++){
			if(x[sa[i]]!=x[sa[i-1]]) gr++;
			rk[sa[i]]=gr;
		}
		vmax=gr;
	}
}

我们还注意到，如果 \(n\) 个字符串都不相同了就可以 break 了，所以可以再加一句 if(vmax==n) break;

发现一句的效果比之前好得多，现在最慢的一个点只跑了 200ms 了：

void get_sa(){
	int vmax=123;
	for(int i=1;i<=n;i++) buc[s[i]]++;
	for(int i=1;i<=vmax;i++) buc[i]+=buc[i-1];
	for(int i=n;i;i--) sa[buc[s[i]]--]=i;
	int gr=0;
	for(int i=1;i<=n;i++){
		if(s[sa[i]]!=s[sa[i-1]]) gr++;
		rk[sa[i]]=gr;
	} vmax=gr;
//	for(int i=1;i<=n;i++) printf("%d%c",rk[i],(i==n)?'\n':' ');
	for(int k=1;k<=n;k<<=1){
		for(int i=1;i<=n;i++){
			if(i+k<=n) x[i]=mp(rk[i],rk[i+k]);
			else x[i]=mp(rk[i],0);
		}
		int num=0;
		for(int i=n-k+1;i<=n;i++) seq[++num]=i;
		for(int i=1;i<=n;i++) if(sa[i]>k) seq[++num]=sa[i]-k;
		memset(buc,0,sizeof(buc));
		for(int i=1;i<=n;i++) buc[x[i].fi]++;
		for(int i=1;i<=vmax;i++) buc[i]+=buc[i-1];
		for(int i=n;i;i--) sa[buc[x[seq[i]].fi]--]=seq[i];
		gr=0;
		for(int i=1;i<=n;i++){
			if(x[sa[i]]!=x[sa[i-1]]) gr++;
			rk[sa[i]]=gr;
		}
		vmax=gr;if(vmax==n) break;
	}
}

至此我们有了常数较小的 SA 的实现。

那么 SA 有什么用呢？SA 一个用途就是与 LCP 结合。

LCP，即 Longest Common Prefix，指两串的最长公共前缀。

我们记 \(LCP(i,j)\) 为后缀 \(sa_i\) 与 \(sa_j\) 的最长公共前缀，注意这边是 \(sa_i\) 与 \(sa_j\) 的最长公共前缀，即排名为 \(i\) 的串与排名为 \(j\) 的串的最长公共前缀。

那么 LCP 有哪些性质呢？

\(LCP(i,j)=LCP(j,i)\)
\(LCP(i,i)=n-sa_i+1\)
\(LCP(i,j)=\min(LCP(i,k),LCP(k,j)) (i<k<j)\)，该定理我们称之为 LCP Lemma
\(LCP(i,j)=\min\limits_{k=i+1}^{j}LCP(k-1,k)(i<j)\)，该定理我们称之为 LCP Theorem

首先结论 \(1,2\) 肯定是显然的，用脚趾头都能想出来。

结论 \(3\) 和 \(4\) 显然是同一类的，证出结论 \(3\) 就能顺带着证出结论 \(4\)。

考虑如何证明结论 \(3\)，设 \(u\) 为后缀 \(sa_i\)，\(v\) 为后缀 \(sa_j\)，\(w\) 为后缀 \(sa_k\)，\(x=\min(LCP(i,k),LCP(k,j))\)

由于 \(x=\min(LCP(i,k),LCP(k,j))\)，\(u,v,w\) 的前 \(x\) 位应当都是相同的，故 \(LCP(i,j)\geq x\)。

我们要证明 \(LCP(i,j)=x\)，即 \(LCP(i,j)>x\) 不可能实现。

考虑反证法，假设 \(LCP(i,j)>x\)，则 \(u\) 的第 \(x+1\) 位等于 \(v\) 的第 \(x+1\) 位。

由于 \(i<k<j\)，故 \(u,w,v\) 在字典序上依次递增，而 \(u,v,w\) 的前 \(x\) 位都相等，故 \(u_{x+1}\leq w_{x+1}\leq v_{x+1}\)。

而 \(x=\min(LCP(i,k),LCP(k,j))\) 说明 \(u_{x+1},v_{x+1},w_{x+1}\) 不可能全相等，就说明 \(u\) 的 \(x+1\) 位不可能等于 \(v\) 的 \(x+1\) 位。

故 \(LCP(i,j)=x\)。

而结论 4 其实是结论 3 的变形。

如果我们记 \(ht_i=LCP(i-1,i)\)，那么 \(LCP(i,j)=\min\limits_{k=i+1}^jht_i\)，这东西显然可以 RMinQ 维护。

现在的问题就转化为怎么求 \(ht_i\)。暴力吗？一氧化氮。

如果我们设 \(h_i\) 为 \(ht_{rk_i}\)，即后缀 \(i\) 与后缀 \(i\) 字典序前一位的 LCP，那么有一个很重要的公式叫做 \(h_i\geq h_{i-1}-1\)

怎么证明这个公式呢？

考虑设 \(j=sa_{rk_i-1}\)，\(k=sa_{rk_{i-1}-1}\)，分两种情况讨论：

若后缀 \(k\) 与后缀 \(i-1\) 第一位不同，则 \(h_{i-1}=0\)，\(h_i\geq -1\) 显然成立。
若后缀 \(k\) 与后缀 \(i-1\) 第一位相同，则刨掉第一位还剩下后缀 \(k+1\) 和后缀 \(i\)，它们的 \(LCP\) 为 \(h_{i-1}-1\)，而由于后缀 \(k\) 字典序排在后缀 \(i-1\) 之前，后缀 \(k\) 与后缀 \(i-1\) 第一位相同，故后面的后缀 \(k+1\) 的字典序小于后缀 \(i\) 的字典序。根据 LCP Theorem，后缀 \(k+1\) 与后缀 \(i\) 的 LCP 为 \(\min\limits_{t=rk_{k+1}+1}^{rk_i}ht_t=h_{i-1}-1\)，故 \(ht_{rk_i}\geq h_{i-1}-1\)，故 \(h_i\geq h_{i-1}-1\)

附图：

有了这个式子，就可以用类似 two pointers 的东西求出 \(h_i\)，也就顺带着求出 \(ht\) 了。

void getht(){
	int k=0;
	for(int i=1;i<=n;i++){//calculate h[i]
		if(rk[i]==1) continue;
		if(k) k--;//h[i]>=h[i-1]-1
		int j=sa[rk[i]-1];
		while(s[i+k]==s[j+k]) k++;//violently match suffix i and suffix sa[rk[i]-1]
		ht[rk[i]]=k;//use formula ht[rk[i]]=h[i] to calculate ht[i]
	}
}

后缀数组求一个子串在原串中的出现次数：

由于这个东西应用比较广泛就在这里一并讲掉了罢。

倘若我们要求 \(s\) 的子串 \(s[l...r]\) 在 \(s\) 中的出现次数，那么等价于求有多少个后缀 \(j\) 满足 \(\text{LCP}(s[j...n],s[l...n])\geq r-l+1\)，因为如果 \(\text{LCP}(s[j...n],s[l...n])\geq r-l+1\)，那么一定有 \(s[j...j+r-l]=s[l...r]\)，也就对 \(s[l...r]\) 在 \(s\) 中的出现次数产生了 \(1\) 的贡献，反之一定有 \(s[j...j+r-l]\neq s[l...r]\)。

于是考虑怎样求有多少个后缀 \(j\) 满足 \(\text{LCP}(s[j...n],s[l...n])\geq r-l+1\)。考虑将这 \(n\) 个后缀放在字典序数组 \(rk_i\) 上。显然对于所有 \(\text{LCP}(s[j...n],s[l...n])\geq r-l+1\)，它们的 \(rk\) 值应当是一段连续的区间 \([L,R]\)。而这个区间又可以通过二分+ST 表求出。具体来说，二分找出最小的满足 \(\min\limits_{i=t+1}^{rk_l}ht_i\ge r-l+1\) 的 \(t\) 即为该区间的左端点 \(L\)；二分找出最大的满足 \(\min\limits_{i=rk_l+1}^{t}ht_i\ge r-l+1\) 的 \(t\) 即为该区间的左端点 \(R\)。然后返回 \(R-L+1\) 即可，正确性显然。

时间复杂度 \(\mathcal O(\log n)\)。

总结

后缀数组的知识点就这么多，学好后缀数组首先需掌握后缀数组的求法以及它与 LCP 相结合的一些定理的应用。

后缀数组最常用的场景是处理与子串有关的问题，因为后缀的前缀即是子串，因此看到“子串”等字眼就要尽量往后缀数组上想。

后缀数组的难题一般难点都不在怎样求后缀数组 \(sa_i\) 及 \(ht_i\)，而是在于求完 \(sa_i\) 和 \(ht_i\) 之后怎样求出答案，由于 LCP Theorem 涉及到区间取 \(\min\) 的问题，所以这部分题也常与数据结构结合，需要较强的数据结构基础。

posted @ 2021-01-10 11:02 tzc_wk 阅读(179) 评论(1) 编辑收藏举报

刷新页面返回顶部

tzc_wk

SA 复习笔记

总结

公告