后缀自动机 Suffix-AutoMaton SAM
更新日志
2025/01/23:开工。2025/01/24:增加map版模板和应用
2025/01/24:增加求后缀数组
前言
本文不会提供很理性细节的证明,主要提供方法、逻辑以及实现,如有证明需要推荐:
本文对第二篇博客有所借鉴,在此鸣谢。
概念
在开始之前,我们先给出一些概念:
- SAM,也就是后缀自动机,是一个储存所有后缀信息(事实上,它储存了所有子串信息)的 DFA(确定性有限状态自动机)
- DFA,你可以视作一个 DAG,状态就是节点,转移就是边。
- 我们用从初始状态 \(root\) 至当前节点 \(u\) 的路径经过的字符组成的字符串代表当前状态。你会发现,按照这个定义,一个节点会代表多个状态,我们后面细说。
- 我们令 \(\mathrm{endpos}(s)\) 表示 \(s\) 子串在整个字符串中出现的结束位置集合(为方便,我们下文中从 \(1\) 开始给下标编号),显然,实际运用中它应该非空。
- 定义 \(\mathrm{endpos}\) 等价类,这是一个字符串集合,其内所有字符串的 \(\mathrm{endpos}\) 相同。我们将使用 \(\mathrm{endpos}(E)\) 表示等价类 \(E\) 的 \(\mathrm{endpos}\) 信息。
前置知识差不多就这些,开始了。
节点信息
我们在概念中提到,一个节点可能储存了多种状态,那究竟是哪些状态?
事实上,一个节点代表着一个 \(\mathrm{endpos}\) 等价类。这样,当我们从初始状态走到当前状态,就获得了当前子串的全部信息,包括出现次数与位置。
构造方案将会在后面提出,在此之前,请先看完所有需要的前置知识。
后缀链接
概念
对于一个节点 \(u\),我们令 \(\mathrm{long}(u)\) 表示 \(u\) 等价类中最长的状态(显然出现位置相同且长度相同的子串唯一),
则后缀链接 \(lnk_u\) 指向 \(\mathrm{long}(u)\) 最长的不在 \(u\) 等价类中的后缀 \(s\) 所属的等价类 \(q\),且 \(\mathrm{long}(q)=s\)。
这有什么意义呢?从 \(u\) 开始一直跳后缀链接直到根节点,你可以获取 \(\mathrm{long}(u)\) 的所有后缀。
性质
我们将简单的给出一些显然性质,并视情况给出感性证明。
- 一个节点中中子串长度是连续的。
显然。一个节点中的等价类相同,也就是结束位置完全一样,那显然内含所有子串长度就是一个长度的连续区间。
- 设 \(\mathrm{short}(u)\) 为 \(u\) 等价类中最短的后缀长度,则 \(\mathrm{short}(u)=\mathrm{long}(lnk_u)+1\)。
显然,它这个等价类的最短长度,就是最长的不是这个等价类的后缀长度再向左扩展一位。我们已经说明子串长度是连续的。
- \(\mathrm{endpos}(u)\) 完全被包含于 \(\mathrm{endpos}(lnk_u)\)。
显然,自己画一下图或者想一下。当前出现的位置其后缀必然会出现,但其后缀出现的位置当前字符串就不一定了。由于二者不在一个等价类内,所以完全包含。
差不多这些就够了。
构造
前置
是时候上手构造 SAM 了!
我们令 \(len_u\) 表示 \(\mathrm{long}(u)\) 的长度,\(lnk_u\) 表示 \(u\) 的后缀链接。
令 \(root=0\) 表示初始状态(根节点),代表空串,令 \(lnk_{root}\) 表示虚拟状态,存为 \(-1\),用于后续判断。
我们的算法是在线的,也就是一个字符一个字符插入,记每次插入到字符为 \(\texttt c\)。
流程
- 记录代表了当前整个字符串的状态 \(lst\),新建状态 \(cur\) 表示插入 \(\texttt c\) 之后新的整个字符串代表的状态,\(len_{cur}\gets len_{lst}+1\)。
- 从 \(lst\) 开始遍历后缀链接,记 \(p\) 为当前节点。如果 \(p\) 没有 \(\texttt c\) 的转移,则新建转移 \(p\overset{\texttt c}{\rightarrow}cur\)。
-
- 如果 \(p=root\),且仍然没有 \(\texttt c\) 的转移,新建转移完毕之后 \(lnk_{cur}\gets root\),跳转第 \(4\) 步。
- 否则,若 \(p\) 存在一条 \(\texttt c\) 标记的转移,记 \(q\) 为这条转移的对象状态。
- 若 \(len_q=len_p+1\),\(lnk_{cur}\gets q\),跳转第 \(4\) 步。
- 否则,我们分裂 \(q\)。
首先,我们复制 \(q\) 至 \(cpy\),\(len_{cpy}=len_p+1\)。
然后,\(lnk_p\gets cpy,lnk_{cur}\gets cpy\)。
接着,从 \(p\) 开始继续遍历后缀链接,若存在转移 \(p\overset{\texttt c}{\rightarrow}q\),则将这个转移重定向为 \(p\overset{\texttt c}{\rightarrow}cpy\)。直到当前 \(p\) 不存在 \(p\overset{\texttt c}{\rightarrow}q\),退出遍历。
- 更新 \(lst\gets cur\),结束。
原理
- 新的整个字符串的状态显然可以由原整个字符串的状态添加一个 \(\texttt c\) 得到。
- 新的状态代表的等价类的唯一出现位置显然就是新字符串的最后一个位置,对于字符串的每一个后缀,首先显然可以在其后面添加新字符串得到新的全字符串的一个后缀,如果原后缀只有原最后一个出现位置,那么新后缀显然也只有新的最后一个位置,所以可以直接指向 \(cur\)。
-
- 所有后缀均位于 \(cur\) 等价类中,不存在 \(lnk\),指向空。
- 否则,则说明 \(\mathrm{long}(p)+\texttt{c}\) 就是 \(\mathrm{long}(lnk_{cur})\)。我们用下面的例子理解一下:
\(\texttt{...abc...abc...ab + c}\)
假如 \(\texttt{...}\) 部分无法匹配,那么显然我们当前遇到的符合条件的节点 \(p=\mathrm{endpos}(\texttt{ab})\),它具有一个转移至 \(q=\mathrm{endpos}(\texttt{abc})\)。显然 \(\mathrm{long}(lnk_{cur})=\texttt{abc}\)。
值得注意的是,我们不能直接 \(lnk_{cur}\gets q\),考虑下面例子:
\(\texttt{...xabc...xabc...yab + c}\)
此时 \(\mathrm{long}(q)=\texttt{xabc}\) 而非 \(\texttt{abc}\),你会发现,加入 \(\texttt{c}\) 之后,\(\texttt{abc}\) 与 \(\texttt{xabc}\) 就不在同一个等价类了,所以我们需要分裂 \(p\) 状态。- 若 \(len_q=len_p+1\),\(lnk_{cur}\gets q\)。什么意思呢,也就是 \(\mathrm{long}(q)\) 恰好就是 \(\mathrm{long}(p)+\texttt{c}\),故而直接建立连接即可。
- 否则,我们分裂 \(q\)。
首先,我们复制 \(q\) 至 \(cpy\),\(len_{cpy}=len_p+1\)。这时候,\(cpy\) 状态内存储了长度小于等于 \(|\mathrm{long}(p)+\texttt{c}|\) 的 \(q\) 的所有后缀。
显然,因为他是从 \(q\) 分裂出来的后缀等价类,故而将 \(q\) 的后缀链接指向它。同时,这个等价类就是用来给 \(cur\) 链接的,所以也链上。
接着,我们需要更新 \(q\) 和 \(cpy\) 在 SAM 上的出入边以及后缀链接。- 首先,关于后缀链接,显然 \(cpy\) 继承了 \(q\) 的,易得。
- 然后,关于出边,不难想到 \(cpy\) 直接继承了 \(q\),它们原先就位于同一个等价类中,在后面加一个字符,能转移到的等价类显然与原先不变。
- 再然后,关于入边,首先我们明确,\(cpy\) 内所有字符串都是 \(\mathrm{long}(p)+\texttt{c}\) 的后缀,所以所有 \(\mathrm{long}(p)\) 的后缀加上 \(\texttt{c}\),如果原先指向了 \(q\) 的话,就应该重定向至 \(cpy\)。如果不再指向 \(q\) 了,就说明指向了别的等价类,那其后缀也不可能指回 \(q\) 了,退出即可。故而遍历 \(p\) 的后缀链接更新即可。
- 更新 \(lst\gets cur\),结束。没啥好说的。
Parent树
虽然整个 SAM 都是基于 DAG 的,但通常我们使用 SAM,都是使用其衍生出的这一棵 Parent 树。
不难发现所有后缀链接将会构成一棵树,这棵树就是 Parent 树。
其具体的运用将会在应用部分举例讲解。
复杂度
空间复杂度
点数最多为 \(2n-1\),故而空间开两倍。
边数最多为 \(3n-4\),可以预估空间复杂度,如果使用 map
的话。但会时间复杂度带一个 \(\log\)。
时间复杂度
不用 map
是 \(O(n)\) 的,否则 \(O(n\log n)\)。
模板
数组版
struct suffix_automaton{
struct node{
int nxt[26];
int len,lnk;
}ns[L<<1];
int siz,lst;
vec<int> g[L<<1];
void init(){
repl(i,0,L<<1)g[i].clear();
siz=lst=0;
ns[0].lnk=-1;
}
void insert(string s){
for(auto ch:s){
int cur=++siz,p=lst,c=ch-'a';
f[cur]=1;
ns[cur].len=ns[lst].len+1;
while(p!=-1&&!ns[p].nxt[c]){
ns[p].nxt[c]=cur;
p=ns[p].lnk;
}
if(p==-1)ns[cur].lnk=0;
else{
int q=ns[p].nxt[c];
if(ns[q].len==ns[p].len+1)ns[cur].lnk=q;
else{
int cpy=++siz;
ns[cpy]=ns[q];
ns[cpy].len=ns[p].len+1;
while(p!=-1&&ns[p].nxt[c]==q){
ns[p].nxt[c]=cpy;
p=ns[p].lnk;
}
ns[q].lnk=ns[cur].lnk=cpy;
}
}
lst=cur;
}
}
void build(){
rep(i,1,siz)g[ns[i].lnk].pub(i);
}
}sam;
map版
inline int get(char c){
if('0'<=c&&c<='9')return c-'0';
if('a'<=c&&c<='z')return 10+c-'a';
if('A'<=c&&c<='Z')return 10+26+c-'A';
}
struct suffix_automaton{
struct node{
map<int,int> nxt;
int len,lnk;
}ns[L<<1];
int siz,lst;
vec<int> g[L<<1];
void init(){
repl(i,0,L<<1)g[i].clear();
siz=lst=0;
ns[0].lnk=-1;
}
void insert(string s){
for(auto ch:s){
int cur=++siz,p=lst,c=get(ch);
ns[cur].len=ns[lst].len+1;
while(p!=-1&&!ns[p].nxt.count(c)){
ns[p].nxt[c]=cur;
p=ns[p].lnk;
}
if(p==-1)ns[cur].lnk=0;
else{
int q=ns[p].nxt[c];
if(ns[q].len==ns[p].len+1)ns[cur].lnk=q;
else{
int cpy=++siz;
ns[cpy]=ns[q];
ns[cpy].len=ns[p].len+1;
while(p!=-1&&ns[p].nxt[c]==q){
ns[p].nxt[c]=cpy;
p=ns[p].lnk;
}
ns[q].lnk=ns[cur].lnk=cpy;
}
}
lst=cur;
}
}
void build(){
rep(i,1,siz)g[ns[i].lnk].pub(i);
}
}sam;
应用
分配结束标记
我们知道,整个 SAM 的原始用途其实是维护每一个后缀,那么我们考虑如何给每一个后缀分配结束标记。
其实很简单,我们从最后的 \(lst\) 开始跳后缀链接,跳到的每一个状态都是整个字符串的后缀。
很显然,不多加讲解。
最长公共后缀
求字符串两个子串的最长公共后缀,不难想到,我们只需要找到这两个子串在 Parent 树上的 LCA,那么其等价类内最长串就是它们的最长公共后缀。
统计出现次数
实际上就是统计 \(\mathrm{endpos}\) 集合的大小。
我们考虑与 SAM 结合,不难发现,其每一次出现,都是作为原串的一个前缀的后缀出现。
我们可以借助后缀链接快速获取后缀,所以我们只需要知道当前节点是哪些前缀的后缀即可。
考虑前缀,我们发现每一条从初始状态开始的路径都代表了一个前缀。
所以,我们给每一个前缀都打上 \(1\) 的价值,最后在 Parent 树上做树型DP即可。
更详细地,我们在新建 \(cur\) 的时候,代表一个新的前缀,给其价值设为 \(1\)。
加入完毕整个字符串后,我们对 Parent 树做树型DP,每个节点的价值都加上其所有子树的价值和。
求后缀数组
也就是把所有后缀按字典序排序,并获取其起始位置。
事实上,我们只需要按字典序深搜 Parent 树,显然就可以得到答案。
所以现在唯一的问题就是如何按照字典序重构树。
为了方便,我们向 SAM 内传入反串,为什么呢?考虑下图,表示 Parent 树中一对父子节点:
显然,红框内的部分完全相同。对于这个父节点的所有子节点,其 \(\texttt{c}\) 必然不同。由于我们传入的是反串,那么我们直接按照 \(\texttt c\) 排序,就是按字典序排序了。
使用基数排序可以减少常数。
考虑 \(\texttt{c}\) 在反串上的位置,记 \(pos_u\) 表示 \(u\) 某一次出现的结束位置,显然:
考虑如何求 \(pos\),如果是新增的节点,那么其 \(\mathrm {endpos}\) 中必然包含 \(|s|\),也就是整个字符串的结尾位置,直接把 \(pos\) 设为当前总长或 \(len_{lst}+1\) 即可。
否则,那么这个节点是复制的,其 \(\mathrm {endpos}\) 中必然包含其复制来的 \(\mathrm{endpos}\),所以直接继承其复制来的节点即可。
最后深搜整棵 Parent 树,每找到一个后缀(预先打好标记)就将其在原串中的起始位置加入 \(sa\),具体的,这个位置为:
代码见例题 \(2\)。
例题
代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef __int128 i128;
typedef double db;
typedef long double ld;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
typedef pair<int,ll> pil;
typedef pair<ll,int> pli;
template <typename Type>
using vec=vector<Type>;
template <typename Type>
using grheap=priority_queue<Type>;
template <typename Type>
using lrheap=priority_queue<Type,vector<Type>,greater<Type> >;
#define fir first
#define sec second
#define pub push_back
#define pob pop_back
#define puf push_front
#define pof pop_front
#define chmax(a,b) a=max(a,b)
#define chmin(a,b) a=min(a,b)
#define rep(i,x,y) for(int i=(x);i<=(y);i++)
#define repl(i,x,y) for(int i=(x);i<(y);i++)
#define per(i,x,y) for(int i=(x);i>=(y);i--)
const int inf=0x3f3f3f3f;
const ll INF=0x3f3f3f3f3f3f3f3f;
const int mod=998244353;
const int L=1e6+5;
ll f[L<<1];
struct suffix_automaton{
struct node{
int nxt[26];
int len,lnk;
void clear(){
len=lnk=0;
memset(nxt,0,sizeof(nxt));
}
}ns[L<<1];
int siz,lst;
vec<int> g[L<<1];
void init(){
repl(i,0,L<<1)g[i].clear();
siz=lst=0;
ns[0].clear();
ns[0].lnk=-1;
}
void insert(string s){
for(auto ch:s){
int cur=++siz,p=lst,c=ch-'a';
f[cur]=1;
ns[cur].clear();
ns[cur].len=ns[lst].len+1;
while(p!=-1&&!ns[p].nxt[c]){
ns[p].nxt[c]=cur;
p=ns[p].lnk;
}
if(p==-1)ns[cur].lnk=0;
else{
int q=ns[p].nxt[c];
if(ns[q].len==ns[p].len+1)ns[cur].lnk=q;
else{
int cpy=++siz;
ns[cpy]=ns[q];
ns[cpy].len=ns[p].len+1;
while(p!=-1&&ns[p].nxt[c]==q){
ns[p].nxt[c]=cpy;
p=ns[p].lnk;
}
ns[q].lnk=ns[cur].lnk=cpy;
}
}
lst=cur;
}
}
void build(){
rep(i,1,siz)g[ns[i].lnk].pub(i);
}
}sam;
string s;
ll ans;
void dfs(int now){
for(auto nxt:sam.g[now]){
dfs(nxt);
f[now]+=f[nxt];
}
if(f[now]>1)chmax(ans,f[now]*sam.ns[now].len);
}
int main(){
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
cin>>s;
sam.init();
sam.insert(s);
sam.build();
dfs(0);
cout<<ans;
return 0;
}
代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef __int128 i128;
typedef double db;
typedef long double ld;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
typedef pair<int,ll> pil;
typedef pair<ll,int> pli;
template <typename Type>
using vec=vector<Type>;
template <typename Type>
using grheap=priority_queue<Type>;
template <typename Type>
using lrheap=priority_queue<Type,vector<Type>,greater<Type> >;
#define fir first
#define sec second
#define pub push_back
#define pob pop_back
#define puf push_front
#define pof pop_front
#define chmax(a,b) a=max(a,b)
#define chmin(a,b) a=min(a,b)
#define rep(i,x,y) for(int i=(x);i<=(y);i++)
#define repl(i,x,y) for(int i=(x);i<(y);i++)
#define per(i,x,y) for(int i=(x);i>=(y);i--)
const int inf=0x3f3f3f3f;
const ll INF=0x3f3f3f3f3f3f3f3f;
const int mod=998244353;
const int L=1e6+5;
inline int get(char c){
if('0'<=c&&c<='9')return c-'0';
if('A'<=c&&c<='Z')return 10+c-'A';
if('a'<=c&&c<='z')return 10+26+c-'a';
}
struct suffix_automaton{
struct node{
map<int,int> nxt;
int len,lnk;
int pos;
bool suf;
}ns[L<<1];
int siz,lst;
vec<int> g[L<<1];
void init(){
repl(i,0,L<<1)g[i].clear();
siz=lst=0;
ns[0].lnk=-1;
}
void insert(string s){
for(auto ch:s){
int cur=++siz,p=lst,c=get(ch);
ns[cur].pos=ns[cur].len=ns[lst].len+1;
ns[cur].suf=1;
while(p!=-1&&!ns[p].nxt.count(c)){
ns[p].nxt[c]=cur;
p=ns[p].lnk;
}
if(p==-1)ns[cur].lnk=0;
else{
int q=ns[p].nxt[c];
if(ns[q].len==ns[p].len+1)ns[cur].lnk=q;
else{
int cpy=++siz;
ns[cpy]=ns[q];ns[cpy].suf=0;
ns[cpy].len=ns[p].len+1;
while(p!=-1&&ns[p].nxt[c]==q){
ns[p].nxt[c]=cpy;
p=ns[p].lnk;
}
ns[q].lnk=ns[cur].lnk=cpy;
}
}
lst=cur;
}
}
int key[L<<1],cnt[64],srt[L<<1];
void build(string s){
rep(i,1,siz)key[i]=get(s[ns[i].pos-ns[ns[i].lnk].len-1]),cnt[key[i]]++;
rep(i,1,62)cnt[i]+=cnt[i-1];
rep(i,1,siz)srt[cnt[key[i]]--]=i;
rep(i,1,siz)g[ns[srt[i]].lnk].pub(srt[i]);
}
}sam;
string s;
int cnt;
int sa[L];
void dfs(int now){
if(sam.ns[now].suf)sa[++cnt]=s.size()-sam.ns[now].len+1;
for(auto nxt:sam.g[now])dfs(nxt);
}
int main(){
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
cin>>s;
reverse(s.begin(),s.end());
sam.init();
sam.insert(s);
sam.build(s);
dfs(0);
rep(i,1,cnt)cout<<sa[i]<<" ";
return 0;
}
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 地球OL攻略 —— 某应届生求职总结
· 周边上新:园子的第一款马克杯温暖上架
· Open-Sora 2.0 重磅开源!
· 提示词工程——AI应用必不可少的技术
· .NET周刊【3月第1期 2025-03-02】