后缀数组 - Suffix Array 防忘笔记

Published on 2016-05-13

记得 VFK 曾经说过:“NOI试机的时候,有很多选手都在照着书打后缀数组,因为不理解。”。
的确,后缀数组真的很容易忘记啊,明明以前搞得很清楚,过一阵子又忘记了,一忘记,又容易写错,一写错,考场上就没有办法了。所以简单的写一下核心内容,方便用时回忆。

简介

KMP 字符串匹配算法,它可以在 O(m+n)O(m + n) 的时间内完成对模版串的匹配。然而,对于多模版串的匹配,每一次都要预处理,十分不高效。这时有 AC 自动机的出现,虽然可以多模匹配,但每次匹配仍然要遍历整个字符串,对于很长很长的字符串,例如 DNA 信息,有上亿的数据,这时候每一次查询遍历就不太科学了。这时就要预处理文本串。
假设文本串为 banana,可以在其末尾加一个字符 $(代表一个没在串中出现过的字符),稍后解释原因。把其所有后缀(banana$, anana$, nana$, ana$, na$, a$, $),放进一棵 Trie 树。

有了这个后缀 Trie,查找一个长度为 mm 的模版时只需要进行一次普通的 Trie 查找,时间复杂度为 O(m)O(m),如上图所示。
而在实际中,我们会把后缀 Trie 中没有分支的链合并在一起,得到所谓的后缀树(Suffix Tree),然而构造起来太麻烦,复杂难懂,而且容易写错,而它的后缀数组则是一个不错的替代品,它时间效率高,而且代码简单,不易写错。空间大概是后缀树的 15\frac{1}{5} ~ 13\frac{1}{3}
或许你们已经注意到了,在绘制上面的后缀 Trie 的时候,我故意把每个节点的所有子节点排好序,字典序小的在左边(规定 $ 最小),然后在叶子节点上标上了后缀首字母的下标。例如后缀 ana$ 所对应的叶子节点标有 33,为了叙述方便,我们把以下标 kk 开头的后缀(包含$),叫做后缀 kk,记作
我们从左到右排列出所有叶子的编号,就可以得到后缀数组 (Suffix Array),banana 的后缀数组为 ,它就是所有后缀按照字典序从小到大排序后的结果。

用后缀数组求解匹配

有了这个后缀数组,我们就可以通过二分搜索(利用后缀数组字典序的单调性)的方法,来在 O(mlogn)O(m\log n) 的时间内完成匹配。

int m, sa[maxn];
int cmp_suffix(const char *pattern, int p) {
    return strncmp(pattern, str + sa[p], m);
}
//在 Suffix Array 数组的(lb, ub) 区间寻找字符串 pattern
void find(const char* pattern, int lb, int ub) {
         m = strlen(pattern);
        if(cmp_suffix(pattern, lb + 1) < 0) return; //如果当前最小的字典序都大于模版串的字典序,肯定无法匹配
        if(cmp_suffix(pattern, ub - 1) > 0) return; //如果当前最大的字典序都小于模版串的字典序,肯定无法匹配
        while(ub - lb > 1) {
            int mid = (lb + ub) / 2;
            int res = cmp_suffix(pattern, mid);
            if(!res) {
                printf("Match! Index:%d\n", sa[mid]);
                find(pattern, lb, mid); //(lb, mid) 可能还有匹配
                find(pattern, mid, ub); //(mid, ub) 可能还有匹配
                return;
            }
            if(res < 0) ub = mid; //如果模版串的字典序比后缀 sa[mid] 小,那么解的范围变为 (lb, mid) 
            else lb = mid;    //如果模版串的字典序比后缀 sa[mid] 小,那么解的范围变为 (mid, ub) 
        }
}

strncmp(s1, s2, size) 的意思是比较是比较 s1,s2s_1, s_2 的前 size 个字符,如果相同返回 0,如果不同,则返回第一个不同位置 kks1[k]s2[k]s_1[k] - s_2[k]

计算 Suffix Array

提到 数组,就不得不说一说 数组。SA[i]\mathrm{SA}[i] 表示的是字典序排名为 ii 的后缀是谁(字典序越小的排名越靠前)。而 rank[i]\mathrm{rank}[i] 表示的是后缀 ii 所对应的排名是多少。事实上,它们是反函数的关系,即 SA[rank[i]]=i\mathrm{SA}[\mathrm{rank}[i]] = i
特别注意,字符串的长度是原串的长度 +1,因为有 $ 这个字符,而在实现的时候我们就用 0 来表示 $,显然这个后缀 $ 排名最靠前。

朴素的排序

那么如何计算 数组呢?根据定义,我们可以通过一次快速排序得到后缀数组,比较次数 O(nlogn)O(n\log n),单次比较时间 O(n)O(n),总体 O(n2logn)O(n^2\log n),时间复杂度太高,难以承受。
不用担心后缀长度不一致的比较问题,因为有 $ 字典序最小,所以总能在不越界的情况下完成比较。
不过这方法几乎不用,非常慢,但是在没有时间的时候还是可以 rush 一发拿点暴力分,毕竟可以音速实现。

倍增思想

下面介绍一种 Manber 和 Myers 发明的倍增算法,它的时间复杂度是 O(nlogn)O(n\log n)

static const int maxNode = 100000 + 3;
int sa[maxNode], rank[maxNode], n;
char str[maxNode];
void buildSA(int m) {
    static int cnt[maxNode], tmpSA[maxNode], rank1[maxNode], rank2[maxNode];
    n = strlen(str) + 1; str[n] = 0; //结尾的 $
    fill(cnt, cnt + m, 0);
    for (int i = 0; i < n; ++i) cnt[(int)str[i]]++;
    for (int i = 1; i < m; ++i) cnt[i] += cnt[i - 1];
    for (int i = 0; i < n; ++i) rank[i] = cnt[(int)str[i]] - 1;
    for (int l = 1; l < n; l <<= 1) {
        for (int i = 0; i < n; ++i)
            rank1[i] = rank[i], rank2[i] = i + l < n ? rank1[i + l] : 0;
        fill(cnt, cnt + n, 0);
        for (int i = 0; i < n; ++i) cnt[rank2[i]]++;
        for (int i = 1; i < n; ++i) cnt[i] += cnt[i - 1];
        for (int i = n - 1; i >= 0; --i) tmpSA[--cnt[rank2[i]]] = i;
        fill(cnt, cnt + n, 0);
        for (int i = 0; i < n; ++i) cnt[rank1[i]]++;
        for (int i = 1; i < n; ++i) cnt[i] += cnt[i - 1];
        for (int i = n - 1; i >= 0; --i) sa[--cnt[rank1[tmpSA[i]]]] = tmpSA[i];
        bool unique = true;
        rank[sa[0]] = 0;
        for (int i = 1; i < n; ++i) {
            rank[sa[i]] = rank[sa[i - 1]];
            if (rank1[sa[i]] == rank1[sa[i - 1]] && rank2[sa[i]] == rank2[sa[i - 1]]) unique = false;
            else rank[sa[i]]++;
        }
        if (unique) break;
    }
}

具体我不细讲,因为这文章只是用来回顾的。
一开始直接利用字符来计算初始的 ,接着开始倍增,以 为第一关键字,也就是对于一个后缀的起点,前 ll 个字符的排名, 为第二关键字,也就是后 ll 个字符的排名,特别注意一旦越界,直接给 00
如何体现第二关键字的作用呢?先计算出 ,字典序越是小的越排在 的前面,用这句话巧妙地实现(注意倒序循环):
for (int i = n - 1; i >= 0; --i) tmpSA[--cnt[rank2[i]]] = i;
再根据第一关键字排序的时候,就利用 的顺序来排序,我们也是要字典序越是小的越排在 的前面,同样这样实现:
for (int i = n - 1; i >= 0; --i) sa[--cnt[rank1[tmpSA[i]]]] = tmpSA[i];
由于 是大到小的,同一个位置越先访问,领取到的排名越靠后。为了保证相同的 下,越靠前的后缀越小,所以我们倒序循环,这也是两次倒序循环的原因。
后面没什么好说的了,注意在 互不相同的时候跳出是一个不错的优化。

DC3 什么的,还是算了吧。巨大的常数。

计算 Height

高度数组: 是一个一维数组,保存了在后缀数组中相邻两个后缀的最长公共前缀(Longest Common Prefix,LCP)长度。

表示存在的最大的 xx,满足对于任何 k[0,x)k\in [0,x),有

除了朴素比较外,可以 O(n)O(n) 求解,代码如下:

void getHeight() {
    int k = 0;
    for (int i = 0; i < n - 1; ++i) {
        if (k) --k;
        int j = sa[rank[i] - 1];
        while (str[i + k] == str[j + k]) k++;
        height[rank[i]] = k;
    }
}

原理是按照 ,也就根据原串的顺序计算。整个求解的核心是这样一个不等式:如果设 h[i]=height[rank[i]]h[i] = \mathrm{height}[\mathrm{rank}[i]] 的话,一定有 h[i]h[i1]1h[i] \ge h[i - 1] - 1,下面来证明一下:
如果设 Suffix(i1)\mathrm{Suffix(i - 1)} 的在后缀数组中的前一个是 Suffix(SA[rank[i1]1])=Suffix(k)\mathrm{Suffix(\mathrm{SA}[\mathrm{rank}[i - 1] - 1])} = \mathrm{Suffix(k)},那么我们已经求出它们的最长公共前缀是 h[i1]h[i - 1],考虑去掉它们的头一位,得到前缀 Suffix(i)\mathrm{Suffix(i)}Suffix(k+1)\mathrm{Suffix(k + 1)},显然它们的最长公共前缀 h[i1]1h[i - 1] - 1,而 Suffix(i)\mathrm{Suffix(i)} 的在后缀数组中的前一个是 Suffix(SA[rank[i]1])=Suffix(p)\mathrm{Suffix(\mathrm{SA}[\mathrm{rank}[i] - 1])} = \mathrm{Suffix(p)},是与 Suffix(i)\mathrm{Suffix}(i) 字典序最为接近的一个(因为 Suffix(k+1)\mathrm{Suffix(k + 1)} 中必定比 Suffix(i)\mathrm{Suffix(i)} 靠前,而 Suffix(p)\mathrm{Suffix(p)} 就是 Suffix(i)\mathrm{Suffix(i)} 之前一个),Suffix(i)\mathrm{Suffix(i)}Suffix(k+1)\mathrm{Suffix(k + 1)} 的最长公共前缀不可能大于 Suffix(i)\mathrm{Suffix(i)}Suffix(p)\mathrm{Suffix(p)} 的最长公共前缀的,所以它们的最长公共前缀至少是 h[i]h[i1]1h[i] \ge h[i - 1] - 1
这样做,就不用从头比较了。

求解 LCP

依照刚刚的理论,对于 Suffix(i)\mathrm{Suffix}(i),只有可能与其在 中相邻的元素有最长的最长公共前缀,如果要求 Suffix(i)\mathrm{Suffix}(i)Suffix(j)\mathrm{Suffix}(j) 的最长公共前缀,且 rank[i]<rank[j]\mathrm{rank}[i] < \mathrm{rank}[j],那么 LCP(Suffix(i),Suffix(j))=minheight[k],k[rank[i]+1,rank[j]]\mathrm{LCP}(\mathrm{Suffix}(i), \mathrm{Suffix}(j)) = \min \mathrm{height}[k], k \in[\mathrm{rank}[i] + 1, \mathrm{rank}[j]]。这是一环扣一环的过程。
问题转化为区间最值查询(Range Minimum/Maximum Query,RMQ)问题,可以使用稀疏表(Sparse Table,ST)算法解决。该算法在 的时间内预处理,并可在 O(1)O(1) 的时间内完成每个询问。

还有单次匹配 O(m+logn)O(m + \log n) 的算法以及一些常用技巧,等我碰到了题再补吧~
鸣谢 Menci 的 SuffixArray 讲解,部分懒得打直接搬过来了。