在字符串处理当中,后缀树和后缀数组都是非常有力的工具,其中后缀树大家了解得比较多,关于后缀数组则很少见于国内的资料。其实后缀数组是后缀树的一个非常精巧的替代品,它比后缀树容易编程实现,能够实现后缀树的很多功能而时间复杂度也不太逊色,并且,它比后缀树所占用的空间小很多。可以说,在信息学竞赛中后缀数组比后缀树要更为实用。因此在本文中笔者想介绍一下后缀数组的基本概念、构造方法,以及配合后缀数组的最长公共前缀数组的构造方法,最后结合一些例子谈谈后缀数组的应用。
基本定义
子串
字符串 S 的子串 r[i..j] , i ≤ j ,表示 r 串中从 i 到 j 这一段,就是顺次排列 r[i],r[i+1],...,r[j] 形成的字符串。
后缀
后缀是指从某个位置 i 开始到整个串末尾结束的一个特殊子串。字符串 r 的从 第 i 个字 符 开 始 的 后 缀 表 示 为 Suffix(i) ,也 就 是Suffix(i)=r[i..len(r)] 。
大小比较
关于字符串的大小比较,是指通常所说的 “ 字典顺序 ” 比较, 也就是对于两个字符串 u 、 v ,令 i 从 1 开始顺次比较 u[i] 和 v[i] ,如果u[i]=v[i] 则令 i 加 1 ,否则若 u[i]<v[i] 则认为 u<v , u[i]>v[i] 则认为 u>v(也就是 v<u ),比较结束。如果 i>len(u) 或者 i>len(v) 仍比较不出结果,那么 若 len(u)<len(v) 则 认 为u<v , 若 len(u)=len(v) 则 认 为 u=v , 若len(u)>len(v) 则 u>v 。从字符串的大小比较的定义来看, S 的两个开头位置不同的后缀 u 和 v 进行比较的结果不可能是相等,因为 u=v 的必要条件 len(u)=len(v) 在这里不可
能满足。
后缀数组
后缀数组 SA 是一个一维数组,它保存 1..n 的某个排列 SA[1 ] ,
SA[2] , …… , SA[n] ,并且保证 Suffix(SA[i]) < Suffix(SA[i+1]) , 1 ≤ i<n 。
也就是将 S 的 n 个后缀从小到大进行排序之后把排好序的后缀的开头位置顺
次放入 SA 中。
名次数组
名次数组 Rank[i] 保存的是 Suffix(i) 在所有后缀中从小到大 排
列的 “ 名次 ” 。
简单的说,后缀数组是 “ 排第几的是谁? ” ,名次数组是 “ 你排第几? ” 。 容
易看出,后缀数组和名次数组为互逆运算。
设字符串的长度为 n 。为了方便比较大小,可以在字符串后面添加一个字 符 ,
这个字符没有在前面的字符中出现过,而且比前面的字符都要小。在求出名次 数
组后,可以仅用 O(1) 的时间比较任意两个后缀的大小。在求出后缀数组或名次
数组中的其中一个以后,便可以用 O(n) 的时间求出另外一个。任意两个后缀如
果直接比较大小,最多需要比较字符 n 次,也就是说最迟在比较第 n 个字符时 一
定能分出 “ 胜负 ” 。
倍增算法
主要思路
如何构造后缀数组呢?最直接最简单的方法当然是把S的后缀都看作一些普通的字符串,按照一般字符串排序的方法对它们从小到大进行排序。
不难看出,这种做法是很笨拙的,因为它没有利用到各个后缀之间的有机联系,所以它的效率不可能很高。即使采用字符串排序中比较高效的Multi-key Quick Sort,最坏情况的时间复杂度仍然是O(n2)的,不能满足我们的需要。
下面介绍倍增算法(Doubling Algorithm),它正是充分利用了各个后缀之间的联系,将构造后缀数组的最坏时间复杂度成功降至O(nlogn)。
对一个字符串u,我们定义u的k-前缀
定义k-前缀比较关系<k、=k和≤k:
设两个字符串u和v,
u<kv 当且仅当 uk<vk
u=kv 当且仅当 uk=vk
u≤kv 当且仅当 uk≤vk
直观地看这些加了一个下标k的比较符号的意义就是对两个字符串的前k个字符进行字典序比较,特别的一点就是在作大于和小于的比较时如果某个字符串的长度不到k也没有关系,只要能够在k个字符比较结束之前得到第一个字符串大于或者小于第二个字符串就可以了。
根据前缀比较符的性质我们可以得到以下的非常重要的性质:
性质1.1 对k≥n,Suffix(i)<kSuffix(j) 等价于 Suffix(i)<Suffix(j)。
性质1.2 Suffix(i)=2kSuffix(j)等价于
Suffix(i)=kSuffix(j) 且 Suffix(i+k)=kSuffix(j+k)。
性质1.3 Suffix(i)<2kSuffix(j) 等价于
Suffix(i)<kS(j) 或 (Suffix(i)=kSuffix(j) 且 Suffix(i+k)<kSuffix(j+k))。
这里有一个问题,当i+k>n或者j+k>n的时候Suffix(i+k)或Suffix(j+k)是无明确定义的表达式,但实际上不需要考虑这个问题,因为此时Suffix(i)或者Suffix(j)的长度不超过k,也就是说它们的k-前缀以'$'结尾,于是k-前缀比较的结果不可能相等,也就是说前k个字符已经能够比出大小,后面的表达式自然可以忽略,这也就看出我们规定S以'$'结尾的特殊用处了。
定义k-后缀数组SAk保存1..n的某个排列SAk[1],SAk[2],…SAk[n]使得Suffix(SAk) ≤kSuffix(SAk[i+1]),1≤i<n。也就是说对所有的后缀在k-前缀比较关系下从小到大排序,并且把排序后的后缀的开头位置顺次放入数组SAk中。
定义k-名次数组Rankk,Rankk代表Suffix(i)在k-前缀关系下从小到大的“名次”,也就是1加上满足Suffix(j)<kSuffix(i)的j的个数。通过SAk很容易在O(n)的时间内求出Rankk。
假设我们已经求出了SAk和Rankk,那么我们可以很方便地求出SA2k和Rank2k,因为根据性质1.2和1.3,2k-前缀比较关系可以由常数个k-前缀比较关系组合起来等价地表达,而Rankk数组实际上给出了在常数时间内进行<k和=k比较的方法,即:
Suffix(i)<kSuffix(j) 当且仅当 Rankk<Rankk[j]
Suffix(i)=kSuffix(j) 当且仅当 Rankk=Rankk[j]
因此,比较Suffix(i)和Suffix(j)在k-前缀比较关系下的大小可以在常数时间内完成,于是对所有的后缀在≤k关系下进行排序也就和一般的排序没有什么区别了,它实际上就相当于每个Suffix(i)有一个主关键字Rankk和一个次关键字Rankk[i+k]。如果采用快速排序之类O(nlogn)的排序,那么从SAk和Rankk构造出SA2k的复杂度就是O(nlogn)。更聪明的方法是采用基数排序,复杂度为O(n)。
求出SA2k之后就可以在O(n)的时间内根据SA2k构造出Rank2k。因此,从SAk和Rankk推出SA2k和Rank2k可以在O(n)时间内完成。
下面只有一个问题需要解决:如何构造出SA1和Rank1。这个问题非常简单:因为<1,=1和≤1这些运算符实际上就是对字符串的第一个字符进行比较,所以只要把每个后缀按照它的第一个字符进行排序就可以求出SA1,不妨就采用快速排序,复杂度为O(nlogn)。
于是,可以在O(nlogn)的时间内求出SA1和Rank1。
求出了SA1和Rank1,我们可以在O(n)的时间内求出SA2和Rank2,同样,我们可以再用O(n)的时间求出SA4和Rank4,这样,我们依次求出:
SA2和Rank2,SA4和Rank4,SA8和Rank8,……直到SAm和Rankm,其中m=2k且m≥n。而根据性质1.1,SAm和SA是等价的。这样一共需要进行logn次O(n)的过程,因此
可以在O(nlogn)的时间内计算出后缀数组SA和名次数组Rank。
具体实现
int wa[maxn],wb[maxn],wv[maxn],ws[maxn];
int cmp(int *r,int a,int b,int l)
{return r[a]==r[b]&&r[a+l]==r[b+l];}
void da(int *r,int *sa,int n,int m)
{
int i,j,p,*x=wa,*y=wb,*t;
for(i=0;i<m;i++) ws[i]=0;
for(i=0;i<n;i++) ws[x[i]=r[i]]++;
for(i=1;i<m;i++) ws[i]+=ws[i-1];
for(i=n-1;i>=0;i--) sa[--ws[x[i]]]=i;
for(j=1,p=1;p<n;j*=2,m=p)
{
for(p=0,i=n-j;i<n;i++) y[p++]=i;
for(i=0;i<n;i++) if(sa[i]>=j) y[p++]=sa[i]-j;
for(i=0;i<n;i++) wv[i]=x[y[i]];
for(i=0;i<m;i++) ws[i]=0;
for(i=0;i<n;i++) ws[wv[i]]++;
for(i=1;i<m;i++) ws[i]+=ws[i-1];
for(i=n-1;i>=0;i--) sa[--ws[wv[i]]]=y[i];
for(t=x,x=y,y=t,p=1,x[sa[0]]=0,i=1;i<n;i++)
x[sa[i]]=cmp(y,sa[i-1],sa[i],j)?p-1:p++;
}
return;
}
待排序的字符串放在 r 数组中 , 从 r[0] 到 r[n-1] , 长度为 n , 且最大值小于 m 。为了函数操作的方便,约定除 r[n-1] 外所有的 r[i] 都大于 0, r[n-1]=0 。函数结束后,结果放在 sa 数组中,从 sa[0] 到 sa[n-1] 。
函数的第一步,要对长度为 1 的字符串进行排序。一般来说,在字符串的 题目中, r 的最大值不会很大,所以这里使用了基数排序。如果 r 的最大值很大,那么把这段代码改成快速排序。代码:
for(i=0;i<m;i++) ws[i]=0;
for(i=0;i<n;i++) ws[x[i]=r[i]]++;
for(i=1;i<m;i++) ws[i]+=ws[i-1];
for(i=n-1;i>=0;i--) sa[--ws[x[i]]]=i;
这里 x 数组保存的值相当于是 rank 值。下面的操作只是用 x 数组来比较字符的大小,所以没有必要求出当前真实的 rank 值。
接下来进行若干次基数排序,在实现的时候,这里有一个小优化。基数排 序要分两次,第一次是对第二关键字排序,第二次是对第一关键字排序。对第二 关键字排序的结果实际上可以利用上一次求得的 sa 直接算出,没有必要再算一次 。代码:
for(p=0,i=n-j;i<n;i++) y[p++]=i;
for(i=0;i<n;i++) if(sa[i]>=j) y[p++]=sa[i]-j;
其中变量 j 是当前字符串的长度,数组 y 保存的是对第二关键字排序的结 果 。然后要对第一关键字进行排序,代码:
for(i=0;i<n;i++) wv[i]=x[y[i]];
for(i=0;i<m;i++) ws[i]=0;
for(i=0;i<n;i++) ws[wv[i]]++;
for(i=1;i<m;i++) ws[i]+=ws[i-1];
for(i=n-1;i>=0;i--) sa[--ws[wv[i]]]=y[i];
这样便求出了新的 sa 值。在求出 sa 后,下一步是计算 rank 值。这里要注意的是,可能有多个字符串的 rank 值是相同的,所以必须比较两个字符串是否完全相同, y 数组的值已经没有必要保存,为了节省空间,这里用 y 数组保存 r ank值。这里又有一个小优化,将 x 和 y 定义为指针类型,复制整个数组的操作可 以用交换指针的值代替,不必将数组中值一个一个的复制。代码:
for(t=x,x=y,y=t,p=1,x[sa[0]]=0,i=1;i<n;i++)
x[sa[i]]=cmp(y,sa[i-1],sa[i],j)?p-1:p++;
其中 cmp 函数的代码是:
int cmp(int *r,int a,int b,int l)
{return r[a]==r[b]&&r[a+l]==r[b+l];}
这里可以看到规定 r[n-1]=0 的好处,如果 r[a]=r[b] ,说明以 r[a] 或 r[b ]开头的长度为 l 的字符串肯定不包括字符 r[n-1] ,所以调用变量 r[a+l] 和 r[b +l]不会导致数组下标越界,这样就不需要做特殊判断。执行完上面的代码后, ran k值保存在 x 数组中,而变量 p 的结果实际上就是不同的字符串的个数。这里可 以加一个小优化,如果 p 等于 n ,那么函数可以结束。因为在当前长度的字符串 中 ,已经没有相同的字符串,接下来的排序不会改变 rank 值。例如图 1 中的第四次排序,实际上是没有必要的。对上面的两段代码,循 环的初始赋值和终止条件 可以这样写:
for(j=1,p=1;p<n;j*=2,m=p) { ………… }
在第一次排序以后, rank 数组中的最大值小于 p ,所以让 m=p 。
整个倍增算法基本写好,代码大约 25 行。
算法分析
倍增算法的时间复杂度比较容易分析。每次基数排序的时间复杂度为 O(n) ,
排序的次数决定于最长公共子串的长度,最坏情况下,排序次数为 logn 次,所
以总的时间复杂度为 O(nlogn) 。
DC3算法
主要思路
DC3 算法分 3 步:
(1)、先将后缀分成两部分,然后对第一部分的后缀排序。
将后缀分成两部分,第一部分是后缀 k ( k 模 3 不等于 0 ),第二部分是后 缀k ( k 模 3 等于 0 )。先对所有起始位置模 3 不等于 0 的后缀进行排序,即对suffix(1), suffix(2), suffix(4), suffix(5), suffix(7) …… 进行排序。做法是将 suffix(1) 和 suffix(2) 连接,如果这两个后缀的长度不是 3 的倍数,那先各自在末尾添 0 使得长度都变成 3 的倍数。然后每 3 个字符为一组,进行基 数排序,将每组字符 “ 合并 ” 成一个新的字符。然后用递归的方法求这个新的字 符串的后缀数组。在得到新的字符串的 sa 后,便可以计算原字符串所有起始位置模 3 不等于 0 的后缀的 sa 。要注意的是,原字符串必须以一个最小的且前面没有出现过的字符结尾,这样才能保证结果正确。
(2)、利用( 1 )的结果,对第二部分的后缀排序。
剩下的后缀是起始位置模 3 等于 0 的后缀,而这些后缀都可以看成是一个 字符加上一个在( 1 )中已经求出 rank 的后缀,所以只要一次基数排序便可以求 出剩下的后缀的 sa 。
(3)、将(1)和(2)的结果合并,即完成对所有后缀排序。
这个合并操作跟合并排序中的合并操作一样。每次需要比较两个后缀的大
小。分两种情况考虑,第一种情况是suffix(3*i)和suffix(3*j+1)的比较,可
以把suffix(3*i)和suffix(3*j+1)表示成:
suffix(3*i) = r[3*i] +suffix(3*i+1)
suffix(3*j+1)= r[3*j+1]+suffix(3*j+2)
其中suffix(3*i+1)和suffix(3*j+2)的比较可以利用(2)的结果快速得到。
第二种情况是suffix(3*i)和suffix(3*j+2)的比较,可以把suffix(3*i)和
suffix(3*j+2)表示成:
suffix(3*i) =r[3*i] +r[3*i+1]+ suffix(3*i+2)
suffix(3*j+2)=r[3*j+2]+r[3*j+3]+ suffix(3*(j+1)+1)
同样的道理,suffix(3*i+2)和suffix(3*(j+1)+1) 的比较可以利用(2)
的结果快速得到。所以每次的比较都可以高效的完成,这也是之前要每3个字符
合并,而不是每2个字符合并的原因。
具体实现
#defineF(x)((x)/3+((x)%3==1?0:tb))
#defineG(x)((x)<tb?(x)*3+1:((x)-tb)*3+2)
intwa[maxn],wb[maxn],wv[maxn],ws[maxn];
intc0(int*r,int a,intb)
{returnr[a]==r[b]&&r[a+1]==r[b+1]&&r[a+2]==r[b+2];}
intc12(intk,int *r,inta,intb)
{if(k==2)return r[a]<r[b]||r[a]==r[b]&&c12(1,r,a+1,b+1);
elsereturnr[a]<r[b]||r[a]==r[b]&&wv[a+1]<wv[b+1];}
voidsort(int *r,int*a,int*b,intn,intm)
{
inti;
for(i=0;i<n;i++)wv[i]=r[a[i]];
for(i=0;i<m;i++)ws[i]=0;
for(i=0;i<n;i++)ws[wv[i]]++;
for(i=1;i<m;i++)ws[i]+=ws[i-1];
for(i=n-1;i>=0;i--)b[--ws[wv[i]]]=a[i];
return;
}
voiddc3(int*r,int *sa,intn,intm)
{
inti,j,*rn=r+n,*san=sa+n,ta=0,tb=(n+1)/3,tbc=0,p;
r[n]=r[n+1]=0;
for(i=0;i<n;i++)if(i%3!=0)wa[tbc++]=i;
sort(r+2,wa,wb,tbc,m);
sort(r+1,wb,wa,tbc,m);
sort(r,wa,wb,tbc,m);
for(p=1,rn[F(wb[0])]=0,i=1;i<tbc;i++)
rn[F(wb[i])]=c0(r,wb[i-1],wb[i])?p-1:p++;
if(p<tbc) dc3(rn,san,tbc,p);
elsefor(i=0;i<tbc;i++)san[rn[i]]=i;
for(i=0;i<tbc;i++)if(san[i]<tb)wb[ta++]=san[i]*3;
if(n%3==1) wb[ta++]=n-1;
sort(r,wb,wa,ta,m);
for(i=0;i<tbc;i++)wv[wb[i]=G(san[i])]=i;
for(i=0,j=0,p=0;i<ta&&j<tbc;p++)
sa[p]=c12(wb[j]%3,r,wa[i],wb[j])?wa[i++]:wb[j++];
for(;i<ta;p++)sa[p]=wa[i++];
for(;j<tbc;p++)sa[p]=wb[j++];
return;
}
各个参数的作用和前面的倍增算法一样,不同的地方是r数组和sa数组的
大小都要是3*n,这为了方便下面的递归处理,不用每次都申请新的内存空间。
函数中用到的变量:
inti,j,*rn=r+n,*san=sa+n,ta=0,tb=(n+1)/3,tbc=0,p;
rn数组保存的是(1)中要递归处理的新字符串,san数组是新字符串的sa。
变量ta表示起始位置模3为0的后缀个数,变量tb表示起始位置模3为1的后
缀个数,已经直接算出。变量tbc表示起始位置模3为1或2的后缀个数。先按
(1)中所说的用基数排序把3个字符“合并”成一个新的字符。为了方便操作,
先将r[n]和r[n+1]赋值为0。
代码:
r[n]=r[n+1]=0;
for(i=0;i<n;i++)if(i%3!=0)wa[tbc++]=i;
sort(r+2,wa,wb,tbc,m);
sort(r+1,wb,wa,tbc,m);
sort(r,wa,wb,tbc,m);
其中sort函数的作用是进行基数排序。代码:
voidsort(int *r,int*a,int*b,intn,intm)
{
inti;
for(i=0;i<n;i++)wv[i]=r[a[i]];
for(i=0;i<m;i++)ws[i]=0;
for(i=0;i<n;i++)ws[wv[i]]++;
for(i=1;i<m;i++)ws[i]+=ws[i-1];
for(i=n-1;i>=0;i--)b[--ws[wv[i]]]=a[i];
return;
}
基数排序结束后,新的字符的排名保存在wb数组中。
跟倍增算法一样,在基数排序以后,求新的字符串时要判断两个字符组是否
完全相同。代码:
for(p=1,rn[F(wb[0])]=0,i=1;i<tbc;i++)
rn[F(wb[i])]=c0(r,wb[i-1],wb[i])?p-1:p++;
其中F(x)是计算出原字符串的suffix(x)在新的字符串中的起始位置,c0
函数是比较是否完全相同,在开头加一段代码:
#defineF(x)((x)/3+((x)%3==1?0:tb))
inlineintc0(int *r,inta,intb)
{returnr[a]==r[b]&&r[a+1]==r[b+1]&&r[a+2]==r[b+2];}
接下来是递归处理新的字符串,这里和倍增算法一样,可以加一个小优化,
如果p等于tbc,那么说明在新的字符串中没有相同的字符,这样可以直接求出
san数组,并不用递归处理。代码:
if(p<tbc)dc3(rn,san,tbc,p);
elsefor(i=0;i<tbc;i++)san[rn[i]]=i;
然后是第(2)步,将所有起始位置模3等于0的后缀进行排序。其中对第
二关键字的排序结果可以由新字符串的sa直接计算得到,没有必要再排一次。
代码:
for(i=0;i<tbc;i++)if(san[i]<tb)wb[ta++]=san[i]*3;
if(n%3==1)wb[ta++]=n-1;
sort(r,wb,wa,ta,m);
for(i=0;i<tbc;i++)wv[wb[i]=G(san[i])]=i;
要注意的是,如果n%3==1,要特殊处理suffix(n-1),因为在san数组里并
没有suffix(n)。G(x)是计算新字符串的suffix(x)在原字符串中的位置,和F(x)
为互逆运算。在开头加一段:
#defineG(x)((x)<tb?(x)*3+1:((x)-tb)*3+2)。
最后是第(3)步,合并所有后缀的排序结果,保存在sa数组中。代码:
for(i=0,j=0,p=0;i<ta&&j<tbc;p++)
sa[p]=c12(wb[j]%3,r,wa[i],wb[j])?wa[i++]:wb[j++];
for(;i<ta;p++)sa[p]=wa[i++];
for(;j<tbc;p++)sa[p]=wb[j++];
其中c12函数是按(3)中所说的比较后缀大小的函数,k=1是第一种情况,
k=2是第二种情况。代码:
intc12(intk,int *r,inta,intb)
{if(k==2)return r[a]<r[b]||r[a]==r[b]&&c12(1,r,a+1,b+1);
elsereturnr[a]<r[b]||r[a]==r[b]&&wv[a+1]<wv[b+1];}
整个DC3算法基本写好,代码大约40行。
算法分析
假设这个算法的时间复杂度为f(n)。容易看出第(1)步排序的时间为O(n)
(一般来说,m比较小,这里忽略不计),新的字符串的长度不超过2n/3,求新
字符串的sa的时间为f(2n/3),第(2)和第(3)步的时间都是O(n)。所以
f(n)=O(n)+f(2n/3)
f(n)≤c×n+f(2n/3)
f(n)≤c×n+c×(2n/3)+c×(4n/9)+c×(8n/27)+……≤3c×n
所以 f(n)=O(n)
由此看出,DC3算法是一个优秀的线性算法。
倍增算法与DC3 算法的比较
从时间复杂度、空间复杂度、编程复杂度和实际效率等方面对倍增算法与DC3算法进行比较。
时间复杂度:
倍增算法的时间复杂度为O(nlogn),DC3算法的时间复杂度为O(n)。从常
数上看,DC3算法的常数要比倍增算法大。
空间复杂度:
倍增算法和DC3算法的空间复杂度都是O(n)。按前面所讲的实现方法,倍
增算法所需数组总大小为6n,DC3算法所需数组总大小为10n。
编程复杂度:
倍增算法的源程序长度为25行,DC3算法的源程序长度为40行。
实际效率:
测试环境:NOI-linux Pentium(R)4CPU2.80GHz
N 倍增算法 DC3算法
200000 192 140
300000 367 244
500000 750 499
1000000 1693 1248
(不包括读入和输出的时间,单位:ms)
从表中可以看出,DC3算法在实际效率上还是有一定优势的。倍增算法容易
实现,DC3算法效率比较高,但是实现起来比倍增算法复杂一些。对于不同的题
目,应当根据数据规模的大小决定使用哪个算法。
最长公共前缀
现在一个字符串S的后缀数组SA可以在O(nlogn)的时间内计算出来。利用SA我们已经可以做很多事情,比如在O(mlogn)的时间内进行模式匹配,其中m,n分别为模式串和待匹配串的长度。但是要想更充分地发挥后缀数组的威力,我们还需要计算一个辅助的工具——最长公共前缀(Longest Common Prefix)。
对两个字符串u,v定义函数lcp(u,v)=max{i|u=iv},也就是从头开始顺次比较u和v的对应字符,对应字符持续相等的最大位置,称为这两个字符串的最长公共前缀。
对正整数i,j定义LCP(i,j)=lcp(Suffix(SA),Suffix(SA[j]),其中i,j均为1至n的整数。LCP(i,j)也就是后缀数组中第i个和第j个后缀的最长公共前缀的长度。
关于LCP有两个显而易见的性质:
性质2.1 LCP(i,j)=LCP(j,i)
性质2.2 LCP(i,i)=len(Suffix(SA))=n-SA+1
这两个性质的用处在于,我们计算LCP(i,j)时只需要考虑i<j的情况,因为i>j时可交换i,j,i=j时可以直接输出结果n-SA+1。
直接根据定义,用顺次比较对应字符的方法来计算LCP(i,j)显然是很低效的,时间复杂度为O(n),所以我们必须进行适当的预处理以降低每次计算LCP的复杂度。
经过仔细分析,我们发现LCP函数有一个非常好的性质:
设i<j,则LCP(i,j)=min{LCP(k-1,k)|i+1≤k≤j} (LCP Theorem)
要证明LCP Theorem,首先证明LCP Lemma:
对任意1≤i<j<k≤n,LCP(i,k)=min{LCP(i,j),LCP(j,k)}
证明:设p=min{LCP(i,j),LCP(j,k)},则有LCP(i,j)≥p,LCP(j,k)≥p。
设Suffix(SA)=u,Suffix(SA[j])=v,Suffix(SA[k])=w。
由u=LCP(i,j)v得u=pv;同理v=pw。
于是Suffix(SA)=pSuffix(SA[k]),即LCP(i,k)≥p。 (1)
又设LCP(i,k)=q>p,则
u[1]=w[1],u[2]=w[2],...u[q]=w[q]。
而min{LCP(i,j),LCP(j,k)}=p说明u[p+1]≠v[p+1]或v[p+1]≠w[q+1],
设u[p+1]=x,v[p+1]=y,w[p+1]=z,显然有x≤y≤z,又由p<q得p+1≤q,应该有x=z,也就是x=y=z,这与u[p+1]≠v[p+1]或v[p+1]≠w[q+1]矛盾。
于是,q>p不成立,即LCP(i,k)≤p。 (2)
综合(1),(2)知 LCP(i,k)=p=min{LCP(i,j),LCP(j,k)},LCP Lemma得证。
于是LCP Theorem可以证明如下:
当j-i=1和j-i=2时,显然成立。
设j-i=m时LCP Theorem成立,当j-i=m+1时,
由LCP Lemma知LCP(i,j)=min{LCP(i,i+1),LCP(i+1,j)},
因j-(i+1)≤m,LCP(i+1,j)=min{LCP(k-1,k)|i+2≤k≤j},故当j-i=m+1时,仍有
LCP(i,j)=min{LCP(i,i+1),min{LCP(k-1,k)|i+2≤k≤j}}=min{LCP(k-1,k}|i+1≤k≤j)
根据数学归纳法,LCP Theorem成立。
根据LCP Theorem得出必然的一个推论:
LCP Corollary 对i≤j<k,LCP(j,k)≥LCP(i,k)。
定义一维数组height,令height=LCP(i-1,i),1<i≤n,并设height[1]=0。
由LCP Theorem,LCP(i,j)=min{height[k]|i+1≤k≤j},也就是说,计算LCP(i,j)等同于询问一维数组height中下标在i+1到j范围内的所有元素的最小值。如果height数组是固定的,这就是非常经典的RMQ(Range Minimum Query)问题。
RMQ问题可以用线段树或静态排序树在O(nlogn)时间内进行预处理,之后每次询问花费时间O(logn),更好的方法是RMQ标准算法,可以在O(n)时间内进行预处理,每次询问可以在常数时间内完成。
对于一个固定的字符串S,其height数组显然是固定的,只要我们能高效地求出height数组,那么运用RMQ方法进行预处理之后,每次计算LCP(i,j)的时间复杂度就是常数级了。于是只有一个问题——如何尽量高效地算出height数组。
根据计算后缀数组的经验,我们不应该把n个后缀看作互不相关的普通字符串,而应该尽量利用它们之间的联系,下面证明一个非常有用的性质:
为了描述方便,设h=height[Rank],即height=h[SA]。h数组满足一个性质:
性质3 对于i>1且Rank>1,一定有h≥h[i-1]-1。
为了证明性质3,我们有必要明确两个事实:
设i<n,j<n,Suffix(i)和Suffix(j)满足lcp(Suffix(i),Suffix(j)>1,则成立以下两点:
Fact 1 Suffix(i)<Suffix(j) 等价于 Suffix(i+1)<Suffix(j+1)。
Fact 2 一定有lcp(Suffix(i+1),Suffix(j+1))=lcp(Suffix(i),Suffix(j))-1。
看起来很神奇,但其实很自然:lcp(Suffix(i),Suffix(j))>1说明Suffix(i)和Suffix(j)的第一个字符是相同的,设它为α,则Suffix(i)相当于α后连接Suffix(i+1),Suffix(j)相当于α后连接Suffix(j+1)。比较Suffix(i)和Suffix(j)时,第一个字符α是一定相等的,于是后面就等价于比较Suffix(i)和Suffix(j),因此Fact 1成立。Fact 2可类似证明。
于是可以证明性质3:
当h[i-1]≤1时,结论显然成立,因h≥0≥h[i-1]-1。
当h[i-1]>1时,也即height[Rank[i-1]]>1,可见Rank[i-1]>1,因height[1]=0。
令j=i-1,k=SA[Rank[j]-1]。显然有Suffix(k)<Suffix(j)。
根据h[i-1]=lcp(Suffix(k),Suffix(j))>1和Suffix(k)<Suffix(j):
由Fact 2知lcp(Suffix(k+1),Suffix(i))=h[i-1]-1。
由Fact 1知Rank[k+1]<Rank,也就是Rank[k+1]≤Rank-1。
于是根据LCP Corollary,有
LCP(Rank-1,Rank)≥LCP(Rank[k+1],Rank)
=lcp(Suffix(k+1),Suffix(i))
=h[i-1]-1
由于h=height[Rank]=LCP(Rank-1,Rank),最终得到 h≥h[i-1]-1。
根据性质3,可以令i从1循环到n按照如下方法依次算出h:
若Rank=1,则h=0。字符比较次数为0。
若i=1或者h[i-1]≤1,则直接将Suffix(i)和Suffix(Rank-1)从第一个字符开始依次比较直到有字符不相同,由此计算出h。字符比较次数为h+1,不超过h-h[i-1]+2。
否则,说明i>1,Rank>1,h[i-1]>1,根据性质3,Suffix(i)和Suffix(Rank-1)至少有前h[i-1]-1个字符是相同的,于是字符比较可以从h[i-1]开始,直到某个字符不相同,由此计算出h。字符比较次数为h-h[i-1]+2。
设SA[1]=p,那么不难看出总的字符比较次数不超过
也就是说,整个算法的复杂度为O(n)。
求出了h数组,根据关系式height=h[SA]可以在O(n)时间内求出height数组,于是
可以在O(n)时间内求出height数组。
结合RMQ方法,在O(n)时间和空间进行预处理之后就能做到在常数时间内计算出对任意(i,j)计算出LCP(i,j)。
因为lcp(Suffix(i),Suffix(j))=LCP(Rank,Rank[j]),所以我们也就可以在常数时间内求出S的任何两个后缀之间的最长公共前缀。这正是后缀数组能强有力地处理很多字符串问题的重要原因之一。
LCS (Longest Common Subsequence) 算法用于找出两个字符串最长公共子串。
算法原理:
(1) 将两个字符串分别以行和列组成矩阵。
(2) 计算每个节点行列字符是否相同,如相同则为 1。
(3) 通过找出值为 1 的最长对角线即可得到最长公共子串。
人 民 共 和 时 代
中 0, 0, 0, 0, 0, 0
华 0, 0, 0, 0, 0, 0
人 1, 0, 0, 0, 0, 0
民 0, 1, 0, 0, 0, 0
共 0, 0, 1, 0, 0, 0
和 0, 0, 0, 1, 0, 0
国 0, 0, 0, 0, 0, 0
为进一步提升该算法,我们可以将字符相同节点(1)的值加上左上角(d[i-1, j-1])的值,这样即可获得最大公用子串的长度。如此一来只需以行号和最大值为条件即可截取最大子串。
人 民 共 和 时 代
中 0, 0, 0, 0, 0, 0
华 0, 0, 0, 0, 0, 0
人 1, 0, 0, 0, 0, 0
民 0, 2, 0, 0, 0, 0
共 0, 0, 3, 0, 0, 0
和 0, 0, 0, 4, 0, 0
国 0, 0, 0, 0, 0, 0
算法实现:
public static string LCS(string s1, string s2)
{
if (s1 == s2)
return s1;
else if (String.IsNullOrEmpty(s1) || String.IsNullOrEmpty(s2))
return null;
var d = new int[s1.Length, s2.Length];
var index = 0;
var length = 0;
for (int i = 0; i < s1.Length; i++)
{
for (int j = 0; j < s2.Length; j++)
{
// 左上角值
var n = i - 1 >= 0 && j - 1 >= 0 ? d[i - 1, j - 1] : 0;
// 当前节点值 = "1 + 左上角值" : "0"
d[i, j] = s1[i] == s2[j] ? 1 + n : 0;
// 如果是最大值,则记录该值和行号
if (d[i, j] > length)
{
length = d[i, j];
index = i;
}
}
}
return s1.Substring(index - length + 1, length);
}
本文来自CSDN博客,转载请标明出处:http://blog.csdn.net/a497785609/archive/2011/04/28/6369707.aspx