Dai Chao

fordring866@gmail.com | Put your faith in the Light!

破解外星代码:详解后缀数组模板

朕听说,后缀数组只需寥寥二十行代码即可求出?

是,陛下。

快,快呈上来!

是,陛下。

int wa[MAXN], wb[MAXN], wv[MAXN], ws[MAXN];
inline bool cmp(int *r, int a, int b, int len) {
	return r[a]==r[b] && r[a+len]==r[b+len];
}
void SA(char *r, int *sa, int n, int m)
{
	int i, j, p, *x = wa, *y = wb, *t;	
	for (i = 0; i < m; i++) ws[i] = 0;
	for (i = 0; i < n; i++) ws[x[i] = r[i]]++;
	for (i = 1; i < m; i++) ws[i] += ws[i - 1];
	for (i = n - 1; i >= 0; i--) sa[--ws[x[i]]] = i;
	for (j = p = 1; p < n; j <<= 1, m = p)
	{
		for (p = 0, i = n - j; i < n; i++) y[p++] = i;
		for (i = 0; i < n; i++) if (sa[i] >= j) y[p++] = sa[i] - j;
		for (i = 0; i < m; i++) ws[i] = 0;
		for (i = 0; i < n; i++) ws[wv[i] = x[y[i]]]++;
		for (i = 1; i < m; i++) ws[i] += ws[i - 1];
		for (i = n - 1; i >= 0; i--) sa[--ws[wv[i]]] = y[i];
		for (t=x,x=y,y=t, x[sa[0]]=0, p=i=1; i<n; i++)
			x[sa[i]] = cmp(y,sa[i-1],sa[i],j) ? p-1 : p++;
	}
}

陛下……陛下!快请御医!

前言

鉴于我也是第一次接触后缀数组,理解程度还很低,再加上类似的博客早已烂大街,因此不求读者能看懂或学到什么东西了。你还可以向我留言提出指导意见,我会非常欢迎。

什么是后缀数组

今天,我们来谈谈这份“后缀数组”的模板。

这份模板,是用来求“后缀数组”的。

后缀很简单,数组也很简单,但是“后缀数组”就需要解释一番了。

一个长为N的字符串,有N个后缀,每个后缀都从不同的位置开始。

我们用suffix(i)表示从第i个字符开始的后缀(从0计)。例如,”abcdefg”的suffix(2)就是”cdefg”。

也就是说,只要确定了字符串,我们就可以用[0, N-1]的数字表示任意一个后缀。

然后,这些后缀作为字符串,是可以按照字典序来排列的。

suffix(0)...suffix(N-1)从小到大排列后,会被打乱为suffix(sa[0])...suffix(sa[i])

这里的sa是0到N-1的一个排列,同时也可以视为数组。

它就是后缀数组(suffix array)sa[i]就是第i小的后缀在原字符串的位置(从0计),而suffix(sa[i])就是第i小的后缀。

那么,我们为什么要求这个数组呢?

不知道,听大佬们说很有用,所以先求出来再说吧。

暴力求法

假如逐一比较两个后缀的字符,就需要比较O(N)次,才能知道两个后缀的相对大小。(因为它们的长度都在1到N之间)

假如用插入排序,上面的比较就需要进行O(N^2)次。

乘起来就是O(N^3),不能忍。

文明求法

目前已经有了好几种O(N)的算法,确实厉害,但我懒得看。

倍增法

前面那段外星代码,使用的就是倍增法。

我先定义一个概念,叫做k-后缀。它也是字符串的后缀,但它长度至多为k——把长于k的部分截断即可。

比如,”abcdefg”的4-后缀分别是”abcd”, “bcde”, “cdef”, “defg”, “efg”, “fg”和”g”。

为了方便表示,我就借用之前的符号,让suffix(i,k)表示从第i个字符开始的k-后缀吧。

相应地,sa(k)就叫做k-后缀数组,suffix(sa(k)[i],k)就是第i小的k-后缀。

至于倍增法,就是先求sa(1),再求sa(2)sa(4),直到k-后缀成为完整的后缀,sa也就出来了。

求sa(1)

sa(1),1-后缀数组,sa(1)[i]就是第i小的1-后缀的位置。

第i小的1-后缀……那不就是第i小的字符么?听上去很简单的样子。

也就是说,我只要知道每个字符是第几小就好了。下面的代码完成了这件事情。

for (i = 0; i < m; i++) ws[i] = 0;
for (i = 0; i < n; i++) ws[x[i] = r[i]]++;
for (i = 1; i < m; i++) ws[i] += ws[i - 1];
for (i = n - 1; i >= 0; i--) sa[--ws[x[i]]] = i;

变量介绍:

for (i = 0; i < m; i++) ws[i] = 0;

清零ws数组。

for (i = 0; i < n; i++) ws[x[i] = r[i]]++;

之后,ws就存储了r(或者x)中每种字符的个数。

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

这句完成后,新的ws[i]成为了原来的ws[0]+ws[1]+...+ws[i],意为x中有多少字符的值不大于i,而不再是等于i

for (i = n - 1; i >= 0; i--) sa[--ws[x[i]]] = i; // 倒序遍历

结合1-后缀数组的定义,这句话的意思是,位置i的字符,是第--ws[x[i]]小的。

短小精悍👍

从sa(k)到sa(2k)

知道了sa(k),怎么求sa(2k)呢?

sa(k)告诉了我们k-后缀的排列顺序,而2k-后缀相当于两个相邻的k-后缀:suffix(i,2k) = suffix(i,k) + suffix(i+k,k)

我们把2k-后缀suffix(i,2k)的前半部分称作第一关键字,后半部分称作第二关键字

因此,对于两个2k后缀,只需先比较第一关键字,再比较第二关键字。

这样的比较并不难,难的是我们需要对第一关键字做基数排序,这样才能保证每次倍增都是O(N)的复杂度。

我们知道,刚才对单个字符做基数排序是非常简单的。基数排序需要一个数组,记录每种待比较元素的个数。待比较元素有多少种,这个数组的尺寸就要有多大。由于字符的大小可以直接用值区分,因此可以直接把字符值作为数组下标。但是多字符组成的k-后缀,是不能直接作为数组下标的。

虽然STL map支持该功能,但其插入复杂度为O(logN),插入N次则为O(NlogN)……你懂我1⃣️4⃣️8⃣️❓

实际上,k-后缀只有n个,因此最多只有n种。有没有一种办法,能够把每个k-后缀用[0,n)的自然数标记呢?

答案是肯定的,这个“标记数组”就是x。

求sa(1)时,x就是字符串本身,x[i]标记的就是r[i]的值。但之后,x[i]将成为小于n的自然数。suffix(i,k)就是第x[i]小的k-后缀。

有了x数组,我们对k-后缀进行基数排序,就非常轻松了。

了解了x数组的意义后,我们就正式进入这层可怕的循环吧。

构造y数组

变量介绍:

for (p = 0, i = n - j; i < n; i++) y[p++] = i;

[n-j, n)依次填入y数组。这些位置的2k-后缀是没有第二关键字的,因此直接排在前面,并且位置从小到大(确保稳定性)。

for (i = 0; i < n; i++) if (sa[i] >= j) y[p++] = sa[i] - j;

这句话想把剩下的2k-后缀,按照第二关键字的大小顺序,依次填入y数组。

遍历sa让我们从小到大地访问了k-后缀的位置,即sa[i]。

把这些k-后缀当作第二关键字,将其位置减去k,就是对应2k-后缀的位置,即sa[i] - j

按此顺序填入2k-后缀的位置,就能将2k-后缀按第二关键字排序。

升级sa数组

变量介绍:

for (i = 0; i < m; i++) ws[i] = 0;
for (i = 0; i < n; i++) ws[wv[i] = x[y[i]]]++;
for (i = 1; i < m; i++) ws[i] += ws[i - 1];
for (i = n - 1; i >= 0; i--) sa[--ws[wv[i]]] = y[i];

这段代码其实与前面求sa(1)是同一个样子,但是第二句话和之前不太一样:遍历名次数组x的顺序不是从0到n-1,而是从y[0]到y[n-1]。但由于y数组本身就是0到n-1的另一个排列,因此ws的计算结果依然正确。同时,又令wv[i]成为x[y[i]],便于后面使用。这一句话相当于做了两句话的事情,非常巧妙。

第四句话也不太一样:把sa[--ws[x[i]]] = i;换成了sa[--ws[x[y[i]]]] = y[i];

这一换,是基数排序的精髓。

如果不换,那么得到的还是稳定排序后的sa(k),如同前面计算sa(1)那样。

但换了之后,k-后缀(2k-后缀的第一关键字)依然是递增的。这是因为,y[i]处的k-后缀越大,wv[i]就越大。由于ws数组递增,ws[wv[i]]也越大,因此该k-后缀的位置在sa的更后面。

当y[i]处的k-后缀相同时,又会如何呢?

可以肯定,它们对应的wv[i]值也相同。

由于ws[wv[i]]会随赋值而递减,因此在这些相同的k-后缀中,先赋值者在sa的下标更大。

由于i递减,因此第二关键字更大的2k-后缀,是先赋值者。

结合上面两句话,对于第一关键字相同的2k-后缀,其第二关键字在sa也是递增的。

于是,sa数组就成为了真正的2k-后缀数组。

服了👍

更新x数组

for (t=x,x=y,y=t, x[sa[0]]=0, p=i=1; i<n; i++)
	x[sa[i]] = cmp(y,sa[i-1],sa[i],j) ? p-1 : p++;
  1. 交换x和y。x即将成为2k-后缀的名次,而y则记录着k-后缀的名次。
  2. x[sa[0]]=0意味着最小的2k-后缀的名次为0,是最小的名次。
  3. 从1遍历到n-1,计算cmp(y,sa[i-1],sa[i],j),也就是比较第i小的2k-后缀,和它前面那个第i-1小的2k-后缀。比较的结果有两种:大于或等于。有了y数组,我们不用逐字符比较。先比较第一关键字,再比较第二关键字即可。
  4. 若两者相等,则cmp返回true,把现有名次p-1赋给x[sa[i]];否则返回false,意味着需要一个更大的名次,于是把p赋给x[sa[i]],再自增。

循环的终止

一次循环结束后,sa数组从sa(k)升级为sa(2k),x数组也随之更新。

此外,代表k值的变量j也要加倍,代表不同2k-后缀个数的m,则成为上一轮通过循环递增算出来的p。

但是循环的终止条件不是k达到n,而是p达到n!

当p达到n时,说明已经有了n个不同的k-后缀,并且都在sa里排好了序。即使k再增大,它们的顺序也会保持不变。因此p达到n时,循环即可终止。

服了👍

后记

之前写代码,命名长得一比,换行频繁,美其名曰“整齐易读”,视压缩代码为毒瘤。

今日见了传说中的后缀数组,才知代码的压缩也是需要水平的。这么多细致的计算,全部约束于几行紧密而规整的for循环中。

下篇应该是LCP数组和后缀树的构造。