后缀数组和高度数组 模板及应用

后缀数组

特别注意: 下面的代码都是来自挑战程序设计竞赛, lcp[i]代表的是后缀Suffix[sa[i]]后缀Suffix[sa[i+1]]的最长公共前缀, 而网上比较流行的那个模板的height[i]则是后缀Suffix[i]与Suffix[sa[i-1]]的最长公共前缀, 挑战上的模板速度慢一些, 但是好写很多, 并且因为用的是快速排序数值范围比较大的时候不需要离散化

原理

为了方便处理一些情况, 我们把空串也当做一个后缀, 所以长度为n的字符串有n+1个后缀.
sa[i] 表示以S[sa[i]]开头的后缀是字符串所有后缀中字典序第i的.
rank[i]表示以S[i]开头的字符串在字符串所有后缀中排第rank[i]位.
用倍增法计算后缀数组.
记rk[i, j]为子串S[i…i+2j-1]在所有长度为2j长度的子串中的字典序排名, 如果i+2j-1超过了字符串的末尾, 则为S[i, n-1].
如果我们知道了所有的rk[i, j], 那么, 所有的rk[i, j+1]可以通过rk[i, j]和rk[i+2j, j]快速计算出来.
而一开始的rk[i, 0]是长度为1的子串, 可以直接通过每个字符的编码值得到(rk[i, 0] = int(S[i]);).

代码

效率O(nlog2n)

int rk[maxn], tmp[maxn], sa[maxn], k, n;bool cmp_sa(int i, int j){    if (rk[i] != rk[j]) return rk[i] < rk[j];    return (i + k <= n ? rk[i + k] : -1) < (j + k <= n ? rk[j + k] : -1);}void construct_sa(char *s){    for (int i = 0; i <= n; ++i)    {        sa[i] = i;        rk[i] = s[i];    }    rk[n] = -1;    for (k = 1; k <= n; k *= 2)//这里的k是全局变量, 不要在这里定义    {        sort(sa, sa + n + 1, cmp_sa);        tmp[sa[0]] = 0;        for (int i = 1; i <= n; ++i)            tmp[sa[i]] = tmp[sa[i - 1]] + (cmp_sa(sa[i - 1], sa[i]) ? 1 : 0);        copy(tmp, tmp + n + 1, rk);    }}

应用

基于后缀数组的字符串匹配

通过预处理目标串的后缀数组, 得到所有后缀的字典序排名, 然后就可以用二分搜索目标串在所有后缀中的位置得到匹配情况.

代码

效率: 预处理O(|S|log2|S|), 匹配O(|T|log|S|)

int contain_l(char *s, char *p){    int l = 0, r = strlen(s), len = strlen(p);    while (l < r)    {        int mid = (l + r) / 2;        if (strncmp(s + sa[mid], p, len) < 0) l = mid + 1;        else r = mid;    }    return strncmp(s+sa[l], p, len)>=0 ? l : l+1;}int contain_r(char *s, char *p){    int l = 0, r = strlen(s), len = strlen(p);    while (l < r)    {        int mid = (l + r) / 2;        if (strncmp(s + sa[mid], p, len) <= 0) l = mid + 1;        else r = mid;    }    return strncmp(s+sa[l], p, len)>0 ? l : l+1;}...int l =contain_l(s, p), r = contain_r(s, p);sa[l], sa[l+1], ..., sa[r-1]为所有匹配的开始位置

高度数组

原理

高度数组是有后缀数组中相邻两个后缀的最长公共前缀的长度组成的数组.
lcp[i]表示后缀S[sa[i]…]和后缀S[sa[i+1]…]的最长公共前缀的长度位lcp[i]
利用类似于尺取法的思想可以快速求出高度数组.
从0开始, 计算后缀S[i…]与后缀S[sa[rank[i]-1]…]的lcp[i], sa[rank[i]-1]是后缀数组中在后缀S[i…]前面一位的后缀, 因为第一个后缀数组中第一个后缀是空串, 所有非空后缀一定有前一个后缀, 所以不用担心sa[rank[i]-1]不存在.
然后求S[i+1…]和S[sa[rank[i+1]-1…]的lcp[i+1]时, 而S[i+1..]相当于S[i…]去掉第一个字母, 而S[sa[rank[i+1]-1…]虽然不一定是S[sa[rank[i]-1]+1…], 但是, 他们前面lcp[i]-1一定是相同的, 所以计算前缀只需要从i+1+lcp[i]-1开始比较.

代码

复杂度: 计算SA数组O(nlog2n)+计算LCP数组O(n)

int lcp[maxn];void construct_lcp(char *s){    for (int i = 0; i <= n; ++i) rk[sa[i]] = i;    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;    }}

应用

任意两个后缀的最长公共前缀

LCP数组中存的是在SA数组相邻后缀的公共前缀, 但有一个规律: 后缀S[i...], S[j...](rank[i]<rank[j])的公共前缀长度为min{lcp[rank[i]], lcp[rank[i]+1], ... , lcp[rank[j]-1]}.
所所以只要利用RMQ算法就可以快速求出任意两个后缀的公共前缀.

最长公共子串(POJ - 2217)

原理

先考虑一个串中出现过两次以上的子串最大长度, 很显然, 出现过两次以上的子串一定是SA数组中相邻两个后缀的前缀, 所以要找最长的一定是高度数组的最大值.

而两个字符串则可以通过拼接成一个字符串来简化成上面一个串的形式, 但是拼接时候要插入在两个字符串中间插入一个不会出现的字符, 这样就能保证我们求出的公共子串一定是只来自于一个串的.

还要判断SA数组中相邻两个后缀是来自两个不同的串的

代码

int main(){    int T;    scanf("%d", &T);    getchar();    while(T--)    {        cin.getline(a, maxn);        cin.getline(b, maxn);        int la = strlen(a), lb = strlen(b);        a[la] = '\0';        copy(b, b + lb, a + la + 1);        n = la + lb + 1;        construct_sa(a, sa);        construct_lcp(a, sa, lcp);        int ans = 0;        for (int i = 0; i < n; ++i)        {            if ((sa[i] < la) != (sa[i + 1] < la))                ans = max(ans, lcp[i]);        }        printf("Nejdelsi spolecny retezec ma delku %d.\n", ans);    }    return 0;}

最长回文串

原理

求回文串有一种更简单的算法: litmxs的博客 - 最长回文串 Manacher

奇数长度的回文串是 从i开始的后缀 与 反转后字符串中从原来的i位子开始的后缀 的公共前缀
偶数长度类似
所以和最长公共子串类似的思路, 将原串S与反转后的S中间加入一个不会出现的字符拼接在一起, 对得到的字符串S’求高度数组
以原来S[i]为中心的奇数长度的回文串长度 = 后缀S’[i…]和后缀 S’[2n-i…]的公共前缀
以原来S[i-1]和S[i]为中心的回文串长度 = 后缀S’[i…]和后缀S’[2n-i+1…]的公共前缀

代码

#include <iostream>#include <cstdio>#include <cstring>#include <algorithm>using namespace std;const int maxn = 2e5 + 100;int a[maxn], n;int mi[maxn][30];void init_rmq(int a[], int n){    for (int i = 0; i < n; ++i) mi[i][0] = a[i];    for (int j = 1; (1 << j) <= n; ++j)    {        for (int i = 0; i + (1 << j) <= n; ++i)        {            mi[i][j] = min(mi[i][j - 1], mi[i + (1 << (j - 1))][j - 1]);        }    }}int rmq(int l, int r){    int k = 0;    while ((1 << (k + 1)) <= r - l + 1) ++k;    return min(mi[l][k], mi[r - (1 << k) + 1][k]);}int rk[maxn], tmp[maxn], sa[maxn], k, lcp[maxn];bool cmp_sa(int i, int j){    if (rk[i] != rk[j]) return rk[i] < rk[j];    return (i + k <= n ? rk[i + k] : -1) < (j + k <= n ? rk[j + k] : -1);}void construct_sa(char *s){    for (int i = 0; i <= n; ++i)    {        sa[i] = i;        rk[i] = s[i];    }    rk[n] = -1;    for (k = 1; k <= n; k *= 2)    {        sort(sa, sa + n + 1, cmp_sa);        tmp[sa[0]] = 0;        for (int i = 1; i <= n; ++i)            tmp[sa[i]] = tmp[sa[i - 1]] + (cmp_sa(sa[i - 1], sa[i]) ? 1 : 0);        copy(tmp, tmp + n + 1, rk);    }}void construct_lcp(char *s){    for (int i = 0; i <= n; ++i) rk[sa[i]] = i;    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;    }}char s[maxn];int main(){    scanf("%s", s);    int len = strlen(s);    s[len] = '$';    reverse_copy(s, s + len, s + len + 1);    n = len * 2 + 1;    s[n] = '\0';    construct_sa(s);    construct_lcp(s);    init_rmq(lcp, n + 1);    int ans = 0;    for (int i = 0; i < len; ++i)    {        int j = len * 2 - i;        int l = rmq(min(rk[i], rk[j]), max(rk[i], rk[j]));        ans = max(ans, 2*l-1);    }    for (int i = 1; i < len; ++i)    {        int j = len * 2 - i + 1;        int l =rmq(min(rk[i], rk[j]), max(rk[i], rk[j]));        ans = max(ans, 2*l);    }    cout << ans << endl;    return 0;}

模板

这个模板要求有效数据后面一个位置必须是0, 所有数字要大于0, height[i]表示sa[i]和sa[i-1]的公共前缀长度

/**suffix array*倍增算法  O(n*logn)*待排序数组长度为n,放在0~n-1中，在最后面补一个0*build_sa( ,n+1, );//注意是n+1;*getHeight(,n);*例如：*n   = 8;*num[]   = { 1, 1, 2, 1, 1, 1, 1, 2, $ };注意num最后一位为0，其他大于0*rank[]  = { 4, 6, 8, 1, 2, 3, 5, 7, 0 };rank[0~n-1]为有效值，rank[n]必定为0无效值*sa[]    = { 8, 3, 4, 5, 0, 6, 1, 7, 2 };sa[1~n]为有效值，sa[0]必定为n是无效值*height[]= { 0, 0, 3, 2, 3, 1, 2, 0, 1 };height[2~n]为有效值**/int sa[MAXN];//SA数组，表示将S的n个后缀从小到大排序后把排好序的             //的后缀的开头位置顺次放入SA中int t1[MAXN],t2[MAXN],c[MAXN];//求SA数组需要的中间变量，不需要赋值int rank[MAXN],height[MAXN];//待排序的字符串放在s数组中，从s[0]到s[n-1],长度为n,且最大值小于m,//除s[n-1]外的所有s[i]都大于0，r[n-1]=0//函数结束以后结果放在sa数组中void build_sa(int s[],int n,int m){    int i,j,p,*x=t1,*y=t2;    //第一轮基数排序，如果s的最大值很大，可改为快速排序    for(i=0;i<m;i++)c[i]=0;    for(i=0;i<n;i++)c[x[i]=s[i]]++;    for(i=1;i<m;i++)c[i]+=c[i-1];    for(i=n-1;i>=0;i--)sa[--c[x[i]]]=i;    for(j=1;j<=n;j<<=1)    {        p=0;        //直接利用sa数组排序第二关键字        for(i=n-j;i<n;i++)y[p++]=i;//后面的j个数第二关键字为空的最小        for(i=0;i<n;i++)if(sa[i]>=j)y[p++]=sa[i]-j;        //这样数组y保存的就是按照第二关键字排序的结果        //基数排序第一关键字        for(i=0;i<m;i++)c[i]=0;        for(i=0;i<n;i++)c[x[y[i]]]++;        for(i=1;i<m;i++)c[i]+=c[i-1];        for(i=n-1;i>=0;i--)sa[--c[x[y[i]]]]=y[i];        //根据sa和x数组计算新的x数组        swap(x,y);        p=1;x[sa[0]]=0;        for(i=1;i<n;i++)            x[sa[i]]=y[sa[i-1]]==y[sa[i]] && y[sa[i-1]+j]==y[sa[i]+j]?p-1:p++;        if(p>=n)break;        m=p;//下次基数排序的最大值    }}void getHeight(int s[],int n){    int i,j,k=0;    for(i=0;i<=n;i++)rank[sa[i]]=i;    for(i=0;i<n;i++)    {        if(k)k--;        j=sa[rank[i]-1];        while(s[i+k]==s[j+k])k++;        height[rank[i]]=k;    }}