后缀数组(倍增)算法原理及应用

后缀数组(倍增)算法原理及应用

概述

后缀数组是指对某一字符串的所有后缀按照字典序排序后的结果,这一信息常用于处理各种字符串问题,维护该信息的算法之一:倍增法的复杂度为O(n \log n)。尽管有更优秀的DC3法,但是倍增法更加通用好写,因此这里只介绍倍增法求后缀数组。

原理与实现

最初的思路

我们知道,如果处理出所有后缀,再使用快速排序等排序方式,复杂度是O(n^2 \log n)的。
我们希望优化复杂度,可以从排序算法上下手。比快排更优秀的排序方法也是有的,比如计数排序。计数排序需要消耗值域规模的空间,如果用在此处空间应该是O(n)的,不成问题。并且,计数排序并非基于比较的排序算法,之前复杂度中计入的O(n)比较字符串字典序现在也不需要了。那么,我们单次排序的复杂度降为O(n)了。
可是,计数排序需要有一个数值来代表每个字符串,基于数值来比较大小。对于代表字符串的数值,我们首先可以想到Hash。不过字符串长了Hash要取模无法比较大小不说,光是计算一遍Hash都需要O(n^2)时间,白白浪费了精力对排序进行优化。

倍增?

接下来是一个很神奇的操作,即倍增法的精髓所在:倍增。
我们可以考虑倍增关键字的长度,采用双关键字排序的方法。比如,当长度为1的时候,两个关键字就是每个后缀的第一个和第二个字符。每次对长度倍增,最多当长度变为 \displaystyle \lfloor \frac{n}{2} \rfloor的时候,所有排序就完成了。这样来看,复杂度就是O(n \log n)了。
接下来,我们要考虑具体的实现问题,比如,如何完成双关键字的计数排序。这一段,不放我们来直接阅读分析实现代码。

代码分析

建议在阅读下面的内容的时候配合后缀数组算法总结 | 远航休息栈一文中的若干示例,动手模拟一下,加深理解。
首先我们给出整个程序的代码。这里用struct简单封装了一下。

struct SuffixArray {
    char s[MAXN];
    int n, m, sa[MAXN], t[MAXN], t1[MAXN], c[MAXN];
    inline void init() {
        memset(this, 0, sizeof(SuffixArray));
    }
    inline void calsa() {
        int *x = t, *y = t1;
        for(int i = 1; i <= m; i++) c[i] = 0;
        for(int i = 1; i <= n; i++) c[x[i] = s[i]]++;
        for(int i = 1; i <= m; i++) c[i] += c[i - 1];
        for(int i = n; i >= 1; i--) sa[c[x[i]]--] = i;
        for(int k = 1; k <= n; k <<= 1) {
            int p = 0;
            for(int i = n - k + 1; i <= n; i++) y[++p] = i;
            for(int i = 1; i <= n; i++) if(sa[i] > k) y[++p] = sa[i] - k;
            for(int i = 1; i <= m; i++) c[i] = 0;
            for(int i = 1; i <= n; i++) c[x[y[i]]]++;
            for(int i = 1; i <= m; i++) c[i] += c[i - 1];
            for(int i = n; i >= 1; i--) sa[c[x[y[i]]]--] = y[i];
            std::swap(x, y);
            p = x[sa[1]] = 1;
            for(int i = 2; i <= n; i++) {
                x[sa[i]] = (y[sa[i - 1]] == y[sa[i]] && y[sa[i - 1] + k] == y[sa[i] + k]) 
                    ? p : ++p;
            }
            if(p >= n) break;
            m = p;
        }
    }
};
// Code by KSkun, 2018/5

变量声明

n代表字符串长度,m代表字符集大小。s是字符串本体,内容从下标1开始。
sa在排序中是一个排名数组,下标为排名,值为该排名对应的后缀首字母位置(后称“后缀编号”)。严格意义上说,这个数组实际上保存了每一次排序的结果,但如遇本次排序仍然存在两个后缀排序关键字相同,则长度大的(编号小的)后缀在前。SA计算完成后,该数组就是SA本体。
tt1是临时数组,用于存放后面的xy的信息。
c是计数排序的计数数组,下标为计数的值,值为该值出现的次数。排序过程中,我们会对其求前缀和以便计算排名。
x是一个中间量数组,意义为得到第一关键字的大小,对于一次排序,下标为代表后缀编号,值为象征对应后缀编号第一关键字大小的值(事实上可以视作排名)。
y是一个中间量数组,意义为第二关键字的排名,对于一次排序,下标为排名,值为第二关键字对应的后缀的编号。

初始化

int *x = t, *y = t1;
for(int i = 1; i <= m; i++) c[i] = 0;
for(int i = 1; i <= n; i++) c[x[i] = s[i]]++;
for(int i = 1; i <= m; i++) c[i] += c[i - 1];
for(int i = n; i >= 1; i--) sa[c[x[i]]--] = i;

第1/8行分配tt1两个内存空间给x数组和y数组。
第2/9行清空用于计数排序的计数数组c
第3/10行,统计该次第一关键字(首字母)的出现次数。
第4/11行,计算前缀和。
第5/12行,计算第一次排序的sa。此次排序,仅区分首字母和后缀编号,首字母小的在前,首字母相同,后缀编号小(后缀长度大)的在前。

倍增:计算y数组

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

第2/15行,有一些后缀拆分后的第二关键字为空,这个循环将这些后缀的排名提到最前,以表示空字符的最高优先级。
第3/16行,对于满的第二关键字排序,利用某一后缀的后缀也是原串的后缀这一特点,由上一次的sa计算而来。考虑那些排名靠前的后缀,如果往前补倍增长度位,那么这个后缀对应的长度前缀就变成了第二关键字。如我们有后缀aaabaa,长度为2,计算abaa的第二关键字的时候,实际上查到的是aa的排名,并且往前补了2位。

倍增:计数排序计算x数组

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

第2/18行,依然是在统计各个量的出现次数。y[i]表示的是第二关键字排名为i的后缀编号,而x[i]拿到的则是第一关键字排序时,后缀编号为i的后缀对应的象征第一关键字大小的值(在初始化过程中,我们将该值设置为首字母的值了,因此,对初始化还有一种理解,即y[i] = i,并没对第二关键字排序)。因此x[y[i]]的值反映的是代表第一关键字大小的值。
第4/20行,对i逆序统计,其用意是:当第一关键字相同的时候,第二关键字排名大的,计算后sa的值也会更大。

倍增:计算下一次的x数组

std::swap(x, y);
p = x[sa[1]] = 1;
for(int i = 2; i <= n; i++) {
    x[sa[i]] = (y[sa[i - 1]] == y[sa[i]] && y[sa[i - 1] + k] == y[sa[i] + k]) ? p : ++p;
}

第1/21行,把旧的x数组值换到y数组上。
第2/22行,规定当前排名第一的后缀x值为1,以后的后缀x值不会小于1。
第4/24行,考虑什么时候一个后缀和前面的后缀下一次的第一关键字(=这一次的第一+第二关键字,因为是倍增)相同,只有这一次的第一第二关键字都相等才可以判定下一次的第一关键字相等,否则应该比上一次的大。

倍增:边界条件

if(p >= n) break;
m = p;

第1/27行,当本次排序中已经产生了不小于n种不同的第一关键字的值,排序便没有必要继续下去了,因为已经区分开了这些后缀。
第2/28行,本次排序更新了关键字的种类数,因此要更新m值。

例题:【P3809】【模板】后缀排序 – 洛谷

把上面的板子打上去就可以了。

// Code by KSkun, 2018/5
#include <cstdio>
#include <cctype>
#include <cstring>

#include <algorithm>

const int MAXN = 1000005;

struct SuffixArray {
    char s[MAXN];
    int n, m, sa[MAXN], t[MAXN], t1[MAXN], c[MAXN];
    inline void init() {
        memset(this, 0, sizeof(SuffixArray));
    }
    inline void calsa() {
        int *x = t, *y = t1;
        for(int i = 1; i <= m; i++) c[i] = 0;
        for(int i = 1; i <= n; i++) c[x[i] = s[i]]++;
        for(int i = 1; i <= m; i++) c[i] += c[i - 1];
        for(int i = n; i >= 1; i--) sa[c[x[i]]--] = i;
        for(int k = 1; k <= n; k <<= 1) {
            int p = 0;
            for(int i = n - k + 1; i <= n; i++) y[++p] = i;
            for(int i = 1; i <= n; i++) if(sa[i] > k) y[++p] = sa[i] - k;
            for(int i = 1; i <= m; i++) c[i] = 0;
            for(int i = 1; i <= n; i++) c[x[y[i]]]++;
            for(int i = 1; i <= m; i++) c[i] += c[i - 1];
            for(int i = n; i >= 1; i--) sa[c[x[y[i]]]--] = y[i];
            std::swap(x, y);
            p = x[sa[1]] = 1;
            for(int i = 2; i <= n; i++) {
                x[sa[i]] = (y[sa[i - 1]] == y[sa[i]] && y[sa[i - 1] + k] == y[sa[i] + k]) 
                    ? p : ++p;
            }
            if(p >= n) break;
            m = p;
        }
    }
} sa;

int main() {
    scanf("%s", sa.s + 1);
    sa.n = strlen(sa.s + 1); sa.m = 256;
    sa.calsa();
    for(int i = 1; i <= sa.n; i++) {
        printf("%d ", sa.sa[i]);
    }
    return 0;
}

应用

height

上文提到过,x数组的值实际上是后缀编号对应的排名,我们将其内容复制一份进入rnk数组待用。
下面给出height数组的定义:代表该后缀与前一排名后缀的最长公共前缀(LCP)长度,下标为后缀编号。以下将该数组简写为hei数组。
我们只需要对每个后缀找到紧接着后面排名的后缀,来统计相同的部分即可。需要注意的是,当我们向后枚举排名的时候,实际上可以看做比较的后缀各删去了第一位,可以用上次的k-1来作为基础,然后再比较前k-1位之后有没有更多的相同的部分,这样可以优化时间。
以下代码就是求hei数组的代码,紧接着求sa的代码。

memcpy(rnk, x, sizeof(rnk));
int k = 0;
for(int i = 1; i <= n; i++) {
    int j = sa[rnk[i] + 1];
    if(k) k--;
    if(!j) continue;
    while(s[i + k] == s[j + k]) k++;
    hei[rnk[i]] = k;
}

求前缀信息

把字符串逆序一下扔进去求个SA就好了。

求前缀的LCS或后缀的LCP

要求两个排名不相邻的后缀的LCP怎么办?我们可以把上面的hei数组建一个ST表,查询排名区间的最小值即可。hei数组本身的含义就是LCP长度,这个性质有一定的传递性。在这里不太能说的明白,倒是可以采用观察一个sa的例子验证一下的做法。
至于前缀的LCS,其实你发现逆序了以后就变成一个后缀LCP问题了。
这个问题的一个经典应用题目是:[NOI2016]优秀的拆分

参考资料



发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注

This site is protected by reCAPTCHA and the Google Privacy Policy and Terms of Service apply.

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据