Magolor 首页文章幻灯片归档标签
求后缀数组的SAM与SA-IS方法
|OI
字数总计: 7.6k|阅读时长: 34 分钟

填一个OI期间一直想学但是没有学的坑: SA-IS求后缀数组。

众所周知后缀数组有$O(n \log n)$的倍增方法以及$O(n)$的DC3方法。而本文将会简单介绍倍增方法并详细讲解黑科技: $O(n)$的SAM方法以及SA-IS方法。

食用本文前建议对后缀数组(SA)和后缀自动机(SAM)有基本了解。

这里我们假设字符集大小$\vert\Sigma\vert$是一个可以接受的比较小的常数(比如字符集是所有小写字母，$\vert\Sigma\vert = 26$)。


倍增方法

基本思路

倍增方法求后缀数组的基本思路是: 如果对每个长度为$j$的子串已经求出一个排名数组，即可以$O(1)$判断两个长度为$2^j$的子串大小，那么就可以通过一个二元组的基数排序得到长度为$2^{j+1}$的子串的排名数组。


基数排序

基数排序:

1
2
3
#define Radix(rk,ps)
    for(memset(book+1,0,m<<2), i = 1; i <= n; ++book[rk], i++);\
    for(i = 1; i <= m; book[i] += book[i-1], i++); for(i = n; i; SA[book[Rank[ps]]--] = ps, i--);

传入的rk表示长度为$j$的子串的排名数组，ps只是为了方便处理下标变换。


倍增

1
2
3
4
5
for(p = L = 1; p < n; m = p, L <<= 1)
{
    for(p = 0, i = n-L+1; i <= n; x[++p] = i++); for(i = 1; i <= n; SA[i]>L?x[++p]=SA[i]-L:0, i++); Radix(Rank[x[i]],x[i]);
    for(memcpy(x+1,Rank+1,n<<2), Rank[SA[1]] = p = 1, i = 2; i <= n; Rank[SA[i]] = p+=x[SA[i]]^x[SA[~-i]]||x[SA[i]+L]^x[SA[~-i]+L], i++);
}

循环倍增步长，每次需要按照(Rank[i],Rank[i+L])顺序排序，基数排序是稳定排序，先按照第二关键字排序再按照第一关键字排序即可。就按照SA[]的顺序预处理x[]是已经按照第二关键字顺序排好序的第一关键字数组，所以再按照第一关键字基数排序即可。最后一行处理一下新的排名(类似于离散化)。


Doubling SA

完整代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
#include <bits/stdc++.h>
using namespace std;
#define MAXN 100000
#define rint register int
inline int rf(){int r;int s=0,c;for(;!isdigit(c=getchar());s=c);for(r=c^48;isdigit(c=getchar());(r*=10)+=c^48);return s^45?r:-r;}
int SA[MAXN+5], Rank[MAXN+5], Height[MAXN+5], n; char s[MAXN+5];
#define Radix(rk,ps)
    for(memset(book+1,0,m<<2), i = 1; i <= n; ++book[rk], i++);\
    for(i = 1; i <= m; book[i] += book[i-1], i++); for(i = n; i; SA[book[Rank[ps]]--] = ps, i--);
void Doubling(rint i=0)
{
    static int book[MAXN+5], x[MAXN+5], m, p, L; m = 27; Radix(Rank[i]=s[i]-'a'+1,i);
    for(p = L = 1; p < n; m = p, L <<= 1)
    {
        for(p = 0, i = n-L+1; i <= n; x[++p] = i++); for(i = 1; i <= n; SA[i]>L?x[++p]=SA[i]-L:0, i++); Radix(Rank[x[i]],x[i]);
        for(memcpy(x+1,Rank+1,n<<2), Rank[SA[1]] = p = 1, i = 2; i <= n; Rank[SA[i]] = p+=x[SA[i]]^x[SA[~-i]]||x[SA[i]+L]^x[SA[~-i]+L], i++);
    }
}
void CalHeight(){for(rint i = 1, L = 0, j; i <= n; Height[Rank[i++]] = L) for(j = SA[~-Rank[i]], L-=L>0; i+L<=n&&j+L<=n&&s[i+L]==s[j+L]; ++L);}
int main()
{
    scanf("%s",s+1), n = strlen(s+1), Doubling(), CalHeight();
    for(rint i = 1; i <= n; printf("%d%c",SA[i],"\n "[i<n]), i++);
    for(rint i = 2; i <= n; printf("%d%c",Height[i],"\n "[i<n]), i++); return 0;
}


小Trick

往往通过使用代码技巧可以降低代码复杂度。

例如上述代码中使用过的宏函数:

1
2
3
#define Radix(rk,ps)
    for(memset(book+1,0,m<<2), i = 1; i <= n; ++book[rk], i++);\
    for(i = 1; i <= m; book[i] += book[i-1], i++); for(i = n; i; SA[book[Rank[ps]]--] = ps, i--);

宏函数会简单地将所有代码段中的rk替换为传入的参数，更夸张的是，宏函数甚至可以:

1
#define Op(a,b) a b

然后使用Op(a,++);，Op(a,--);甚至Op(int x,[3]);等写法。在后面的SA-IS方法代码中会使用类似于这种神奇的宏函数。

另外，求后缀数组过程中往往用到很多数组之间来回折腾的操作，这时候使用C++自带函数会好很多: memcpy(或copy)，memset(或fill)来代替复制和清零；partial_sum(a.begin(),a.end(),b.begin())求出a[]的前缀和b[]；transform(a.begin(),a.end(),b.begin(),function)将a[]中每个元素调用function的返回值存入b[](配合C++11 lambda函数使用更方便)。

例如上述代码中的for(i = 1; i <= m; book[i] += book[i-1], i++);就可以替换为partial_sum(book,book+m+1,book);。后面的SA-IS算法中可以看到更多的使用范例。


SAM方法

基本思路

SAM方法求后缀数组的基本思路是: 后缀树的DFS序就是后缀数组，非常方便，我们之所以不这样做是因为构建后缀树的Ukkonen算法很糟糕我太菜了(被大佬打脸QwQ)。然而根据后缀自动机的性质，有后缀自动机的parent树(我称为link树)就是反串后缀的树。

UPD2019.08.18: 已更新Ukkonen算法，见Ukkonen算法与后缀树构造后缀数组。


SAM

请先行自学SAM基础知识。

这里只贴出一个完整后缀自动机节点功能的代码:

1
2
3
4
5
6
7
8
9
struct Node{int c[26],maxlen,link,mx,ch,end; ll num0,num1; bool clone,accept;}t[MAXN+5];
int tot, last, book[(MAXN>>1)+5], o[MAXN+5], in[MAXN+5], q[MAXN+5], head, tail, n; char s[(MAXN>>1)+5];
#define New(p,w,x,i) (++tot, maxlen(tot) = p?maxlen(p)+1:0, mx(tot) = w, ch(tot) = x, end(tot) = i, tot)
inline void Append(int x, int i)
{
    rint p = last, np = last = New(p,1,x,i), q, nq; for(; p && !c(p,x); c(p,x) = np, p = link(p));
    if(!p) link(np) = 1; else if(maxlen(q=c(p,x))==maxlen(p)+1) link(np) = q; else
    {for(nq = New(p,0,x,end(q)), memcpy(t[nq].c,t[q].c,104); p && c(p,x)==q; c(p,x) = nq, p = link(p)); clone(nq) = true, link(nq) = link(q), link(q) = link(np) = nq;}
}

注意常用的maxlen，link(parent)，end和mx。

注意字符集很大的时候需要用std::map。这一部分会成为时间复杂度瓶颈，使用手写哈希表或者unordered_map可以做到$O(1)$。但即使是这样，在处理SAMSA的时候需要排序，也不能逃离一个至少$O(n \log n)$(对于过大字符集离散化或暴力排序)或$O(n + \vert\Sigma\vert)$(对于中等大小的字符集基数排序)的复杂度。事实上，任何后缀数组在处理任意大字符集情况下复杂度都至少是$O(\min(n+\vert\Sigma\vert,n \log n))$，这个复杂度以排序问题作为下界。

这些信息的具体维护方式如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
last = New(0,0,0,0), scanf("%s",s+1), n = strlen(s+1); for(rint i = 1; i <= n; Append(s[i]-'a',i), i++);

for(rint p = last; p; accept(p) = true, p = link(p));

for(rint i = 1; i <= tot; ++book[maxlen(i)], i++); for(rint i = 1; i <= n; book[i] += book[i-1], i++);
for(rint i = tot; i; o[book[maxlen(i)]--] = i, i--); for(rint i = tot; i; mx(link(o[i])) += mx(o[i]),i--); mx(1) = 0;

for(rint p = 1, x; p <= tot; p++) for(x = 0; x < 26; ++in[c(p,x)], x++);
for(rint i = 1; i <= tot; !in[i] ? q[++tail] = i : 0, i++);
for(rint p, x; head<=tail; )
    for(p = q[++head], x = 0; x < 26; x++)
        !--in[c(p,x)] ? q[++tail] = c(p,x) : 0;
for(rint i = tot, p, x; i; i--)
    for(p = q[i], num0(p) = p>1, num1(p) = mx(p), x = 0; x < 26; x++)
        num0(p) += num0(c(p,x)), num1(p) += num1(c(p,x));

完整代码以及更多SAM、SAMSA相关可以看我的博客OI资源分享中的Algorithm Template。这里有快速访问的文档SAM。

请先行自学SAM基础知识。


SAM构建后缀树

对于构建后缀树的Ukkonen算法，由于我不会，这里只先给出@xehoth大佬的代码: BZOJ4516-SuffixTree.cpp。感谢@xehoth大佬。

后缀自动机link树与反串后缀树同构。问题就是如何将自动机的link树转化为后缀树，因为link边的字典序关系是不确定的。注意到link边就已经确定了$LCP$，例如所有连向2的子树(1-2-3串和1-2-7串)都有$LCP =$ a。又由于后缀树是由Trie压缩而来，同一LCP位置(深度相同的位置)，转移边字符一定不同，所以后缀树每条边开头的第一个字符作为唯一关键字就可以完成后缀边排序了。

所以问题的关键就是确定每条link边代表的字符串的首字母。上图就是原串dbabbaa逆序建后缀自动机: 对于每个节点记录这个状态最后一次出现的时间(注意对于复制状态也要更新时间)，因为逆序建了后缀自动机，所以这实际上找到了当前状态第一次出现的开头位置(例如ab的反串ba的最后一位a就是原串ab的最前一位)，所以记为beg(i)。要找到其状态的$LCP$，就是其前驱link状态的最长长度，自然就是maxlen(i)。因此求出$LCP$，这个节点应该分叉的位置就是beg(i)+maxlen(link(i))。

也就是说，每条后缀树的边(i到i的父亲link(i))的开头字符就是s[beg(i)+maxlen(link(i))]。

所以最后用后缀自动机构建后缀数组的方法就有了:从$n$到$1$建后缀自动机并维护beg(i)，将每个节点按照s[beg(i)+maxlen(link(i))]桶排，从小到大加边(如果是vector从小到大；如果是邻接表，边后进先出，从大到小)，然后DFS就可以了。DFS的时候只考虑非复制节点，将其beg赋为当前的SA就可以了。


UPD2019.08.18 Ukkonen算法构建后缀树

学习了一下Ukkonen算法，于是在博客上补充了新的内容: Ukkonen算法与后缀树构造后缀数组。


SAMSA

完整代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#include <bits/stdc++.h>
using namespace std;
#define MAXN 200000
#define c(x,y) t[x].c[y]
#define maxlen(x) t[x].maxlen
#define link(x) t[x].link
#define mx(x) t[x].mx
#define end(x) t[x].end
#define rint register int
inline int rf(){int r;int s=0,c;for(;!isdigit(c=getchar());s=c);for(r=c^48;isdigit(c=getchar());(r*=10)+=c^48);return s^45?r:-r;}
struct Node{int c[26],maxlen,link,end; bool mx;}t[MAXN+5]; int first[MAXN+5], start, tot, last, n, rk;
#define New(p,w,i) (++tot, first[tot] = -1, maxlen(tot) = p?maxlen(p)+1:0, mx(tot) = w, end(tot) = i, tot)
struct Edge{int to,nex; Edge(){} Edge(int _to, int _nex):to(_to),nex(_nex){}}e[MAXN+5]; int etot;
inline void Add(int a, int b){e[etot] = Edge(b,first[a]), first[a] = etot++;}
char s[(MAXN>>1)+5]; int book[26], w[MAXN+5], o[MAXN+5], SA[(MAXN>>1)+5], Rank[(MAXN>>1)+5], Height[(MAXN>>1)+5];
inline void Append(int x, int i)
{
    rint p = last, np = last = New(p,1,i), q, nq; for(; p && !c(p,x); c(p,x) = np, p = link(p));
    if(!p) link(np) = start; else if(maxlen(q=c(p,x))==maxlen(p)+1) link(np) = q;
    else
    {
        for(nq = New(p,0,i), memcpy(t[nq].c,t[q].c,104); p && c(p,x)==q; c(p,x) = nq, p = link(p));
        link(nq) = link(q), link(q) = link(np) = nq;
    }
}
void DFS(int p){mx(p) ? SA[Rank[end(p)]=++rk] = end(p) : 0; for(rint u = first[p], v; ~u; DFS(e[u].to), u = e[u].nex);}
void CalHeight(){for(rint i = 1, L = 0, j; i <= n; Height[Rank[i++]] = L) for(j = SA[~-Rank[i]], L-=L>0; i+L<=n&&j+L<=n&&s[i+L]==s[j+L]; ++L);}
int main()
{
    start = 1, New(0,0,0), scanf("%s",s+1), n = strlen(s+1); last = 1; for(rint i = n; i; Append(s[i]-'a',i), i--); book[0] = 1;
    for(rint i = 2; i <= tot; ++book[w[i]=s[end(i)+maxlen(link(i))]-'a'], i++); for(rint i = 1; i < 26; book[i] += book[i-1], i++);
    for(rint i = tot; i; o[book[w[i]]--] = i, i--); for(rint i = tot; i>1; Add(link(o[i]),o[i]), i--); DFS(1); CalHeight();
    for(rint i = 1; i <= n; printf("%d%c",SA[i],"\n "[i<n]), i++); for(rint i = 2; i <= n; printf("%d%c",Height[i],"\n "[i<n]), i++); return 0;
}


SA-IS方法

(敲黑板)重点重点！

概念

首先定义一些概念。

如果一个后缀s[i,n]比s[i+1,n]字典序要小，那么记为$S$型后缀(Smaller)，否则记为$L$型后缀(Larger)，显然后缀不会相等。

当两个后缀第一个字符相同时只需要递归比较既可以判断类型，于是很显然地有:

1
t[n] = 0; for(rint i = n-1; i; t[i] = s[i]^s[i+1]?s[i]>s[i+1]:t[i+1], i--);

这里默认在最后一个位置加入一个字典序最小的特殊字符，且规定最后一位是$S$型后缀。上述代码中记$S$型后缀为$0$，$L$型后缀为$1$。

每个左侧为$L$型后缀的$S$型后缀是一个$LMS$型后缀(Left-Most-S)。例如:

1
2
3
4
1 2 3 4 5 6 7 8 9 A B C D
a a b a a b b a c b a b #
S S L S S L L S L L S L S
      *       *     *   *

所有星号标出的位置就是$LMS$型后缀。而$LMS$型后缀切割出的每一段都是一个$LMS$子串: aab, aabb, acb, ab 和 #。

然后我们为方便起见认为字符串是横着的，后缀数组是竖着的，即用一个后缀的”左”和”右”来称呼其最前面加上或者减去一个字符的后缀，用”上”和”下”来称呼比其字典序小和大的后缀。


基本思路

SA-IS方法求后缀数组的基本思路是:

如果得到了$LMS$型后缀的排名，又因为通过$S$和$L$的定义天然可以知道相邻两个后缀的大小关系，那么获得所有后缀的排名就是一个多路归并排序: 每个$LMS$型左侧的若干个$L$根据定义字典序递增，多个$LMS$型左侧的递增序列归并；右侧的$S$递减序列同理归并。这种方法称为诱导排序。

而获得$LMS$型后缀的排名的方法是将$LMS$子串”离散化”，例如之前的串aabaabbacbab#划分为aab, aabb, acb, ab 和 #，就根据相互的大小关系分别可以离散化为b, c, e, d 和 a。那么求$LMS$型后缀排名就是求bceda的后缀数组，就变成了一个递归的子问题。

而离散化过程需要求出$LMS$子串的排名，做法也是一个诱导排序。

这样实现一个SA-IS算法: 其中s[]表示原串，t[]表示每个位置后缀的类型，tmp[]存储所有$LMS$型后缀的下标，n表示原串长度，m表示原串字符集数量:

1
void SAIS(int s[], bool t[], int tmp[], int n, int m)

首先要预处理出后缀类型并标注出$LMS$后缀:

1
2
t[n] = 0; for(rint i = n-1; i; t[i] = s[i]^s[i+1]?s[i]>s[i+1]:t[i+1], i--);
for(rint i = 2; i <= n; Rank[i] = t[i-1]&&!t[i]?tmp[++n1]=i,n1:0, i++);

tmp[i]表示从左向右第i个$LMS$型后缀的下标，Rank[i]表示位置i是从左向右第几个$LMS$型后缀(如果位置i根本不是$LMS$型后缀那么就是$0$)。

接着我们调用诱导排序获得$LMS$子串的排名，得到一个初步排序的排名数组。注意这个类似于后缀数组，但是其中不保证后缀有序，而是只对于所有$LMS$型后缀位置满足其排名一定符合$LMS$子串的相对大小关系。假设已经诱导排序得到了这个数组，接下来是离散化部分: 

1
2
3
4
5
6
for(rint i = 1, x, y = 0; i <= n; i++) if(x=Rank[SA[i]]){
    if(m1 <= 1 || tmp[x+1]-tmp[x]!=tmp[y+1]-tmp[y]) ++m1;
    else for(rint a = tmp[x], b = tmp[y]; a <= tmp[x+1]; a++, b++)
            if((s[a]<<1|t[a])^(s[b]<<1|t[b])){++m1; break;}
    s1[y=x] = m1;
}

假设诱导排序得到的SA[]中存储了一个排名，其中名次越大的$LMS$型后缀下标对应的$LMS$子串字典序越大(而不是对应的后缀字典序越大)。那么随i递增枚举Rank[SA[i]]，如果不为$0$，就是按字典序增大枚举了$LMS$子串，如果当前枚举到x=Rank[SA[i]]，就对应$LMS$子串s[tmp[x],tmp[x+1])。

那么这段代码就是: 首先判断当前$LMS$子串和之前的$LMS$子串长度不等一定是新的$LMS$子串，要离散化为一个新字符；否则按位比较，如果有一位字符不一样或后缀类型不一样都离散化为一个新字符。

离散化完得到了一个新的字符串s1[]，那么递归即可:

1
2
if(m1 < n1) SAIS(s1,t+n,tmp+n1,n1,m1); else for(rint i = 1; i <= n1; SA[s1[i]] = i, i++);
for(rint i = 1; i <= n1; s1[i] = tmp[SA[i]], i++);

这里如果m1 < n1不成立，即m1 == n1，就说明新字符串每个字符都不一样，那么求其后缀数组就是一个很trivial的问题，基数排序即可，否则递归。

已经求出了s1的后缀数组SA[]，那么就可以构造出$LMS$型后缀的排名关系，这里临时使用s1[]作为新的tmp[]来按照一个新的顺序存储各个$LMS$后缀的顺序。

那么现在已经得到了$LMS$后缀的顺序表，通过诱导排序求出整个后缀顺序表就做完了。将诱导排序记为IS()，代码就是:

1
2
3
4
5
6
7
8
9
10
11
12
13
void SAIS(int s[], bool t[], int tmp[], int n, int m){
    int n1 = 0, m1 = Rank[1] = 0, *s1 = s+n; t[n] = 0;
    for(rint i = n-1; i; t[i] = s[i]^s[i+1]?s[i]>s[i+1]:t[i+1], i--);
    for(rint i = 2; i <= n; Rank[i] = t[i-1]&&!t[i]?tmp[++n1]=i,n1:0, i++); IS(tmp);
    for(rint i = 1, x, y = 0; i <= n; i++) if(x=Rank[SA[i]]){
        if(m1 <= 1 || tmp[x+1]-tmp[x]!=tmp[y+1]-tmp[y]) ++m1;
        else for(rint a = tmp[x], b = tmp[y]; a <= tmp[x+1]; a++, b++)
                if((s[a]<<1|t[a])^(s[b]<<1|t[b])){++m1; break;}
        s1[y=x] = m1;
    }
    if(m1 < n1) SAIS(s1,t+n,tmp+n1,n1,m1); else for(rint i = 1; i <= n1; SA[s1[i]] = i, i++);
    for(rint i = 1; i <= n1; s1[i] = tmp[SA[i]], i++); IS(s1);
}

这里s[]和t[]都被复用了，这样可以避免每次递归动态开空间造成效率损失。


IS

所以遗留下来的问题就是两次诱导排序: 第一次是已知按照乱序(初始顺序)存放的所有$LMS$型后缀位置(存储在tmp[]里)，要求所有$LMS$子串的排名(一个不完全的子串SA[])；第二次是已知按照$LMS$型后缀字典序递增存放的所有$LMS$型后缀(存储在s1[]里)，要求所有后缀的排名(SA[])。

先考虑第二次诱导排序是怎么做的，首先已知$LMS$型后缀的相对大小关系:

1
2
memset(SA+1,0,n<<2); memset(c+1,0,m<<2); for(rint i = 1; i <= n; ++c[s[i++]]);\
partial_sum(c+1,c+m+1,c+1); memcpy(p+1,c+1,m<<2); for(rint i = n1; i; Ar(s1[i],--), i--);\

这就是先用基数排序按照首字母分组(后面会用到)，然后插入所有的$LMS$型后缀到同一组的末尾。

接下来要对所有的$L$型后缀确定相对大小关系。现在每个$LMS$型后缀左侧的连续$L$型序列已经形成了一个字典序递增序列，接下来要归并这些序列，例如:

1
2
3
4
1 2 3 4 5 6 7 8 9 A B C D
a a b a a b b a c b a b #
S S L S S L L S L L S L S
      *       *     *   *

即需要归并suf[4],suf[3]，suf[8],suf[7],suf[6]，suf[B],suf[A],suf[9]，suf[D],suf[C]得到一个递减序列。做法是按着SA[]顺序扫一遍:

1
transform(c,c+m,p+1,[](int a){return++a;}); for(rint i = 1; i <= n; SA[i]>1&&t[SA[i]-1]?Ar(SA[i]-1,++):0, i++);\

因为本身SA[]维护的就是顺序，所以扫描就是不断找到这若干个序列的最小值，将其左侧的后缀加入到其下方某个位置。这里的下方某个位置取法是将后缀按照首字母分组，组内顺序从上向下存放。为什么这样是可以的？

也就是需要理解，为什么要比较suf[SA[i]-1]和suf[SA[j]-1]，只需要比较suf[SA[i]]和suf[SA[j]]。因为之前提到过，类型的比较是递归的。如果suf[SA[i]-1]和suf[SA[j]-1]的首字母不一样，那么根据首字母分组的做法他们一定被分到不同组，那么相对顺序一定是正确的；如果首字母一样，那么类型比较递归就是比较suf[SA[i]]和suf[SA[j]]，所以我们只要保证按照i递增即SA[i]字典序递增的顺序枚举，最后得到的相对关系就是合法的。

例如字符串mmiissiissiippii和它的后缀数组:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
m  m  i  i  s  s  i  i  s  s  i  i  p  p  i  i  #
L  L  S  S  L  L  S  S  L  L  S  S  L  L  L  L  S
      *           *           *                 *

#                   S   *
i#                  L   
ii#                 L   
iippii#             S   *
iissiippii#         S   *
iissiissiippii#     S   *
ippii#              S   
issiippii#          S   
issiissiippii#      S   
miissiissiippii#    L   
mmiissiissiippii#   L   
pii#                L   
ppii#               L   
siippii#            L   
siissiippii#        L   
ssiippii#           L   
ssiissiippii#       L

在诱导排序的时候首先奇数排序按照首字母分组，将所有$LMS$型后缀插入到同一组的末尾得到右侧数组:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
m  m  i  i  s  s  i  i  s  s  i  i  p  p  i  i  #
L  L  S  S  L  L  S  S  L  L  S  S  L  L  L  L  S
      *           *           *                 *

#                   S   *   |   #                   S   *
i#                  L       |   -
ii#                 L       |   -
iippii#             S   *   |   -
iissiippii#         S   *   |   -
iissiissiippii#     S   *   |   -
ippii#              S       |   iippii#             S   *
issiippii#          S       |   iissiippii#         S   *
issiissiippii#      S       |   iissiissiippii#     S   *
miissiissiippii#    L       |   -
mmiissiissiippii#   L       |   -
pii#                L       |   -
ppii#               L       |   -
siippii#            L       |   -
siissiippii#        L       |   -
ssiippii#           L       |   -
ssiissiippii#       L       |   -

然后按照这个数组的顺序向下归并$L$型后缀，先处理了#向下归并产生i#，而i#产生ii#，ii#产生pii#:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
m  m  i  i  s  s  i  i  s  s  i  i  p  p  i  i  #
L  L  S  S  L  L  S  S  L  L  S  S  L  L  L  L  S
      *           *           *                 *

#                   S   *   |   #                   S   *
i#                  L       |   i#                  L   
ii#                 L       |   ii#                 L   
iippii#             S   *   |   -
iissiippii#         S   *   |   -
iissiissiippii#     S   *   |   -
ippii#              S       |   iippii#             S   *
issiippii#          S       |   iissiippii#         S   *
issiissiippii#      S       |   iissiissiippii#     S   *
miissiissiippii#    L       |   -
mmiissiissiippii#   L       |   -
pii#                L       |   pii#                L   
ppii#               L       |   -
siippii#            L       |   -
siissiippii#        L       |   -
ssiippii#           L       |   -
ssiissiippii#       L       |   -

接下来处理iippii#向下归并产生siippii#，iissiippii#向下归并产生siissiippii#，iissiissiippii#向下归并产生miissiissiippii#:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
m  m  i  i  s  s  i  i  s  s  i  i  p  p  i  i  #
L  L  S  S  L  L  S  S  L  L  S  S  L  L  L  L  S
      *           *           *                 *

#                   S   *   |   #                   S   *
i#                  L       |   i#                  L   
ii#                 L       |   ii#                 L   
iippii#             S   *   |   -
iissiippii#         S   *   |   -
iissiissiippii#     S   *   |   -
ippii#              S       |   iippii#             S   *
issiippii#          S       |   iissiippii#         S   *
issiissiippii#      S       |   iissiissiippii#     S   *
miissiissiippii#    L       |   miissiissiippii#    L   
mmiissiissiippii#   L       |   -
pii#                L       |   pii#                L   
ppii#               L       |   -
siippii#            L       |   siippii#            L   
siissiippii#        L       |   siissiippii#        L   
ssiippii#           L       |   -
ssiissiippii#       L       |   -

接下来类似地依次处理miissiissiippii#，mmiissiissiippii#，pii#，ppii#，siippii#，siissiippii#，而ipii#，ssiippii#和ssiissiippii#产生的是$S$子串先不管，得到:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
m  m  i  i  s  s  i  i  s  s  i  i  p  p  i  i  #
L  L  S  S  L  L  S  S  L  L  S  S  L  L  L  L  S
      *           *           *                 *

#                   S   *   |   #                   S   *
i#                  L       |   i#                  L   
ii#                 L       |   ii#                 L   
iippii#             S   *   |   -
iissiippii#         S   *   |   -
iissiissiippii#     S   *   |   -
ippii#              S       |   iippii#             S   *
issiippii#          S       |   iissiippii#         S   *
issiissiippii#      S       |   iissiissiippii#     S   *
miissiissiippii#    L       |   miissiissiippii#    L   
mmiissiissiippii#   L       |   mmiissiissiippii#   L
pii#                L       |   pii#                L   
ppii#               L       |   ppii#               L   
siippii#            L       |   siippii#            L   
siissiippii#        L       |   siissiippii#        L   
ssiippii#           L       |   ssiippii#           L   
ssiissiippii#       L       |   ssiissiippii#       L

然后倒着扫一遍用归并$L$型的方法来归并$S$型即可。注意由于一开始$LMS$型的位置是随便猜的(同组首字母末尾)，这里$LMS$型要作为$S$型重新排一遍。因为$S$型是字典序递减的，因此从下往上扫。处理完就得到了正确的后缀数组。

诱导排序完整代码:

1
2
3
4
5
6
#define Ar(x,a) SA[p[s[x]]a]=x
#define IS(s1)\
    memset(SA+1,0,n<<2); memset(c+1,0,m<<2); for(rint i = 1; i <= n; ++c[s[i++]]);\
    partial_sum(c+1,c+m+1,c+1); memcpy(p+1,c+1,m<<2); for(rint i = n1; i; Ar(s1[i],--), i--);\
    transform(c,c+m,p+1,[](int a){return++a;}); for(rint i = 1; i <= n; SA[i]>1&&t[SA[i]-1]?Ar(SA[i]-1,++):0, i++);\
    memcpy(p+1,c+1,m<<2); for(rint i = n; i; SA[i]>1&&!t[SA[i]-1]?Ar(SA[i]-1,--):0, i--);

回忆一下这样完成了第二次诱导排序: 已知按照$LMS$型后缀字典序递增存放的所有$LMS$型后缀(存储在s1[]里)，要求所有后缀的排名(SA[])。

那么第一次呢？第一次是已知按照乱序(初始顺序)存放的所有$LMS$型后缀位置(存储在tmp[]里)，要求所有$LMS$子串的排名(一个不完全的子串SA[])。

用tmp[]代替s1[]调用同一个算法就完成了。

按照之前算法的流程，我们会将$LMS$型后缀乱序放入同首字母末尾，然后归并$L$，最后重新归并$S$(这时候$LMS$子串被重新归并成正确顺序)。

为什么？直观理解一下，还是已之前的字符串为例，假设现在乱序排序$LMS$型后缀，然后归并$L$型:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
m  m  i  i  s  s  i  i  s  s  i  i  p  p  i  i  #
L  L  S  S  L  L  S  S  L  L  S  S  L  L  L  L  S
      *           *           *                 *

#                   S   *   |   #                   S   *
i#                  L       |   i#                  L   
ii#                 L       |   ii#                 L   
iippii#             S   *   |   -
iissiippii#         S   *   |   -
iissiissiippii#     S   *   |   -
ippii#              S       |   iissiippii#         S   *
issiippii#          S       |   iissiissiippii#     S   *
issiissiippii#      S       |   iippii#             S   *
miissiissiippii#    L       |   miissiissiippii#    L   
mmiissiissiippii#   L       |   mmiissiissiippii#   L   
pii#                L       |   pii#                L   
ppii#               L       |   ppii#               L   
siippii#            L       |   siissiippii#        L   
siissiippii#        L       |   siippii#            L   
ssiippii#           L       |   ssiissiippii#       L   
ssiissiippii#       L       |   ssiippii#           L

这里刻意打乱了初始$LMS$型后缀的插入顺序，这样最后几个$L$型的排序似乎是错误的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
m  m  i  i  s  s  i  i  s  s  i  i  p  p  i  i  #
L  L  S  S  L  L  S  S  L  L  S  S  L  L  L  L  S
      *           *           *                 *

#                   S   *   |   #                   S   *
i#                  L       |   i#                  L   
ii#                 L       |   ii#                 L   
iippii#             S   *   |   iippii#             S   *
iissiippii#         S   *   |   iissiissiippii#     S   *
iissiissiippii#     S   *   |   iissiippii#         S   *
ippii#              S       |   ippii#              S   
issiippii#          S       |   issiissiippii#      S   
issiissiippii#      S       |   issiippii#          S   
miissiissiippii#    L       |   miissiissiippii#    L   
mmiissiissiippii#   L       |   mmiissiissiippii#   L   
pii#                L       |   pii#                L   
ppii#               L       |   ppii#               L   
siippii#            L       |   siissiippii#        L   
siissiippii#        L       |   siippii#            L   
ssiippii#           L       |   ssiissiippii#       L   
ssiissiippii#       L       |   ssiippii#           L

我们发现这个错误没有影响最后的子串排序，注意我们要排序的不是iippii#，iissiippii#和iissiissiippii#，而是子串iippii，iiss和iiss。

会发现如果在$LMS$子串中有不同(例如iiss和iippii)，就会反映在下方的某个首字母不同上，pii#开始归并比siippii#和siissiippii#都要晚。

严谨论证可以看riteme的博客: 诱导排序与 SA-IS 算法。


SA-IS

最后注意下SAIS()处理出的SA[]数组里第一个是#要去掉，而且没有处理好Rank[]数组，需要手写:

1
SAIS(s,t,tmp,n,27); --n; for(rint i = 1; i <= n; Rank[SA[i]=SA[i+1]]=i, i++); CalHeight();

这里传入的27是字符集上界(SA-IS字符集只支持正整数，且需要保证#是最小的1号)，如果字符集过大或含有$0$或负数需要离散化(基数排序)，假设原串是str[]:

1
2
3
4
str[++n] = 0; m = *max_element(str+1,str+n+1);
for(rint i = 1; i <= n; Rank[str[i]] = 1, i++);
partial_sum(Rank,Rank+m+1,Rank);
for(rint i = 1; i <= n; s[i] = Rank[str[i]], i++);

如果原串字符集范围过大可能需要平移Rank[]下标或者将它开为std::map类型(然而这又会成为复杂度瓶颈2333)。

于是给出完整代码(没有离散化):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#include <bits/stdc++.h>
using namespace std;
#define MAXN 100000
#define rint register int
inline int rf(){int r;int s=0,c;for(;!isdigit(c=getchar());s=c);for(r=c^48;isdigit(c=getchar());(r*=10)+=c^48);return s^45?r:-r;}
int s[(MAXN<<1)+5], SA[MAXN+5], Rank[MAXN+5], Height[MAXN+5], c[MAXN+5], p[MAXN+5], tmp[MAXN+5], n, m; char str[MAXN+5]; bool t[(MAXN<<1)+5];
#define Ar(x,a) SA[p[s[x]]a]=x
#define IS(s1)\
    memset(SA+1,0,n<<2); memset(c+1,0,m<<2); for(rint i = 1; i <= n; ++c[s[i++]]);\
    partial_sum(c+1,c+m+1,c+1); memcpy(p+1,c+1,m<<2); for(rint i = n1; i; Ar(s1[i],--), i--);\
    transform(c,c+m,p+1,[](int a){return++a;}); for(rint i = 1; i <= n; SA[i]>1&&t[SA[i]-1]?Ar(SA[i]-1,++):0, i++);\
    memcpy(p+1,c+1,m<<2); for(rint i = n; i; SA[i]>1&&!t[SA[i]-1]?Ar(SA[i]-1,--):0, i--);
void SAIS(int s[], bool t[], int tmp[], int n, int m){
    int n1 = 0, m1 = Rank[1] = 0, *s1 = s+n; t[n] = 0;
    for(rint i = n-1; i; t[i] = s[i]^s[i+1]?s[i]>s[i+1]:t[i+1], i--);
    for(rint i = 2; i <= n; Rank[i] = t[i-1]&&!t[i]?tmp[++n1]=i,n1:0, i++); IS(tmp);
    for(rint i = 1, x, y = 0; i <= n; i++) if(x=Rank[SA[i]]){
        if(m1 <= 1 || tmp[x+1]-tmp[x]!=tmp[y+1]-tmp[y]) ++m1;
        else for(rint a = tmp[x], b = tmp[y]; a <= tmp[x+1]; a++, b++)
                if((s[a]<<1|t[a])^(s[b]<<1|t[b])){++m1; break;}
        s1[y=x] = m1;
    }
    if(m1 < n1) SAIS(s1,t+n,tmp+n1,n1,m1); else for(rint i = 1; i <= n1; SA[s1[i]] = i, i++);
    for(rint i = 1; i <= n1; s1[i] = tmp[SA[i]], i++); IS(s1);
}
void CalHeight(){for(rint i = 1, L = 0, j; i <= n; Height[Rank[i++]] = L) for(j = SA[~-Rank[i]], L-=L>0; i+L<=n&&j+L<=n&&s[i+L]==s[j+L]; ++L);}
int main(){
    scanf("%s",str+1); n = strlen(str+1); for(rint i = 1; i <= n; s[i] = str[i]-'a'+2, i++);
    s[++n] = 1; SAIS(s,t,tmp,n,27); --n; for(rint i = 1; i <= n; Rank[SA[i]=SA[i+1]]=i, i++); CalHeight();
    for(rint i = 1; i <= n; printf("%d%c",SA[i],"\n "[i<n]), i++);
    for(rint i = 2; i <= n; printf("%d%c",Height[i],"\n "[i<n]), i++); return 0;
}


比较

上述三份代码都可以直接通过UOJ#35 后缀排序。

在附加相同的读入优化且不使用输出优化情况下: 倍增方法约为$350\;\text{ms}$；SAMSA方法约为$330\;\text{ms}$；SA-IS方法约为$250\;\text{ms}$。

在附加相同的读入输出优化情况下: 倍增方法约为$210\;\text{ms}$；SAMSA方法约为$250\;\text{ms}$(!!!!!!??????)；SA-IS方法约为$100\;\text{ms}$。

在稳定的评测鸭上测试同一套数据得到的数据(这里应该是最慢数据点): 倍增方法$28.535\;\text{ms}$；SAMSA方法$46.304\;\text{ms}$(!!!!!!??????)；SA-IS方法约为$9.68\;\text{ms}$。

可以看到SA-IS明显表现优秀；SAMSA表现有点奇怪，而且不适用于大字符集。

但是在掌握SAM的基础上，SAMSA依然是我认为最简单最好理解的后缀数组算法。


扫描二维码即可在手机上查看这篇文章，或者转发二维码来分享这篇文章:


文章作者: Magolor
文章链接: https://magolor.cn/2019/08/13/2019-08-13-blog-01/
版权声明: 本博客所有文章除特别声明外，均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 Magolor！
OI资源黑科技后缀数组后缀树
扫描二维码在手机上查看或转发二维码以分享Magolor的博客