还在写倍增后缀数组? SA-IS算法了解一下~

shadowice1984

2019-02-27 07:25:39

算法·理论

引子

众所周知,写串串题是绕不开后缀数据结构的(有些时候还要上回文数据结构)

说到后缀数据结构,我就想起了后缀数组和后缀自动机,对于后缀自动机来讲,大家常用的构建算法是O(n\Sigma)的,如果使用map存储转移边,那么复杂度可以降到O(nlog\Sigma),使用动态扩张的hash表,理论上我们可以做到O(n)的期望复杂度,但是实际应用的时候相当吃hash函数而且常数也不小

但是对于后缀数组来讲,最常见的写法是大常数O(nlogn)的倍增法,这导致后缀数组经常被人扣上复杂度高的帽子,甚至有毒瘤出题人恶意开大数据范围来卡后缀数组……

我们也有一个DC3算法可以在O(n)的时间内构造后缀数组,但是想必大家刚学后缀数组的时候就被告知这个算法常数极大实现复杂不太适合当板子背

那么有没有一个实现简单,效率高常数小的后缀数组构建算法可以帮后缀数组洗掉复杂度高常数大的帽子呢?,答案是有的,就是sais算法

SA-IS

这个算法的名字大概是Suffix-Array-Induce-Sort的意思,induce sort指的是一种被称之为诱导排序的算法也是这个算法的精髓

sais算法的精髓在于它比较后缀的方式,它根据两个后缀首字母是否相同来比较两个后缀

如果两个后缀(i,j)的首字母不同,那么两个后缀可以直接比较,否则我们可以将两个后缀删掉他们的首字母递归比较(i+1,j+1)这两个后缀的大小

在接下来的叙述当中我们认为我们处理的字符串是字符数组a,a[i]表示字符串a的第i个字符

S型后缀和L型后缀

为了方便我们推导各种各样的性质我们先把一个字符串的后缀分成两类,S型后缀和L型后缀,并且为了规避各种各样的分情况讨论,我们在字符串的后面需要加一个哨兵字符“#”

S型后缀和L型后缀的定义也很简单

"#"是一个S型后缀

如果后缀ii+1i就是S型后缀,否则就是L型后缀

那么这个定义有什么用呢?根据这个定义我们可以得到两个神奇的结论

结论1:若a[i]和a[i+1]相同,则后缀i和后缀i+1的类型相同

证明很简单,两个后缀(i,i+1)的首字母相同,显然要去递归的比较(i+1,i+2),那么(i,i+1)(i+1,i+2)的大小关系是相同的,自然类型也就一致了

根据这个结论我们可以倒着扫一遍字符串确定所有后缀的类型

结论2:若a[i]和a[j]相同且i是S型后缀j是L型后缀,则i比j大

还是那个思想,首字母相同了要去递归的比较(i+1,j+1)

由于我们知到i是一个S型后缀而j是一个L型后缀,那么我们可以很快的推出a[i+1] \geq a[i],a[j+1] \leq a[j],所以说除非a[i+1]=a[i]=a[j+1]=a[j],我们都能得到i比j大的结论,

a[i+1]=a[j+1]的时候我们可以接着首字母删去递归下去,这样我们的结论就成立了

诱导排序(induce sort)

sais的核心是诱导排序,但是诱导排序的正确性是基于归纳法证明的.所以这个算法不太好理解,它很简单,但是你可能理解不了为什么它是对的

根据我们推出来的结论我们可以知道,在后缀数组当中首字母相同的后缀排成了连续的一段;而首字母相同的后缀里面,前面全部是L型的后缀后面全部是S型的后缀

那么我们想个办法,先把所有L型后缀排好,再把所有的S型后缀排好,这样就可以把所有的后缀排好序了

那么我们来考虑一下怎么把所有的L型后缀排好序,那么我们还是采用这样的方式来比较大小

如果两个L型后缀(i,j)的首字母不一样那么可以直接比较,否则我们递归的比较(i+1,j+1)

那么我们倒着想,假设我们已经有了一部分后缀的大小关系,然后我们每次取出手头已有的最小后缀i,执行这样一个操作,如果i-1是一个L型后缀,那么把i-1加入我们的集合,重复这个过程若干次,我们就能得到一部分L型后缀的大小关系

如果用堆来维护刚才的过程,复杂度带log无法接受,但是我们有一个巧妙的做法可以完成上述过程

我们直接维护sa数组,先处理出一个cnt数组表示首字母为i的后缀从sa数组的哪里结束(和你倍增法基数排序的时候原理很像),之后把sa数组全部初始化为-1

然后把我们倒着扫一遍手头里已经有的后缀全部插入到sa数组对应的位置(操作方式和倍增法的基排是一样的)

接下来我们从左向右扫一遍sa数组,如果第i位比1大并且sa(i-1)是一个L型后缀,我们就把L型后缀插入到对应的位置上

(你可以把插入的过程理解为将sa数组按照首字母划分为若干个桶,然后将这个后缀push到对应的桶里去,代码上写起来和倍增法的基排没啥区别)

如此这般我们就可以在O(n)的时间内完成上述算法

如果我们一开始掌握的信息足够多,比如说我们知道了全部的S型后缀排序结果,那么我们就可以用上面的算法还原出所有L型后缀的排序结果,原理很简单,在刚才的算法当中,每一个L型后缀i都是被后缀i+1加入的,因为是L型,所以在sa数组当中i肯定在i+1后面,所以没有一个L型后缀会被漏掉

接下来我们需要证明的就是对于任意的两个L型后缀ij,如果ij小那么i肯定比j先加入

证明如下:如果ij的首字母不同那么ij被分在了不同的桶里,大小关系一目了然,否则两个后缀的大小关系依赖于被塞到桶里的时间,也就是依赖于(i+1,j+1)的大小关系,符合我们递归比较的原理

同理我们知道了全部的L型后缀的排序结果,我们只需要把刚才的算法中的正序扫描改成倒序扫描就可以还原S型后缀的全部信息

这样只要我们一开始掌握的后缀的大小信息无误并且足够的多,我们就能正确的还原出一个有序的sa数组,具体点,我们知道了关于S型的大小信息可以生成L型后缀的信息,知道了L型后缀的信息也可以还原出S型后缀的所有信息

递归,lms子串与lms后缀

诱导排序算法允许我们从L型后缀诱导到S型后缀上,也允许我们从S型后缀诱导到L型后缀上,但是我们从S型后缀诱导到L型后缀的时候并不需要全部的S型后缀的大小信息,相反,我们仅仅需要一部分S型后缀的信息照样可以还原出L型后缀的信息

把诱导排序算法的正确性证明念一遍,我们会发现关键在于递归比较(i+1,j+1)的大小关系这一段,如果ij都是L型后缀,我们不能保证(i+1,j+1)都是L型后缀

但是倒过来看,我们也仅仅需要那些左边就是L型后缀的S型后缀的大小信息了

我们在这里称一个后缀为lms(Left-Most-Suffix)后缀当且仅当这个后缀左边是L型后缀而它本身是一个S型后缀,根据这个定义,还原所有的L型后缀的排序关系仅仅需要所有lms后缀的大小关系

那么我们可以把刚才的诱导排序算法升级一下,我们先搞出来所有lms后缀的大小关系,然后当做"种子"生成L型后缀的大小关系,然后从L型后缀的大小关系反过来生成S型后缀的大小关系,如此这般我们就从lms后缀的排序结果诱导出了我们想要的后缀数组

那么一个显然的事实是长度为n的字符串最多有\frac{n}{2}个lms后缀,如果我们可以想办法将这个问题递归下去,我们就可以得到一个复杂度满足以下递归式的算法

T(n)=T(\frac{n}{2})+O(n)=O(n)

那么问题来了我们怎么递归问题呢?答案是用lms子串

我们定义lms子串表示从左向右数第i个lms后缀和第i+1个lms后缀之间的字符串

给个例子:

Num 1  2  3  4  5  6  7  8  9  10 11 12 13 14 15 16 17
Str m  m  i  i  s  s  i  i  s  s  i  i  p  p  i  i  #
Typ L  L  S  S  L  L  S  S  L  L  S  S  L  L  L  L  S
Lms       *           *           *

对于字符串“mmiissiissppii“来讲,lms子串指的就是[3,7],[7,11],[11,17]这三个子串

注意相邻的两个lms子串会重叠一个字符

那么我们的递归方法是,将lms子串离散化之后按顺序拼起来形成一个新的字符串,比如说在上面的例子当中我们会得到一个字符串“221”,

那么我们这里有一个结论是对于原串的两个lms后缀p1,p2,如果

看起来挺显然的,不就是把一块字符缩成了一个字符嘛…… 但是其实并没有我们想的那么显然,我们接下来会严谨的证明为什么刚才的结论是对的 首先先明确一下比较lms子串的规则,我们比较两个lms子串的时候对于每一个字符同时比较字符的字典序大小和所在位置后缀的SL类型,举个例子,我们认为(a,S)和(a,L)是不相等的,尽管字符一样但是所在位置的S和L类型却不一样,同时我们认为S类型大于L类型 假设我们有两个lms后缀$p1,p2$,现在我们想要比较两个后缀的大小,我们可以根据两个后缀在新字符串上第一位的字符(不妨设为第$i$位和第$j$位)是否相等来进行分情况讨论 ### case1:第一位的字符是相等的 我们发现新字符串上两个字符相等,等价于原来字符串上两个lms子串相等,那么既然两个lms子串相等了那么这两个子串的长度$len$必然是一样的,因此我们可以将两个后缀同时去掉长度为$len-1$的部分递归比较剩下的部分,在新字符串上也就是递归到了第$i+1$位和第$j+1$位,符合我们字符串比较的逻辑 ### case2:第一位的字符是不等的 那么你可能会说:既然新字符串的两个字符都不一样了,就说明在原来的字符串上肯定有不一样的位,所以这种情况是显然的 那好像我们忽视了一种情况哎,两个lms子串代表的字符不一样是因为其中一个是另一个的前缀 这种情况下我们的算法就会锅掉,因为两个后缀明明没有分出大小我们就强行钦定一个大于另一个了 事实果真如此吗? 答案是这种情况不存在,我们比较lms子串的规则是同时比较字符和SL类型的后缀,在这种情况下不存在一个lms子串是另一个的前缀,因为一个lms子串的SL类型必须是形如SSSSLLLLS这样的,前面是一串S中间是一串L,最后是恰好一个S,那么我们发现如果一个是另一个的前缀,两个串必然会在短的字符串的结尾处SL类型不匹配,从而不存在这种情况 ------ 那么此时我们就可以放心大胆把所有lms子串进行离散化,然后递归进行sais算法了,看起来我们得到了一个优秀的$O(n)$算法? ### 离散化lms子串与诱导排序 现在我们只差最后一步就能实现完整的sais算法了,也就是离散所有的lms子串的过程 如果我们实现这一步使用的是**基数排序**那么我们将会得到传说中的KA-algorithm,但是这东西很慢啊……我们想想有没有什么优秀的排序做法 答案还是诱导排序,传统的诱导排序要求我们传一个有序的lms后缀数组进去,但是现在我们可以乱序传入lms后缀数组,当然诱导排序会返回一个错误的结果, 但是我们断言,在返回的后缀数组当中关于lms后缀的**子序列**是按照每个lms后缀**开头的lms子串**为关键字进行排序的 具体点来讲依然使用归纳法证明算法正确性 我们知道初始的时候的lms数组是乱的,因此归纳的初始条件直接被打破了 但是如果我们换个命题进行归纳呢? 我们一开始传进去的lms数组,如果只看开头的第一个字母,在后缀数组当中是有序的 那么在第一轮诱导过后,对于每个L型后缀,如果只看这个后缀的开头到第一个S型字符位置的位置,他们也是有序的 在第二轮诱导过后,对于每个S型后缀,它们有序的部分恰好是这个后缀的开头到第一个lms字符所在的位置,如果这个后缀是一个lms型后缀那么这个后缀有序的部分是这个后缀开头的lms字符 所以我们证明了只需要再来一轮诱导排序就可以给所有lms子串排好序,然后我们就可以很轻松的完成离散化啦~ 至此我们完成了SA-IS算法的最后一块拼图,此时我们已经可以尝试实现sais算法了~ ## 实现 由于sais算法的实现难度不小(主要是如何避免冗余代码上),这里给出一个范例代码和简略的注释 实际上实现精妙的话sais算法的代码量相当小,甚至和倍增法差不多 (基本是抄的uoj的[这份提交记录](http://uoj.ac/submission/309676)) ```Cpp /****************************************** 这里的诱导排序函数由于需要传非常多的参数,所以使用了宏的方式来实现了一个函数 具体的实现方式可以看下面的inds(lms)这个宏 由于多行宏旁边写不了注释我会在行上面写注释 sais的算法流程 : 1.确定每个后缀的类型 2.确定第i个lms后缀的编号以及编号为i的后缀是第几个lms后缀 3.进行一次诱导排序给lms子串排序 4.扫描一遍得到的sa数组,生成用于递归的字符串 5.如果新字符串中的每一个字符都不同,后缀数组可以直接计算,此处不递归 6.否则递归进行一次sais 7.根据之前记录了第i个lms后缀的编号和返回的后缀数组计算lms后缀的顺序 8.进行一次诱导排序并返回 ******************************************/ #include <cstdio> #include <algorithm> using namespace std; const int N = 1e6 + 10; typedef long long ll; template <typename T> inline void gm(T*& bas, int siz, T*& op) { op = bas; bas += siz;// 这里把内存提前申请好然后用指针分配内存 } //这里sum表示桶的位置,cur是sum的一份copy,避免每次都要做前缀和 // S型后缀的桶是倒着开的也要倒着扫,L型后缀的桶的是正着开的也要正着扫 #define pus(x) (sa[cur[a[x]]--] = x) // 插入S型后缀 #define pul(x) (sa[cur[a[x]]++] = x) //插入诱导排序 //诱导排序的宏,进行了两轮诱导过程:分别是从lms诱导到L和从L诱导到S #define inds(lms) \ for (int i = 1; i <= n; i++) sa[i] = -1; \ for (int i = 1; i <= n; i++) sum[i] = 0; \ for (int i = 1; i <= n; i++) sum[a[i]]++; \ for (int i = 1; i <= n; i++) sum[i] += sum[i - 1]; \ for (int i = 1; i <= n; i++) cur[i] = sum[i]; \ for (int i = m; i >= 1; i--) pus(lms[i]); \ for (int i = 1; i <= n; i++) cur[i] = sum[i - 1] + 1; \ for (int i = 1; i <= n; i++) \ if (sa[i] > 1 && !tp[sa[i] - 1]) \ pul(sa[i] - 1); \ for (int i = 1; i <= n; i++) cur[i] = sum[i]; \ for (int i = n; i >= 1; i--) \ if (sa[i] > 1 && tp[sa[i] - 1]) \ pus(sa[i] - 1); int sa[N]; int sum[N]; int cur[N]; int rk[N]; int A_bas[N << 4]; int* A_t; inline void sais(int n, int* a) { int* tp; //现开数组现分配内存 gm(A_t, n + 1, tp); int* p; gm(A_t, n + 2, p); tp[n] = 1; //倒着扫一遍确定后缀类型 for (int i = n - 1; i >= 1; i--) tp[i] = (a[i] == a[i + 1]) ? tp[i + 1] : (a[i] < a[i + 1]); //确定第i个lms后缀的编号和编号为i的后缀是第几个lms后缀 int m = 0; for (int i = 1; i <= n; i++) rk[i] = (tp[i] && !tp[i - 1]) ? (p[++m] = i, m) : -1; //进行一轮诱导排序,生成新的字符串 inds(p); int tot = 0; int* a1; gm(A_t, m + 1, a1); p[m + 1] = n; //扫一遍对lms子串进行离散化 for (int i = 1, x, y; i <= n; i++) if ((x = rk[sa[i]]) != -1) { if (tot == 0 || p[x + 1] - p[x] != p[y + 1] - p[y]) //如果长度不一致显然不等 tot++; //否则暴力for一遍比较两个lms串是否相等 else for (int p1 = p[x], p2 = p[y]; p2 <= p[y + 1]; p1++, p2++) if ((a[p1] << 1 | tp[p1]) != (a[p2] << 1 | tp[p2])) { tot++; break; } a1[y = x] = tot; } if (tot == m) for (int i = 1; i <= m; i++) sa[a1[i]] = i; //如果字符互不相同就直接计算后缀数组,否则递归 else sais(m, a1); for (int i = 1; i <= m; i++) a1[i] = p[sa[i]]; //还原lms子串的顺序,进行诱导排序 inds(a1); } char mde[N]; int n; int a[N]; int tr[300]; char buf[20]; int cnt; int main() { A_t = A_bas; scanf("%s", mde + 1); while (mde[n + 1] != '\0') n++; for (int i = 1; i <= n; i++) tr[mde[i]] = 1; for (int i = 1; i < 300; i++) tr[i] += tr[i - 1]; for (int i = 1; i <= n; i++) a[i] = tr[mde[i]] + 1; a[++n] = 1; sais(n, a); for (int i = 2; i <= n; i++) { int tmp = sa[i]; while (tmp) buf[++cnt] = tmp % 10 + 48, tmp /= 10; while (cnt) putchar(buf[cnt--]); putchar(' '); } return 0; //拜拜程序~ } ``` ## 例题 大部分时候后缀数组都需要套数据结构所以用到$O(n)$后缀数组的时候不多…… 但是不能说没有这种题啊,比如说这道题 ## [codeforcesgym102114 J Just So You know](https://codeforces.com/gym/102114/problem/J) 一句话题意是求一个字符串所有子串的哈夫曼树代价和 那么我们可以借助SAIS算法$O(n)$构建后缀树,统计出每个子串出现了多少次,由于出现次数小于n的节点的权值只有$O(n)$种,而权值大于$n$的节点的数目不超过$O(n)$,那么对于前面一种节点我们拿数组存,后面一种节点拿队列维护一下就可以做到线性的复杂度了~ 这里是代码~ ```Cpp #include<cstdio> #include<algorithm> using namespace std;const int N=1e6+10;typedef long long ll; template <typename T>inline void gm(T*& bas,int siz,T*& op){op=bas;bas+=siz;} inline ll gcd(ll a,ll b){if(a<b)swap(a,b);while(b){ll c=a%b;a=b;b=c;}return a;} int st[N];int tp;int siz[N];ll cnt[N];int tr[300]; inline void push(int len,int si) { int lst=tp;int psiz=0;while(tp&&st[tp]>len)tp--; if(tp!=lst) { for(int j=lst;j>=tp+2;j--) {siz[j-1]+=siz[j],cnt[siz[j]]+=st[j]-st[j-1];} psiz=siz[tp+1];cnt[siz[tp+1]]+=st[tp+1]-len; }if(len==st[tp]){siz[tp]+=psiz+si;}else {st[++tp]=len,siz[tp]=psiz+si;} } namespace suffixarray { # define pus(x) (sa[cur[a[x]]--]=x) # define pul(x) (sa[cur[a[x]]++]=x) # define inds(lms)\ for(int i=1;i<=n;i++)sa[i]=-1;for(int i=1;i<=n;i++)sum[i]=0;\ for(int i=1;i<=n;i++)sum[a[i]]++;for(int i=1;i<=n;i++)sum[i]+=sum[i-1];\ for(int i=1;i<=n;i++)cur[i]=sum[i];for(int i=m;i>=1;i--)pus(lms[i]);\ for(int i=1;i<=n;i++)cur[i]=sum[i-1]+1;\ for(int i=1;i<=n;i++)if(sa[i]>1&&!tp[sa[i]-1])pul(sa[i]-1);\ for(int i=1;i<=n;i++)cur[i]=sum[i];\ for(int i=n;i>=1;i--)if(sa[i]>1&&tp[sa[i]-1])pus(sa[i]-1); int sum[N];int cur[N];int A_bas[N<<4];int* A_t;int sa[N];int rk[N];int h[N]; inline void sais(int n,int* a) { int* tp;gm(A_t,n+1,tp);int* p;gm(A_t,n+2,p);tp[n]=1; for(int i=n-1;i>=1;i--)tp[i]=(a[i]==a[i+1])?tp[i+1]:(a[i]<a[i+1]); int m=0;for(int i=1;i<=n;i++)rk[i]=(tp[i]&&!tp[i-1])?(p[++m]=i,m):-1; inds(p);int* a1;gm(A_t,m+1,a1);p[m+1]=n;int tot=0; for(int i=1,x,y;i<=n;i++) if((x=rk[sa[i]])!=-1) { if(tot==0||p[x+1]-p[x]!=p[y+1]-p[y])tot++; else for(int p1=p[x],p2=p[y];p2<=p[y+1];p1++,p2++) if((a[p1]<<1|tp[p1])!=(a[p2]<<1|tp[p2])){tot++;break;} a1[y=x]=tot; }if(tot==m){for(int i=1;i<=m;i++)sa[a1[i]]=i;}else sais(m,a1); for(int i=1;i<=m;i++)a1[i]=p[sa[i]];inds(a1); } inline void create_sa(int n,int* a) { for(int i=0;i<150;i++)tr[i]=0;for(int i=1;i<=n;i++)tr[a[i]]=1; for(int i=1;i<150;i++)tr[i]+=tr[i-1]; for(int i=1;i<=n;i++)a[i]=tr[a[i]]+1;a[++n]=1;sais(n,a);n--; for(int i=1;i<=n;i++)sa[i]=sa[i+1];for(int i=1;i<=n;i++)rk[sa[i]]=i; for(int i=1,j=0,k=0;i<=n;h[rk[i++]]=k) for(k=k?k-1:k,j=sa[rk[i]-1];a[i+k]==a[j+k];k++); for(int i=1;i<=n;i++) {if(i-1)push(h[i],0);push(n-sa[i]+1,1);} while(tp)siz[tp-1]+=siz[tp],cnt[siz[tp]]+=st[tp]-st[tp-1],tp--; } inline void clr(){++A_t;for(int* p=A_bas;p!=A_t;++p)*p=0;A_t=A_bas;} } struct data{ll we;int ct;}q[N<<3];int hd;int tl;int n;int a[N];ll ans; inline void mg(int w1,int w2,ll ct) { ans+=(w1+w2)*(ll)ct;cnt[w1]-=ct;cnt[w2]-=ct; if(w1+w2<=n)cnt[w1+w2]+=ct;else q[++tl]=(data){w1+w2,ct}; } inline void mg1(ll we,int ct){ans+=(ll)we*ct;q[++tl]=(data){we,ct};} inline void solve() { scanf("%d",&n);for(int i=1;i<=n;i++)scanf("%d",&a[i]); ll crc=0;suffixarray::create_sa(n,a);hd=2;tl=1; for(int i=1;i<=n;i++)crc+=(ll)i*cnt[i]; for(int i=1,lst=0;i<=n;i++) { if(cnt[i]==0)continue;if(cnt[lst])mg(lst,i,1); if(cnt[i]/2)mg(i,i,cnt[i]/2);lst=i; }for(int i=1;i<=n;i++) if(cnt[i]==1) { if(hd<=tl)mg1(i+q[hd].we,1);q[hd].ct--; if(q[hd].ct==0)hd++; } while(tl-hd>0) { if(q[hd].ct>1)mg1(q[hd].we*2,q[hd].ct/2),q[hd].ct&=1; else mg1(q[hd].we+q[hd+1].we,1),q[++hd].ct--; if(q[hd].ct==0)hd++; }ll ansq=(ll)n*(n+1)/2;ll d=gcd(ansq,ans);ansq/=d;ans/=d; printf("%I64d",ans);if(ansq!=1)printf("/%I64d",ansq);printf("\n"); } inline void clear(){for(int i=1;i<=n;i++)cnt[i]=0;ans=0;suffixarray::clr();} int main() { suffixarray::A_t=suffixarray::A_bas; int T;scanf("%d",&T);for(int z=1;z<=T;z++)solve(),clear();return 0; } ``` 还有一个例题就是[P5112](https://www.luogu.org/problemnew/show/P5112) 虽然说这题可以拿hash或者后缀自动机搞过去,但是使用SAIS的后缀数组应该是这题最快的稳定做法了(毕竟hash是一个不稳定做法,很容易被卡住),想练手的同学可以尝试着写一写这题