最长公共子串问题的后缀数组解法

計算機科學 5 Comments »1,796 views

[最长公共子串]

最长公共子串(Longest Common Substring ,简称LCS)问题,是指求给定的一组字符串长度最大的共有的子串的问题。例如字符串”abcb”,”bca”,”acbc”的LCS就是”bc”。

求多串的LCS,显然穷举法是极端低效的算法。改进一些的算法是用一个串的每个后缀对其他所有串进行部分匹配,用KMP算法,时间复杂度为O(N*L^2),其中N为字符串个数,L为每个串的长度。更优秀的有广义后缀树的方法,时间可以达到 O(N*L)。本文介绍一种基于后缀数组的LCS解法,利用二分查找技术,时间复杂度可以达到O(N*L*logL)。

[最长公共子串问题的后缀数组解法]

关于后缀数组的构建方法以及Height数组的性质,本文不再具体介绍,可以参阅IOI国家集训队2004年论文《后缀数组》(许智磊)和IOI国家集训队2009年论文《后缀数组——处理字符串的有力工具》(罗穗骞)

回顾一下后缀数组,SA[i]表示排名第i的后缀的位置,Height[i]表示后缀SA[i]和SA[i-1]的最长公共前缀(Longest Common Prefix,LCP),简记为Height[i]=LCP(SA[i],SA[i-1])。连续的一段后缀SA[i..j]的最长公共前缀,就是H[i-1..j]的最小值,即LCP(SA[i..j])=Min(H[i-1..j])。

求N个串的最长公共子串,可以转化为求一些后缀的最长公共前缀的最大值,这些后缀应分属于N个串。具体方法如下:

设N个串分别为S1,S2,S3,…,SN,首先建立一个串S,把这N个串用不同的分隔符连接起来。S=S1[P1]S2[P2]S3…SN-1[PN-1]SN,P1,P2,…PN-1应为不同的N-1个不在字符集中的字符,作为分隔符(后面会解释为什么)。

接下来,求出字符串S的后缀数组和Height数组,可以用倍增算法,或DC3算法。

然后二分枚举答案A,假设N个串可以有长度为A的公共字串,并对A的可行性进行验证。如果验证A可行,A’(A’<A)也一定可行,尝试增大A,反之尝试缩小A。最终可以取得A的最大可行值,就是这N个串的最长公共子串的长度。可以证明,尝试次数是O(logL)的。

于是问题就集中到了,如何验证给定的长度A是否为可行解。方法是,找出在Height数组中找出连续的一段Height[i..j],使得i<=k<=j均满足Height[k]>=A,并且i-1<=k<=j中,SA[k]分属于原有N个串S1..SN。如果能找到这样的一段,那么A就是可行解,否则A不是可行解。

具体查找i..j时,可以先从前到后枚举i的位置,如果发现Height[i]>=A,则开始从i向后枚举j的位置,直到找到了Height[j+1]<A,判断[i..j]这个区间内SA是否分属于S1..SN。如果满足,则A为可行解,然后直接返回,否则令i=j+1继续向后枚举。S中每个字符被访问了O(1)次,S的长度为N*L+N-1,所以验证的时间复杂度为O(N*L)。

到这里,我们就可以理解为什么分隔符P1..PN-1必须是不同的N-1个不在字符集中的字符了,因为这样才能保证S的后缀的公共前缀不会跨出一个原有串的范围。

后缀数组是一种处理字符串的强大的数据结构,配合LCP函数与Height数组的性质,后缀数组更是如虎添翼。利用后缀数组,容易地求出了多个串的LCS,而且时空复杂度也相当优秀了。虽然比起后缀树的解法有所不如,但其简明的思路和容易编程的特点却在实际的应用中并不输于后缀树。

附:POI 2000 Repetitions 最长公共子串

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
/* 
 * Problem: POI2000 pow
 * Author: Guo Jiabao
 * Time: 2009.4.16 21:00
 * State: Solved
 * Memo: 多串最长公共子串 后缀数组 二分答案
*/
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
const int MAXL=10011,MAXN=6;
using namespace std;
struct SuffixArray
{
	struct RadixElement
	{
		int id,k[2];
	}RE[MAXL],RT[MAXL];
	int N,A[MAXL],SA[MAXL],Rank[MAXL],Height[MAXL],C[MAXL];
	void RadixSort()
	{
		int i,y;
		for (y=1;y>=0;y--)
		{
			memset(C,0,sizeof(C));
			for (i=1;i<=N;i++) C[RE[i].k[y]]++;
			for (i=1;i<MAXL;i++) C[i]+=C[i-1];
			for (i=N;i>=1;i--) RT[C[RE[i].k[y]]--]=RE[i];
			for (i=1;i<=N;i++) RE[i]=RT[i];
		}
		for (i=1;i<=N;i++)
		{
			Rank[ RE[i].id ]=Rank[ RE[i-1].id ];
			if (RE[i].k[0]!=RE[i-1].k[0] || RE[i].k[1]!=RE[i-1].k[1])
				Rank[ RE[i].id ]++;
		}
	}
	void CalcSA()
	{
		int i,k;
		RE[0].k[0]=-1;
		for (i=1;i<=N;i++)
			RE[i].id=i,RE[i].k[0]=A[i],RE[i].k[1]=0;
		RadixSort();
		for (k=1;k+1<=N;k*=2)
		{
			for (i=1;i<=N;i++)
				RE[i].id=i,RE[i].k[0]=Rank[i],RE[i].k[1]=i+k<=N?Rank[i+k]:0;
			RadixSort();
		}
		for (i=1;i<=N;i++)
			SA[ Rank[i] ]=i;
	}
	void CalcHeight()
	{
		int i,k,h=0;
		for (i=1;i<=N;i++)
		{
			if (Rank[i]==1)
				h=0;
			else
			{
				k=SA[Rank[i]-1];
				if (--h<0) h=0;
				for (;A[i+h]==A[k+h];h++);
			}
			Height[Rank[i]]=h;
		}
	}
}SA;
int N,Ans,Bel[MAXL];
char S[MAXL];
void init()
{
	int i;
	freopen("pow.in","r",stdin);
	freopen("pow.out","w",stdout);
	scanf("%d",&N);
	SA.N=0;
	for (i=1;i<=N;i++)
	{
		scanf("%s",S);
		for (char *p=S;*p;p++)
		{
			SA.A[++SA.N]=*p-'a'+1;
			Bel[SA.N]=i;
		}
		if (i<N)
			SA.A[++SA.N]=30+i;
	}
}
bool check(int A)
{
	int i,j,k;
	bool ba[MAXN];
	for (i=1;i<=SA.N;i++)
	{
		if (SA.Height[i]>=A)
		{
			for (j=i;SA.Height[j]>=A && j<=SA.N;j++);
			j--;
			memset(ba,0,sizeof(ba));
			for (k=i-1;k<=j;k++)
				ba[Bel[SA.SA[k]]]=true;
			for (k=1;ba[k] && k<=N;k++);
			if (k==N+1)
				return true;
			i=j;
		}
	}
	return false;
}
void solve()
{
	int a,b,m;
	SA.CalcSA();
	SA.CalcHeight();
	a=0;b=SA.N;
	while (a+1<b)
	{
		m=(a+b)/2;
		if (check(m))
			a=m;
		else
			b=m-1;
	}
	if (check(b))
		a=b;
	Ans=a;
}
int main()
{
	init();
	solve();
	printf("%dn",Ans);
	return 0;
}
标签:, , ,

省队培训日志6

競賽歷程 3 Comments »334 views

今天上午做了一个BT的交互式题,我用一个很烂的贪心策略,居然过了6组数据,而崔万云大牛的筛选法正解只是30分。还有一个就是很经典的问题,求N的字符串的最长公共子串(LSC问题)。有人用朴素的算法,直接调用C标准库中的strstr函数,居然拿了满分。我写了后缀字符串树,但是建树写的不太对,所以超时了一部分,得了70分。

(据说gcc的strstr使用BM写的,C语言果然有无比的优越性啊)

马浩和秦凌聪带过来了笔记本电脑,因为电脑太有问题了,不是不能联网,就是不能用USB,还不能用光驱,甚至不时还短路。昨天我把主板上的电池换了,开机好开多了,不用在反复开关20多次了。

郑州的夏天真热。昨天中午去的早了,机房没开门,我在楼梯上站了2个小时,汗水湿透了全身。后来张华翼也来陪我站了,有人一起说话也就不感觉太热了。张华翼只比我大1个月零6天,却已经该上高三。原来他初中参加的是新乡一中的“少年班”,5年学完初中和高中的全部知识。只是他后来走了,为了OI,去了河南师范大学附中。听他说他以前“少年班”的同学有很多要复读高三了,看来是“少年班”不行啊。张华翼两年前作出的决策是十分正确的。

在OIBH上听说明天鹰牛要来铁一中看看,无限期待中。

谁会用Cena评测交互式和提交答案的题,教教我啊,今天评测交互式题很痛苦。

标签:, , , , , , ,
25 queries. 0.992 seconds. Designed by NattyWP .
Images by desEXign.