Suffix Array

简介

后缀数组（或者叫后缀排序算法），能够将一个字符串每一个后缀按照字典序排序。

暴力的 \(\mathcal{O}(n^{2}\log n)\) 不用细讲，使用哈希优化后的 \(\mathcal{O}(n\log^{2}n)\) 也不讲。

\(\mathcal{O}(n \log^{2} n)\) 做法

一些定义：

\(sa_{i}\) 表示后缀排序后，排名为 \(i\) 的后缀起点。
\(rk_{i}\) 表示起点为 \(i\) 的后缀排序后的排名。

显然有 \(sa_{rk_{i}} = i\)。

暴力做法满的地方在于，如果直接按照字典序比较的方式，一位一位地进行比较，其单次比较的时间复杂度是 \(\mathcal{O}(|S|)\)（\(S\) 为字符串）的。

如何优化这个比较过程？答案是使用倍增，我也不知道倍增这个方法是怎么 yy 出来的，但是不难理解就不去深究了。

假设我们要比较两个长度为 \(L\) 的字符串大小（如果有一个字符串长度不足就在后面补上空字符），那么可以分成两部分来比较，前面一部分长 \(\left\lfloor\dfrac{L}{2}\right\rfloor\) 后一部分长 \(\left\lceil\dfrac{L}{2}\right\rceil\)（长度反过来一样可以）。

先比较长为 \(\left\lfloor\dfrac{L}{2}\right\rfloor\) 的前半部分，再比较长为 \(\left\lceil\dfrac{L}{2}\right\rceil\) 的后半部分即可，前半部分和后半部分互不影响，只是前半部分的优先级高于后半部分而已，这就是倍增处理后缀数组的核心思想。

首先我们只看每个后缀的首字母来排序，初步得到一个 \(rk\) 数组。

接下来假设我们已经处理好了每个后缀按其长为 \(L\) 的前缀来排序的后缀数组。

对于一个长为 \(2L\) 的子串（一个子串是某个后缀的前缀），可以按照上述方式分成两个长为 \(L\) 的部分进行比较。

而在比较前半部分时，我们不需要重新比较一遍，因为我们已经处理好了每个长为 \(L\) 的子串的排名（临时储存在 \(rk\) 数组中），只需查询其排名就可以知道它们的大小关系。后半部分同理。

初始时我们已经处理好了 \(L = 1\) 的情况了，接下来的过程每次 \(L\) 会翻倍，所以倍增过程是 \(\mathcal{O}(n \log n)\) 的，再算上快速排序就是 \(\mathcal{O}(n \log^{2} n)\) 的。

注意每轮倍增结束后要更新一遍 \(rk\)。

代码实现：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int n, sa[1000005], rk[2000005], ork[2000005];
string s;
int main() {
	ios::sync_with_stdio(0);
	cin.tie(0), cout.tie(0);
	cin >> s;
	n = s.length();
	s = "V" + s;
	for(int i = 1; i <= n; ++i) sa[i] = i, rk[i] = s[i];
	for(int w = 1; w < n; w <<= 1) {
		stable_sort(sa + 1, sa + 1 + n, [&w](int p1, int p2) -> bool {return make_pair(rk[p1], rk[p1 + w]) < make_pair(rk[p2], rk[p2 + w]);});
		copy(rk + 1, rk + 1 + n, ork + 1);
		for(int i = 1, now = 0; i <= n; ++i) {
			if(ork[sa[i]] != ork[sa[i - 1]] || ork[sa[i] + w] != ork[sa[i - 1] + w]) ++now;
			rk[sa[i]] = now;
		}
	}
	for(int i = 1; i <= n; ++i) cout << sa[i] << " ";
	return 0;
}

\(\mathcal{O}(n \log n)\) 做法

上面的后缀排序是 \(\mathcal{O}(n \log^{2} n)\) 的，但是它的常数实在是比哈希大，只有降低复杂度才有可能让它有不被哈希替代的充分的理由。

注意到快速排序的过程中，值域是 \(\mathcal{O}(n)\) 的，并且是一个双关键字排序，所以可以使用基数排序和计数排序进行优化，将复杂度降至 \(\mathcal{O}(n \log n)\)。

基数排序和计数排序不讲了。

代码就暂时不放了。

常数优化

为什么不放代码？因为它的常数还是太大了。

如果就拿没有卡过常的代码交题的话， @Super_Cube 会拿着他的小常数哈希来嘲讽你，这会让你非常难受，同时在数据范围稍大的题中（为了卡双 \(\log\) 做法之类的）你也有可能被卡常。

（就像多项式牛顿迭代的 \(\mathcal{O}(n \log n)\) 卡不过半在线卷积的 \(\mathcal{O}(n \log^{2} n)\)，悲）

第二关键字无需计数排序

~~照抄 OI wiki 吗？除了这个小标题我想不到别的了。~~

基数排序中我们会先按第二关键字进行排序，这个时候其实可以特殊讨论一下来减少常数。

我们会发现第二关键字排序过程中有两类元素，一类没有超出范围的，它们会按照 \(rk\) 进行排序，一类超出范围的，它们会被统一当作空字符处理。

这个时候先把超出范围的部分放到前面，再把没有超出范围的正常排序接到后面即可。

注意到我们在处理这一部分时是使用的已经处理好的 \(rk\) 和 \(sa\)，此时无需再进行多余的排序，直接按照原顺序接到后面即可。

优化计数排序值域范围

每次更新完 \(rk\) 时会遗留下一个变量 \(now\)，它是 \(rk\) 的最大值也是 \(rk\) 的范围，每次将值域改为 \(now\) 即可。

排名均不相同时算法可直接结束

就像冒泡排序在数组有序时可以退出一样，所有排名均不相同时代表着每个后缀长为 \(L\) 的前缀均不相同，此时再进行后续的比较也不会有改变了（第一关键字的优先级高于第二关键字。）

减少不连续的内存访问

点名 st 表，倍增求树上 LCA 被卡常的时候。

特别是数据范围比较大时，不连续的内存访问真的能要命。

可以将排序过程中会重复用到的，或者可能会造成严重的内存不连续访问的变量存下来。

优化后的代码：

有点轻微压行，不过我觉得这样压了之后还挺便于记忆的。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int n, sa[1000005], rk[2000005], ork[2000005], id[1000005], ky[1000005], cnt[1000005];
string s;
int main() {
	ios::sync_with_stdio(0);
	cin.tie(0), cout.tie(0);
	cin >> s;
	n = s.length();
	s = "V" + s;
	for(int i = 1; i <= n; ++i) ++cnt[rk[i] = s[i]];
	for(int i = 1; i <= 256; ++i) cnt[i] += cnt[i - 1];
	for(int i = n; i >= 1; --i) sa[cnt[rk[i]]--] = i;
	for(int w = 1, limit = max(256, n + 1), now = 0; limit != n; w <<= 1, limit = now, now = 0) {
		for(int i = n; i > n - w; --i) id[++now] = i;
		for(int i = 1; i <= n; ++i) {
			if(sa[i] > w) id[++now] = sa[i] - w;
		}
		fill(cnt, cnt + 1 + limit, 0);
		for(int i = 1; i <= n; ++i) ++cnt[ky[i] = rk[id[i]]];
		for(int i = 1; i <= limit; ++i) cnt[i] += cnt[i - 1];
		for(int i = n; i >= 1; --i) sa[cnt[ky[i]]--] = id[i];
		copy(rk + 1, rk + 1 + n, ork + 1);
		now = 0;
		for(int i = 1; i <= n; ++i) {
			if(ork[sa[i]] != ork[sa[i - 1]] || ork[sa[i] + w] != ork[sa[i - 1] + w]) ++now;
			rk[sa[i]] = now;
		}
	}
	for(int i = 1; i <= n; ++i) cout << sa[i] << " ";
	return 0;
}

（说是便于记忆是指中间排序的三行，那里的变量使用有点混乱，但是三行长度相同这一点可以帮助检查和记忆。码风不一样的话就别管了。）

这份代码在 ber 中机子上跑出的战绩如下：范围较小时耗时为 @Super_Cube 的小常数哈希的 \(\dfrac{1}{2}\) 左右，范围较大时耗时在 @Super_Cube 的小常数哈希的 \(\dfrac{1}{5}\) 左右。

模板题除外，只快了 100ms。

反正不会被 @Super_Cube 嘲讽了就是

\(height\) 数组

SA 还有 \(O(n)\) 做法？学它干嘛，反正小常数单 \(\log\) 够用了，而且大部分题目除了后缀数组的部分还是有 \(\log\) 的。

接下来就是 \(height\) 数组的介绍和求法了，大部分时候后缀数组的题目都会用到 \(height\) 数组。

求法

令 \(\operatorname{LCP}(i, j)\) 表示以 \(i\) 为起点的后缀和以 \(j\) 为起点的后缀的最长公共前缀。

那么有 \(height_{i} = \operatorname{LCP}(sa_{i - 1}, sa_{i})\)。

哈希是双 \(\log\) 的，我要单 \(\log\) 甚至线性怎么办？

Lemma：\(height_{rk_{i}} \geqslant height_{rk_{i - 1}} - 1\)。

证明如下：

\(height_{rk_{i - 1}} \leqslant 1\) 时，结论显然成立。
\(height_{rk_{i - 1}} \geqslant 2\) 时，有：

\(height_{rk_{i - 1}} = \operatorname{LCP}(sa_{rk_{i - 1}}, sa_{rk_{i - 1} - 1}) = \operatorname{LCP}(i - 1, sa_{rk_{i - 1} - 1})\)。（\(sa_{rk_{i}} = i\) 在开头提到了。）

这说明分别从 \(i - 1\)，\(sa_{rk_{i - 1} - 1}\) 开始，它们有长度大于等于 \(2\) 的公共前缀，那么可以令以 \(i - 1\) 开始的后缀为 \(xAB\)，以 \(sa_{rk_{i - 1} - 1}\) 开始的后缀为 \(xAC\)（\(x\) 为一个字符，\(A, B, C\) 都为字符串，\(A\) 非空）。

而 \(height_{rk_{i}} = \operatorname{LCP}(sa_{rk_{i}}, sa_{rk_{i} - 1}) = \operatorname{LCP}(i, sa_{rk_{i} - 1})\)。

根据上面的约定，从 \(i\) 开始的后缀就为 \(AB\)。

注意到 \(rk_{i} - 1 < rk_{i}\)，所以以 \(sa_{rk_{i} - 1}\) 为起点的后缀是小于以 \(sa_{rk_{i}} = i\) 为起点的后缀的，并且恰好小一位。

所以以 \(sa_{rk_{i} - 1}\) 为起点的后缀小于 \(AB\)，并且 \(AC < AB\)。

因为 \(sa_{rk_{i} - 1}\) 比 \(AB\) 恰好小一位，所以 \(sa_{rk_{i} - 1} \geqslant AC\)，这一点是比较难理解的，可以类比：\(a - b = 1, c < a \Rightarrow c \leqslant b < a (a, b, c \in \mathbb{Z})\)。

因为 \(AC \leqslant \operatorname{Suffix}(sa_{rk_{i} - 1}) < AB\)，所以后缀 \(i\) 和后缀 \(sa_{rk_{i} - 1}\) 至少都有 \(A\) 这个前缀，也就是 \(height_{rk_{i}} \geqslant |A| = height_{rk_{i - 1}} - 1\)。

有了这个引理之后就很好做了，每次将 \(height_{rk_{i}}\) 设为 \(height_{rk_{i - 1}} - 1\)，在重复尝试扩大，均摊下来是 \(\mathcal{O}(n)\) 的（最多 \(n\) 次减法和 \(2n\) 次加法）。

代码：

for(int i = 1, k = 0; i <= n; ++i) {
	if(!rk[i]) continue;
	k -= (bool)k;
	while(s[i + k] == s[sa[rk[i] - 1] + k]) ++k;
	height[rk[i]] = k;
}

作用

大部分时候是来求任意两个后缀的 \(\operatorname{LCP}\) 的。

有 \(\operatorname{LCP}(i, j) = \operatorname{LCP}(rk_{i}, rk_{i} + 1, \cdots, rk_{j}) = \min\limits_{rk_{i} < k \leqslant rk_{j}}height_{k}(rk_{i} < rk_{j})\)。

感性理解的话，就是后缀数组中相邻的两个后缀的 \(\operatorname{LCP}\) 一定比不相邻的长，任意多个连续后缀之间的 \(height\) 就代表着它们的 \(\operatorname{LCP}\)，可以利用相等关系的传递性加深理解。

例题

鸽子薇薇当然会把例题鸽到最后一点一点来写啦。

B. 优秀的拆分

求连续相同子串的经典题目。

拆分的本质是把两个 \(AA\) 式的子串拼在一起，所以对于前面和后面的两个 \(AA\) 式子串分别求出其个数在乘上求和即可。

如何求一个 \(AA\) 式子串？考虑枚举 \(|A|\)，但是显然不能把整个字符串遍历一遍，尝试隔 \(|A|\) 个字符枚举一次（或者说以 \(|A|\) 为块长将字符串分块，取每一块的块首）。

任意一个 \(AA\) 式子串一定会跨过两个选取点的。

再求出以每两个相邻选取点为起点/终点的最长公共后缀/最长公共前缀，这里可以用后缀数组做。

令最长公共前缀的长度为 \(lcp\)，最长公共后缀的长度为 \(lcs\)，如果 \(lcp + lcs < |A|\) 那一定不可能有 \(AA\) 式子串跨过这两个点。

否则 \(lcp\) 和 \(lcs\) 将有一部分重叠，令重叠长度为 \(tmp\)，如下图：

\[\underbrace{\fcolorbox{000000}{33667f}{\kern{20pt}}}_{lcs}\fcolorbox{000000}{66ccff}{\kern{20pt}}\overbrace{\fcolorbox{000000}{19333f}{\kern{16pt}}}^{tmp}\fcolorbox{000000}{33667f}{\kern{10pt}}\underbrace{\fcolorbox{000000}{66ccff}{\kern{38pt}}}_{lcp} \]

手玩一下可得：在这一段区间中，一共会有 \(tmp\) 个 \(AA\) 式子串出现，并且出现的位置是连续的。

那么用一个差分数组记录一下出现的位置即可，最后统计答案。

Code：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define lg __lg
struct Suffix_Array {
	string s;
	int n, sa[30005], rk[60005], ork[60005], id[30005], ky[30005], cnt[30005], height[30005], st[16][30005];
	void build(const string& str) {
		s = "V" + str + "V";
		n = str.length();
		rk[n + 1] = 0, ork[n + 1] = 0;
		fill(cnt + 1, cnt + 256, 0);
		for(int i = 1; i <= n; ++i) ++cnt[rk[i] = s[i]];
		for(int i = 1; i <= 256; ++i) cnt[i] += cnt[i - 1];
		for(int i = n; i >= 1; --i) sa[cnt[s[i]]--] = i;
		for(int w = 1, limit = max(256, n + 1), now = 0; limit != n; w <<= 1, limit = now, now = 0) {
			for(int i = n; i > n - w; --i) id[++now] = i;
			for(int i = 1; i <= n; ++i) {
				if(sa[i] > w) id[++now] = sa[i] - w;
			}
			fill(cnt + 1, cnt + limit + 1, 0);
			for(int i = 1; i <= n; ++i) ++cnt[ky[i] = rk[id[i]]];
			for(int i = 1; i <= limit; ++i) cnt[i] += cnt[i - 1];
			for(int i = n; i >= 1; --i) sa[cnt[ky[i]]--] = id[i];
			copy(rk + 1, rk + 1 + n, ork + 1);
			now = 0;
			for(int i = 1; i <= n; ++i) {
				if(ork[sa[i]] != ork[sa[i - 1]] || ork[sa[i] + w] != ork[sa[i - 1] + w]) ++now;
				rk[sa[i]] = now;
			}
		}
		for(int i = 1, k = 0; i <= n; ++i) {
			height[rk[i]] = 0;
			if(!rk[i]) continue;
			k -= (bool)k;
			while(s[i + k] == s[sa[rk[i] - 1] + k]) ++k;
			st[0][rk[i]] = height[rk[i]] = k;
		}
		for(int i = 1; i <= lg(n); ++i) {
			for(int j = 1; j <= n - (1 << i) + 1; ++j) {
				st[i][j] = min(st[i - 1][j], st[i - 1][j + (1 << (i - 1))]);
			}
		}
	}
	#define k lg(r - l + 1)
	int ask(int l, int r) {
		l = rk[l], r = rk[r];
		if(l > r) swap(l, r);
		++l;
		return min(st[k][l], st[k][r - (1 << k) + 1]);
	}
	#undef k
} sa, rsa;
int t, n, lcp, lcs, l[30005], r[30005];
ll ans;
string s;
int main() {
	ios::sync_with_stdio(0);
	cin.tie(0), cout.tie(0);
	cin >> t;
	while(t--) {
		cin >> s;
		n = s.length();
		sa.build(s);
		reverse(s.begin(), s.end());
		rsa.build(s);
		ans = 0;
		fill(l, l + 1 + n, 0);
		fill(r, r + 1 + n, 0);
		for(int len = 1; len <= (n >> 1); ++len) {
			for(int p1 = len, p2 = (len << 1); p2 <= n; p1 += len, p2 += len) {
				lcp = min(sa.ask(p1, p2), len), lcs = min(rsa.ask(n - p1 + 2, n - p2 + 2), len - 1);
				if(lcp + lcs >= len) {
					++r[p1 - lcs], --r[p1 + lcp - len + 1];
					++l[p2 - lcs + len - 1], --l[p2 + lcp];
				}
			}
		}
		for(int i = 1; i <= n; ++i) l[i] += l[i - 1], r[i] += r[i - 1];
		for(int i = 2; i <= n; ++i) ans += (ll)l[i - 1] * r[i];
		cout << ans << '\n';
	}
	return 0;
}

posted @ 2024-04-11 20:10 A_box_of_yogurt 阅读(12) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

一只下饭菜的blog

简简单单，可可爱爱