zoukankan      html  css  js  c++  java
  • 后缀数组(SA)及height数组

    最近感觉自己越来越蒟蒻了……后缀数组不会,费用流不会……
    看着别人切一道又一道的题,我真是很无奈啊……
    然后,我花了好长时间,终于弄懂了后缀数组。

    后缀数组是什么?

    后缀SASA数组

    给你一个字符串,让你将每个后缀排序,就是一个后缀数组
    比如,字符串为ababa,就会搞出一个这样的东西:

    a
    aba
    ababa
    ba
    baba
    SA={4,2,0,3,1};
    

    其中,每个后缀用开始的位置来表示。

    rankrank数组

    相当于逆着的SASArank[sa[i]]=irank[sa[i]]=i
    #后缀数组怎么求?

    方法一:O(N3)O(N^3)暴力

    打个选择排序,每次比较用O(N)O(N)的方法。
    当然,这样的暴力出不了奇迹。

    方法二:O(N2lgN)O(N^2lg N)快排

    仅仅是将选择排序变成快排罢了。
    ##方法三:倍增

    倍增1.0

    这就是本文的重点了!!!
    想想其它的倍增是怎么做的,再想想字符串怎么倍增。
    首先,给每个字符赋一个排名,像这样:

    'a'->1
    'b'->2
    

    现在rank={1,2,1,2,1}
    然后像下面的这幅图一样搞一搞。
    倍增
    倍增算法的精髓就在这幅图中。
    枚举一个ii,对于每个位置jj,每次将jjj+2ij+2^i合并到一起(不够补0),然后排名(可以当做是离散化)。
    这样就可以求出rankrank,然后求出SASA
    具体没什么好讲的,只要看懂这幅图,就懂了倍增算法的思想。
    这样就可以优化到O(Nlg2N)O(Nlg^2N)了。

    倍增2.0

    然而,这不是倍增的极限。
    可以用基数排序进一步地优化!!!
    什么是基数排序,基数排序怎么打,请百度一下(基数排序可以用一维的数组来打,具体看代码实践)。
    在此,我有意提醒的是,你可以将数看做一个m进制,在倍增合并两个数时,可以将其看做m进制的两位数,然后对它进行基数排序。
    这样就有O(NlgN)O(Nlg N)的时间复杂度了。
    具体见代码。

    DC3算法

    笔者暂时不会……


    后缀数组怎么打?

    后缀数组其实是比较好理解的,但是,为了追求完美,我们不应光靠自己的理解打模板。
    因为自己打的有时会非常丑陋……
    看了,理解的网上的标,综合我自己的风格,就打出了个这样的标:

    int y[2000003],ws[2000003],wv[2000003];
    void getSA(char s[],int rank[],int sa[],int n,int m)
    {
    	memset(ws,0,sizeof(int)*m);
    	memset(y,255,sizeof y);
    	memset(rank,255,sizeof rank);
    	for (int i=0;i<n;++i)
    		++ws[rank[i]=s[i]];
    	for (int i=1;i<m;++i)
    		ws[i]+=ws[i-1];
    	for (int i=n-1;i>=0;--i)
    		sa[--ws[rank[i]]]=i;
    	for (int i=1;p<n;i<<=1,m=p)
    	{
    		p=0;
    		for (int j=n-i;j<n;++j)
    			y[p++]=j;
    		for (int j=0;j<n;++j)
    			if (sa[j]>=i)
    				y[p++]=sa[j]-i;
    		for (int j=0;j<n;++j)
    			wv[j]=rank[y[j]];
    		memset(ws,0,sizeof(int)*m);
    		for (int j=0;j<n;++j)
    			++ws[wv[j]];
    		for (int j=1;j<m;++j)
    			ws[j]+=ws[j-1];
    		for (int j=n-1;j>=0;--j)
    			sa[--ws[wv[j]]]=y[j];
    		swap(rank,y);
    		p=1;
    		rank[sa[0]]=0;
    		for (int j=1;j<n;++j)
    			rank[sa[j]]=(y[sa[j-1]]==y[sa[j]] && y[sa[j-1]+i]==y[sa[j]+i]?p-1:p++);
    	}
    

    这就是网上通常的打法,当然,风格会有些不一样……
    是不是看了后,一头雾水?
    别急,慢慢解释。
    首先说一下,在这个程序中,rankrank取值是在[0,n)[0,n)范围内的,和上面那张图不一样!

    	memset(ws,0,sizeof(int)*m);
    	memset(y,255,sizeof y);
    	memset(rank,255,sizeof rank);
    	for (int i=0;i<n;++i)
    		++ws[rank[i]=s[i]];
    	for (int i=1;i<m;++i)
    		ws[i]+=ws[i-1];
    	for (int i=n-1;i>=0;--i)
    		sa[--ws[rank[i]]]=i;
    

    前面三行赋初值。
    这是处理最开始的rankranksasa,也就是还没有合并时。
    wsws数组是一个桶,用于辅助基数排序。
    注意第三行rank[i]=s[i]。我们在实践的时候一开始不需要将真正的排名弄出来,我们只需知道它们的相对大小。而sis_i作为单个字符,是可以表示它们的相对大小的。
    其它就没什么了,要理解好一维的基数排序

    for (int i=1,p=1;p<n;i<<=1,m=p)
    

    ii表示的是对于一个位置jj,在这一轮中要用jjj+ij+i合并。
    pp表示不同的字符串的个数,初值设为11是为了循环条件p&lt;np&lt;n,显然1n1geq n时就没必要做了。
    为什么循环条件是p&lt;np&lt;n呢?因为我们发现,最后的rankrank数组一定是一个范围在[0,n)[0,n)的排列
    所以pp顶多为nn,想想,p=np=n时,那么其实已经排好序了,没必要再做下去,比如,可以看看上面那张图,可以发现最后一轮是没有必要的。
    mm表示的也是不同字符串的个数,只是因为在下面pp要被用作计数器罢了。

    		p=0;
    		for (int j=n-i;j<n;++j)
    			y[p++]=j;
    		for (int j=0;j<n;++j)
    			if (sa[j]>=i)
    				y[p++]=sa[j]-i;
    

    当初我看得最久的是这一段……
    这其实是一个小优化。
    yy是一个临时的rankrank数组。
    在合并后,其实第二关键字可以通过上一次的sasa数组求出
    先看看二、三行。显然,[ni,n)[n-i,n)这段区间内,如果要和后面的合并,只能00,应该说是补1-1,因为rankrank数组在这个程序中的取值是[0,n)[0,n)
    1-1一定是最小的,所以先把它们排在前面。
    然后看倒数三行,这个就比较难理解了。
    对于位置jj,在这一轮中会对jij-i有影响,所以说,ji0j-i geq 0
    因为sasa是有序的,所以我们顺序枚举sajsa_j,将其中满足以上条件的加入yy中。

    		for (int j=0;j<n;++j)
    			wv[j]=rank[y[j]];
    		memset(ws,0,sizeof(int)*m);
    		for (int j=0;j<n;++j)
    			++ws[wv[j]];
    		for (int j=1;j<m;++j)
    			ws[j]+=ws[j-1];
    		for (int j=n-1;j>=0;--j)
    			sa[--ws[wv[j]]]=y[j];
    

    这一段的作用就是以第一关键字来进行一次基数排序,和上面的那个一样的道理。

    		swap(rank,y);
    		p=1;
    		rank[sa[0]]=0;
    		for (int j=1;j<n;++j)
    			rank[sa[j]]=(y[sa[j-1]]==y[sa[j]] && y[sa[j-1]+i]==y[sa[j]+i]?p-1:p++);
    

    pp起到计数器的作用。
    重点是最后一行y[sa[j-1]]==y[sa[j]] && y[sa[j-1]+i]==y[sa[j]+i]
    这是在比较saj1sa_{j-1}sajsa_j是否相等。如果相等,那么rankrank值应该要一样(不过注意,到最后时rankrank值一定是不同的!)
    网上的标这样比较,就不怕爆掉吗?对此,我很不理解,只是开了两倍的数组来解决这个问题。


    关于LCP

    一些概念

    LCP(i,j)LCP(i,j)表示suffix(i)suffix(i)suffix(j)suffix(j)的公共最长前缀。
    height(i)=LCP(SA[i1],SA[i])height(i)=LCP(SA[i-1],SA[i])
    很明显,若ranki&lt;rankjrank_i&lt;rank_j,则
    LCP(i,j)=mink(ranki,rankj]heightkLCP(i,j)=min_{kin left(rank_i,rank_j ight]}{height_k}

    如何求heightheight?

    首先,我们要知道一个性质:
    hi=heightrankih_i=height_{rank_i},也就是suffix(i)suffix(i)与它前一名的最长公共前缀。
    那么hihi11h_igeq h_{i-1}-1
    证明:
    suffix(k)suffix(k)表示suffix(i1)suffix(i-1)前一名的后缀,hi1h_{i-1}即是它们的最长公共前缀。
    hi11h_{i-1} leq 1时,显然等式成立。
    hi1&gt;1h_{i-1} &gt;1时,可以发现suffix(k+1)suffix(k+1)suffix(i)suffix(i)的公共后缀至少为hi11h_{i-1}-1。可以画张图理解一下。
    所以,综上,hihi11h_igeq h_{i-1}-1
    利用这个性质,我们可以在O(N)O(N)的时间内求出heightheight数组

    代码

    void getheight(char s[],int rank[],int sa[],int height[])
    {
    	for (int i=0,k=0;i<n;++i)
    		if (rank[i])
    		{
    			if (k)
    				--k;
    			int j=sa[rank[i]-1];
    			while (s[i+k]==s[j+k])
    				++k;
    			height[rank[i]]=k;
    		}
    }
    

    其实不必真正地构出个hh数组。


    其它

    建议数组从11开始,或者将rankrank及辅助数组yy初值设为1-1,因为在比较时,后面的要补00(或1-1),我就因为这样调了很久……(被罗穗骞大佬的论文坑了)

  • 相关阅读:
    kaggle CTR预估
    基于大规模语料的新词发现算法【转自matix67】
    vim E437: terminal capability "cm" required
    makefile 中的符号替换($@、$^、$<、$?)
    【转】Makefile 中:= ?= += =的区别
    python urljoin问题
    python 写文件刷新缓存
    python Popen卡死问题
    nohup 日志切割
    换行和回车野史
  • 原文地址:https://www.cnblogs.com/jz-597/p/11145287.html
Copyright © 2011-2022 走看看