本文介绍后缀数组的构造和简单应用。
大部分内容可以在罗穗骞2009年集训队中看到,本文进行了部分精简。
(刷新以获取数学公式)
定义
1、后缀(\(Suffix\)):对于字符串\(S[1..n]\),用\(Suffix(i)\)表示子串\(S[i..n]\)
2、后缀数组(\(Sa\)):后缀数组 \(Sa\) 是一个一维数组,它保存 \(1..n\) 的某个排列 \(SA[1], SA[2],……,SA[n],\)并且保证 \(Suffix(SA[i])\) < \(Suffix(SA[i+1])\),\(1≤i\)<\(n\)。
3、名次数组(\(Rank\)):名次数组 \(Rank[i]\) 保存的是 \(Suffix(i)\) 在所有后缀中从小到大排列的“名次”。
4、最长公共前缀(\(lcp\)):对于字符串\(A\),\(B\),定义\(lcp(A,B)\)为字符串\(A\),\(B\)的最长公共前缀。对于整数\(i\),\(j\),定义\(lcp(i,j)=lcp(Suffix(i),Suffix(j))\)
5、\(height\)数组:\(height[i]=lcp(sa[i-1],sa[i])\),即排名为\(i\)和\(i-1\)的两个后缀的最长公共前缀。
由\(height\)数组的定义可以得到它的一个性质:对于\(i\)和\(j\)不妨设\(rank[i]\)<\(rank[j]\),那么\(lcp(i,j)=min_{i+1≤k≤j}\{height[k]\}\)
倍增构造后缀数组
基数排序
对于\(n\)个双关键字数据,数据范围在\([1..m]\)之间,已知两个数组:\(x[i]\)表示第\(i\)个数第一关键字的值,\(y[i]\)表示第二关键字排名为\(i\)的数的位置(这里\(y[i]\)是两两不同的,如果两数两个关键字都相同,我们可以用序号比较出他们的第二关键字)我们可以用基数排序在\(O(n+m)\)的复杂度内完成对他们的排序。
对第一关键字开\(m\)个桶,然后按第二关键字从大到小的顺序将他们扔进桶中即可完成排序。
for (int i = 1; i <= m; ++i) cnt[i] = 0;
for (int i = 1; i <= n; ++i) ++cnt[x[i]];
for (int i = 2; i <= m; ++i) cnt[i] += cnt[i-1];
for (int i = n; i >= 1; --i) sa[cnt[x[y[i]]]--] = y[i];
倍增
我们进行\(O(logn)\)次排序.第\(x\)次排序以\(S[i..i+2^x)\)为第一关键字,以\(S[i+2^x..i+2^{x+1})\)为第二关键字.不过最开始要以第二关键字是序号(即\(y[i]=i\))做一次.
#include<iostream>
#include<cstdio>
#include<cstring>
#define rint register int
#define inv inline void
#define ini inline int
#define maxn 1000050
using namespace std;
char s[maxn];
int y[maxn],x[maxn],c[maxn],sa[maxn],rk[maxn],height[maxn],wt[30];
int n,m;
liti
inv putout(int x)
{
if(!x) {putchar(48);return;}
rint l=0;
while(x) wt[++l]=x%10,x/=10;
while(l) putchar(wt[l--]+48);
}
inv get_SA()
{
for (rint i=1;i<=n;++i) ++c[x[i]=s[i]];
//c数组是桶
//x[i]是第i个元素的第一关键字
for (rint i=2;i<=m;++i) c[i]+=c[i-1];
//做c的前缀和,我们就可以得出每个关键字最多是在第几名
for (rint i=n;i>=1;--i) sa[c[x[i]]--]=i;
for (rint k=1;k<=n;k<<=1)
{
rint num=0;
for (rint i=n-k+1;i<=n;++i) y[++num]=i;
//y[i]表示第二关键字排名为i的数,第一关键字的位置
//第n-k+1到第n位是没有第二关键字的 所以排名在最前面
for (rint i=1;i<=n;++i) if (sa[i]>k) y[++num]=sa[i]-k;
//排名为i的数 在数组中是否在第k位以后
//如果满足(sa[i]>k) 那么它可以作为别人的第二关键字,就把它的第一关键字的位置添加进y就行了
//所以i枚举的是第二关键字的排名,第二关键字靠前的先入队
for (rint i=1;i<=m;++i) c[i]=0;
//初始化c桶
for (rint i=1;i<=n;++i) ++c[x[i]];
//因为上一次循环已经算出了这次的第一关键字 所以直接加就行了
for (rint i=2;i<=m;++i) c[i]+=c[i-1];//第一关键字排名为1~i的数有多少个
for (rint i=n;i>=1;--i) sa[c[x[y[i]]]--]=y[i],y[i]=0;
//因为y的顺序是按照第二关键字的顺序来排的
//第二关键字靠后的,在同一个第一关键字桶中排名越靠后
//基数排序
swap(x,y);
//这里不用想太多,因为要生成新的x时要用到旧的,就把旧的复制下来,没别的意思
x[sa[1]]=1;num=1;
for (rint i=2;i<=n;++i)
x[sa[i]]=(y[sa[i]]==y[sa[i-1]] && y[sa[i]+k]==y[sa[i-1]+k]) ? num : ++num;
//因为sa[i]已经排好序了,所以可以按排名枚举,生成下一次的第一关键字
if (num==n) break;
m=num;
//这里就不用那个122了,因为都有新的编号了
}
for (rint i=1;i<=n;++i) putout(sa[i]),putchar(' ');
}
int main()
{
gets(s+1);
n=strlen(s+1);m=122;
//我们第一次读入字符直接不用转化,按原来的ascll码来就可以了
get_SA();
}
\(height\)数组的构造
性质
定义\(h[i]=height[rank[i]]\),也就是 \(suffix(i)\)和在它前一名的后缀的最长公共前缀。
\(h\) 数组有以下性质: \(h[i]≥h[i-1]-1\)
当\(h[i-1]=1\)时显然成立.当\(h[i-1]>1\)时,设 \(suffix(k)\) 是排在 \(suffix(i-1)\)前一名的后缀,则它们的最长公共前缀是\(h[i-1]\)。那么\(suffix(k+1)\)将排在\(suffix(i)\)的前面并且 \(suffix(k+1)\)和 \(suffix(i)\)的最长公共前缀是\(h[i-1]-1\),所以 \(suffix(i)\)和在它前一名的后缀的最长公共前缀至少是 \(h[i-1]-1\)。
构造
利用这个\(h\)数组的性质,按照 \(h[1],h[2],……,h[n]\)的顺序可以在线性时间内处理出\(height\)数组.
void calc_height() {
for (int i = 1; i <= n; ++i) rnk[sa[i]] = i;
for (int i = 1; i <= n; ++i) {
int p = sa[rnk[i] - 1], k = max(0, hgt[rnk[i - 1]] - 1);
while (s[p + k] == s[i + k]) ++k;
hgt[rnk[i]] = k;
}
}
例题
可重叠的 k 次最长重复子串
来源:POJ3261
题意
给定一个字符串,求至少出现 \(k\) 次的最长重复子串,这 \(k\) 个子串可以重叠。
保证\(s.length ≤ 100000\)
题解
后缀数组的一个技巧:用\(height\)数组将后缀分组。
题目等价于:求一个最长的子串,它是至少 \(k\) 个后缀的前缀。
先二分一个答案 \(mid\) ,如果 \(height[i] ≥ mid\), 将 \(sa[i]\) 和 \(sa[i - 1]\) 分在一组。只需记录有没有大小超过 \(k\) 的组即可。
#include<cstdio>
#include<iostream>
using namespace std;
const int N = 1e6 + 10;
int n, k;
int a[N], x[N], y[N], sa[N], rnk[N], cnt[N], hgt[N];
void calc_sa() {
int m = 1e6 + 1;
for (int i = 1; i <= n; ++i) ++cnt[x[i] = a[i]];
for (int i = 2; i <= m; ++i) cnt[i] += cnt[i - 1];
for (int i = n; i >= 1; --i) sa[cnt[x[i]]--] = i;
for (int w = 1, p; w <= n; w <<= 1, m = p) {
p = 0;
for (int i = n - w + 1; i <= n; ++i) y[++p] = i;
for (int i = 1; i <= n; ++i) if (sa[i] > w) y[++p] = sa[i] - w;
for (int i = 1; i <= m; ++i) cnt[i] = 0;
for (int i = 1; i <= n; ++i) ++cnt[x[i]];
for (int i = 1; i <= m; ++i) cnt[i] += cnt[i - 1];
for (int i = n; i >= 1; --i) sa[cnt[x[y[i]]]--] = y[i];
swap(x, y);
x[sa[1]] = 1; p = 1;
for (int i = 2; i <= n; ++i) {
int f = (y[sa[i]] == y[sa[i - 1]] && y[sa[i] + w] == y[sa[i - 1] + w]);
x[sa[i]] = f ? p : ++p;
}
if (p == n) break;
}
}
void calc_height() {
for (int i = 1; i <= n; ++i) rnk[sa[i]] = i;
for (int i = 1; i <= n; ++i) {
int p = sa[rnk[i] - 1], k = max(0, hgt[rnk[i - 1]] - 1);
while (a[p + k] == a[i + k]) ++k;
hgt[rnk[i]] = k;
}
}
int check(int mid) {
int cnt = 0;
for (int i = 2; i <= n; ++i) {
if (hgt[i] < mid) cnt = 1;
else if (++cnt >= k) return 1;
}
return 0;
}
int main() {
cin >> n >> k;
for (int i = 1; i <= n; ++i) {
scanf("%d", &a[i]);
++a[i];
}
calc_sa();
calc_height();
int l = 1, r = n, ans = 0;
while (l <= r) {
int mid = (l + r) >> 1;
if (check(mid)) {
l = mid + 1;
ans = mid;
}
else r = mid - 1;
}
cout << ans << endl;
return 0;
}
本质不同子串个数
来源:HDU4622
题意
给定字符串\(S\), \(q\)组询问,每次询问\(S[l..r]\)本质不同子串个数。
保证 \(s.length≤2000\), \(q ≤ 10000\)
题解
先考虑计算\(S[1..n]\)本质不同子串个数。
字符串\(suffix(i)\)的前缀有\(n-sa[i]+1\)个,其中\(height[i]\)个和\(suffix(i-1)\)重复,所以每次将 \(n-sa[i]+1-lcp(i, i - 1)\)累加进答案即可。由于保证了字典序,很显然这种做法是正确的。
而在询问一个区间时,每个后缀被删掉了一个后缀。如果\(suffix(sa[i]) < suffix(sa[i-1])\),一定有\(suffix(sa[i])\)变为了\(suffix(sa[i-1])\)的前缀,那么我们直接无视掉\(suffix(sa[i])\)。对于剩下的串,按顺序将\(n-sa[i]+1-lcp(i, lcp)\)累加进答案即可。
还有一种做法是对于每个前缀求一遍后缀数组,然后从而将问题转化为询问\(S[l..n]\)本质不同子串个数。
void calc_sa();
void calc_height();
void Solve() {
scanf("%s", s + 1);
n = strlen(s + 1);
calc_sa();
calc_height();
cin >> q;
while (q--) {
int l, r;
scanf("%d%d", &l, &r);
int ans = 0, lcp = 0;
for (int i = 1; i <= n; ++i) {
lcp = min(lcp, hgt[i]);
if (l <= sa[i] && sa[i] <= r) {
if (r - sa[i] + 1 - lcp > 0) {
ans += max(0, r - sa[i] + 1 - lcp);
lcp = r - sa[i] + 1;
}
}
}
printf("%d\n", ans);
}
}
最长回文子串
题意
给定字符串\(S\),求最长回文子串。
保证\(s.length≤100000\)
题解
manacher模板题,也可以用后缀数组做到线性复杂度。
将字符串翻转加在原来的字符串后面,中间用一个特殊符号连接。
比如:\(S="aabaaaab"\),\(S'="aabaaaab?baaaabaa"\)
原来的一个回文串就可以转化为\(S'\)中两个后缀的最长公共前缀。
比如以\(S[3]\)为中心的回文串,对应\("aaaab?baaaabaa"\)和\("aa"\)的公共前缀。
以\(S[3]和S[4]\)为中心的回文串,对应\("aaaab?baaaabaa"\)和\("baa"\)的公共前缀。
用DC3求后缀数组,\(O(N)\)预处理RMQ就可以做到线性复杂度了。
大概是没有人愿意打的。
最长公共子串
来源:POJ2774
题意
给定字符串\(A\),\(B\),求最长公共子串。
保证\(a,b.length≤100000\)
题解
将字符串拼成 \(A?B\),找在\('?'\)两边的最大\(height\)即可。
void calc_sa();
void calc_height();
int main() {
scanf("%s", s + 1); n = strlen(s + 1);
scanf("%s", t + 1); m = strlen(t + 1);
s[n + 1] = 'A';
for (int i = 1; i <= m; ++i) s[n + i + 1] = t[i];
calc_sa(n + m + 1);
calc_height(n + m + 1);
int ans = 0;
for (int i = 1; i <= n + m + 1; ++i) {
if (sa[i] == n + 1) continue;
if (sa[i - 1] == n + 1) continue;
int flag = (sa[i] <= n) ^ (sa[i - 1] <= n);
if (flag) ans = max(ans, hgt[i]);
}
cout << ans << endl;
return 0;
}
长度不小于 k 的公共子串的个数
来源:POJ3415
题意
给定字符串\(A\),\(B\),求三元组\((i,j,k)\)的数目,使得\(A[i..i+k-1]=B[j..j+k-1],k>=K\)
保证\(a,b.length≤100000\)
题解
将字符串拼成 \(A?B\)。加入\(suffix(sa[i])\)时,对于\(j\)<\(i\),产生的贡献是\(([sa_i≤n] xor [sa_j≤n])*min_{j+1≤k≤i}\{height_k\}\)
对\(sa_i≤n\)和\(sa_i>n\)分别维护单调栈。每次将单调栈中所有大于\(height_{i+1}\)的元素更改为\(height_{i+1}\)。注意要把相同的元素压在一起,这样每弹出一次栈中元素就少一个,每次只加入一个元素,就保证了这部分复杂度线性。
void calc_sa();
void calc_height();
int main() {
while (scanf("%d", &k) && k) {
scanf("%s", s + 1); n = strlen(s + 1);
scanf("%s", t + 1); m = strlen(t + 1);
s[n + 1] = '!';
for (int i = 1; i <= m; ++i) s[n + i + 1] = t[i];
calc_sa(n + m + 1);
calc_height(n + m + 1);
stack<pair<int, int> > st[2];
long long ans = 0, val[2] = {0, 0};
for (int i = 1; i <= n + m + 1; ++i) {
if (sa[i] == n + 1) continue;
int p = sa[i] <= n, v = hgt[i + 1], cnt[2] = {0, 0};
ans += val[p];
++cnt[p ^ 1];
val[p ^ 1] += max(0, v - k + 1);
for (int j = 0; j < 2; ++j) {
while (!st[j].empty() && st[j].top().first > v) {
pair<int, int> tmp = st[j].top();
st[j].pop();
val[j] -= 1LL * tmp.second * max(0, tmp.first - k + 1);
val[j] += 1LL * tmp.second * max(0, v - k + 1);
cnt[j] += tmp.second;
}
if (cnt[j]) st[j].push(make_pair(v, cnt[j]));
}
}
printf("%lld\n", ans);
}
return 0;
}