「算法笔记」后缀数组

一、定义

后缀数组（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 以及 论文。