zoukankan      html  css  js  c++  java
  • SA 复习笔记

    大家好,由于蒟蒻 tzc 最近被动态点分治这个学也学不会的毒瘤玩意儿虐得不轻,所以就准备换换脑筋来 Van 同样学也学不会的后缀数组了。

    考虑一个非常经典的问题:【模板】后缀排序

    一些定义(very important,不要搞混):

    • (sa_i) 表示字典序排在第 (i) 的字符串是谁,也就是一个排名 ( o) 下标的映射
    • (rk_i) 表示后缀 (i)((注:在本文中,后缀 (i) 就是开头位置为 (i) 的后缀)的排名是啥,也就是一个下标 ( o) 排名的映射

    字符串比较是 (mathcal O(n)) 的,所以暴力排序复杂度是 (n^2log n) 的,显然不是我们所满足的。

    考虑模拟暴力排序的过程,我们以这 (n) 个后缀的第一位为第一关键字,再以它们的第二位为第二关键字,再以它们的第三位为第三关键字……如果后缀长度不足 (n) 就在它后面补上 (c_0),其中 (c_0) 为某个小于 (s) 中所有字符的字符,比如说 (!)

    不难发现整个过程可以分“阶段”进行,我们记第 (k) 阶段为比较完前 (k) 个字符后的结果,设 (rk_i(k)) 表示将所有后缀的长度为 (k) 的前缀排序后,后缀 (i) 的排名。如果有相同那么它们的 (rk) 也应该相同

    很明显暴力排序的过程就是每次从第 (k) 阶段转移到第 (k+1) 阶段,具体方法是对每个 (i)((rk_i(k),s_{i+k})) 当作一个 pair 进行排序,如果 (i+k>n) 那么 (s_{i+k}=0),即上文中所说的 (c_0)

    那么我们考虑,每次能不能不从第 (k) 阶段转移到第 (k+1) 阶段,而是从(k) 阶段转移到第 (2k) 阶段呢

    答案是肯定的。

    我们从第 (k) 阶段转移到第 (k+1) 阶段的思路是:一个长度为 (k+1) 个字符串可以断成一个长度为 (k) 的子串和一个字符,而这个长度为 (k) 的字符串的大小关系我们是知道的,单个字符的大小关系我们也是知道的,所以可以转移。

    同样道理一个长度为 (2k) 个字符串可以断成两个长度为 (k) 的子串,而这两个长度为 (k) 的字符串的大小关系我们也是知道的,所以也可以转移。

    具体来说,你将 ((rk_i(k),rk_{i+k}(k))) 当作一个 pair 进行排序,若 (i+k>n)(rk_{i+k}(k)=0),就可以得到 (rk_{i}(2k))

    算下时间复杂度:倍增 1 log,sort n log,总复杂度 (nlog^2n)

    for(int i=1;i<=n;i++) x[i]=mp(mp(s[i],0),i);
    sort(x+1,x+n+1);int gr=0;//stage 1
    for(int i=1;i<=n;i++){
    	if(x[i].fi!=x[i-1].fi) gr++;//note that if s[i]=s[j],rk(1,i)=rk(1,j), so we have to do a process similar to discretion
    	rk[x[i].se]=gr; 
    }
    for(int k=1;k<n;k<<=1){//from stage k to stage 2k
    	for(int i=1;i<=n;i++){
    		if(i+k<=n) x[i]=mp(mp(rk[i],rk[i+k]),i);//sort the pairs of (rk(i,k),rk(i+k,k))
    		else x[i]=mp(mp(rk[i],0),i);//i+k>n
    	}
    	sort(x+1,x+n+1);gr=0;
    	for(int i=1;i<=n;i++){
    		if(x[i].fi!=x[i-1].fi) gr++;
    		rk[x[i].se]=gr; 
    	}
    }
    

    然而这个程序是 TLE 90 的,这说明有比 (nlog^2n) 更优的算法。那么怎么进一步优化呢?

    不难发现在我们所有 pair 中的两个数最大不过 (n),于是我们回忆起一个东西叫做 gay♂ 基数排序。

    基数排序是一种能够在 (mathcal O(n+m)) 的时间内对某个数组进行排序的算法,其中 (m)数组的值域

    感觉 SA 鸽了这么久就是因为不会这玩意儿啊 qwq。

    先考虑一维的基数排序,我们及一个桶 (c_i) 表示值为 (i)(a_i) 有多少个。然后对 (c_i) 做一遍前缀和。

    显然,如果 (a_i) 互不相同排名为 (c_{a_i}) 的数就是 (a_i),因为 (c_{a_i}) 就是比 (a_i) 小的个数。

    那如果有 (a_i) 相同怎么办呢?那显然排名在 ((c_{a_i-1},c_{a_i}]) 之间的数都等于 (a_i)

    考虑从后到前扫一遍,每一次将 (a_i) 的排名设为 (c_{a_i}) 并令 (c_{a_i})(1),这样你发现所有相同的 (a_i) 都被赋上了一个不同的排名,并且对于相同的 (a_i),越靠后的排名越靠后,也就是我们对 ((a_i,i)) 这个 pair 排了一遍序。

    回过头来,如果是二维的,例如在后缀排序中遇到的这些 pair 怎么办呢?

    考虑我们在刚才的基数排序中为什么要从后到前扫,因为这样一来如果有 (a_i) 相同的那么靠后的排名永远比考前的排名高,也就是按照第二维降序的顺序更新桶。那这样一来就好办了,先将所有 (y_i) 排个序,然后按照 (y_i) 从大到小的顺序更新桶。就可以实现 (mathcal O(n)) 排序了。

    注意到如果有相同的 pair 那么它们的 (rk) 也应该相同,所以排完序后还要做一遍类似于离散化的操作。

    const int M=123;//the ASCII of 'z' is 122
    pii x[MAXN+5];
    int n,sa[MAXN+5],rk[MAXN+5],buc[MAXN+5],seq[MAXN+5];
    void get_sa(){
    	for(int i=1;i<=n;i++) buc[s[i]]++;
    	for(int i=1;i<=M;i++) buc[i]+=buc[i-1];
    	for(int i=n;i;i--) sa[buc[s[i]]--]=i;
    	int gr=0;
    	for(int i=1;i<=n;i++){
    		if(s[sa[i]]!=s[sa[i-1]]) gr++;
    		rk[sa[i]]=gr;
    	}//calculate the ranks of stage 1
    	for(int k=1;k<n;k<<=1){//from stage k to stage 2k
    		for(int i=1;i<=n;i++){
    			if(i+k<=n) x[i]=mp(rk[i],rk[i+k]);//sort the pairs of (rk(i,k),rk(i+k,k))
    			else x[i]=mp(rk[i],0);
    		}
    		memset(buc,0,sizeof(buc));//clear the bucket
    		for(int i=1;i<=n;i++) buc[x[i].se]++;//sort the values of the second dimension
    		for(int i=1;i<=n;i++) buc[i]+=buc[i-1];
    		for(int i=n;i;i--) seq[buc[x[i].se]--]=i;//seq[i] means the i-smallest index sorted by the second dimension
    		memset(buc,0,sizeof(buc));
    		for(int i=1;i<=n;i++) buc[x[i].fi]++;//sort the values of the first dimension
    		for(int i=1;i<=n;i++) buc[i]+=buc[i-1];
    		for(int i=n;i;i--) sa[buc[x[seq[i]].fi]--]=seq[i];//update the bucket in descending order of the second dimension
    		gr=0;
    		for(int i=1;i<=n;i++){
    			if(x[sa[i]]!=x[sa[i-1]]) gr++;//note that if x[i]=x[j],rk(i)=rk(j), so we have to do a process similar to discretion
    			rk[sa[i]]=gr;
    		}
    	}
    }
    

    这样过是过了,但是我们发现最慢的一个点跑了整整 1s!显然常数上天了,我们考虑进行常数优化。

    首先我们发现其实不用用桶第二维进行排序。对于 (i=n-k+1sim n) 显然它们的第二维都是 (0),应该排在最前面;对于 (ileq n-k) 我们可以用之前求得的 (k) 阶段的排名算出,所以 2 遍循环就行了。

    现在最慢的一个点跑了 750ms,大约省去了 (1/4) 的常数

    void get_sa(){
    	for(int i=1;i<=n;i++) buc[s[i]]++;
    	for(int i=1;i<=M;i++) buc[i]+=buc[i-1];
    	for(int i=n;i;i--) sa[buc[s[i]]--]=i;
    	int gr=0;
    	for(int i=1;i<=n;i++){
    		if(s[sa[i]]!=s[sa[i-1]]) gr++;
    		rk[sa[i]]=gr;
    	}
    //	for(int i=1;i<=n;i++) printf("%d%c",rk[i],(i==n)?'
    ':' ');
    	for(int k=1;k<=n;k<<=1){
    		for(int i=1;i<=n;i++){
    			if(i+k<=n) x[i]=mp(rk[i],rk[i+k]);
    			else x[i]=mp(rk[i],0);
    		}
    		int num=0;
    		for(int i=n-k+1;i<=n;i++) seq[++num]=i;
    		for(int i=1;i<=n;i++) if(sa[i]>k) seq[++num]=sa[i]-k;
    		memset(buc,0,sizeof(buc));
    		for(int i=1;i<=n;i++) buc[x[i].fi]++;
    		for(int i=1;i<=n;i++) buc[i]+=buc[i-1];
    		for(int i=n;i;i--) sa[buc[x[seq[i]].fi]--]=seq[i];
    		gr=0;
    		for(int i=1;i<=n;i++){
    			if(x[sa[i]]!=x[sa[i-1]]) gr++;
    			rk[sa[i]]=gr;
    		}
    	}
    }
    

    继续卡!我们发现,并不是每次的值域都能达到 (n) 那么大,所以可以记录一个 (vmax) 表示当前阶段的值域,更新桶的时候只用对 ([1,vmax]) 做前缀和。

    现在最慢的一个点跑了 600ms:

    void get_sa(){
    	int vmax=123;
    	for(int i=1;i<=n;i++) buc[s[i]]++;
    	for(int i=1;i<=vmax;i++) buc[i]+=buc[i-1];
    	for(int i=n;i;i--) sa[buc[s[i]]--]=i;
    	int gr=0;
    	for(int i=1;i<=n;i++){
    		if(s[sa[i]]!=s[sa[i-1]]) gr++;
    		rk[sa[i]]=gr;
    	} vmax=gr;
    //	for(int i=1;i<=n;i++) printf("%d%c",rk[i],(i==n)?'
    ':' ');
    	for(int k=1;k<=n;k<<=1){
    		for(int i=1;i<=n;i++){
    			if(i+k<=n) x[i]=mp(rk[i],rk[i+k]);
    			else x[i]=mp(rk[i],0);
    		}
    		int num=0;
    		for(int i=n-k+1;i<=n;i++) seq[++num]=i;
    		for(int i=1;i<=n;i++) if(sa[i]>k) seq[++num]=sa[i]-k;
    		memset(buc,0,sizeof(buc));
    		for(int i=1;i<=n;i++) buc[x[i].fi]++;
    		for(int i=1;i<=vmax;i++) buc[i]+=buc[i-1];
    		for(int i=n;i;i--) sa[buc[x[seq[i]].fi]--]=seq[i];
    		gr=0;
    		for(int i=1;i<=n;i++){
    			if(x[sa[i]]!=x[sa[i-1]]) gr++;
    			rk[sa[i]]=gr;
    		}
    		vmax=gr;
    	}
    }
    

    我们还注意到,如果 (n) 个字符串都不相同了就可以 break 了,所以可以再加一句 if(vmax==n) break;

    发现一句的效果比之前好得多,现在最慢的一个点只跑了 200ms 了:

    void get_sa(){
    	int vmax=123;
    	for(int i=1;i<=n;i++) buc[s[i]]++;
    	for(int i=1;i<=vmax;i++) buc[i]+=buc[i-1];
    	for(int i=n;i;i--) sa[buc[s[i]]--]=i;
    	int gr=0;
    	for(int i=1;i<=n;i++){
    		if(s[sa[i]]!=s[sa[i-1]]) gr++;
    		rk[sa[i]]=gr;
    	} vmax=gr;
    //	for(int i=1;i<=n;i++) printf("%d%c",rk[i],(i==n)?'
    ':' ');
    	for(int k=1;k<=n;k<<=1){
    		for(int i=1;i<=n;i++){
    			if(i+k<=n) x[i]=mp(rk[i],rk[i+k]);
    			else x[i]=mp(rk[i],0);
    		}
    		int num=0;
    		for(int i=n-k+1;i<=n;i++) seq[++num]=i;
    		for(int i=1;i<=n;i++) if(sa[i]>k) seq[++num]=sa[i]-k;
    		memset(buc,0,sizeof(buc));
    		for(int i=1;i<=n;i++) buc[x[i].fi]++;
    		for(int i=1;i<=vmax;i++) buc[i]+=buc[i-1];
    		for(int i=n;i;i--) sa[buc[x[seq[i]].fi]--]=seq[i];
    		gr=0;
    		for(int i=1;i<=n;i++){
    			if(x[sa[i]]!=x[sa[i-1]]) gr++;
    			rk[sa[i]]=gr;
    		}
    		vmax=gr;if(vmax==n) break;
    	}
    }
    

    至此我们有了常数较小的 SA 的实现。

    那么 SA 有什么用呢?SA 一个用途就是与 LCP 结合。

    LCP,即 Longest Common Prefix,指两串的最长公共前缀。

    我们记 (LCP(i,j)) 为后缀 (sa_i)(sa_j) 的最长公共前缀,注意这边是 (sa_i)(sa_j) 的最长公共前缀,即排名为 (i) 的串与排名为 (j) 的串的最长公共前缀

    那么 LCP 有哪些性质呢?

    1. (LCP(i,j)=LCP(j,i))

    2. (LCP(i,i)=n-sa_i+1)

    3. (LCP(i,j)=min(LCP(i,k),LCP(k,j)) (i<k<j)),该定理我们称之为 LCP Lemma

    4. (LCP(i,j)=minlimits_{k=i+1}^{j}LCP(k-1,k)(i<j)),该定理我们称之为 LCP Theorem

    首先结论 (1,2) 肯定是显然的,用脚趾头都能想出来。

    结论 3 和 4 显然是同一类的,证出结论 3 就能顺带着证出结论 4.

    考虑如何证明结论 (3),设 (u) 为后缀 (sa_i)(v) 为后缀 (sa_j)(w) 为后缀 (sa_k)(x=min(LCP(i,k),LCP(k,j)))

    由于 (x=min(LCP(i,k),LCP(k,j)))(u,v,w) 的前 (x) 位应当都是相同的,故 (LCP(i,j)geq x)

    我们要证明 (LCP(i,j)=x),即 (LCP(i,j)>x) 不可能实现。

    考虑反证法,假设 (LCP(i,j)>x),则 (u) 的第 (x+1) 位等于 (v) 的第 (x+1) 位。

    由于 (i<k<j),故 (u,w,v) 在字典序上依次递增,而 (u,v,w) 的前 (x) 位都相等,故 (u_{x+1}leq w_{x+1}leq v_{x+1})

    (x=min(LCP(i,k),LCP(k,j))) 说明 (u_{x+1},v_{x+1},w_{x+1}) 不可能全相等,就说明 (u)(x+1) 位不可能等于 (v)(x+1) 位。

    (LCP(i,j)=x)

    而结论 4 其实是结论 3 的变形。

    如果我们记 (ht_i=LCP(i-1,i)),那么 (LCP(i,j)=minlimits_{k=i+1}^jht_i),这东西显然可以 RMinQ 维护。

    现在的问题就转化为怎么求 (ht_i)。暴力吗?一氧化氮。

    如果我们设 (h_i)(ht_{rk_i})即后缀 (i) 与后缀 (i) 字典序前一位的 LCP,那么有一个很重要的公式叫做 (h_igeq h_{i-1}-1)

    怎么证明这个公式呢?

    考虑设 (j=sa_{rk_i-1})(k=sa_{rk_{i-1}-1}),分两种情况讨论:

    1. 若后缀 (k) 与后缀 (i-1) 第一位不同,则 (h_{i-1}=0)(h_igeq -1) 显然成立。
    2. 若后缀 (k) 与后缀 (i-1) 第一位相同,则刨掉第一位还剩下后缀 (k+1) 和后缀 (i),它们的 (LCP)(h_{i-1}-1),而由于后缀 (k) 字典序排在后缀 (i-1) 之前,后缀 (k) 与后缀 (i-1) 第一位相同,故后面的后缀 (k+1) 的字典序小于后缀 (i) 的字典序。根据 LCP Theorem,后缀 (k+1) 与后缀 (i) 的 LCP 为 (minlimits_{t=rk_{k+1}+1}^{rk_i}ht_t=h_{i-1}-1),故 (ht_{rk_i}geq h_{i-1}-1),故 (h_igeq h_{i-1}-1)

    附图:

    有了这个式子,就可以用类似 two pointers 的东西求出 (h_i),也就顺带着求出 (ht) 了。

    void getht(){
    	int k=0;
    	for(int i=1;i<=n;i++){//calculate h[i]
    		if(rk[i]==1) continue;
    		if(k) k--;//h[i]>=h[i-1]-1
    		int j=sa[rk[i]-1];
    		while(s[i+k]==s[j+k]) k++;//violently match suffix i and suffix sa[rk[i]-1]
    		ht[rk[i]]=k;//use formula ht[rk[i]]=h[i] to calculate ht[i]
    	}
    }
    

    例题:

    1. P2408 不同子串个数

    u1s1 这题作为 SA 的第一题还是特别合适的。

    我们先求出 (ht_i)。考虑按照字典序依次加入这 (n) 个后缀。

    考虑加入排名为 (i) 的后缀的时候,会有多少个子串被重复计算。

    显然,对于该后缀的每个长度为 (k) 的前缀,如果存在某个已经加入的后缀,其排名为 (j)(LCP(i,j)geq k),那么这个长度为 (k) 的前缀就会被重复计算。

    而对于 (j<i)(LCP(i,j)=(minlimits_{k=j+1}^iht_k)geq ht_i),也就是说 (LCP(i,j),j<i) 的最大值为 (ht_i)

    故总共会有 (ht_i) 个被重复计算。

    答案即为 (dfrac{n(n+1)}{2}-sumlimits_{i=2}^nht_i)

    2. P2852 [USACO06DEC]Milk Patterns G

    原题等价于,选出原数组的 (k) 个互不相同的后缀,求这 (k) 个后缀的 LCP 的最大值。

    老套路,求出 (ht_i),假设这 (k) 个后缀中排名最靠前的后缀的排名为 (i),最靠后的为 (j)

    那么 (LCP=minlimits_{t=i+1}^jht_t)

    按照贪心的思想,我们肯定希望这 (k) 个后缀在字典序上是相邻的,即 (j=i+k-1)

    故需求出 (minlimits_{t=i+1}^{i+k-1}ht_t)

    单调队列一下就完事了

    3. P4248 [AHOI2013]差异

    首先,把 (sum) 拆开,拆成 (sumlimits_{1leq i<jleq n}len(T_i)+len(T_j)-2 imessumlimits_{1leq i<jleq n}lcp(T_i,T_j))

    左边的 (sum) 很好算,每个后缀都会产生 (n-1) 次贡献,总贡献为 (dfrac{n(n+1)(n-1)}{2})

    重点是后面的 (sum),也就是任意两个后缀的 LCP 的和。

    还是求出 (ht_i),如果我们用 (ht_i) 的视角考虑这个问题,那么任意两个后缀的 LCP 之和就是 (sumlimits_{1leq l<rleq n}minlimits_{i=l+1}^rht_i)

    到这里问题就变得异常套路了。考虑每个 (ht_i) 会对多少个区间产生贡献,设 (l_i)(i) 前面第一个 (ht_j>ht_i)(j)(r_i)(i) 前面第一个 (ht_j>ht_i)(j),那么 (ht_i) 会对所有 (l_i<lleq i,ileq r<r_i) 的区间 ([l,r]) 产生贡献,单调栈一下就可以了。

    4. P3181 [HAOI2016]找相同字符

    ……有什么意义吗?——ymx

    这次终于轮到我说这句话了。

    考虑枚举两个子串的开头位置,易知答案为 (sumlimits_{i=1}^nsumlimits_{j=1}^mlcp(s[i...n],t[j...m]))

    (u=s+'|'+t),对 (u) 执行一遍上一题的操作,算得的总贡献可以分为三部分:

    1. 两个后缀全在 (s)
    2. 两个后缀全在 (t)
    3. 一个后缀在 (s) 中,一个后缀在 (t) 中。

    显然我们只要第三部分的答案,故再对 (s,t) 分别执行一遍上一题的操作,拿总贡献减去这两次算得的贡献就行了。

    5. P2743 [USACO5.1]乐曲主题Musical Themes

    考虑这个“转调”是甚么东西。两个数组可以通过转调互相转化,当且仅当它们的差分数组相同。

    故我们只需求出原数组的差分数组中相同的没有公共部分的子串的最大长度。

    那么这个答案怎么求呢?显然答案具有可二分性,我们二分一个 (mid)

    长度 (mid) 满足条件当且仅当存在两个后缀 (i,j) 满足它们的 LCP (geq mid) 并且这两个后缀长度为 (mid) 的前缀没有公共部分。

    考虑求出高度数组 (ht_i),我们根据 LCP Theorem 将 (ht_i) 划分成若干个极长的段,使得每段中任意两个后缀的 LCP 都 (geq mid),i.e.,对于段 ([l,r])(ht_{l+1},ht_{l+2},dots,ht_r)(geq mid),显然不在同一段中的两个后缀不可能满足条件。

    然后我们对于每一段,找出是否存在某两个长度为 (mid) 的子串,它们没有公共部分。也就是这段中 (sa_i) 的最大值 (-) 最小值 (>k)。如果存在证明 (mid) 满足条件。

    6. P2463 [SDOI2008]Sandy的卡片

    按照上一题的套路,还是求出每个数组的差分数组。考虑将全部 (n) 个差分数组粘在一起,中间用某个别的数隔开。

    搞出后缀数组 (ht_i)。原题等价于在 (n) 个数组中各选择一个后缀,求这 (n) 个后缀的 LCP 的最大值。

    考虑在 (sa) 数组上扫一遍,固定住左端点,用 two pointer 求出右端点至少在什么位置才能包含全部 (n) 个数组。

    单调队列维护最小值,时间复杂度 (mathcal O(nlog n))。除去后缀数组,剩余部分复杂度 (mathcal O(n))

    7. P2178 [NOI2015] 品酒大会

    首先考虑对于某个特定的 (r) 怎么计算答案。搞出后缀数组,按照第 5 题的套路,将 (ht_i) 切成若干个极长段,使得每段中任意两个后缀的 LCP 都 (geq r)

    对于第一问,显然每段中任意两个后缀都是 (r) 相似的,任意两个不在同一段中的后缀都不是 (r) 相似的,故直接把每个段的贡献加起来就行了。一段长度为 (l) 的贡献为 (inom{l}{2})

    对于第二问,显然就是每一段中,任选两个不同后缀,所得的 (a_p imes a_q) 的最大值,注意到这题是可以小于 (0) 的,所以用个什么东西随便记录一下最大值、次大值、最小值、次小值就行了。

    考虑怎样对于每个 (r) 维护这两个值。显然当我们的 (r)(r+1) 变为 (r) 的时候,对于 (ht_{i}=r)(i)(sa_i)(sa_{i-1}) 所在的两段会合并成一段,这很显然可以用冰茶姬维护。

    时间复杂度 (mathcal O(nlog n))。除去后缀数组,剩余部分复杂度 (mathcal O(n))

    8. P4051 [JSOI2007]字符加密

    ……这题真的只是来调节一下期末考试自闭的心情的……

    断环成链,求一遍 SA 不就行了吗。。。。。。

  • 相关阅读:
    vue父组件促发子组件中的方法
    油猴脚本:油猴脚本自动点击 | 自动检测元素并点击、休眠、顺序执行、单页面也适用
    油猴脚本:使用layer.js mobx lodash jquery
    vue项目统计src目录下代码行数
    常用mobx响应新值变化函数autorun和observe
    uni app使用mobx | uni app状态管理mobx
    File and Code Templates | webstorm代码文件模板 vue typescript
    javascript立即执行函数简单介绍
    VSCode 安装GitLens插件不生效问题
    常用的浅拷贝实现方法
  • 原文地址:https://www.cnblogs.com/ET2006/p/sa.html
Copyright © 2011-2022 走看看