「算法笔记」后缀数组

一、定义

后缀数组（Suffix Array），简称 SA。

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

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

例如 $S=\text{banana\$}$ 有后缀（假设 $\text{\$}$ 的字典序小于任何字母）：

「算法笔记」后缀数组

对所有后缀按字典序升序排序后：（可得 $A=[7,6,4,2,1,5,3]$）

「算法笔记」后缀数组

二、倍增算法求解

1. 一些定义

方便起见，我们定义：

$\text{suffix}(i)$ 表示第 $i$ 个位置开始的后缀，即 $S[i\sim n]$。
$sa(i)$ 表示将所有后缀排序后排名为 $i$ 的后缀的起始下标。即排名为 $i$ 的后缀为 $\text{suffix}(sa(i))$。
$rk(i)$ 表示以起始下标为 $i$ 的后缀的排名，也就是 $\text{suffix}(i)$ 的排名。显然有 $sa(rk(i))=i,rk(sa(i))=i$。

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

2. 主要思路

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

如何比较两个长度为 $2^k$ 的串？先将所有的 $S[i]$ 进行排序，然后每次通过 $S[i\sim i+2^{k-1}-1]$ 的大小关系求出 $S[i\sim 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}(n\log^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?'\n':' ');
    return 0;
}

在刚刚的 $\mathcal{O}(n\log^2 n)$ 做法中，用 sort 单次排序是 $\mathcal{O}(n\log n)$ 的，若能 $\mathcal{O}(n)$ 排序，就能 $\mathcal{O}(n\log 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}(n\log 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?'\n':' ');
    return 0;
}

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

Tips：

//swap(t,rk) 有时非常慢！必要时 for(int i=1;i<=n;i++) swap(t[i],rk[i]);
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);    //sa[i]+k<=n&&sa[i-1]+k<=n，多组数据的时候要注意 
//if(p==n) break; //一个小优化

Update：感觉基数排序的部分讲的还不太清楚，来补充一下（参考 tls 博客）。

先考虑一维的基数排序（假如我们要对 $a$ 排序），我们记一个桶 $c_i$ 表示值为 $i$ 的 $a_i$ 有多少个，然后对 $c_i$ 做一遍前缀和，得到 $\leq i$ 的数的个数。那么如果 $a_i$ 互不相同的话，排名为 $c_{a_i}$ 的数就是 $a_i$。而如果有 $a_i$ 相同，则显然排名在 $(c_{a_i-1},c_{a_i}]$ 之间的数都等于 $a_i$。考虑从后往前扫一遍，每次将 $a_i$ 的排名设为 $c_{a_i}$ 并令 $c_{a_i}$ 减 $1$，这样所有相同的 $a_i$ 都被赋上了一个不同的排名，并且对于相同的 $a_i$，位置越靠后的排名越靠后，也就是我们对 $(a_i,i)$ 这个二元组排了一遍序。

考虑我们在刚才的基数排序中为什么要从后往前扫，因为这样一来如果有 $a_i$ 相同的那么位置靠后的排名永远比位置靠前的排名高，也就是 按第二维降序 的顺序更新桶。那么二维也是一样的。

具体来说，我们要对 $rk_i$ 按第二关键字排序，记 $c_i$ 表示值为 $i$ 的 $rk_i$ 有多少个，然后对 $c_i$ 做一遍前缀和，再记一个 $t_i$ 表示按第二关键字排序后排名为 $i$ 的后缀的起始位置，从后往前扫一遍，将普通情况下的“设 $i$ 的排名为 $cnt_{rk_i}$，再令 $cnt_{rk_i}$ 减 $1$”，改为“设 $t_i$ 的排名为 $cnt_{rk_{t_i}}$，再令 $cnt_{rk_{t_i}}$ 减 $1$”，也就是将 $i$ 换成 $t_i$。

三、求最长公共前缀

求原串任意两个后缀的 最长公共前缀（LCP）。

1. 一些定义

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

接下来再定义一些值：

$\text{LCP}(i,j)=\text{lcp}(sa(i),sa(j))$，即排名为 $i,j$ 的后缀的最长公共前缀。
$height(i)=\text{LCP}(i-1,i)$，即排名为 $i-1,i$ 的后缀的最长公共前缀。显然有 $height(1)=0$（排名为 $1$ 的字符串和空串的 $\text{lcp}$ 显然为 $0$）。
$h(i)=height(rk(i))$，即 $\text{suffix}(i)$ 与它前一个排名的后缀的最长公共前缀。

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

两个基本的等式：$\text{LCP}(i,j)=\text{LCP}(j,i),\text{LCP}(i,i)=n-sa(i)+1$。

2. LCP 的性质

在证明接下来的内容之前，先证明一个引理（性质 $1$，称为 LCP Lemma）

$\forall 1\leq i<k<j\leq n,\text{LCP}(i,j)=\min(\text{LCP}(i,k),\text{LCP}(k,j))$

证明：设 $p=\min(\text{LCP}(i,k),\text{LCP}(k,j))$，则有 $\text{LCP}(i,k)\geq p,\text{LCP}(k,j)\geq p$。

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

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

推论：根据上面的引理，可以得到一个推论（性质 $2$，称为 LCP Theorem）：

$\forall 1\leq i<j\leq n,\text{LCP}(i,j)=\min\limits_{i<k\leq j}\text{LCP}(k-1,k)$

证明：结合引理可得：$\text{LCP}(i,j)=\min(\text{LCP}(i,i+1),\text{LCP}(i+1,j))$。

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

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

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

4. 最重要的性质

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

$h(i) \geq h(i-1)-1$

证明：

若 $h(i-1)=0$，结论显然成立。
否则在 $\text{suffix}(i-1)$ 与它前一个排名的后缀的最长公共前缀的基础上，去掉 第一个字符，就找到了一个后缀与 $\text{suffix}(i)$ 的最长公共前缀 $\geq h(i-1)-1$。进而有 $height(rk(i))\geq height(rk(i-1))-1$，即 $h(i)\geq h(i-1)-1$。

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)$。

后缀数组求 $\text{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\n",query(x,y));
    }
    return 0;
}

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

四、其他应用

先讲一下后缀数组求一个子串在原串中的出现次数：

求 $s$ 的子串 $s[l\sim r]$ 在 $s$ 中的出现次数，就等价于求有多少个 $\text{suffix}(j)$ 满足 $\text{lcp}(j,l)\geq r-l+1$。

考虑将 $n$ 个后缀放在字典序排名数组 $rk_i$ 上，显然对于所有 $\text{lcp}(j,l)\geq r-l+1$，它们的 $rk$ 值应道是一段连续的区间 $[L,R]$，而这个区间又可以通过二分 + ST 表求出。具体来说，二分找到最小的满足 $\min\limits_{i=t+1}^{rk_l} height_i\geq r-l+1$ 的 $t$ 记为该区间的左端点 $L$；同理找到最大的满足 $\min\limits_{i=rk_l+1}^t height_i\geq r-l+1$ 的 $t$ 即为该区间的左端点 $R$。然后返回 $R-L+1$ 即可，正确性显然。

时间复杂度 $\mathcal O(\log n)$（不考虑预处理的话）。

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

Upd on 2021.7.27：更新了一些应用，详见「笔记」2021.7.27 后缀数组，密码可以找我要鸭。