【算法学习】后缀数组

一个字符串的题，有姿势水平的OIers的脑中应该要浮现出许多算法……

但是我没有姿势，也没有水平，除了KMP和trie树，什么也想不起来。

直到我学了它——后缀数组！

多亏这玩意儿，我现在什么都想不起来了。

---

后缀数组干嘛用的？

主要处理同一个字符串中的重复子串问题。

如何实现？

注意到每一个子串，都是一个后缀的某个前缀，这个后缀和前缀都是唯一确定的。

而后缀相同的前缀，和他们的字典序有密切联系。你有没有想过，字典中的相邻单词，他们的公共前缀总是很长。

一个字符串的任意后缀，都能用它的起始位置的下标唯一确定。

一个字符串的后缀数组(SuffixArray/SA)，就是它的所有后缀字典序排序后的结果，当然不是存储了所有的后缀，仅仅是存储了它们的起始位置罢了。

这篇文章并不是要详细说明后缀数组能干嘛的，怎么实现。它只是让你们对后缀数组的实现有更直观的理解，弄清楚代码的每一行都是干嘛的。

这篇文章适合用作补充材料。如果你是初学者，想要对理论有着更基础的认识，建议百度一下，或者看看这篇blog或者这篇blog。

因为博主学习的时候，也弄不太清楚SA到底怎么求出的，那个神奇的计数排序如何实现等等……

下面贴出一个注释详尽的模板，来为大家理清思路。

#include<cstdio>
#include<cstring>
#define MAXN 1000005
int len, m=127;
//len : 字符串长度 , m : 不同排名的个数 / 字符集大小(初始时)
char str[MAXN];
//字符串 , 从 1 开始编号
int rank[MAXN], SA2[MAXN];
//rank : 第 i 次排序后的 rank / 第 i+1 次排序前的第一关键字
//SA2[k] : 第 i+1 次排序前的第二关键字第 k 小的后缀下标
int cnt[MAXN], tmp[MAXN];
//cnt : 桶排序用的"桶" , tmp : 重新计算 rank 的辅助数组 , 这两个数组可以合并
int SA[MAXN];
//后缀数组 , SA[i] : 排名为 i 的后缀的起始位置
int Height[MAXN];
//Height[i] = LCP(Suffix(SA[i]), Suffix(SA[i-1]))
void getHeight(){
    //计算方法 : Height[rank[i]] >= Height[rank[i-1]] - 1
    //即意 : "长度为 i 的后缀和它排名前一名的 LCP 长度" , 不比 "长度为 i+1 的后缀与其排名前一名的 LCP 长度" 减一来得小
    //所以从最长的后缀(整个字符串,即rank[1],对应Height[rank[1]])开始计算 , 然后缩短长度直至 1
    //保证复杂度 O(len)
    //PS : Height[1]=0
    int k=0;
    for(int i=1;i<=len;++i){
        if(rank[i]==1) {k=Height[1]=0; continue;}
        if(k) --k;
        int j=SA[rank[i]-1];
        while(i+k<=len&&j+k<=len&&str[i+k]==str[j+k]) ++k;
        Height[rank[i]]=k;
    }
//    for(int i=1;i<=len;++i) printf("%d ",Height[i]); puts("");
}
inline void Rsort() {
    //通过 rank , SA2 两个数组 , 确定新的 SA
    for(int i=1; i<=m; ++i) cnt[i]=0;
    //这里用 memset 也可以
    for(int i=1; i<=len; ++i) ++cnt[rank[i]];
    for(int i=1; i<=m; ++i) cnt[i]+=cnt[i-1];
    //计数排序的 nb 操作
    for(int i=len; i>=1; --i) SA[cnt[rank[SA2[i]]]--]=SA2[i];
    //更新 SA
}
void getSA() {
    for(int i=1; i<=len; ++i) rank[i]=str[i], SA2[i]=i;
    Rsort();
//    printf("rk : "); for(int i=1;i<=len;++i) printf("%d ",rank[i]); puts("");
//    printf("SA : "); for(int i=1;i<=len;++i) printf("%d ",SA[i]); puts("");
    for(int j=1; j<=len; j<<=1) {
        int p=0;
        for(int i=len-j+1; i<=len; ++i) SA2[++p]=i;
        for(int i=1; i<=len; ++i) if(SA[i]>j) SA2[++p]=SA[i]-j;
        //计算 SA2
        Rsort();
        //计算新的 SA
        tmp[SA[1]]=p=1;
        for(int i=2; i<=len; ++i) {
            if(rank[SA[i]]!=rank[SA[i-1]]||rank[SA[i]+j]!=rank[SA[i-1]+j]) ++p;
            tmp[SA[i]]=p;
        }
        for(int i=1; i<=len; ++i) rank[i]=tmp[i];
        //根据 SA 更新新的 rank
        m=p;
        if(m==len) break;
        //不同排名的个数 m 更新成 p , 当所有排名都不同时 , 就可以退出了
        //实践证明这一句话加上会快(数据较随机) , 但是理论上可以被卡掉
//        printf("rk : "); for(int i=1;i<=len;++i) printf("%d ",rank[i]); puts("");
//        printf("SA : "); for(int i=1;i<=len;++i) printf("%d ",SA[i]); puts("");
//        printf("SA2: "); for(int i=1;i<=len;++i) printf("%d ",SA2[i]); puts("");
        //输出每步结果
    }
    getHeight();
}
int main() {
    scanf("%s",str+1);
    len=strlen(str+1);
    getSA();
    return 0;
}

接下来解释一下代码中难懂的部分，也是众多blog没有说清楚的地方：计数排序。

代码37行：

for(int i=1; i<=len; ++i) ++cnt[rank[i]];

rank[i]表示第一关键字，cnt就是桶。
这只是在统计桶中元素罢了。

代码38行：

for(int i=1; i<=m; ++i) cnt[i]+=cnt[i-1];

它把桶做了前缀和，这样的目的在于，我们知道了对于某一个第一关键字，它所对应的区间：( cnt[i-1] , cnt[i] ]，这是一个左开右闭区间。
仔细思考一下是不是这样。

代码40行：

for(int i=len; i>=1; --i) SA[cnt[rank[SA2[i]]]--]=SA2[i];

对于确定的第一关键字，我们知道对应区间，但是在同一个区间中，第二关键字还要保证升序啊。
看看SA2数组，它是第二关键字从小到大排序对应的后缀下标，那么在 i 从大到小的循环过程中，SA2[i]就是满足第二关键字从大到小的后缀下标。
那么再在外面套一个rank[SA2[i]]呢？如果有两个 i 的rank[SA2[i]]是相同的，即它们的第一关键字相同，那么我们说，是第二关键字大的先访问到。
那么在SA数组中，这就表示成，在同一个第一关键字下，按照着第二关键字从大到小的顺序，数组从后往前填充着后缀下标。
那就是说，在满足了第一关键字有序的情况下，从前往后，第二关键字也是从小往大的。

这就是计数排序的原理，还是非常巧妙的。

---

新模板：

#include <cstdio>
#include <cstring>
#include <algorithm>

const int MN = 100005;

int N, M;
char str[MN];

struct SuffixArray {
	int N, Str[MN], SA[MN], rk[MN], SA2[MN];
	int Sig, Buk[MN], Tmp[MN];
	int Height[MN];
	
	inline void RSort() {
		for (int i = 1; i <= Sig; ++i) Buk[i] = 0;
		for (int i = 1; i <= N; ++i) ++Buk[rk[i]];
		for (int i = 1; i <= Sig; ++i) Buk[i] += Buk[i - 1];
		for (int i = N; i >= 1; --i) SA[Buk[rk[SA2[i]]]--] = SA2[i];
	}
	
	inline void BuildSA() {
		/* Init Str */
		for (int i = 1; i <= N; ++i) rk[i] = Str[i], SA2[i] = i;
		rk[N + 1] = 0;
		Sig = 26 /* maximum letter in Str */, RSort();
		for (int j = 1; j < N; j <<= 1) {
			int P = 0;
			for (int i = 1; i <= j; ++i) SA2[++P] = N - j + i;
			for (int i = 1; i <= N; ++i) if (SA[i] > j) SA2[++P] = SA[i] - j;
			RSort();
			Tmp[SA[1]] = P = 1;
			for (int i = 2; i <= N; ++i) {
				if (rk[SA[i]] != rk[SA[i - 1]] || rk[SA[i] + j] != rk[SA[i - 1] + j]) ++P;
				Tmp[SA[i]] = P;
			}
			for (int i = 1; i <= N; ++i) rk[i] = Tmp[i];
			Sig = P;
			if (Sig == N) break;
		}
	}
	
	inline void GetHeight() {
		int k = 0;
		for (int i = 1; i <= N; ++i) {
			if (rk[i] == 1) k = Height[1] = 0;
			else {
				if (k) --k;
				int j = SA[rk[i] - 1];
				while (i + k <= N && j + k <= N && Str[i + k] == Str[j + k]) ++k;
				Height[rk[i]] = k;
			}
		}
	}
	
	int BLCP[MN][21], Bt;
	
	inline void InitST() {
		for (int i = 1; i <= N; ++i) BLCP[i][0] = Height[i];
		while (2 << Bt <= N) ++Bt;
		for (int j = 1; j <= Bt; ++j)
			for (int i = 1 << j; i <= N; ++i)
				BLCP[i][j] = std::min(BLCP[i][j - 1], BLCP[i - (1 << (j - 1))][j - 1]);
	}
	
	inline int LCP(int x, int y) {
		if (x == y) return N + 1;
		x = rk[x], y = rk[y];
		if (x > y) std::swap(x, y);
		int j = 31 - __builtin_clz(y - x);
		return std::min(BLCP[y][j], BLCP[x + (1 << j)][j]);
	}
} C1, C2;

int main() {
	scanf("%d%d", &N, &M);
	scanf("%s", str + 1);
	for (int i = 1; i <= N; ++i) C2.Str[N - i + 1] = C1.Str[i] = str[i] - 'a' + 1;
	C2.N = C1.N = N, C1.BuildSA(), C2.BuildSA(), C1.GetHeight(), C2.GetHeight();
	for (int i = 1; i <= N; ++i) printf("%3d : %s\n", i, str + C1.SA[i]);
	return 0;
}

ideone 链接：https://ideone.com/RHDZdU。