SAM初探

SAM,即Suffix Automaton,后缀自动机。


关于字符串有很多玩法,有很多算法都是围绕字符串展开的。为什么?我的理解是:相较于数字组成的序列,字母组成的序列中每个单位上元素的个数是有限的。对于有限的东西,相较于无限的东西就会具有一些奇妙的性质。最简单的,就是序列扩展成的每个节点的儿子数是有限的。所以根据这个,从字符串Hash,到KMP,再到Suffix Array,Suffix Automaton,纷纷诞生。


后缀数组在处理字符串上相当于一把好钢,他能应付在字符串的大多数问题。那么为什么还要学SAM,很简单。SAM很有用,而且很短。优雅而有用,这就是要学的原因。



如果不说没用的话就是这玩意不容易写崩,而且可以直接构造后缀数组。


具体的细节不深究,搞OI的讲实用,所以这里面值谈构造应用,以题为主。




 


构造


众所周知,后缀数组的编写很糟心,常数也比较大。但是SAM好就好在他可以在线构造


先放代码:

struct SAM{
    int pre[MAXN],son[MAXN][26],cnt,len,now,step[MAXN],np,nq,p,q;
    SAM(){
        memset(pre,0,sizeof(pre));
        memset(son,0,sizeof(son));
        cnt=1;len=0;now=1;
    }
    void Extend(int nxt){
        p=now;np=++cnt;
        step[np]=step[now]+1;
        now=np;
        while(p&&!son[p][nxt]){
            son[p][nxt]=np;
            p=pre[p];
        }
        if(!p)pre[np]=1;
        else{
            q=son[p][nxt];
            if(step[q]==step[p]+1)pre[np]=q;
            else{
                step[(nq=++cnt)]=step[p]+1;
                memcpy(son[nq],son[q],sizeof(son[q]));
                pre[nq]=pre[q];
                pre[np]=pre[q]=nq;
                while(p&&son[p][nxt]==q){
                    son[p][nxt]=nq;
                    p=pre[p];
                }
            }
        }
    }
    int Walk(int nxt){
        while(!son[now][nxt]&&pre[now]){
            now=pre[now];
            Len=step[now];
        }
        if(!son[now][nxt])return 0;
        now=son[now][nxt];Len++;
        return Len;
    }
    void Build(){
        scanf("%s",s+1);
        int len=strlen(s+1);
        up(i,1,len)Extend(s[i]-'a');
    }
}sam;

关于构造我简单介绍一下。


首先,后缀自动机的结构是DAG。其共用了一些状态节点,同时和AC自动机类似,具有失配边。同时节点数为$O(2N)$。


后缀自动机的构造以在线添加的形式体现。就是sam.Extrend(nxt)。每调用一次最多添加$2$个节点。


每个节点有son和pre,pre指向的即为失配节点,son指向的即为儿子,同时有一个性质就是pre指向的点是可接受态的节点。


可接受态节点上代表的状态就是可以接受字母状态的节点。就是可以扩展后缀。


同时每个节点还有一个状态为len,就是当前节点到root的最长距离。


 


沿用构造里的变量,假设当前状态为now,我们要添加一个字符s[i],先把他转化成数字nxt然后调用Extend。


对于当前节点now,赋给p。同时新建一个节点np表示新增的一个状态同时把np的len置为len[p]+1。然后p不断沿着pre向上走,沿途把不存在nxt儿子的节点的nxt儿子赋为np,如果一直到了根节点, 说明这个字母之前字符串没有出现过。把pre[np]指向root即可。


如果走到了一个节点,他的son[p][nxt]已经被占用了,不慌,分两类情况讨论:


1.如果满足len[son[p][nxt]]==len[p]+1,很简单,直接把pre[np]指向son[p][nxt]即可。


2.如果不满足,就比较麻烦了。


先把son[p][nxt]赋给q,然后新建一个节点nq,使得len[nq]=len[q]+1,


然后把q的pre和son全部赋给nq。接着把q和np的pre指向nq。然后p沿着pre向上走,把所有满足son[p][nxt]==q的置为nq。


结束。


一些小细节建议参考代码。为什么这样构造我搞懂了会补上QAQ。


一些性质


1.$max(pre_s)=min(s)-1$


2.从$root$开始的任意一条路径都代表一个子串。


3.$right(trans(s,nxt)) \subset right(trans(pre_s,nxt))$


4.$|right(fa)|=\sum |right(son)|$


题目


SPOJ1811 LCS


很简单,直接暴力走就OK了。

//SPOJ1811
//by Cydiater
//2016.12.19
#include <iostream>
#include <queue>
#include <map>
#include <ctime>
#include <cstring>
#include <string>
#include <iomanip>
#include <cmath>
#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <bitset>
#include <set>
using namespace std;
#define ll long long
#define up(i,j,n)       for(int i=j;i<=n;i++)
#define down(i,j,n)     for(int i=j;i>=n;i--)
#define cmax(a,b)       a=max(a,b)
#define cmin(a,b)       a=min(a,b)
const int MAXN=5e5+5;
const int oo=0x3f3f3f3f;
char s[MAXN],t[MAXN];
int len,ans=0,Len=0;
struct SAM{
    int pre[MAXN],son[MAXN][26],cnt,len,now,step[MAXN],np,nq,p,q;
    SAM(){
        memset(pre,0,sizeof(pre));
        memset(son,0,sizeof(son));
        cnt=1;len=0;now=1;
    }
    void Extend(int nxt){
        p=now;np=++cnt;
        step[np]=step[now]+1;
        now=np;
        while(p&&!son[p][nxt]){
            son[p][nxt]=np;
            p=pre[p];
        }
        if(!p)pre[np]=1;
        else{
            q=son[p][nxt];
            if(step[q]==step[p]+1)pre[np]=q;
            else{
                step[(nq=++cnt)]=step[p]+1;
                memcpy(son[nq],son[q],sizeof(son[q]));
                pre[nq]=pre[q];
                pre[np]=pre[q]=nq;
                while(p&&son[p][nxt]==q){
                    son[p][nxt]=nq;
                    p=pre[p];
                }
            }
        }
    }
    int Walk(int nxt){
        while(!son[now][nxt]&&pre[now]){
            now=pre[now];
            Len=step[now];
        }
        if(!son[now][nxt])return 0;
        now=son[now][nxt];Len++;
        return Len;
    }
    void Build(){
        scanf("%s",s+1);
        int len=strlen(s+1);
        up(i,1,len)Extend(s[i]-'a');
    }
}sam;
namespace solution{
    void Prepare(){
        sam.Build();
        scanf("%s",t+1);
        len=strlen(t+1);
    }
    void Slove(){
        sam.now=1;
        up(i,1,len)cmax(ans,sam.Walk(t[i]-'a'));
        cout<<ans<<endl;
    }
}
int main(){
    //freopen("input.in","r",stdin);
    using namespace solution;
    Prepare();
    Solve();
    return 0;
}

 SPOJ1812 LCS II


相对来说麻烦一些。


还是随便挑了一个串建SAM,然后剩下的串在SAM上跑,对于每个串,标记在每个节点的最大值。对于所有串,取每个节点的最小值。跑完一次后,按照拓扑序更新父亲。

//SPOJ1812
//by Cydiater
//2018.12.23
#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <cmath>
#include <ctime>
#include <cstdlib>
#include <iomanip>
#include <algorithm>
#include <queue>
#include <map>
#include <bitset>
#include <set>
#include <vector>
using namespace std;
#define ll long long
#define up(i,j,n)   for(int i=j;i<=n;i++)
#define down(i,j,n) for(int i=j;i>=n;i--)
#define cmax(a,b)   a=max(a,b)
#define cmin(a,b)   a=min(a,b)
const int MAXN=2e5+5;
const int oo=0x3f3f3f3f;
inline int read(){
    char ch=getchar();int x=0,f=1;
    while(ch>'9'||ch<'0'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}
char s[MAXN],t[MAXN];
int flag[MAXN][15],ID=0,LEN[MAXN],ans=0,ml[MAXN],mat[MAXN],Len=0;
struct SAM{
    int son[MAXN][26],pre[MAXN],step[MAXN],p,q,np,nq,cnt,now,ma[MAXN],mi[MAXN];
    SAM(){
        cnt=now=1;
        memset(son,0,sizeof(son));
        memset(pre,0,sizeof(pre));
        memset(step,0,sizeof(step));
    }
    void Extend(int nxt){
        p=now;np=++cnt;
        step[np]=step[p]+1;
        now=np;
        while(p&&!son[p][nxt]){
            son[p][nxt]=np;
            p=pre[p];
        }
        if(!p)pre[np]=1;
        else{
            q=son[p][nxt];
            if(step[q]==step[p]+1)pre[np]=q;
            else{
                nq=++cnt;step[nq]=step[p]+1;
                memcpy(son[nq],son[q],sizeof(son[q]));
                pre[nq]=pre[q];
                pre[q]=pre[np]=nq;
                while(q&&son[p][nxt]==q){
                    son[p][nxt]=nq;
                    p=pre[p];
                }
            }
        }
    }
    void Build(){
        scanf("%s",s+1);
        int len=strlen(s+1);
        up(i,1,len)Extend(s[i]-'a');
        up(i,1,cnt)ml[step[i]]++;
        up(i,1,len)ml[i]+=ml[i-1];
        up(i,1,cnt)mat[ml[step[i]]--]=i;
        up(i,1,cnt)mi[i]=step[i];
    }
    void Walk(int nxt){
        while(!son[now][nxt]&&now){
            now=pre[now];
            Len=step[now];
        }
        if(!now){
            now=1;Len=0;
            return;
        }else{
            Len++;now=son[now][nxt];
            cmax(ma[now],Len);
        }
    }
}sam;
namespace solution{
    void Prepare(){
        sam.Build();
    }
    void Slove(){
        while(scanf("%s",t+1)!=EOF){
            int len=strlen(t+1);sam.now=1;Len=0;
            up(i,1,len)sam.Walk(t[i]-'a');
            down(i,sam.cnt,1){
                int t=mat[i];
                cmin(sam.mi[t],sam.ma[t]);
                if(sam.ma[t]&&sam.pre[t])sam.ma[sam.pre[t]]=sam.step[sam.pre[t]];
                sam.ma[t]=0;
            }
        }
        up(i,1,sam.cnt)cmax(ans,sam.mi[i]);
        cout<<ans<<endl;
    }
}
int main(){
    //freopen("input.in","r",stdin);
    using namespace solution;
    Prepare();
    Slove();
    return 0;
}

 SPOJ8222


后缀自动机主要还是要靠DP。


分为树边DP和转移边DP,现在这两个具体怎么应用还不是很清楚,多写写题吧。


这道题就是树边DP,因为这个相当于求$right$集大小。而父亲节点的$right$集等于儿子的并集(而且儿子的$right$集不存在交集)。所以topsort后DP即可。

//SPOJ 8222
//by Cydiater
//2017.1.19
#include <iostream>
#include <queue>
#include <map>
#include <ctime>
#include <cmath>
#include <cstring>
#include <string>
#include <algorithm>
#include <cstdlib>
#include <cstdio>
#include <iomanip>
#include <bitset>
#include <set>
#include <vector>
using namespace std;
#define ll long long
#define up(i,j,n)   for(int i=j;i<=n;i++)
#define down(i,j,n) for(int i=j;i>=n;i--)
#define cmax(a,b)   a=max(a,b)
#define cmin(a,b)   a=min(a,b)
const int MAXN=1e6+5;
const int oo=0x3f3f3f3f;
char s[MAXN];
int N,cnt[MAXN],rnk[MAXN],f[MAXN],g[MAXN];
struct SAM{
    int son[MAXN][26],now,step[MAXN],cnt,pre[MAXN];
    SAM(){cnt=now=1;}
    void Extend(int nxt){
        int p=now,np=++cnt;now=np;
        step[np]=step[p]+1;g[np]++;
        for(;(!son[p][nxt])&&p;p=pre[p])son[p][nxt]=np;
        if(!p)pre[np]=1;
        else{
            int q=son[p][nxt];
            if(step[q]==step[p]+1)pre[np]=q;
            else{
                int nq=++cnt;
                step[nq]=step[p]+1;
                memcpy(son[nq],son[q],sizeof(son[q]));
                pre[nq]=pre[q];
                pre[q]=pre[np]=nq;
                for(;son[p][nxt]==q;p=pre[p])son[p][nxt]=nq;
            }
        }
    }
}sam;
namespace solution{
    void Prepare(){
        scanf("%s",s+1);
        N=strlen(s+1);
        up(i,1,N)sam.Extend(s[i]-'a');
    }
    void Solve(){
        up(i,2,sam.cnt)cnt[sam.step[i]]++;
        up(i,1,N)cnt[i]+=cnt[i-1];
        up(i,2,sam.cnt)rnk[cnt[sam.step[i]]--]=i;
        down(i,sam.cnt,2)g[sam.pre[rnk[i]]]+=g[rnk[i]];
        up(i,2,sam.cnt)cmax(f[sam.step[i]],g[i]);
        up(i,1,N)printf("%d\n",f[i]);
    }
}
int main(){
    //freopen("input.in","r",stdin);
    using namespace solution;
    Prepare();
    Solve();
    return 0;
}

 SPOJ7258


对SAM理解还是不够透彻。


根据转移边DP。设$f[i]$表示从状态$i$发出可以到达的子串数量。可以得到


$f[i]=\sum f[child] +1$


然后对于每个询问走一遍即可。

//SPOJ 7258
//by Cydiater
//2017.1.20
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <cstring>
#include <string>
#include <iomanip>
#include <algorithm>
#include <queue>
#include <map>
#include <bitset>
#include <set>
#include <vector>
using namespace std;
#define ll long long
#define up(i,j,n)   for(int i=j;i<=n;i++)
#define down(i,j,n) for(int i=j;i>=n;i--)
#define cmax(a,b)   a=max(a,b)
#define cmin(a,b)   a=min(a,b)
const int MAXN=2e5+5;
const int oo=0x3f3f3f3f;
inline int read(){
    char ch=getchar();int x=0,f=1;
    while(ch>'9'||ch<'0'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}
int rnk[MAXN],N,M,label[MAXN],f[MAXN];
char s[MAXN];
struct SAM{
    int step[MAXN],son[MAXN][26],pre[MAXN],now,cnt;
    SAM(){now=cnt=1;}
    void Extend(int nxt){
        int p=now,np=++cnt;now=np;
        step[np]=step[p]+1;
        for(;p&&!son[p][nxt];p=pre[p])son[p][nxt]=np;
        if(!p)pre[np]=1;
        else{
            int q=son[p][nxt],nq;
            if(step[q]==step[p]+1)pre[np]=q;
            else{
                step[(nq=++cnt)]=step[p]+1;
                memcpy(son[nq],son[q],sizeof(son[q]));
                pre[nq]=pre[q];
                pre[np]=pre[q]=nq;
                for(;son[p][nxt]==q;p=pre[p])son[p][nxt]=nq;
            }
        }
    }
    void Topsort(){
        up(i,1,cnt)label[step[i]]++;
        up(i,1,N)label[i]+=label[i-1];
        up(i,1,cnt)rnk[label[step[i]]--]=i;
    }
}sam;
namespace solution{
    void Prepare(){
        scanf("%s",s+1);
        N=strlen(s+1);
        up(i,1,N)sam.Extend(s[i]-'a');
    }
    void Solve(){
        sam.Topsort();
        down(i,sam.cnt,1){
            int node=rnk[i];
            f[node]=1;
            up(j,0,25)f[node]+=f[sam.son[node][j]];
        }
        M=read();
        while(M--){
            int len=read(),now=1;
            while(len){
                up(i,0,25)if(sam.son[now][i]){
                    if(f[sam.son[now][i]]>=len){
                        len--;
                        now=sam.son[now][i];
                        putchar('a'+i);
                        break;
                    }else len-=f[sam.son[now][i]];
                }
            }
            printf("\n");
        }
    }
}
int main(){
    //freopen("input.in","r",stdin);
    using namespace solution;
    Prepare();
    Solve();
    return 0;
}

 POJ1743


充分利用SAM的性质。


我们可以知道,每个节点的${right}$是他的出现位置,而$step_i$代表的是每个节点的长度,如果可以重复的话,我们对每个$|right| \geq 2$的集合统计$step$,取max即可。


但是这个要求不可重复,所以只需要在取max前统计出每个点$right$的最大值和最小值,然后做差和$step$取min。这样才能满足不能重叠和最大两个要求。

//POJ1743
//by Cydiater
//2017.1.21
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include <cstring>
#include <string>
#include <cmath>
#include <ctime>
#include <algorithm>
#include <map>
#include <queue>
#include <bitset>
#include <set>
#include <vector>
using namespace std;
#define ll long long
#define up(i,j,n)   for(int i=j;i<=n;i++)
#define down(i,j,n) for(int i=j;i>=n;i--)
#define cmax(a,b)   a=max(a,b)
#define cmin(a,b)   a=min(a,b)
const int MAXN=40005;
const int oo=0x3f3f3f3f;
inline int read(){
    char ch=getchar();int x=0,f=1;
    while(ch>'9'||ch<'0'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}
int N,arr[MAXN>>1],L[MAXN],R[MAXN],rnk[MAXN],label[MAXN>>1],ans=0;
struct SAM{
    int pre[MAXN],son[MAXN][180],now,cnt,step[MAXN];
    SAM(){now=cnt=1;}
    void Clear(){
        memset(pre,0,sizeof(pre));
        memset(son,0,sizeof(son));
        memset(step,0,sizeof(step));
        now=cnt=1;
    }
    void Extend(int nxt){
        int p=now,np=++cnt;now=np;
        step[np]=step[p]+1;L[np]=R[np]=step[np];
        for(;p&&!son[p][nxt];p=pre[p])son[p][nxt]=np;
        if(!p)pre[np]=1;
        else{
            int q=son[p][nxt],nq;
            if(step[q]==step[p]+1)pre[np]=q;
            else{
                step[(nq=++cnt)]=step[p]+1;
                memcpy(son[nq],son[q],sizeof(son[q]));
                pre[nq]=pre[q];
                pre[np]=pre[q]=nq;
                for(;son[p][nxt]==q;p=pre[p])son[p][nxt]=nq;
            }
        }
    }
    void Topsort(){
        up(i,1,cnt)label[step[i]]++;
        up(i,1,N-1)label[i]+=label[i-1];
        up(i,1,cnt)rnk[label[step[i]]--]=i;
    }
}sam;
namespace solution{
    void Prepare(){
        N=read();if(N==0)exit(0);
        memset(label,0,sizeof(label));
        memset(L,10,sizeof(L));
        memset(R,0,sizeof(R));
        sam.Clear();ans=0;
        up(i,1,N)arr[i]=read();
        up(i,1,N-1)arr[i]=arr[i+1]-arr[i]+88;
        up(i,1,N-1)sam.Extend(arr[i]);
        sam.Topsort();
    }
    void Solve(){
        down(i,sam.cnt,1){
            int node=rnk[i];
            cmin(L[sam.pre[node]],L[node]);
            cmax(R[sam.pre[node]],R[node]);
        }
        up(i,1,sam.cnt){
            int node=rnk[i];
            cmax(ans,min(R[node]-L[node],sam.step[node]));
        }
        printf("%d\n",ans<4?0:ans+1);
    }
}
int main(){
    //freopen("input.in","r",stdin);
    using namespace solution;
    while(true){
        Prepare();
        Solve();
    }
    return 0;
}

 小结

后缀自动机学了大概一个半月(主要是因为中间停了一段时间取搞文化课),学的越深越能感到后缀自动机的牛逼。


还是那句话:出现次数向父亲传递,接收串数从儿子获取。具体的构造什么的推荐去看CLJ的ppt,那个讲的最全面最透彻,可能不是很容易理解但是建议一定要多看几遍肛出来。


如果具体的概念没有完全理解后面的题基本上无从下手。


后缀自动机很重要的两个东西就是${right}$和$Step$,而这两个东西都可以在结合基数排序的基础上在$O(N)$的时间复杂度内求出。


其中$right$代表了每个串的出现位置,$Step$代表最长的长度。根据这个可以做一些简单的DP之类的。


另外有一个性质就是$Step_{fa}+1=Min$,而且不同的长度不存在交集,而且长度是连续的,这样子也可以方便我们搞令一些事情。


另外,后缀自动机的Parent树是逆序的后缀树,所以我们也就可以根据这个求出后缀树,在后缀树上做一些数据结构的事之类的。


总之,应用非常丰富。