后缀数组

 

参考了:victorique - 后缀数组 最详细讲解

比较基础的知识点,虽然有不少玩法

快速过掉后认真学SAM

(mrfz周年池沉了,难受)

(UPD:我出货了,我又好了)

 


 

~ 后缀数组(SA)的定义 ~

 

给定一个长度为$n$的字符串$s$,将其所有后缀$s[i...n]$按照字典序大小进行排序

我们的目标是求两个数组$sa,rk$

$sa[i]$表示,字典序从小到大第$i$个字符串的起点(那么该字符串为$s[sa[i]...n]$)

$rk[i]$表示,以$i$为起点的字符串$s[i...n]$在所有后缀中的排名(rank)

显然若已知$rk[i]$,那么$sa[i]$就能直接得到,故真正进行求解的是$rk$数组

 


 

~ 求rk数组 ~

 

显然不能使用常规的字符串排序来求$sa$,否则时间复杂度为$O(n^2logn)$

所以需要用到两个trick来将复杂度降到$O(nlogn)$

 

Part 1 - 倍增

对于所有子串暴力比较是不可取的,我们需要利用一下特殊性质:后缀

这有什么用处呢?我们可以像Manacher、EXKMP一样,利用【已经求得】的信息来简化计算

在对两个字符串进行字典序比较的过程中,我们需要从前向后依次比较;而在这里也是一样的,我们需要先对所有后缀的第一个字符进行比较

例如,令$s$=abacab,那么可以给所有$s[i]$进行排名(大小相同则取同排名)

i

1 2 3 4 5 6
s[i] a b a c a b
排名 1 2 1 3 1 2

但是一遍排名过后依然有排名相同的元素;而$s$的所有后缀均不相同(毕竟长度不同),显然排名应该互不相同;所以需要将所有子串的第二个字符纳入考虑,再次进行排名

这一次排名,我们可以利用上一次求得的排名进行计算

比如,$i=1,s[1]=a$的排名是$1$,$i=2,s[2]=b$的排名是$2$,那么我们可以用一个二元组$1/2$来表示$s[1,2]$的大小;对二元组进行比较时,先比较第一维度,再比较第二维度

对于$i=6,s[6]=b$,它的第二维度应当设定为$0$,因为长度小的子串字典序更小

i 1 2 3 4 5 6
s[i] a b a c a b
二元组 1/2 2/1 1/3 3/1 1/2 2/0
排名 1 4 2 5 1 3

进行一遍排名后,可以看出排名相同的元素减少了

我们在下一次求排名的时候,就不是将所有子串的第三个字符纳入考虑了,而是将第$3,4$位;因为在第二轮排名结束后,求得的排名是以$i$为起点的、长度为$2$的子串的排名,所以我们在下一轮利用这个排名的时候,相当于直接获得了两个字符的大小信息

相应地,这次构造二元组的时候,并不是将下一个位置的排名作为第二维度,而是将下下个位置

比如$i=1,s[1,2]=ab$的排名是$1$,$i=3,s[3,4]=ac$的排名是$2$,那么构造的二元组就是$1/2$;这是将$s[1,2]$与$s[3,4]$拼起来、以表示$s[1...4]=abac$的大小

对于$i=5,i=6$,其第二维度为$0$

i 1 2 3 4 5 6
s[i] a b a c a b
二元组 1/2 4/5 2/1 5/3 1/0 3/0
排名 2 5 3 6 1 4

发现排名互不相同,于是可以结束了

如果字符串的长度更大,就是不断地将倍增进行下去

第$i$轮是将从每个位置开始、长度为$2^i$的子串进行排名;而通过我们的倍增转化,其实只需要对$n$个二元组进行比较

总复杂度为$O(n\cdot (logn)^2)$,因为需要比较$logn$轮,而每轮需要对$n$个二元组进行排序

 

Part 2 - 桶排序优化

在上面的进行排名的过程中,我们仅需要比较$n$个二元组

那么就可以利用桶排序降低时间复杂度;共有$n+1$个桶表示排名

接着就是标准的桶排序操作了:先按照第二维度将二元组放到桶中,再从小到大地从桶中取出二元组得到一个顺序,最后按顺序依照第一维度将二元组放到桶中,即可完成排序

于是二元组排序的时间复杂度变成了$O(n)$,那么总的时间复杂度就降到了$O(nlogn)$

 


 

~ 代码实现 ~

 

虽然上面说的逻辑比较清晰,但代码实现会混乱一些(为了简短)

在具体实现的过程中,我们并不是单独求$rk$数组(否则需要开数组表示二元组及对应关系),而是在$rk,sa$之间反复横跳

先介绍一下需要用的数组

int sa[N],rk[N];
int tmp[N<<1],top[N];

$rk,sa$:和定义中的意义相同

$top$:桶排时用到的桶(更准确地说是用来占位的,下面会具体说明)

$tmp$:表示二元组的第二维度中 从小到大第$i$个元素在字符串中的位置(可以理解成对于第二维度的$sa$);开成两倍是为了防止在后面的代码中数组越界(为了简化一个判断条件,可能需要访问$2N$以内的$tmp[i]$)

 

先康康桶排是如何实现的:

//n: 待排元素数  m: 桶数
void quicksort(int n,int m)
{
    for(int i=1;i<=m;i++)
        top[i]=0;
    for(int i=1;i<=n;i++)
        top[rk[i]]++;
    for(int i=1;i<=m;i++)
        top[i]=top[i-1]+top[i];
    for(int i=n;i>=1;i--)
        sa[top[rk[tmp[i]]]--]=tmp[i];
}

这个实现和桶排的一般实现不太一样,但思想是相同的

先看看$top$数组做了什么事:先统计第一维度($rk[i]$)出现了多少次,然后做了一次前缀和;这样相当于第一维度为$i$的元素占据了 排序后$n$个位置中 的一段连续区间,区间的右端点是$top[i]$

这样一来,如果我们对于第一维度相同的元素,将第二维度按从大到小的次序依次加入,并让$top[i]--$,就能恰好将二元组排好序

具体看最后一个for循环:$tmp[i]$是第二维度的$sa$,那么从$n$循环到$1$就是 按照第二维度从大到小 依次获得二元组的下标;$top[rk[tmp[i]]]$就是这个二元组 第一维度所在区间 的右端点

 

剩下的部分也需要慢慢想一下:

void getsa(char *s,int n)
{
    for(int i=1;i<=n;i++)
        rk[i]=s[i],tmp[i]=i;
    quicksort(n,200);
    
    //开始倍增
    for(int i=1;i<n;i<<=1)
    {
        int cnt=0;
        for(int j=n-i+1;j<=n;j++)
            tmp[++cnt]=j;
        for(int j=1;j<=n;j++)
            if(sa[j]>i)
                tmp[++cnt]=sa[j]-i;
        
        quicksort(n,max(cnt,200));
        
        for(int j=1;j<=n;j++)
            swap(rk[j],tmp[j]);
        
        rk[sa[1]]=cnt=1;
        for(int j=2;j<=n;j++)
        {
            if(tmp[sa[j]]!=tmp[sa[j-1]] || tmp[sa[j]+i]!=tmp[sa[j-1]+i])
                cnt++;
            rk[sa[j]]=cnt;
        }
        
        if(cnt==n)
            break;
    }
}

一开始对$rk,tmp$做初始化:由于我们在第一遍排序只关注每个子串首字符的大小,所以$tmp$取任意一个$1$到$n$的排列即可

排完序以后,我们得到了正确的$rk,sa$(仅需要关注$rk$之间的相对大小,并不需要在意其具体值;在这里$rk$只是不从$1$开始,而相对大小是正确的),可以开始倍增了

在每一轮倍增中,我们都是利用上一轮的$rk,sa$构造二元组;每个二元组的第一维度是$rk[j]$,第二维度是$rk[j+i]$:

若$j+i>n$,那么该二元组所对应字符串长度小于$i$(即不存在后半部分),在第二维度的排序中就具有最高优先度,应放在$tmp$数组中的前面;否则就顺着$sa$数组从小到大放,此时需要注意放在$tmp$中的是$sa[j]-i$,即二元组所对应字符串的起点(枚举的$j$表示字符串后半部分的起点)

这样已经足够我们进行排序了,桶排一波得到正确的$sa$;而此时的$rk$仍然是上一轮的,所以我们需要参照$sa$重新计算一次$rk$

$tmp$已经完成了它的使命,所以不妨将上一轮的$rk$放到$tmp$中,并利用它计算当前轮的$rk$

我们求得的$sa$虽然是正确的,但是对于相等的二元组 相当于强行给它们定了一个顺序,而它们的$rk$应该是相等的;于是依次判断是否相等即可(这里在判断第二维度是否相等时,下标可能会超出$n$,故$tmp$需要开成两倍空间)

倍增过程中的$cnt$表示桶的数量,即第一维度的数量;当$cnt=n$时,$rk$已经互不相同了,即可结束倍增

 

模板题:Luogu P3809 (【模板】后缀排序)

若不让$cnt=n$时break,就会超时

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

const int N=1000005;

int sa[N],rk[N];
int tmp[N<<1],top[N];

void quicksort(int n,int m)
{
    for(int i=1;i<=m;i++)
        top[i]=0;
    for(int i=1;i<=n;i++)
        top[rk[i]]++;
    for(int i=1;i<=m;i++)
        top[i]=top[i-1]+top[i];
    for(int i=n;i>=1;i--)
        sa[top[rk[tmp[i]]]--]=tmp[i];
}

void getsa(char *s,int n)
{
    for(int i=1;i<=n;i++)
        rk[i]=s[i],tmp[i]=i;
    quicksort(n,200);
    
    for(int i=1;i<n;i<<=1)
    {
        int cnt=0;
        for(int j=n-i+1;j<=n;j++)
            tmp[++cnt]=j;
        for(int j=1;j<=n;j++)
            if(sa[j]>i)
                tmp[++cnt]=sa[j]-i;
        
        quicksort(n,max(cnt,200));
        
        for(int j=1;j<=n;j++)
            swap(rk[j],tmp[j]);
        
        rk[sa[1]]=cnt=1;
        for(int j=2;j<=n;j++)
        {
            if(tmp[sa[j]]!=tmp[sa[j-1]] || tmp[sa[j]+i]!=tmp[sa[j-1]+i])
                cnt++;
            rk[sa[j]]=cnt;
        }
        
        if(cnt==n)
            break;
    }
}

int n;
char s[N];

int main()
{
    scanf("%s",s+1);
    n=strlen(s+1);
    
    getsa(s,n);
    for(int i=1;i<=n;i++)
        printf("%d ",sa[i]);
    return 0;
}
View Code

 


 

~ 拓展:求LCP ~

 

$LCP(i,j)$指的是两个后缀$s[sa[i]...n],s[sa[j]...n]$的最长公共前缀长度

 

性质1:$LCP(i,k)=min\{LCP(j,j+1)\}$,其中$i\leq j<k$

继续以之前的字符串$s=$abacab为例

i sa[i] s[sa[i]...n]
1 5 a
2 1 abacab
3 3 acab
4 6 b
5 2 bacab
6 4 cab

可以通过观察发现,$LCP(i,k)$会随着$k$的增大而单调不增(推论1

那么最小的$LCP(j,j+1)$相当于一个分界点:在它之前的字符串与$s[sa[i]...n]$拥有更多的公共前缀,而在它之后字符串的至少在第$LCP(j,j+1)+1$位与$s[sa[i]...n]$不同,那么$LCP(i,k)$显然不会超过$LCP(j,j+1)$;而$LCP(i,k)$不小于$min\{LCP(j,j+1)\}$是显然的

 

性质1把求$LCP(i,k)$转化成了求$min\{LCP(j,j+1)\}$,那么我们显然需要提前将$LCP(i,i+1),1\leq i<n$预处理出来

记$height[i]=LCP(i,i+1)$,$h[i]=height[rk[i]]$

可以这样理解这两个数组:$height[i]$对应的是子串$s[sa[i]...n]$与$s[sa[i+1]...n]$的最长公共前缀长度,而$h[i]$对应的是子串$s[i...n]$与$s[sa[rk[i]+1]...n]$的公共前缀长度;即前者的下标指的是在排序后第$i$小的子串,而后者的下标指的就是子串$s[i...n]$

性质2:$h[i+1]\geq h[i]-1$(证明有点绕,但最好还是要搞懂)

可以这样证明:假设将所有后缀进行排序后,$s[i...n]$的后面是$s[j...n]$,那么有$h[i]=height[rk[i]]=LCP(s[i...n],s[j...n])=L$;而考虑到$s[i+1...n]$在排序后的子串中应该也很接近$s[j+1...n]$,因为前$L-1$位是相等的;那么若在排序后的后缀中$s[i+1...n]$的后面是$s[j+1...n]$,则$h[i+1]=height[rk[i+1]]=LCP(s[i+1...n],s[j+1...n])=L-1$

根据推论1,$LCP(i,k)$随着$k$的增大而单调不减;那么若在排序后的后缀中$s[i+1...n]$的后面不是$s[j+1...n]$,即$rk[i+1]+1<rk[j+1]$,就有$h[i+1]=LCP(s[i+1...n],s[sa[rk[i+1]+1]...n)\geq LCP(s[i+1...n],s[j+1...n])=L-1$

综上,有$h[i+1]\geq h[i]-1$

 

有了强力的性质2,我们可以$O(n)$地预处理出$height$数组($h$数组除了推出一个性质2,实际上用不到)

考虑从$1$到$n$依次算出$height[rk[i]]$(即$h[i]$)

根据$height[rk[i+1]]\geq height[rk[i]]-1$(即$h[i+1]\geq h[i]-1$),能够发现$\sum_{i=1}^{n} (height[rk[i+1]]-height[rk[i]])\leq n+n$,那么每次均暴力进行子串比较就行了

int height[N];

void getheight(char *s,int n)
{
    int h=0;
    for(int i=1;i<=n;i++)
    {
        if(rk[i]==n)//height[n]=0
            continue;
        h=(h>0?h-1:h);
        
        int j=sa[rk[i]+1];
        while(j+h<=n && i+h<=n && s[i+h]==s[j+h])
            h++;
        height[rk[i]]=h;
    }
}

 


 

~ 后缀数组的应用 ~

 

上面写了那么多其实都是模板,想要图一乐还得看具体题目

 

Luogu SP705  ($New\ Distinct\ Substrings$)

经典应用——求一个字符串中本质不同的子串的数量

字符串中的每一个子串$s[l...r]$都是某个后缀$s[l...n]$的前缀,按照这个思路尝试考虑

$s[l...n]$这个后缀共有$n-l+1$个前缀,在这些前缀中,若有一些之前被统计过了,那么它们对答案的贡献就是$0$;如何判断是否被统计过呢?可以利用$LCP(i,i+1)$

我们顺着$sa[i]$进行考虑,这样一来所有的后缀都是按照字典序从小到大排列;假设$LCP(i,i+1)=x$,那么$s[sa[i+1]...n]$的前$x$个前缀都在之前被统计过了(至少在$s[sa[i]...n]$中已被统计),而其余的前缀都未被统计过(因为均与之前其他后缀的前$x+1$个字符不同)

于是,答案就是$\frac{n(n+1)}{2}-\sum_{i=1}^{n} height[i]$

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

typedef long long ll;
const int N=50005;

int sa[N],rk[N];
int tmp[N<<1],top[N];

void quicksort(int n,int m)
{
    for(int i=1;i<=m;i++)
        top[i]=0;
    for(int i=1;i<=n;i++)
        top[rk[i]]++;
    for(int i=1;i<=m;i++)
        top[i]=top[i-1]+top[i];
    for(int i=n;i>=1;i--)
        sa[top[rk[tmp[i]]]--]=tmp[i];
}

void getsa(char *s,int n)
{
    for(int i=1;i<=n;i++)
        rk[i]=s[i],tmp[i]=i;
    quicksort(n,200);//可能需要改成最大字符值 
    
    for(int i=1;i<n;i<<=1)
    {
        int cnt=0;
        for(int j=n-i+1;j<=n;j++)
            tmp[++cnt]=j;
        for(int j=1;j<=n;j++)
            if(sa[j]>i)
                tmp[++cnt]=sa[j]-i;
        
        quicksort(n,max(cnt,200));
        
        for(int j=1;j<=n;j++)
            swap(rk[j],tmp[j]);
        
        rk[sa[1]]=cnt=1;
        for(int j=2;j<=n;j++)
        {
            if(tmp[sa[j]]!=tmp[sa[j-1]] || tmp[sa[j]+i]!=tmp[sa[j-1]+i])
                cnt++;
            rk[sa[j]]=cnt;
        }
        
        if(cnt==n)
            break;
    }
}

int height[N];

void getheight(char *s,int n)
{
    int h=0;
    for(int i=1;i<=n;i++)
    {
        if(rk[i]==n)//height[n]=0
            continue;
        h=(h>0?h-1:h);
        
        int j=sa[rk[i]+1];
        while(j+h<=n && i+h<=n && s[i+h]==s[j+h])
            h++;
        height[rk[i]]=h;
    }
}

int n;
char s[N];

int main()
{
    int T;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%s",s+1);
        n=strlen(s+1);
        
        getsa(s,n);
        getheight(s,n);
        
        ll ans=1LL*n*(n+1)/2;
        for(int i=1;i<n;i++)
            ans-=height[i];
        printf("%lld\n",ans);
    }
    return 0;
}
View Code

 

(待续)

posted @ 2020-05-04 23:23  LiuRunky  阅读(210)  评论(0编辑  收藏  举报