后缀数组(SA),高度数组(LCP)及其应用

先导知识

字符串的大小比较

字符串的大小比较,通常是按照字典序来比较的,正如其名,按照字符在字母表中出现的顺序,从左右到右依次比较。如果比较到任意一个字符串的末尾仍然未出现较大者或者较小者,则字符长度短的小
例如,$a<b,ab<ac,ab<abc$

子串

字符串 $S$ 的子串 $r[i..j]$,$i<j$,表示 $r$ 串中从 $i$ 到 $j$ 这一段, 也就是顺次排列$ r[i]$,$r[i+1]$......$r[j]$形成的字符串。

后缀

后缀是指从某个位置$ i $开始到整个串末尾结束的一个特殊子串,字符串$ r $的从 第$ i $个字符开始的后缀表示为 $Suffix(i)$ ,也就是 $Suffix(i)=r[i..len(r)]$

简而言之,字串是起点和终点均在字符串中的一段连续的更小的字符串,而后缀是起点在字符串中,终点在字符串的末尾的特殊的字串。

后缀数组(Suffix Array):

对一个字符串的所有的后缀按字典序进行排序,所得到的数组叫做后缀数组SA(suffix array)

例如,字符串$ABAADCB$的后缀数组为

AADCB
ABAADCB 
ADCB
B
BAADCB
CB
DCB

Rank数组

和字符数组相反,给定后缀的下标,返回其字典序,$RANK(SA[i])=i$
例如:

sa[2]=5;//字典序第二小的字串是从下标为2处的字符开始的。
rk[5]=2;//下标从5开始的字串的字典序为2

$rk$数组和$sa$数组互为逆运算。

后缀数组(SA)的构造

一,朴素法

朴素法构造SA的思路很简单,将字符串的所有字串放入一个字符串数组,然后对其进行排序即可。然而由于C++两个字符串的比较是每个字符从左到右依次比较的,所以其本身会耗费$O(N)$的复杂度,然后再加上排序的复杂度,即使使用快排,总复杂度也高达$O(n^2lgn)$,如此高的复杂度,很明显是无法满足我们解题的要求的。故此在一般的算法竞赛中我们一般使用倍增法构造后缀数组。

struct comp{
    bool operator()(const pair<int,string> &a,const pair<int,string> &b)//由于使用了pair所以需要自定义comp谓词用于排序
    {
        return a.second<b.second;//朴素的comp很简单,使用pair中的两个字符串进行比较就行了,string的比较已在string类中实现了运算符重载
    }
};
string str;
vector<pair<int,string> >sa//构造一个pair用于存放下标和字串
cin>>str;
for(int i=0;i<str.length();i++)
{
    sa.push_back(make_pair(i,str.substr(i)));//获取所有字串
}
sort(sa.begin(),sa.end(),comp());

二,倍增法

倍增法一般常指二倍增算法,除此之外还有三倍增算法,由于不常使用,故此不在本文讨论范围内。
倍增算法的主要思路是:用倍增的方法对每个字符开始的长度为$2k$的子字
符串排序,求出排名,即 $rk$ 值。然后再使用上一轮的$rk$值,进行新一轮排序,长度变为$2k$。$k$ 从$ 0$ 开始,每轮$+1$,当 $2k$大于$ n $以 后,每个字符开始的长度为 $2k$的子字符串便相当于所有的后缀。

例如下图为字符串$aabaaaab$的$SA$的构造过程;
第$0$次排序,初始化$SA$中的值为字符串的每一个字符。比较大小按照其字典序列进行排序,将其排名存放在$rk$数组中。此时$K=0$
第$1$次排序,此时$K++$,这个时候从字符串的每个字符开始,向后扩大两倍。那么我们该如何确定此时的字典序呢。如果每个字符逐位比较那么则和朴素法无异,此时由于我们已经有了上一轮的排名,即首字符的排名,那么我们便可以先根据上一轮的$rk$比较第一个字符,如果不等,那么直接出结果不用再比较,如果相等我们再比较第二个字符。此时$K=1$;

然后将其排序更新$rk$。由此往复,直到$2K>=N$,停止倍增。此时的$SA$即为后缀数组。
需要注意的是$SA$中存储的一直都是每个字符的下标,图中所示只是为了方便理解。

在代码实现中,我们使用$STL$库中封装的$sort$,与朴素方法不同的是我们需要在自定义$COMP$谓词并且将倍增的思想体现在其中,下面我们通过代码来分析。

#include <iostream>
#include <algorithm>
#include <string>
#include <vector>
using namespace std;

vector<int>rk;
vector<int>tmp;
int n, k;//n:字符串字符数,k:循环变量
struct comp {
    bool operator()(const int& i, const int& j)//用于比较的二元谓词
    {
        if (rk[i] != rk[j])return rk[i] < rk[j];//第一个字典序如果不一样则直接比较
        else//一样则判断第二个首字符
        {
            int ni = i + k <= n ? rk[i + k] : -1;//如果不够倍增,则长度短的子串必然字典序小,所以规定为-1
            int nj = j + k <= n ? rk[j + k] : -1;
            return ni < nj;
        }
    }
};

void construct_sa(const string& s, vector<int>& sa)//创建sa
{
    n = s.size();//初始化n的值为字符串的长度
    for (int i = 0; i <= n; i++)//初始化SA
    {
        sa[i] = i;
        rk[i] = i < n ? s[i] : -1;//-1代表空串
    }
    for (k = 1; k <= n; k *= 2)//注意k要<=n
    {
        sort(sa.begin(), sa.end(), comp());//使用comp进行排序
        tmp[sa[0]] = 0;//使用tmp临时储存新rank
        for (int i = 1; i <= n; i++)tmp[sa[i]] = tmp[sa[i - 1]] + (comp()(sa[i - 1], sa[i]) ? 1 : 0);//如果两个字符串相同则字典序不增加
        for (int i = 0; i <= n; i++)rk[i] = tmp[i];//更新rank
    }
}

int main()
{
    string s = "abracadabra";
    string target = "dabr";
    rk.resize(s.size() + 1);//根据str的大小为rk分配空间
    tmp.resize(s.size() + 1);
    vector<int>sa(s.size() + 1);
    construct_sa(s, sa);
    //for_each(sa.begin(), sa.end(), print<int>());
    cout << bs(sa, target, s, 0, s.size() - 1);
    return 0;
}

倍增算法的时间复杂度,决定于最长公共子串的长度,最坏情况下,排序次数为 $logn$ 次,所以总的时间复杂度为 $O(nlogn)$。

高度数组(LCP Array,Longest Common Prefix Array)

高度数组是指后缀数组中相邻两个后缀的最长公共前缀的长度组成的数组,例如后缀$S[sa[i]..]$和后缀$S[sa[i+1]..]$的最长公共前缀的长度记作为$lcp[i]$,求$LCP$的思路非常简单,直接两两比较即可,但仔细思索一下这种情况求$LCP$的复杂度为$O(N^2)$。显然是无法接受的,下面介绍一个$O(N)$的方法,不过在此之前我们先来了解以下字符串和后缀数组的一些特点。
对于字符串S,有两个后缀,起始下标分别为$j$和$k$,设$rk[j]<rk[k]$,即后缀$J$的字典序小于$K$,则有以下性质:
$suffix(j)$和$suffix(k)$的最长公共前缀为$LCP[RK[J]+1]$,$LCP[RK[j]+2]$,$LCP[RK[J]+3]$...$LCP[RK[K]]$中的最小值
例如字符串为$aabaaaab$,后缀$abaaaab$和后缀$aaab$的最长公共前缀如图所示:

相信大家不难理解为什么,因为字典序的比较都是从左到右的比较,由图我们可以看出,其关系类似于一种多重集合取交集的模式,所以$i$和$j$的$LCP$必然会大于最小值,在此处也就是$1$。


有了上面的铺垫,我便可以开始介绍如何高效的求出高度数组了。
首先我们定义$LCP[i]=LCP[rank[i]]$为$suffix(i)$和在它前一名的后缀的最长公共前缀,那么则有$LCP[i]>=LCP[i-1]-1$
以下为证明:
设$suffix(k)$是排在$surffix(i)$前一名的后缀,则他们的最长公共前缀是$LCP[i]$(这里的$LCP[i]>1$,因为如果$LCP[i]<=1$,那么$LCP[i]>=0$必然成立),那么$suffix(k+1)$将会排在$suffix(i+1)$的前面(因为这两个后缀有公共前缀,后移相当于删去首位字符,删去公共前缀并不会影响他们的字典序),所以$suffix(k+1)$和$suffix(i+1)$的最长公共前缀是$LCP[i-1]-1$,上面讲到两个后缀的最长共前缀必然大于在这两个后缀字典序之间的$LCP$的最小值,那么换句话说,两个字典序之间的$LCP$的最小值必然大于该两个后缀的$LCP$.所以$suffix(i)$和他前一名的后缀的公共前缀至少是$LCP[i-1]-1$,所以由此按照$LCP[1]$,$LCP[2]$...$LCP[N]$的顺序计算,并且利用这一性质,时间复杂度可以降为$O(n)$,代码如下:

为了方便,在代码中$LCP[i]$为后缀$i$和比后缀$i$大$1$的后缀的最长公共前缀

void construct_lcp(string s, vector<int>& lcp, const vector<int>& sa)
{
    n = s.length();//初始化n的值为字符串的长度
    int h = 0;//记录上一个子串的LCP,h一直增加,因为下一个后缀lcp必然大于等于前者
    lcp[0] = 0;//初始化特殊情况
    for (int i = 0; i < n; i++)//从字符串的首位后缀开始依次计算。
    {
        int j = sa[rk[i] - 1];//j是比以i开始的后缀字典序少1的后缀,(即SA中的前一位)
        if (h > 0)h--;//应用上面讲的特性
        for (; j + h < n && i + h < n; h++)//直接从h开始比较,因为前面都一样
        {
            if (s[j + h] != s[i + h])break;//如果出现不同则跳出循环
        }
        lcp[rk[i] - 1] = h;//记录lcp值
    }
}

由于$SA$和$LCP$的构造较为繁琐并且比较难以理解,故此建议理解后背模板

后缀数组的应用

1.使用SA和二分查找实现$O(lgn)$字符串匹配

首先在讲这个应用前先讲一个公理吧

子串一定是某个后缀的前缀

简单的来说,就是一个字符串的任意一段子串都会是该字符串的某一个后缀的前缀,这很好理解,实在不理解的去举几个例子也就明白了,这里不再多讲。
由于SA的后缀是有序的,通过比较待查找字符串和SA中中位的字典序的大小,便可以判断待查找字符串处于拿一侧,再通过递归不断缩小这个范围,最后检索出这个子串所在的前缀。对于二分法还存在问题的童鞋可以去看看我之前写的关于递归的文章。

需要注意的是,如果后缀长于待搜索子串,比较的时候需要截取后缀长度和子串长度相等!例如,$aab和aa$不会匹配,即使aab为子串aa存在的后缀。

int bs(const vector<int>& sa, const string& t, const string& s, int sp, int ep)
{

    // 递归形式
    if (sp > ep)return -1;//找不到返回-1
    int mid = sp + ((ep - sp) >> 2);//通过位运算代替除法计算中间值
     if (s.compare(sa[mid], t.size(), t) < 0)return bs(sa, t, s, mid + 1, ep);//调用string的compare比较字符串
    else if (s.compare(sa[mid], t.size(), t) > 0)return bs(sa, t, s, sp, mid - 1);
    else return mid;

}

int bs(const vector<int>& sa, const string& t, const string& s)//t为待搜索子串,s为完整的字符串
{
    //循环形式
    int st = 0, ed = s.size();//st为开始下标,ed为结束下标
    while (ed - st> 1)
    {
        int mid = (st + ed) >> 1;//使用左移代替除以二
        int res = s.compare(sa[mid], t.size(), t);
        //使用string的compare方法比较截取后的后缀(s中起始下标为SA[MID]长度为待查找子串的长度)和待查找子串,compare返回三种情况,若待查找子串字典序小于后缀,则返回值大于0,如果匹配则返回值等于0,如果待查找子串字典序大于后缀,则返回值小于0.
        if (res < 0)//待比较字符串大于后缀,在后面
            st  = mid;
        else if (res == 0)//查找成功
        {
            cout << "match at: " << sa[mid] << endl;
            return sa[mid];
        }
        else rhs = mid;//待比较字符串小于后缀,在前面
    }
}

2.最长公共子串问题

让我们来看看如何应用$SA$解决最长公共子串问题。

最长公共子串(POJ 2217)
给定两个字符串$S$和$T$,请计算两个字符串最长的公共字符串子串的长度。
限制条件: $1\leq |S|,|T|\leq 10000$
样例输入:ABRACADABRA ECADADABRBCRDAR
样例输出:ADABR 5

注意区别子序列和子串的区别,子串要求是连续的,而子序列则不需要,只需要保证不调换顺序即可

思路

这道题可以简化为求一个字符串中出现两次的子串,那么很明显后缀数组中相邻的高度数组不为0的两个后缀的公共前缀都符合条件,不过题目要求最大,那么很明显,满足题意的是高度数组中的最大值。至于题目给的是两个字符串,很简单,将两个字符串合并即可。不过需要注意的是合并的时候需要加上一个分割符,因为这两个子串必须来自两个字符串,在这里我使用'0'

#include <iostream>
#include <string>
#include <vector>
#include <algorithm>
using namespace std;
int n, k;
vector<int>rk;
vector<int>tmp;
vector<int>sa;
vector<int>lcp;
string s;
struct comp {
    bool operator()(const int& i, const int& j)
    {
        if (rk[i] != rk[j])return rk[i] < rk[j];
        else
        {
            int ni = i + k <= n ? rk[i + k] : -1;
            int nj = j + k <= n ? rk[j + k] : -1;
            return ni < nj;
        }
    }
};
void construct_sa()
{
    for (int i = 0; i <= n; i++)
    {
        sa[i] = i;
        rk[i] = i < n ? s[i] : -1;
    }
    for (k = 1; k <= n; k *= 2)
    {
        sort(sa.begin(), sa.end(), comp());
        tmp[sa[0]] = 0;
        for (int i = 1; i <= n; i++)tmp[sa[i]] = tmp[sa[i - 1]] + (comp()(sa[i - 1], sa[i]) ? 1 : 0);
        for (int i = 0; i <= n; i++)rk[i] = tmp[i];
    }
}
void construct_lcp()
{
    int h = 0;
    lcp[0] = 0;
    for (int i = 0; i < n; i++)
    {
        int j = sa[rk[i] - 1];
        if (h > 0)h--;
        for (; j + h < n && i + h < n; h++)
        {
            if (s[j + h] != s[i + h])break;
        }
        lcp[rk[i] - 1] = h;
    }
}
int main()
{
    string s1, t1;
    cin >> s1 >> t1;
    s = s1 +'\0' t1;
    n = s.length();
    rk.resize(n + 1);
    tmp.resize(n + 1);
    sa.resize(n + 1);
    lcp.resize(n + 1);
    construct_sa();
    construct_lcp();
    int p = 0;
    for (int i = 0; i < lcp.size(); i++)
    {
        if (lcp[i] > p)
        {
            int a = sa[i], b = sa[i + 1];
            int fg = s1.size() + 1;
            if((a<fg&&b>fg)||(a>fg&&b<fg))p = i;
        }
    }
    cout << s.substr(sa[p], lcp[p])<< "  "<<lcp[p];

    return 0;
}

参考文献

IOI2009国家集训队论文 - 罗穗骞《后缀数组——处理字符串的有力工具》

Last modification:September 27th, 2019 at 05:18 pm
如果觉得我的文章对你有用,请随意赞赏

Leave a Comment