后缀数组(SA)习记
前言
在学习 \(SA\) 之前,有必要复习一下什么是 计数排序、基数排序、桶排序 以及三者之间的区别。
然后细说一下 后缀数组 的三种构造方式 倍增、DC3、SA_IS
部分代码和内容转自:
如侵权可联系删除。
目录:
计数排序
计数排序(Counting sort)是一种线性时间的排序算法。
算法流程
- 计算每个数出现了几次;
- 求出每个数出现次数的 前缀和;
- 利用出现次数的前缀和,从右至左计算每个数的排名。
算法分析
- 是一种稳定的排序算法。
- 时间复杂度为 \(O(n + w)\), w 为值域大小
缺点:
- \(O(w)>O(n*log(n))\) 时,不如 \(O(nlogn)\) 的排序算法。
代码实现
#include<bits/stdc++.h>
using namespace std;
const int N = 100010;
const int W = 100010;
int n, w, a[N], cnt[W], b[N];
void counting_sort() {
memset(cnt, 0, sizeof(cnt));
for (int i = 1; i <= n; ++i) ++cnt[a[i]];
for (int i = 1; i <= w; ++i) cnt[i] += cnt[i - 1];
// 倒序是因为要保证原数组里面相同项相对位置不变,原来在后面的还在后面
for (int i = n; i >= 1; --i) b[cnt[a[i]]--] = a[i];
}
int main(){
cin >> n;
for (int i = 1; i <= n; i++) {
cin >> a[i];
w = max(w, a[i]);
}
counting_sort();
for (int i = 1; i <= n; i++) {
cout << b[i] << " ";
}
return 0;
}
基数排序
基数排序(Radix sort)是一种非比较型的排序算法,最早用于解决卡片排序的问题
算法流程
将待排序的元素拆分为 \(k\) 个关键字。先对第 \(k\) 个关键字进行稳定排序,然后对第 \(k-1\) 个关键字稳定排序,直到对第 \(1\) 个关键字排序完成。
- 基数排序主要是一种思想,对某个关键字的内部排序是依靠其他排序算法来完成。如计数排序。
算法分析
- 是一种稳定的排序算法。
- 时间复杂度: 一般来说,如果每个关键字的值域都不大,就可以使用 计数排序 作为内层排序,此时的复杂度为 \(O(k*n + \Sigma_{i=1}{k}w_i)\) ,其中 \(w_i\) 为第 \(i\) 关键字的值域大小。如果关键字值域很大,就可以直接使用基于比较的 \(O(k*n*logn)\) 排序而无需使用基数排序了。
- 空间复杂度: \(O(k + n)\);
代码实现
#include <bits/stdc++.h>
using namespace std;
const int N = 100010;
const int W = 100010;
const int K = 100;
int n, w[K], k, cnt[W];
struct Element {
int key[K];
bool operator<(const Element &y) const {
// 两个元素的比较流程
for (int i = 1; i <= k; ++i) {
if (key[i] == y.key[i])
continue;
return key[i] < y.key[i];
}
return false;
}
} a[N], b[N];
void counting_sort(int p) { // 内层计数排序
memset(cnt, 0, sizeof(cnt));
for (int i = 1; i <= n; ++i)
++cnt[a[i].key[p]];
for (int i = 1; i <= w[p]; ++i)
cnt[i] += cnt[i - 1];
// 为保证排序的稳定性,此处循环i应从n到1
// 即当两元素关键字的值相同时,原先排在后面的元素在排序后仍应排在后面
for (int i = n; i >= 1; --i)
b[cnt[a[i].key[p]]--] = a[i];
memcpy(a, b, sizeof(a));
}
void radix_sort() {
for (int i = k; i >= 1; --i) {
// 借助计数排序完成对关键字的排序
counting_sort(i);
}
}
桶排序
桶排序(Bucket sort)是排序算法的一种,适用于待排序数据值域较大但分布比较均匀的情况。
算法流程
- 将值域分块,设置一个定量的数组当作空桶,一个捅对应一个块的元素
- 遍历序列,并将每个块中元素一个个放到对应的桶中;
- 对每个不是空的桶进行排序,通常是插入排序;
- 从不是空的桶里把元素再放回原来的序列中。
算法分析
- 稳定性:
- 如果使用稳定的内层排序,并且将元素插入桶中时不改变元素间的相对顺序,那么桶排序就是一种稳定的排序算法。
- 由于每块元素不多,一般使用插入排序。此时桶排序是一种稳定的排序算法。
- 时间复杂度:
- 桶排序的平均时间复杂度为 \(O(n + n^2/k + k)\) (将值域平均分成 \(n\) 块 + 排序 + 重新合并元素),当 \(k\approx n\) 时为 \(O(n)\).
- 桶排序的最坏时间复杂度为 。
代码实现
#include<bits/stdc++.h>
const int N = 100010;
int n, w, a[N];
vector<int> bucket[N];
void insertion_sort(vector<int> & A) {
for (int i = 1; i < A.size(); ++i)
{
int key = A[i];
int j = i - 1;
while (j >= 0 && A[j] > key)
{
A[j + 1] = A[j];
--j;
}
A[j + 1] = key;
}
}
void bucket_sort() {
int bucket_size = w / n + 1;
for (int i = 0; i < n; ++i) {
bucket[i].clear();
}
for (int i = 1; i <= n; ++i) {
bucket[a[i] / bucket_size].push_back(a[i]);
}
int p = 0;
for (int i = 0; i < n; ++i) {
insertion_sort(bucket[i]);
for (int j = 0; j < bucket[i].size(); ++j) {
a[++p] = bucket[i][j];
}
}
}
SA基本概念及性质
SA基本概念
- 后缀:\(S[i]=S[i, |S|]\)
- 字典序:从左往右找两个字符串第一个不同字母,空字符设为最小。
- 后缀排序:将所有后缀 \(S[i]\) 看作独立的串,放在一起按照字典序进行升序排序。
- 后缀排名 \(rk[i]\):\(rk[i]\) 表示后缀 \(S[i]\) 在后缀排序中的排名,即他是第几小的后缀。
- 后缀数组 \(sa[i]\):\(sa[i]\) 表示排名第 \(i\) 小的后缀。
- LCP: Longest Common Prefix, 最长公共前缀。
一个重要的等式:rk[sa[i]] = sa[rk[i]] i
。
Naive 求法。
对所有后缀字符串进行 std::sort
,用 哈希 + 二分 重写 cmp()
比较两个后缀的 LCP 字典序大小, 时间复杂度为 \(O(n*log^2n)\), 同时哈希检测次数达到了 \(n*log^2n\),非常容易冲突。 显然不是理想的算法。
LCP 最长公共前缀
问:设有一组排序过的字符串 \(A = [A_1, A_2, · · · , A_n]\)。如何快速的求任意 \(A_i\) 与 \(A_j\) 的 LCP?
需要一个关于 LCP 的"区间可加性":
对于任意的 \(k \in [i, j]\),
故有:
证明:
- 令 \(X= A_i[LCP_{ik} + 1], Y = A_k[LCP_{ik} + 1], Z = A_j[LCP_{jk} + 1]\)
- \(LCP(A_i,A_k) \neq LCP(A_k,A_j)\) 时
- 由于 \(X\neq Y, Y=Z\),所以 \(X\neq Z\)。
- 则 \(LCP(A_i,A_j)=LCP(A_i,A_k)=min(LCP(A_i,A_k),LCP(A_k,A_j))\)
- \(LCP(A_i,A_k) = LCP(A_k,A_j)\) 时
- 已知 \(X\neq Y \& Y\neq Z\), 且字典序 \(A_i<A_k<A_j\),所以 \(X<Y<Z\),所以 \(X\neq Z\),结论同样成立
Height 数组
SA 中非常重要的数组
定义: \(height[i] = LCP(sa[i], sa[i - 1])\) , 排名为 \(i\) 的后缀与排名为 \(i-1\) 的后缀的 LCP 长度, 特别地,height[1] = 0
显然有了 Height 数组,刚才的问题就变成了 区间最小值查询 啦!
那么如何求 Height ?
引理:\(height[rk[i]] \geq height[rk[i-1]] - 1\)
展开引理:\(LCP(sa[rk[i]], sa[rk[i]-1]) >= LCP(sa[rk[i-1]],sa[rk[i-1]]) - 1\) ,省略 (S[])。
等价于: \(LCP(i, sa[rk[i] - 1]) >= LCP(i-1, sa[rk[i-1]]) -1\)
因此,不妨设 \(H[i] = LCP(S[i], S[sa[rk[i-1]]])\) ,表示后缀 \(i\) 与排名比他小 \(1\) 的后缀的 LCP
即证: \(H[i] \geq H[i-1] - 1\)
令 \(K1 = sa[rk[i-1]-1]\) , \(K2 = sa[rk[i]-1]\) 。
则 \(H[i-1]=LCP(S[i-1],S[K1])\), \(H[i] = LCP(S[i],S[K2])\)
- 当 \(H[i-1]\leq 1\) 时
- \(H[i]\geq 0\) 显然成立。
- 当 \(H[i-1]> 1\) 时
- 对于 \(H[i-1]\)\(S[K1], S[i-1]\) 去首字母后变为 \(S[K1+1],S[i]\) ,字典序关系是不变的。\(S[K1+1] < S[i]\), 且 \(LCP(S[K1+1],S[i]) = H[i-1] - 1\)
- 而 \(S[K2] < S[i]\), 且中间没有其他后缀,所以有 \(S[K1+1]\leq S[K2]\) 。
- 故 \(S[K1+1]\leq S[K2] \leq S[i]\), 由“区间可加性”得和上式得
\[\begin{aligned} H[i-1] - 1&=LCP(S[K1+1],S[i]) \\ &=min(LCP(S[K1+1],S[K2]),LCP(S[K2],S[i])) \\ &=min(LCP(S[K1+1],S[K2]),H[i]) \end{aligned} \]- 结论依然成立,则 \(height[rk[i]]\geq height[rk[i-1]] - 1\) 成立
代码实现
for (i = 1, k = 0; i <= n; ++i) {
if (rk[i] == 0) continue;
if (k) --k;
while (s[i + k] == s[sa[rk[i] - 1] + k]) ++k;
height[rk[i]] = k;
}
倍增法构造SA
思路
将 比较字典序的二分求 LCP 转化为倍增求 LCP。
首先等效的认为在字符串的末尾增添无限个空字符 \0
按照通常的倍增思路:
- 定义 \(S(i, k) = S[i, i + 2^k+1]\),即以 \(i\) 位置开头,长度为 \(2^k\) 的子串。
- 后缀 \(S[i]\) 与 \(S[j]\) 的字典序关系等价于 \(S(i, ∞)\) 与 \(S(j, ∞)\) 的字典序关系。
- 事实上,只需要将 \(S(i, ⌈log2n⌉)\),\(i = 1, 2, · · · , n\) 排序即可。
然后可以倍增的进行排序
- 假设当前已经得到了 \(S(i, k)\) 的排序结果,
\(rk[S(i, k)]\) 与 \(sa[S(i, k)]\) ,思考如何利用它们排序 \(S(i, k + 1)\)。 - 由于 \(S(i, k + 1)\) 是由 \(S(i, k)\) 和 \(S(i + 2^k, k)\) 前后拼接而成。
- 因此比较 \(S(i, k + 1)\) 与 \(S(j, k + 1)\) 字典序可以转化为先比较 \(S(i, k)\) 与 \(S(j, k)\),
- 再比较 \(S(i + 2^k, k)\) 与 \(S(j + 2^k, k)\)。
- 因此可以将 \(S(i, k + 1)\) 看作一个两位数,高位是 \(rk[S(i, k)]\),低位是 \(rk[S(i + 2^k, k)]\)。
显然有 \(O(n*log^2n)\) 的做法
#include<bits/stdc++.h>
using namespace std;
const int N = 1000010;
char s[N];
int n, sa[N], rk[N << 1], oldrk[N << 1];
// 为了防止访问 rk[i+w] 导致数组越界,开两倍数组。
// 当然也可以在访问前判断是否越界,但直接开两倍数组方便一些。
int main() {
int p;
scanf("%s", s + 1);
n = strlen(s + 1);
for (int i = 1; i <= n; ++i)
sa[i] = i, rk[i] = s[i];
for (int w = 1; w < n; w <<= 1) {
sort(sa + 1, sa + n + 1, [](int x, int y)
{ return rk[x] == rk[y] ? rk[x + w] < rk[y + w] : rk[x] < rk[y]; }); // 这里用到了 lambda
memcpy(oldrk, rk, sizeof(rk));
// 由于计算 rk 的时候原来的 rk 会被覆盖,要先复制一份
for (p = 0, i = 1; i <= n; ++i) {
if (oldrk[sa[i]] == oldrk[sa[i - 1]] &&
oldrk[sa[i] + w] == oldrk[sa[i - 1] + w]) {
rk[sa[i]] = p;
}
else {
rk[sa[i]] = ++p;
} // 若两个子串相同,它们对应的 rk 也需要相同,所以要去重
}
}
for (int i = 1; i <= n; ++i)
printf("%d ", sa[i]);
return 0;
}
而对于两位数的排序,我们有基数排序!
- 将他们排序时,需要先按照高位排序,高位相同时,按照低位排序。此过程为基数排序
- 实际代码运行时,先进行的是低位的排序。
此时借用一下,葫芦爷的课件!其实一直都在借用
算法分析
时间复杂度:
- 总共需要运行 \(logn\) 轮,每轮使用基数排序,复杂度为 \(O(n)\), 整体复杂度为 \(O(n*logn)\)
代码实现
直接干!
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
const int N = 1000010;
char s[N];
int n, sa[N], rk[N << 1], oldrk[N << 1], id[N], cnt[N];
int main() {
int m, p;
scanf("%s", s + 1);
n = strlen(s + 1);
m = max(n, 300);
// 先对长度为 1 的子串进行计数排序
for (int i = 1; i <= n; ++i) ++cnt[rk[i] = s[i]];
for (int i = 1; i <= m; ++i) cnt[i] += cnt[i - 1];
for (int i = n; i >= 1; --i) sa[cnt[rk[i]]--] = i;
for (int w = 1; w < n; w <<= 1) {
// 对第二关键字 rk[id[i] + w] 进行计数排序, id[i] 作为 sa[i] 的备份
memset(cnt, 0, sizeof(cnt));
for (int i = 1; i <= n; ++i) id[i] = sa[i];
for (int i = 1; i <= n; ++i) ++cnt[rk[id[i] + w]];
for (int i = 1; i <= m; ++i) cnt[i] += cnt[i - 1];
for (int i = n; i >= 1; --i) sa[cnt[rk[id[i] + w]]--] = id[i];
memset(cnt, 0, sizeof(cnt));
// 对第一关键字 rk[id[i]] 进行计数排序
for (int i = 1; i <= n; ++i) id[i] = sa[i];
for (int i = 1; i <= n; ++i) ++cnt[rk[id[i]]];
for (int i = 1; i <= m; ++i) cnt[i] += cnt[i - 1];
for (int i = n; i >= 1; --i) sa[cnt[rk[id[i]]]--] = id[i];
memcpy(oldrk, rk, sizeof(rk));
for (int p = 0, i = 1; i <= n; ++i) {
if (oldrk[sa[i]] == oldrk[sa[i - 1]] &&
oldrk[sa[i] + w] == oldrk[sa[i - 1] + w]) {
rk[sa[i]] = p;
}
else {
rk[sa[i]] = ++p;
}
}
}
for (int i = 1; i <= n; ++i)
printf("%d ", sa[i]);
return 0;
}
上述代码可以进行一些常数优化,即:
第二关键字无需计数排序
for (int i = n; i > n - w; --i) // 第二关键字无穷小先放进去
id[++p] = i;
for (int i = 1; i <= n; ++i)
if (sa[i] > w) id[++p] = sa[i] - w; // 顺次放入 s[sa[i]-w] 的第二关键字排名
优化计数排序值域
- 每次对 \(rk\) 进行去重之后,我们都计算了一个 \(p\) ,这个 \(p\) 即是 \(rk\) 的值域,将值域赋值为 \(p\)
- 将 rk[id[i]] 存下来,减少不连续内存访问, 这个优化在数据范围较大时效果非常明显。这个模板
若排名都不相同直接生成后缀数组
考虑新的 \(rk\) 数组,若其值域为 \([1,n]\) 那么每个排名都不同,此时无需再排序。
常数优化版
/*height[i] = lcp(S[sa[i]],S[sa[i-1]]), h[i]=height[rk[i]], h[i]>=h[i-1]-1, lcp(s[i],s[j])=min(height[rk[i]+1],...,height[rk[j]])*/
int n, sa[maxn], rk[maxn], id[maxn], cnt[maxn], height[maxn], px[maxn];
void get_sa(const char* s, int _n) { // get sa and height
n = _n;
int m = 300, p = 0; // m 是值域, 初始化为字符集大小
for (int i = 0; i <= m; i++) cnt[i] = 0;
for (int i = 1; i <= n; ++i) cnt[rk[i] = (int)s[i]] ++; // 先对1个字符大小的子串进行计数排序
for (int i = 1; i <= m; ++i) cnt[i] += cnt[i - 1];
for (int i = n; i >= 1; --i) sa[cnt[rk[i]]--] = i;
for (int w = 1; w <= n; w <<= 1, m = p, p = 0) { // m=p 就是优化计数排序值域
for (int i = n - w + 1; i <= n; ++i) // 第二关键字无穷小先放进去
id[++p] = i;
for (int i = 1; i <= n; ++i)
if (sa[i] > w) id[++p] = sa[i] - w; // 顺次放入 s[sa[i]-w] 的第二关键字排名
for (int i = 0; i <= m; ++i) cnt[i] = 0;
for (int i = 1; i <= n; ++i) ++cnt[rk[i]], px[i] = rk[id[i]];
for (int i = 1; i <= m; ++i) cnt[i] += cnt[i - 1];
for (int i = n; i >= 1; --i) sa[cnt[px[i]]--] = id[i];
for (int i = 1; i <= n; ++i) swap(rk[i], id[i]);
rk[sa[1]] = p = 1;
for (int i = 2; i <= n; ++i) {
rk[sa[i]] = (id[sa[i]] == id[sa[i - 1]] && id[sa[i] + w] == id[sa[i - 1] + w] ? p : ++p);
}
if (p >= n) { // 排名已经更新出来了
break;
}
}
}
void get_height(const char* s){
for (int i = 1, k = 0; i <= n; ++i) { // 获取 height数组
if (k) --k;
int j = sa[rk[i] - 1];
while (s[i + k] == s[j + k]) ++k;
height[rk[i]] = k;
}
#ifdef _DEBUG
for (int i = 1; i <= n; ++i)
cout<<"height["<<i<<"] = "<<height[i]<<endl;
#endif
}
DC3 构造SA
待填
SA_IS
待填
应用
SA 数组性质应用
实现字符串最小表示
将字符串 \(S\) 复制一份变成 \(SS\),找长度大于 \(|S|\) 字典序最小的的后缀对应位置就是最小表示位置
字符串中查找子串
- 在线地在主串 \(T\) 中寻找模式串 \(S\)
- 子串一定是是某个后缀的前缀。先跑一次 SA, 可以在 \(|S|\) 个后缀中按照和 \(T\) 的字典序关系进行二分查找,就可以找到是否出现。
- 如果子串出现多次,如果多次出现可以再次二分查找,两次二分查找可以分别找出在 SA 中的最左位置和最右位置,并且可以依次得到出现位置。
从字符串首尾取字符最小化字典序
- 优化暴力,每次从正尾取是看正串和反串谁小,那么我们就把字符串拼成正串+反串,跑一次 SA。
- 记录首尾位置,比较对应的“后缀”字典序即可!
例题-[USACO007DEC]Best Cow Line G
Height 数组性质应用
求最长公共子串
-
求 \(s,t\) 两串最长公共子串,先将两串用分割符拼接起来。
-
遍历 \(t\) 在 \(sa[]\) 中的所有位置,找到第一个小于 i 和第一个大于 i 的 s 的两个后缀,对应和 t[i] 取最长 LCP 即可,即用数据结构来查最值即可$。
-
求本质不同公共子串个数
- 类似上面的求法,只不过在每次计算的时候需要减去 T 的后缀 \(i\) 与上一个 T 的后缀的 LCP,然后取与 0 取max
比较字符串两个子串大小关系
- 若比较 \(A=S[a,...,b],\; B=S[c...d]\) 大小关系
- 如果 \(LCP(a,c)\geq min(|A|,|B|)\),则
A<B <-> |A| < |B|
- 其中一个是另一个的前缀,长度小的字典序一定不会大于长度大的
- 否则
A < B <-> rk[a] < rk[c]
- 从两串下标为 \([LCP(a,c) + 1]\) 的位置开始看,谁小谁字典序更小
求本质不同子串数目
- 按字典序从小到大枚举所有后缀,统计有多少个新出现的前缀即可。
- 对于排名第 i 的后缀 \(S[sa[i], n]\),共有 \(n - sa[i] + 1\) 个前缀,其中有 \(Height[i]\)
个前缀同时出现在前一个排名的后缀 \(S[sa[i-1], n]\) 中,因此减掉即可。 - 上述证明是不完整的,还需要证明所有在 S[sa[i], n] 中出现,但没有在
\(S[sa[i-1], n]\) 中出现的前缀,他们在所有更小排名的后缀串也都没有出
现。 - 证明就是,如果出现过会破坏我们求 \(lcp\) 时的性质。
- 求本质不同同构子串数目(字符集大小为 \(3\))
- 枚举所有字符集转换,共有 \(3!\) 种。将形成的 \(6\) 种字符串通过不同的分隔符分割,大串跑 SA
- 统计大串 S 的本质不同子串数目,观察发现子串中有 \(2\) 种不同字符时会在 S 中出现 \(6\) 次总次数为 \(cntA\),如果是只有一种字符只会出现 \(3\) 次总次数为 \(cntB\),因此计算时需要将 \(cntB\) 补到 6 次进行计算,统计全 \(a/b/c\) 串只需要一个 for 循环记录最长连续相同子串即可记为 \(single\),因此答案为 \((cntA+cntB+single*3) / 6\)。
- 计算答案前还需要减去包含分隔符的子串,类似双指针的思想
ans -= (n + 1)*(len - i + 1)
查找出现 \(K\) 次的子串的最大长度。
- 一个字符串出现 \(K\) 次表明在相邻的 \(k-1\) 个 height 数组中均出现。
- 那么只需要查找相邻 \(k-1\) 个 height 数组的最小值的最大值即可,单调队列能 \(O(n)\), 也可以用区间最值查询。
洛谷P2852 [USACO06DEC]Milk Patterns G
是否有某字符串在文本串中至少不重叠地出现了两次
可以二分目标串的长度 \(|s|\),将 \(height\) 数组划分成若干个连续 LCP 大于等于 \(|s|\) 的段,利用 RMQ 对每个段求其中出现的数中最大和最小的下标,若这两个下标的距离满足条件,则一定有长度为 \(|s|\) 的字符串不重叠地出现了两次。
连续的若干个相同子串
- 询问 \(S\) 多少个子串满足优秀拆分,即拆分为 \(AABB\) 形式的串。
- 按贡献考虑,观察 AA 和 BB 交界,记录 \(a[i],b[i]\) 分别代表 i 为 AA 串结尾, i 为 BB(与AA同理)串开头。
- 那么最终答案是 \(\Sigma_{i=1}^{n-1}a[i]*b[i+1]\)
- \(O(n^2)\) 加哈希可以拿 \(95\) 分,SA 正解参见 AcFunction's blog
- 第一次做感觉边界没有搞得太明白。。同时注意下多组数据清空的问题。
所有后缀之间 \(lcp\) 的加和。
- 即求
height[]
数组所有区间的最小值加和。 - 定义状态
f[i]
表示前缀 \(i\) 的所有后缀区间的最小值之和。 - 考虑
height[i]
与height[i-1]
的关系:- 如果 height[i] >= height[i - 1],f[i] 能取到 f[i - 1] 的所有值
- 反之,height[i] 不能取到 f[i - 1] 的值,要往前继续找 小于等于 height[i] 的下标。
- 上点可以用单调栈来维护左边第一个小于 height[i] 的下标。
询问不超过 k 次修改单个字符的连续子串匹配个数
- 给定两个字符串 \(S,T\),询问对于 \(S\) 中所有子串,有多少个长度为 \(|T|\) 的联系子串满足不超过 \(K\) 次修改,能变成 \(T\)
- 将 T 接在 S 的后面,跑一次 SA,然后用最多进行 \(k\) 次求 LCP 模拟匹配。