zoukankan      html  css  js  c++  java
  • 「算法笔记」后缀数组

    一、定义

    后缀数组(Suffix Array),简称 SA。

    后缀数组 是一个通过对字符串的所有后缀经过排序后得到的数组。形式化的定义如下:

    • 令字符串 (S=S_1 S_2 cdots S_n)(S[isim j]) 表示 (S) 的子串,下标从 (i)(j)(S)后缀数组 (A) 被定义为一个数组,内容是 (S) 的所有后缀经过字典排序后的起始下标。即,(forall 1<ileq n),有 (S[A_{i-1}sim n]<S[A_isim n]) 成立。

    例如 (S= ext{banana$}) 有后缀(假设 ( ext{$}) 的字典序小于任何字母):

    对所有后缀按字典序升序排序后:(可得 (A=[7,6,4,2,1,5,3])

    二、倍增算法求解

    1. 一些定义

    方便起见,我们定义:

    • ( ext{suffix}(i)) 表示第 (i) 个位置开始的后缀,即 (S[isim n])

    • (sa(i)) 表示将所有后缀排序后排名为 (i) 的后缀的起始下标。即排名为 (i) 的后缀为 ( ext{suffix}(sa(i)))

    • (rk(i)) 表示以起始下标为 (i) 的后缀的排名,也就是 ( ext{suffix}(i)) 的排名。显然有 (sa(rk(i))=i,rk(sa(i))=i)

    朴素的求后缀数组的方法是,直接对 (n) 个后缀进行排序,由于比较两个字符串是 (mathcal{O}(n)) 的,所以排序是 (mathcal{O}(n^2log n)) 的,这些不再赘述。

    2. 主要思路

    倍增算法的主要思路是:对所有长度为 (2^k) 的字符串进行排序,并求出它们的排名(即 (rk(i)) 数组)。当 (2^kgeq n) 后,每个位置开始的长度为 (2^k) 的子串就相当于所有的后缀了,此时的 (rk(i)) 数组就是答案(此时这些子串一定都已经比较出大小,即 (rk(i)) 数组中没有相同的值)。

    如何比较两个长度为 (2^k) 的串?先将所有的 (S[i]) 进行排序,然后每次通过 (S[isim i+2^{k-1}-1]) 的大小关系求出 (S[isim i+2^k-1]) 的大小关系(串的右端点省略了和 (n) 取最小值的操作)。

    具体来说,假如当前我们已经把长度为 (2^{k-1}) 的子串排好序并求出了 ({rk'}(i))。那么长度为 (2^k) 的子串可以用两个长度为 (2^{k-1}) 的子串的排名作为关键字表示,记为二元组 (({rk}'(i),{rk}'(i+2^{k-1}))),以该二元组的第一部分为第一关键字,第二部分为第二关键字;当 (i+2^{k-1}>n) 时,则令 (rk'(i+2^{k-1})) 的值为 (0),因为其对应子串为空串,字典序最小。

    得到这些二元组后,就可以将长度为 (2^k) 的子串进行排序,得到它们的 (rk(i)) 值。

    3. 具体实现

    用 sort 进行排序:(mathcal{O}(nlog^2 n))

    #include<bits/stdc++.h>
    #define int long long
    using namespace std;
    const int N=1e6+5;
    int n,m,sa[N],rk[N],t[N],p; 
    char s[N];
    bool cmp(int x,int y){
        return rk[x]==rk[y]?rk[x+m]<rk[y+m]:rk[x]<rk[y]; 
    }
    signed main(){
        scanf("%s",s+1),n=strlen(s+1);
        for(int i=1;i<=n;i++) rk[i]=s[i],sa[i]=i;
        sort(sa+1,sa+1+n,cmp);
        for(int k=1;k<=n;k<<=1){
            m=k,p=0,sort(sa+1,sa+1+n,cmp),memcpy(t,rk,sizeof(rk));    //t 是原先的 rk 数组。 
            for(int i=1;i<=n;i++)
                rk[sa[i]]=(t[sa[i]]==t[sa[i-1]]&&t[sa[i]+k]==t[sa[i-1]+k]?p:++p);    //去重操作。若两个子串相同,它们对应的 rk 也需要相同。
        }
        for(int i=1;i<=n;i++)
            printf("%lld%c",sa[i],i==n?'
    ':' ');
        return 0;
    } 

    在刚刚的 (mathcal{O}(nlog^2 n)) 做法中,用 sort 单次排序是 (mathcal{O}(nlog n)) 的,若能 (mathcal{O}(n)) 排序,就能 (mathcal{O}(nlog n)) 计算后缀数组了。

    考虑使用 基数排序。那么这个基数排序具体怎么实现呢?

    现在我们已经得到了长度为 (2^{k-1}) 的子串的排名,需要计算长度为 (2^k) 的子串的排名。

    定义 (t(i)) 表示 按第二关键字 排序后排名为 (i) 的后缀的起始位置,即 (t(i)) 满足以 (t(i)+2^{k-1}) 开头的长度为 (2^{k-1}) 的子串按字典序排列。(cnt(i)) 表示基数排序需要的桶。

    假设我们已经求出了 (t(i))。我们知道,(t(i)) 是已经以第二关键字排好序的,所以我们要做的就是,以第一关键字排序,并保持第二关键字的相对顺序不变,并把排序的结果存在 (sa(i)) 中。

    实现二元组排序的方式是,把第一关键字放在桶里,然后按第二关键字的顺序拿出来。具体地:

    • 首先用桶(即 (cnt(i)) 数组)统计第一关键字(也就是所有长度为 (2^{k-1}) 的子串)的每个排名的出现次数,并对其做前缀和,得到小于等于当前排名的子串个数。
    • 同一个桶内的元素的顺序由第二关键字就决定。倒序枚举第二关键字的排名 (i),它对应的第一关键字的排名为 (rk(t(i))),其对应的桶为 (cnt(rk(t(i)))),那么 (sa(cnt(rk(t(i))))=t(i)),然后 (cnt(rk(t(i)))=cnt(rk(t(i)))-1)。我们倒着枚举倒着拿走桶中的元素,这样就能够维护第二关键字的相对顺序了。

    时间复杂度:(mathcal{O}(nlog n))

    //Luogu P3809
    #include<bits/stdc++.h>
    #define int long long
    using namespace std;
    const int N=1e6+5;
    int n,m,sa[N],rk[N],t[N],cnt[N],p,tot; 
    char s[N];
    void rsort(){    //基数排序,解释过了 
        for(int i=1;i<=m;i++) cnt[i]=0;
        for(int i=1;i<=n;i++) cnt[rk[i]]++;
        for(int i=1;i<=m;i++) cnt[i]+=cnt[i-1]; 
        for(int i=n;i>=1;i--) sa[cnt[rk[t[i]]]--]=t[i];
    }
    signed main(){
        scanf("%s",s+1),n=strlen(s+1),m=max(n,(int)'z');    //m: 桶的大小 
        for(int i=1;i<=n;i++) rk[i]=s[i],t[i]=i;
        rsort();
        for(int k=1;k<=n;k<<=1){
            tot=p=0;
            for(int i=n-k+1;i<=n;i++) t[++tot]=i;    //这些后缀的第二关键字为空串,因此排在最前面
            for(int i=1;i<=n;i++)
                if(sa[i]-k>0) t[++tot]=sa[i]-k;    //sa[i]-k>0,它可以作为别的子串第二关键字
            rsort(),swap(t,rk);    //swap(t,rk): 此时的 t 已经没用了,而我们还要更新 rk,所以 t 变为上一次的 rk 用来备份 
            for(int i=1;i<=n;i++)
                rk[sa[i]]=(t[sa[i]]==t[sa[i-1]]&&t[sa[i]+k]==t[sa[i-1]+k]?p:++p);
        }
        for(int i=1;i<=n;i++)
            printf("%lld%c",sa[i],i==n?'
    ':' ');
        return 0;
    } 

    事实上,还有不常见的 DC3 算法,时间复杂度为线性,但常数比较大,实际运行时间与倍增法相当。

    三、求最长公共前缀

    求原串任意两个后缀的 最长公共前缀(LCP)

    1. 一些定义

    对于后缀 ( ext{suffix}(i))( ext{suffix}(j)),我们定义 ( ext{lcp}(i,j)=k),其中 (k) 为满足 (forall 1leq pleq k, ext{suffix}(i)_p= ext{suffix}(j)_p) 的最大值。

    接下来再定义一些值:

    • ( ext{LCP}(i,j)= ext{lcp}(sa(i),sa(j))),即排名为 (i,j) 的后缀的最长公共前缀。

    • (height(i)= ext{LCP}(i-1,i)),即排名为 (i-1,i) 的后缀的最长公共前缀。显然有 (height(1)=0)(排名为 (1) 的字符串和空串的 ( ext{lcp}) 显然为 (0))。

    • (h(i)=height(rk(i))),即 ( ext{suffix}(i)) 与它前一个排名的后缀的最长公共前缀。

    我们的最终目的是求解 (height(i)) 数组,原因会在后面说明。

    两个基本的等式:( ext{LCP}(i,j)= ext{LCP}(j,i), ext{LCP}(i,i)=n-sa(i)+1)

    2. LCP 的性质

    在证明接下来的内容之前,先证明一个引理(性质 (1)

    (forall 1leq i<k<jleq n, ext{LCP}(i,j)=min( ext{LCP}(i,k), ext{LCP}(k,j)))

    证明:(p=min( ext{LCP}(i,k), ext{LCP}(k,j))),则有 ( ext{LCP}(i,k)geq p, ext{LCP}(k,j)geq p)

    ( ext{suffix}(sa(i))=u, ext{suffix}(sa(k))=v, ext{suffix}(sa(j))=w),则 (u,v) 的前 (p) 个字符相等,(v,w) 的前 (p) 个字符相等。这样得到 (u,w) 的前 (p) 个字符相等,即 ( ext{LCP}(i,j)geq p)

    我们考虑反证法。假如 ( ext{LCP}(i,j)geq p+1),则有 (u_{p+1}=w_{p+1}),由于排名 (i<k<j),那么 (u_{p+1}=v_{p+1}=w_{p+1}),这与条件矛盾。故 ( ext{LCP}(i,j)=p)

    推论:根据上面的引理,可以得到一个推论(性质 (2)):

    (forall 1leq i<jleq n, ext{LCP}(i,j)=minlimits_{i<kleq j} ext{LCP}(k-1,k))

    证明:结合引理可得:( ext{LCP}(i,j)=min( ext{LCP}(i,i+1), ext{LCP}(i+1,j)))

    可以将 ( ext{LCP}(i+1,j)) 可以继续拆下去,正确性显然。

    根据这个推论,由于 (height(i)= ext{LCP}(i-1,i)),所以 ( ext{LCP}(i,j)=minlimits_{i<kleq j}height(k))。如果能求出 (height) 数组,那么 ( ext{lcp}) 问题就转化为了一个 区间最小值问题,显然可以通过 ST 表 解决。

    这也是求解 (height(i)) 数组的原因。

    4. 最重要的性质

    通过上述分析,我们还是没法求 (height) 数组。于是我们要证明一个 最重要的性质:

    (h(i) geq h(i-1)-1)

    证明:

    • (h(i-1)=0),结论显然成立。

    • 否则在 ( ext{suffix}(i-1)) 与它前一个排名的后缀的最长公共前缀的基础上,去掉第一个字符,就找到了一个后缀与 ( ext{suffix}(i)) 的最长公共前缀 (geq h(i-1)-1)。进而有 (height(rk(i))geq height(rk(i-1))-1),即 (h(i)geq h(i-1)-1)

    较严谨的证明(考虑到公式加载慢的问题,把它删掉了 QAQ)

    5. 代码实现

    (height) 数组的代码:

    for(int i=1,k=0;i<=n;i++){
        if(k) k--;    //h[i]>=h[i-1]-1 
        int j=sa[rk[i]-1];    //h[i]=height[rk[i]]=LCP(rk[i]-1,rk[i])=lcp(sa[rk[i]-1],sa[rk[i]])=lcp(sa[rk[i]-1],i)
        while(i+k<=n&&j+k<=n&&s[i+k]==s[j+k]) ++k;    //因为 h[i]>=h[i-1]-1,所以可以直接从 k 开始算 
        ht[rk[i]]=k;    //h[i]=height[rk[i]], 这里 height 简写成了 ht 
    }

    (k) 不会超过 (n),最多减 (n) 次,所以最多加 (2n) 次,总复杂度就是 (mathcal{O}(n))

    后缀数组求  ( ext{lcp}) 完整代码:

    #include<bits/stdc++.h>
    #define int long long
    using namespace std;
    const int N=1e5+5;
    int n,m,q,sa[N],rk[N],t[N],cnt[N],p,tot,ht[N],f[N][30],x,y;
    char s[N];
    void rsort(){
        for(int i=1;i<=m;i++) cnt[i]=0;
        for(int i=1;i<=n;i++) cnt[rk[i]]++;
        for(int i=1;i<=m;i++) cnt[i]+=cnt[i-1];
        for(int i=n;i>=1;i--) sa[cnt[rk[t[i]]]--]=t[i];
    }
    void getht(){    //求 height 数组 
        for(int i=1,k=0;i<=n;i++){
            if(k) k--;
            int j=sa[rk[i]-1];
            while(i+k<=n&&j+k<=n&&s[i+k]==s[j+k]) ++k;
            ht[rk[i]]=k;
        }
    } 
    void getst(){
        for(int i=1;i<=n;i++) f[i][0]=ht[i];
        for(int j=1;j<=log2(n);j++)
            for(int i=1;i+(1<<j)-1<=n;i++)
                f[i][j]=min(f[i][j-1],f[i+(1<<(j-1))][j-1]);
    }
    int query(int x,int y){
        if(x==y) return n-x+1;
        x=rk[x],y=rk[y];    //将 lcp 变为 LCP 
        if(x>y) swap(x,y);x++;    //LCP(i,j)=height(k) (i<k<=j),k 取不到 i 
        int k=log2(y-x+1);
        return min(f[x][k],f[y-(1<<k)+1][k]);
    }
    signed main(){
        scanf("%s",s+1),n=strlen(s+1),m=max(n,(int)'z');
        for(int i=1;i<=n;i++) rk[i]=s[i],t[i]=i;
        rsort();
        for(int k=1;k<=n;k<<=1){
            tot=p=0;
            for(int i=n-k+1;i<=n;i++) t[++tot]=i;
            for(int i=1;i<=n;i++)
                if(sa[i]-k>0) t[++tot]=sa[i]-k;
            rsort(),swap(t,rk);
            for(int i=1;i<=n;i++)
                rk[sa[i]]=(t[sa[i]]==t[sa[i-1]]&&t[sa[i]+k]==t[sa[i-1]+k]?p:++p);
        }
        getht(),getst(),scanf("%lld",&q);
        while(q--){
            scanf("%lld%lld",&x,&y);
            printf("%lld
    ",query(x,y));
        }
        return 0;
    } 

     (附写到一半然后被删掉的 height 数组的应用,建议跳过这个直接去看论文 qwq)

    四、其他应用

    可参考 后缀数组 (SA) -  OI Wiki 以及 论文

  • 相关阅读:
    idea gson使用
    spring对象工具类
    java反射判断对象空字段
    说说沟通乱这件事
    RandomShuffleQueue
    学习材料
    python异常处理
    tensorboard基础使用
    深度学习——特殊应用:人脸识别和神经风格转换[13]
    深度学习——目标检测[12]
  • 原文地址:https://www.cnblogs.com/maoyiting/p/14187838.html
Copyright © 2011-2022 走看看