后缀数组学习笔记

后缀数组，即SA，是一种简单实用的字符串算法。当我们要处理一些关于子串的问题时，我们发现字符串中任意的一个子串都可以表示成该字符串一个后缀的前缀，因此我们只需要维护所有后缀的信息即可。基于这个思想，人们发明了后缀数组。

把所有的后缀按字典序排序

这是后缀数组算法的第一个任务。

我们定义：

• ( ext{suffix}(i))：以位置(i)开头的后缀；

• (sa[i])：排名为(i)的后缀的开头位置；

• (rk[i])：以位置(i)开头的后缀的排名。

当我们在对所有后缀按字典序排序时，我们实际上以每个后缀的第一个字母为第一关键字，第二个字母为第二关键字......对所有后缀进行比较。

因此我们先把所有的后缀按其第一个字母排序，当然，这其中肯定有首字母相同的情况，对于这种情况，我们暂时以后缀的开头位置为第二关键字。这样我们就得到了一个初步的排名。

接下来该处理每个后缀的第二个字母了，我们发现( ext{suffix}(i))的第二个字母就是( ext{suffix}(i+1))的第一个字母。即：( ext{suffix}(i))的第二关键字就是( ext{suffix}(i+1))在上一轮排序中的第一关键字。于是我们在这一轮排序中把( ext{suffix}(i))和( ext{suffix}(i+1))的排名分别作为第一和第二关键字合并起来，作为这一轮中( ext{suffix}(i))的排名。当(i+1>n)时，我们补(0)来替代。感性地理解，相当于现在每一个( ext{suffix})的排名都是一个两位数，十位上是原本( ext{suffix}(i))的排名，个位上是( ext{suffix}(i+1))的排名。当然，得到的这些“两位数”不一定连续，我们“离散化”一下就得到了这一轮的排名。

现在我们有了前两轮的排名，即按每个( ext{suffix}(i))的前两个字母做排名的结果。上一轮中，我们将后缀的排名一个一个合并。仿照上面的过程，现在我们将后缀的排名两个一组进行合并。想象一下：因为现在第(i)位上的关键字是( ext{suffix}(i))的前两个字符的排名，第(i+2)位置上的关键字是( ext{suffix}(i+2))的前两个字符的排名，这两个一合并，不就是( ext{suffix}(i))的前四个字符的排名吗？

同理，接下来我们每4个一组合并，然后8个一组合并。这就是一个倍增的过程！什么时候停止呢？显然是当所有的后缀排名各不相同的时候停止了。

考虑到我们总是合并“两位数”，因此我们可以用基数排序将排序的复杂度优化到(O(n))，总的复杂度变成优秀的(O(nlog n))。

参考代码：

struct Hou{
    int m,x[MAXN],y[MAXN],c[MAXN],sa[MAXN],rk[MAXN],ht[MAXN],wt[30];
    void init() {
        m=(int)'z';
        for(int i=1;i<=n;++i) c[x[i]=s[i]]++;
        for(int i=1;i<=m;++i) c[i]+=c[i-1];
        for(int i=n;i>=1;--i) sa[c[x[i]]--]=i;
        /*以每个后缀的开头字母为第一关键字,
        开头位置为临时第二关键字,得到一个初步的排名*/
        for(int k=1;k<=n;k<<=1) {/*倍增法开始*/
            int num=0;
            for(int i=n-k+1;i<=n;++i) y[++num]=i;
            /*y[i]表示第二关键字排名为i的数,第一关键字的位置
            第n-k+1到第n位第二关键字可以视为0,因此排在最前面*/
            for(int i=1;i<=n;++i) if(sa[i]>k) y[++num]=sa[i]-k;
            /*sa[i]>k,它可以作为别的第二关键字*/
            for(int i=1;i<=m;++i) c[i]=0;
            for(int i=1;i<=n;++i) c[x[i]]++;
            for(int i=1;i<=m;++i) c[i]+=c[i-1];
            for(int i=n;i>=1;--i) sa[c[x[y[i]]]--]=y[i],y[i]=0;
            /*体会这里对第一、第二关键字的基数排序*/
            for(int i=1;i<=n;++i) swap(x[i],y[i]);

            x[sa[1]]=1; num=1;
            for(int i=2;i<=n;++i) {
                x[sa[i]]=(y[sa[i]]==y[sa[i-1]] && y[sa[i]+k]==y[sa[i-1]+k])?num:++num;
            }
            if(num==n) break;
            m=num;
        }
        /*注意：x[]是第一关键字,是可以有重复的
        sa[]是当前的排名,不允许重复,强行构造了先后关系*/
    }
}SA;

有两个小细节值得讨论。一是如果有多测，那么清空时，for(int i=1;i<=max(n,(int)'z');++i)c[i]=0;。也就是c数组要清空的范围是[1,max(n,'z')]，只写n或'z'都是错的。二是y[n+1]这个位置也可能被访问到。所以要保证y[n+1]=0，多测时尤其要注意。

最基本的应用：求LCP

后缀数组最基本的应用就是求LCP，也即最长公共前缀。

我们定义( ext{lcp}(i,j))为( ext{suffix}(sa[i]))和( ext{suffix}(sa[j]))的最长公共前缀，即排名为(i)的后缀和排名为(j)的后缀的最长公共前缀。

我们来研究一下LCP有什么性质。显然的两条是：

1.( ext{lcp}(i,j)= ext{lcp}(j,i))

2.( ext{lcp}(i,i)=n-sa[i]+1)

有了这两条性质，我们接下来只需要讨论(i<j)的情况即可。

尝试使用我们上面辛辛苦苦求出的关于字典序的排名的特殊性质。仔细思考不难发现：

3.( ext{lcp}(i,k)=min( ext{lcp}(i,j), ext{lcp}(j,k))quad (1leq ileq jleq kleq n))

根据性质3我们可以有一个推论：( ext{lcp}(i,k)=min{ ext{lcp}(j,j-1)}quad (j<ileq k))

于是我们设：( ext{height}[i]= ext{lcp}(i,i-1))，显然( ext{height}[1]=0)。

如果我们能正确地求出( ext{height})数组，那么( ext{lcp})问题就转化为了一个区间最小值问题，显然可以通过ST表的预处理做到(O(1))回答。

怎么求( ext{height})？还是从其特殊性质入手。经过仔细观察，合理猜想，大胆假设，小心求证，我们发现：

[ ext{height}[rk[i]]geq ext{height}[rk[i-1]]-1 ]

正确性请读者自行意会。对证明感兴趣的请戳这里。

然后我们就可以线性推出( ext{height})啦！

参考代码：

void LCPinit() {
    /*预处理出ht数组*/
    for(int i=1;i<=n;++i) rk[sa[i]]=i;
    int k=0;
    for(int i=1;i<=n;++i) {
        if(rk[i]==1) continue;
        if(k) --k;/*ht[rk[i]]>=ht[rk[i-1]]-1*/
        int j=sa[rk[i]-1];
        while(j+k<=n && i+k<=n && s[i+k]==s[j+k]) ++k;
        ht[rk[i]]=k;
    }
}

简单的应用

例题1：POJ1743 Musical Theme

根据题意，我们把给出的数列两两做差，就把问题转化为求不可重叠的最长重复子串。

我们二分一个长度(k)，则要求是否存在(|x-y|geq k)使得( ext{suffix}(x))和( ext{suffix}(y))的最长公共前缀(geq k)，也即：( ext{lcp}(rk[x],rk[y])geq k)。

我们根据LCP的性质，( ext{lcp}(rk[x],rk[y])=min{ ext{height}[i]}(x<ileq y))
。

于是我们现在要找一段连续的( ext{height})使得这段( ext{height})的最小值(geq k)且这段中所有对应的位置的最大值和最小值之差（即(y-x)）(geq k)。

具体实现时，只要( ext{height}[i]geq k)，我们就把(sa[i])和(sa[i-1])所在的集合合并，然后看是否存在一个集合使(max{sa[i]}-min{sa[i]}geq k)，只要存在就说明(k)可行，否则不可行。找到最大的可行的(k)。

参考代码：

//POJ1743
#include <iostream>
#define ll long long
#define pb push_back
#define mk make_pair
#define pii pair<int,int>
#define fst first
#define scd second
using namespace std;
inline ll read() {
	ll f=1,x=0;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
	while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}
	return x*f;
}
const int MAXN=2e4+5;
const int INF=1e9;
int n,s[MAXN];
struct houzhuishuzu {
	int m,x[MAXN],y[MAXN],c[MAXN],sa[MAXN],rk[MAXN],ht[MAXN];
	void init() {
		m=200;
		for(int i=1;i<=m;++i) c[i]=0;
		for(int i=1;i<=n;++i) c[x[i]=s[i]]++;
		for(int i=1;i<=m;++i) c[i]+=c[i-1];
		for(int i=n;i>=1;--i) sa[c[x[i]]--]=i;
		for(int k=1;k<=n;k<<=1) {
			int num=0;
			for(int i=n-k+1;i<=n;++i) y[++num]=i;
			for(int i=1;i<=n;++i) if(sa[i]>k) y[++num]=sa[i]-k;
			for(int i=1;i<=m;++i) c[i]=0;
			for(int i=1;i<=n;++i) c[x[i]]++;
			for(int i=1;i<=m;++i) c[i]+=c[i-1];
			for(int i=n;i>=1;--i) sa[c[x[y[i]]]--]=y[i],y[i]=0;
			for(int i=1;i<=n;++i) swap(x[i],y[i]);
			x[sa[1]]=1; num=1;
			for(int i=2;i<=n;++i) {
				x[sa[i]]=(y[sa[i]]==y[sa[i-1]] && y[sa[i]+k]==y[sa[i-1]+k])?num:++num;
			}
			if(num==n) break;
			m=num;
		}
		for(int i=1;i<=n;++i) rk[sa[i]]=i;
		int k=0;
		for(int i=1;i<=n;++i) {
			if(rk[i]==1) continue;
			if(k) --k;
			int j=sa[rk[i]-1];
			while(j+k<=n&&i+k<=n&&s[i+k]==s[j+k])++k;
			ht[rk[i]]=k;
		}
	}
	bool check(int k) {
		int mx=-INF,mn=INF;
		for(int i=2;i<=n;++i) {
			if(ht[i]>=k) {
				mx=max(mx,max(sa[i],sa[i-1]));
				mn=min(mn,min(sa[i],sa[i-1]));
				if(mx-mn>k) return true;
			} else {
				mx=-INF;mn=INF;
			}
		}
		return false;
	}
}SA;
int main() {
//	freopen("","r",stdin);
//	freopen("","w",stdout);
	ios::sync_with_stdio(0);cin.tie(0);/*syn加速*/
	while(true) {
		n=read();if(!n)return 0;
		for(int i=1;i<=n;++i) s[i]=read();
		for(int i=1;i<n;++i) s[i]=s[i+1]-s[i]+100;
		s[n]=0;
		SA.init();
		int l=0,r=n/2;
		while(l<r) {
			int mid=(l+r+1)>>1;
			if(SA.check(mid))l=mid;
			else r=mid-1;
		}
		if(l+1>=5) cout<<l+1<<endl;
		else cout<<0<<endl;
	}
	return 0;
}

高级的应用：SA+主席树

例题2：洛谷 P4094 [HEOI2016/TJOI2016]字符串

先二分答案。考虑check长度(K)是否能达到，即( ext{suffix}(a)... ext{suffix}(b-K+1))中是否存在一个后缀，和( ext{suffix}(c))的最长公共前缀大于等于(K)。

我们发现，如果把所有后缀按字典序排好，那么和( ext{suffix}(c))的最长公共前缀大于等于(K)的后缀一定在一段连续的区间里。即([rk[c]-a...rk[c]+b])的区间内。我们通过预处理ST表可以(O(1))询问两个后缀的LCP。再通过二分把([rk[c]-a...rk[c]+b])这个区间确定下来。

然后只要查询( ext{suffix}(a)... ext{suffix}(b-K+1))中是否有后缀的(rk)恰在我们二分出的范围内。这样一个“二维”的区间查询，我们可以用主席树解决。

总复杂度(O(nlog^2n))

参考代码：

//P4094
#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define mk make_pair
#define pii pair<int,int>
#define fst first
#define scd second
using namespace std;
inline ll read() {
    ll f=1,x=0;char ch=getchar();
    while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
    while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}
void print(int x)
{
    if(x<0)//负数
    {
        putchar('-');
        x=-x;
    }
    if(x>9) 
        print(x/10);
    putchar(x%10+'0'); 
}
const int MAXN=1e5+5;
const int INF=2e9;
int n,m;
char s[MAXN];
struct Hou {
	int m,x[MAXN],y[MAXN],c[MAXN],sa[MAXN],rk[MAXN],ht[MAXN],st[30][MAXN],LOG[MAXN];
	void init() {
		m=(int)'z';
		for(int i=1;i<=n;++i) c[x[i]=s[i]]++;
		for(int i=1;i<=m;++i) c[i]+=c[i-1];
		for(int i=n;i>=1;--i) sa[c[x[i]]--]=i;
		for(int k=1;k<=n;k<<=1) {
			int num=0;
			for(int i=n-k+1;i<=n;++i) y[++num]=i;
			for(int i=1;i<=n;++i) if(sa[i]>k) y[++num]=sa[i]-k;
			for(int i=1;i<=m;++i) c[i]=0;
			for(int i=1;i<=n;++i) c[x[i]]++;
			for(int i=1;i<=m;++i) c[i]+=c[i-1];
			for(int i=n;i>=1;--i) sa[c[x[y[i]]]--]=y[i],y[i]=0;
			for(int i=1;i<=n;++i) swap(x[i],y[i]);//swap(x,y);
			x[sa[1]]=1; num=1;
			for(int i=2;i<=n;++i) {
				x[sa[i]]=(y[sa[i]]==y[sa[i-1]]&&y[sa[i]+k]==y[sa[i-1]+k])?num:++num;
			}
			if(num==n) break;
			m=num;
		}
		for(int i=1;i<=n;++i) rk[sa[i]]=i;
		int k=0;
		for(int i=1;i<=n;++i) {
			if(rk[i]==1) continue;
			if(k)--k;//ht[rk[i]>=ht[rk[i-1]]-1;
			int j=sa[rk[i]-1];
			while(i+k<=n && j+k<=n && s[i+k]==s[j+k]) ++k;
			ht[rk[i]]=k;
		}
		for(int i=1;i<=n;++i) st[0][i]=ht[i];
		for(int j=1;(1<<j)<=n;++j) {
			for(int i=1;i+(1<<(j-1))-1<=n;++i) {
				st[j][i]=min(st[j-1][i],st[j-1][i+(1<<(j-1))]);
			}
		}
		
		LOG[0]=-1;
		for(int i=1;i<=n;++i) LOG[i]=LOG[i>>1]+1;
	}
	inline int lcp(int i,int j) {
		if(i>j) return INF;
		int lg=LOG[j-i+1];
		return min(st[lg][i],st[lg][j-(1<<lg)+1]);
	}
}SA;
struct Zhu {
	int root[MAXN],cnt;
	struct node {
		int ls,rs,sum;
	}t[MAXN*40];
	void updata(int &x,int y,int l,int r,int pos) {
		x=++cnt; t[x]=t[y]; t[x].sum++;
		if(l==r) return ;
		int mid=(l+r)>>1;
		if(pos<=mid) updata(t[x].ls,t[y].ls,l,mid,pos);
		else updata(t[x].rs,t[y].rs,mid+1,r,pos);
	}
	int query(int y,int x,int tl,int tr,int ql,int qr) {
		if(ql>qr || t[x].sum-t[y].sum==0) return 0;
		if(ql<=tl&&qr>=tr) return t[x].sum-t[y].sum;
		int mid=(tl+tr)>>1,res=0;
		if(ql<=mid) res+=query(t[y].ls,t[x].ls,tl,mid,ql,qr);
		if(qr>mid) res+=query(t[y].rs,t[x].rs,mid+1,tr,ql,qr);
		return res;
	}
}T;
int main() {
//	freopen("","r",stdin);
//	freopen("","w",stdout);
	ios::sync_with_stdio(0);cin.tie(0);/*syn加速*/
	n=read();m=read();
	scanf("%s",s+1);
	SA.init();
	for(int i=1;i<=n;++i) {
		T.updata(T.root[i],T.root[i-1],1,n,SA.rk[i]);
	}
	while(m--) {
		int a=read(),b=read(),c=read(),d=read();
		int L=0,R=min(b-a+1,d-c+1);
		while(L<R) {
			int M=(L+R+1)>>1;
			int lrange,rrange;
			int l=1,r=SA.rk[c];
			while(l<r) {
				int mid=(l+r)>>1;
				if(SA.lcp(mid+1,SA.rk[c])>=M) r=mid;
				else l=mid+1;
			}
			lrange=l;
			l=SA.rk[c];r=n;
			while(l<r) {
				int mid=(l+r+1)>>1;
				if(SA.lcp(SA.rk[c]+1,mid)>=M) l=mid;
				else r=mid-1;
			}
			rrange=l;
			if(T.query(T.root[a-1],T.root[b-M+1],1,n,lrange,rrange)) L=M;
			else R=M-1;
		}
		print(L);puts("");
	}
	return 0;
}