求若干个字符串的最长公共子串,是后缀数组一个简单的运用,为说明简单起见,设要求的字符串个数是两个。

 

  首先我们先来看一个定义,最长公共前驱,就是常见的LCP,定义lcp(u,v)=max{i|u=v} 也就是从头开始比较u和v的对应字符持续相等的最远值,定义LCP(i,j)=lcp(Suffix(SA[i]),Suffix(SA[j])),也就是SA数组中第i个和第j个后缀的最长公共前缀,关于LCP有一个定理,LCP Theorem:对任何1≤i<j≤n LCP(i,j)=min{LCP(k-1,k) | i+1≤k≤j},设height[i]=LCP(i-1,i) 根据LCP Theorem, LCP(i,j)=min{height[k] | i+1≤k≤j},这样只要求出height数组,求最长公共前驱就是RMQ问题,运用经典的ST算法,O(nlogn)时间预处理,常数时间查询。

    那么,怎么求height数组呢,定义h[i]=height[rank[i]],也就是suffix(i)和在它前一名的后缀的最长公共前缀。h 数组有以下性质:h[i]≥h[i-1]-1,证明见老罗的论文,我就不详细说了,按照h[1],h[2],……,h[n]的顺序计算,并利用h 数组的性质,时间复杂度可以降为O(n)。h数组求出以后就可以用O(n)的时间求出height数组了,代码如下:

//这段代码并没有先求h数组然后转换成height数组,它只是按照h数组从1到n的顺序求height数组,省去了转换的时间。
int   height[max],rank[max];
char s[max];
void get_height(int *sa,int n)
{
     int i,j,k=0;
     //n表示sa数组的长度,一般编程时都是0到n-1
     for (i = 1; i < n; i++) rank[sa[i]] = i;
     //一般在求sa的时候在串最后加的一个0会排在sa的第零位,我们用不到,所以从1开始
   for (i = 0; i < n-1; height[rank[i++]] = k)
          for (k?k--:0,j = sa[rank[i]-1]; s[i+k] == s[j+k]; k++) ;
}

  求出了height数组,与多串的公共子串有什么关系呢,设串u,v,要求他们俩的最长公共子串,先把他们合并在一起,中间用一个两个串中都从未出现过的字符隔开,求新串的sa,height,

不难发现,要求的答案就是新串的分属于两个原串的两个后缀的最长公共前缀,也就是LCP。当suffix(sa[i-1]) 和suffix(sa[i])不是同一个字符串中的两个后缀时,height[i]的最大值就是我们要求的答案了。同理可得多串最长公共子串的算法。算法时间复杂度0(|s1|+|s2|+|s3|+....+|sn|),|sn|表示串sn的长度。

pku 2774就是这类题的典型,代码如下:

#include<stdio.h>
#include<string.h>
#include<math.h>
#define MAX 200010
int wx[MAX],wy[MAX],bar[MAX],nrank[MAX];
int sa[MAX],height[MAX];
int query[MAX][30];
int cmp(int *in,int a,int b,int l)
{
	return in[a] == in[b]&&in[a+l] == in[b+l];
}
int MIN(int x,int y)
{
	return x > y ? y : x;
}
/*void get_query(int *ai,int len)
{
	int i,k;
	for (i = 0; i< len; i++)
		query[i][0] = ai[i];
	for (k = 1;(int)pow(2.0,k) - 1 < len; k++)
		for (i = 0; i < len; i++)
			if (i + (int)pow(2.0,k) - 1 < len)
				query[i][k] = MIN(query[i][k-1],query[i+(int)pow(2.0,k-1)][k-1]);
			else break;
}*/
void get_sa(char *s,int *sa)
{
	int len = strlen(s),*rank = wx,*sa_y = wy,*t,i,m = 255;
	for (i = 0; i<= m; i++) bar[i] = 0;
	for (i = 0; i< len; i++) bar[rank[i] = s[i]]++;
	for (i = 0; i< m; i++) bar[i+1] += bar[i];
	for (i = len-1; i>= 0; i--)	sa[--bar[rank[i]]] = i;
	//sa数组保存的是某个名词的后缀的起始位置,对k=1时的排序用的是桶排序,先一个一个进桶,然后再倒出来
	for (int p = 1,k = 1; p < len; k *= 2,m = p)
	{
		for (p = 0,i = len-k; i < len; i++) sa_y[p++] = i;
		for (i = 0; i< len; i++) if (sa[i] >= k) sa_y[p++] = sa[i] - k;
		//要求2k-前缀的sa,需用k-前缀倍增得到,用y数组保存2k-前缀时第二关键字的排序结果,即k-前缀时某位置+k的sa排名
		//容易得到k-前缀时某位置+k有一些已经超出了s的范围,还有一些位置不能作为第二关键字起始位置,可以看出,他们的数目是相等的
		for (i = 0; i< len; i++) nrank[i] = rank[sa_y[i]];
		for (i = 0; i<= m; i++) bar[i] = 0;
		for (i = 0; i< len; i++) bar[nrank[i]] ++;
		for (i = 0; i< m; i++) bar[i+1] += bar[i];
		for (i = len-1; i>= 0; i--) sa[--bar[nrank[i]]] = sa_y[i];
		//对关于第二关键字排好序的序列再按第一关键字排序,用的还是桶排序,nrank[i]纪录的是按第二关键字排好序以后的序列的第一关键字名次序列
		for (t = sa_y,sa_y = rank,rank = t,p = 1,rank[sa[0]] = 0,i = 1; i < len; i++)
			rank[sa[i]] = cmp(sa_y,sa[i-1],sa[i],k)?p-1:p++;
		//这段代码求新的rank数组,原来的rank用sa_y来保存,比较排名相邻的两个后缀的第一和第二关键字是否完全一样,确定不同后缀的数目,用p纪录
	}
}
void get_height(char *s,int *sa,int *height)
{
	int len = strlen(s),*rank = wx,i,j,k = 0;
	for (i = 0; i < len; i++) rank[sa[i]] = i;
	for (i = 0; i < len-1; height[rank[i++]] = k)
		for (k?k--:0,j = sa[rank[i] - 1]; s[i+k] == s[j+k]; k++) ;
	//计算height用到了辅助数组h,这段代码中没有保留h,只是按照h从0到len-2的顺序计算的,这段代码非常精简,for循环用的出神入化
	//h数组的定义为h[i] = height[rank[i]],所以代码中用后者直接代替了前者,h数组有一个规律 h[i] >= h[i-1] - 1,证明论文里有
	//根据规律,就可以按照顺序计算h数组,计算方法分三种情况,老罗的论文里有
	return;
}
int main()
{
	
	char str[MAX]; int len;
	memset(str,0,sizeof(str));
	scanf ("%s",&str);
	len = strlen(str);
	str[len++] = '$';
	int start = 0,end = len;
	scanf ("%s",str+len);
	len = strlen(str);
	str[len] = '#';

	get_sa(str,sa);
	get_height(str,sa,height);
	//get_query(height,len);

    int min = 0,i,j;
	for (i = 1; i< len; i++)
	{
		if ((sa[i] >= end && sa[i-1] + height[i] < end) || (sa[i-1] >= end && sa[i] + height[i] < end))
			min = min>height[i]?min:height[i];
	}
	printf ("%d\n",min);
	return 0;
}

 

posted on 2010-07-18 15:37  looker  阅读(5053)  评论(0编辑  收藏  举报