后缀数组求最长回文串
[cpp] view plaincopy
1. #include<cstdio>
2. #include<cmath>
3. #include<cstring>
4. #include<algorithm>
5. using namespace std;
6. #define MAXN 2010
7. int dp[2010][30];
8. char r[MAXN],rr[MAXN];
9. int sa[MAXN];
10. int wa[MAXN],wb[MAXN],wv[MAXN],ws[MAXN];
11. int height[MAXN],rk[MAXN];
12. inline bool cmp(int *r,int a,int b,int len){
13. return r[a]==r[b]&&r[a+len]==r[b+len];
14. }
15. void SA(int n,int m){
16. int i,j,p,*x=wa,*y=wb,*t;
17. for(i=0;i<m;i++)
18. ws[i]=0;
19. for(i=0;i<n;i++)
20. ws[x[i]=r[i]]++;
21. for(i=1;i<m;i++)
22. ws[i]+=ws[i-1];
23. for(i=n-1;i>=0;i--)
24. sa[--ws[x[i]]]=i;
25. for(j=p=1;p<n;j<<=1,m=p){
26. for(p=0,i=n-j;i<n;i++)
27. y[p++]=i;
28. for(i=0;i<n;i++){
29. if(sa[i]>=j)
30. y[p++]=sa[i]-j;
31. }
32. for(i=0;i<m;i++)
33. ws[i]=0;
34. for(i=0;i<n;i++)
35. ws[wv[i]=x[y[i]]]++;
36. for(i=1;i<m;i++)
37. ws[i]+=ws[i-1];
38. for(i=n-1;i>=0;i--)
39. sa[--ws[wv[i]]]=y[i];
40. for(t=x,x=y,y=t,x[sa[0]]=0,p=i=1;i<n;i++)
41. x[sa[i]]=cmp(y,sa[i-1],sa[i],j)?p-1:p++;
42. }
43. }
44. void Height(int n){
45. int i,j,k=0;
46. for(i=0;i<=n;i++) //这里sa[0]为‘\0’开始的子串
47. rk[sa[i]]=i;
48. for(i=0;i<n;height[rk[i++]]=k)
49. for(k?k--:0,j=sa[rk[i]-1];r[i+k]==r[j+k];k++);
50. }
51. int init(){ //把串复制一遍倒着连在rr[]后面
52. int i,j,len;
53. memset(height,0,sizeof(height));
54. len=strlen(rr);
55. for(i=0;i<len;i++)
56. r[i]=rr[i];
57. r[i]='$';
58. for(j=0;j<len;j++)
59. r[j+len+1]=r[len-1-j];
60. r[j+len+1]='\0';
61. return len*2+1;
62. }
63. void st(int n){
64. int i,j,p,q;
65. for(i=1;i<=n;i++)
66. dp[i][0]=i;
67. for(j=1;j<=(int)(log((double)n)/log(2.0));j++)
68. for(i=1;i+(1<<j)-1<=n;i++){
69. p=height[dp[i][j-1]];
70. q=height[dp[i+(1<<(j-1))][j-1]];
71. if(p>q)
72. dp[i][j]=dp[i+(1<<(j-1))][j-1];
73. else
74. dp[i][j]=dp[i][j-1];
75. }
76.
77. }
78. int RMQ_MIN(int i,int j){
79. int tem;
80. if(i>j){
81. tem=i;
82. i=j;
83. j=tem;
84. }
85. i++; //交换后小的要加一
86. int k=(int)(log((double)(j-i+1))/log(2.0));
87. return min(height[dp[i][k]],height[dp[j-(1<<k)+1][k]]);
88. }
89. void solve(int n){
90. int i,j,ans=0,s;
91. st(n);
92. for(i=0;i<n/2;i++){
93. j=RMQ_MIN(rk[i],rk[n-1-i]);
94. if(j*2-1>ans){
95. ans=j*2-1;
96. s=i-j+1;
97. }
98. j=RMQ_MIN(rk[i],rk[n-i]);
99. if(j*2>ans){
100. ans=j*2;
101. s=i-j;
102. }
103. }
104. for(i=s;i<s+ans;i++)
105. printf("%c",rr[i]);
106. printf("\n");
107. }
108. int main(){
109. int i,j,n;
110. scanf("%s",rr);
111. n=init();
112. SA(n+1,130);
113. Height(n);
114. solve(n);
115. }
题意：给定两个字符串 A和 B，求最长公共子串。
思路：后缀数组。（摘自罗穗骞的国家集训队论文）字符串的任何一个子串都是这个字符串的某个后缀的前缀。求 A和 B 的最长公共子串等价于求 A 的后缀和 B 的后缀的最长公共前缀的最大值。如果枚举A和 B的所有的后缀，那么这样做显然效率低下。由于要计算 A的后缀和 B 的后缀的最长公共前缀，所以先将第二个字符串写在第一个字符串后面，中间用一个没有出现过的字符隔开，再求这个新的字符串的后缀数组。观察一下，看看能不能从这个新的字符串的后缀数组中找到一些规律。以 A=“ aaaba ”，B=“ abaa ”为例，如图 8 所示。
那么是不是所有的 height值中的最大值就是答案呢？不一定！有可能这两个后缀是在同一个字符串中的，所以实际上只有当suffix(sa[i-1])和suffix(sa[i])不是同一个字符串中的两个后缀时，height[i]才是满足条件的。而这其中的最大值就是答案。记字符串 A和字符串 B 的长度分别为|A|和|B|。求新的字符串的后缀数组和 height数组的时间是 O(|A|+|B|)，然后求排名相邻但原来不在同一个字符串中的两个后缀的height值的最大值，时间也是O(|A|+|B|)，所以整个做法的时间复杂度为 O(|A|+|B|) 。时间复杂度已经取到下限，由此看出，这是一个非常优秀的算法。
源代码：(5676K，313MS)
#include<iostream>
#include<cstring>
#definemax(a,b) ((a)>(b)?(a):(b))
usingnamespace std;
const intMAX = 200050;
intsa[MAX], rank[MAX], height[MAX];
intwa[MAX], wb[MAX], wv[MAX], wd[MAX];
intcmp(int *r, int a, int b, int l)
{
returnr[a] == r[b] && r[a+l] == r[b+l];
}
voidda(int *r, int n, int m) // 倍增算法0(nlgn)。
{
inti, j, p, *x = wa, *y = wb, *t;
for(i= 0; i < m; i ++) wd[i] = 0;
for(i= 0; i < n; i ++) wd[x[i]=r[i]] ++;
for(i= 1; i < m; i ++) wd[i] += wd[i-1];
for(i= n-1; i >= 0; i --) sa[-- wd[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 ++) wd[i] = 0;
for(i= 0; i < n; i ++) wd[wv[i]] ++;
for(i= 1; i < m; i ++) wd[i] += wd[i-1];
for(i= n-1; i >= 0; i --) sa[-- wd[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 ++;
}
}
}
voidcalHeight(int *r, int n) // 求height数组。
{
inti, j, k = 0;
for(i= 1; i <= n; i ++) rank[sa[i]] = i;
for(i= 0; i < n; height[rank[i ++]] = k)
{
for(k? k -- : 0, j = sa[rank[i]-1]; r[i+k] == r[j+k]; k ++);
}
}
int main()
{
inti, k, l1, l2, num[MAX];
charstr[MAX];
scanf("%s",str);
l1= strlen(str);
for(k= i = 0; i < l1; i ++)
{
num[k++] = str[i] - 'a' + 2;
}
num[k++] = 1; // 相当于论文中的'$'分隔符。
scanf("%s",str);
l2= strlen(str);
for(i= 0; i < l2; i ++)
{
num[k++] = str[i] - 'a' + 2;
}
intn = l1 + l2;
da(num,n + 1, 30);
calHeight(num,n);
intans = 0;
for(i= 2; i < k; i ++)
if((sa[i]< l1 && sa[i-1] > l1) || (sa[i-1] < l1 && sa[i] >l1))
{
ans= max(ans, height[i]);
}
printf("%d\n",ans);
return0;
}
例2：最长重复子串(不重叠)(poj1743)
单个字符串的相关问题：
这类问题的一个常用做法是先求后缀数组和height数组，然后利用height
数组进行求解。
【题意】
一篇n个音符的乐谱，求他最长的音乐主题
乐谱由n个小于88的正整数表示
所谓的音乐主题就是两端变化相同的长度大于五的不重合乐谱
比如1 2 3 4 5 6 7 89 10
1 2 3 4 5的变化与6 7 8 9 10就相同，且两端不重合
【输入】
多组数据
每组数据第一行为n，表示音符的个数
接下来n个数字表示各个音符
输入以n=0结束
【输出】
对于每组数据输出一个数表示最长的音乐主题，音乐主题不存在输出0
后缀数组的典型应用
09年论文《后缀数组——处理字符串的有力工具》上的例题
大意是先对变化求height，然后二分答案，将问题转化为证明是否存在长度为k的不重合序列
根据后缀数组的性质可知height数组是呈波峰状的，所以按height大于k的条件将其分组，若有一组中最小的sa和最大的sa之差大于k，则存在。
需要注意的是，n=1的情况。
容易看出，有希望成为最长公共前缀不小于k的两个后缀一定在同一组。然
后对于每组后缀，只须判断每个后缀的sa值的最大值和最小值之差是否不小于
k。如果有一组满足，则说明存在，否则不存在。整个做法的时间复杂度为
O(nlogn)。本题中利用height值对后缀进行分组的方法很常用，请读者认真体
会。
参考：http://blog.csdn.net/weixinding/article/details/7219236
参考：http://bbezxcy.iteye.com/blog/1407395
待排序的字符串放在r数组中，从r[0]到r[n-1]，长度为n，且最大值小
于m。为了函数操作的方便，约定除r[n-1]外所有的r[i]都大于0, r[n-1]=0。
函数结束后，结果放在sa数组中，从sa[0]到sa[n-1]。
height 数组：定义height[i]=suffix(sa[i-1])和suffix(sa[i])的最长公
共前缀，也就是排名相邻的两个后缀的最长公共前缀。
#include <cstring>
#include <cstdio>
using namespace std;
const int N = 200001;
int num[N];
int sa[N],rank[N],height[N];
int wa[N],wb[N],wv[N],wd[N];
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 n,int m){//倍增算法
int i, j, p,*x = wa,*y = wb,*t;
for(i = 0; i < m; i++) wd[i] = 0;
for(i = 0; i < n; i++) wd[x[i]=r[i]]++;
for(i = 1; i < m; i++) wd[i]+=wd[i-1];
for(i = n-1; i >= 0; i--) sa[--wd[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++) wd[i] =0;
for(i = 0; i < n; i++)wd[wv[i]]++;
for(i = 1; i < m; i++)wd[i]+=wd[i-1];
for(i = n-1; i >= 0; i--) sa[--wd[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++;
}
}
}
void calheight(int *r,int n){
int i,j,k = 0;
for(i = 1; i <= n; i++) rank[sa[i]] = i;
for(i = 0; i < n; height[rank[i++]] = k){
for(k ? k-- : 0, j = sa[rank[i]-1];r[i+k] == r[j+k]; k++);
}
}
int loc[N],m;
char str[N],res[N];
bool vis[1004];
int abs(int a) {return a >= 0 ? a : -a; }
bool check(int mid,int len){
int i,j;
for(i = 2; i <= len; i++){
if(height[i] < mid) continue;
for(j = i - 1; j >= 2; j--){
if(abs(sa[i]-sa[j])>=mid)
return 1;
if(height[j]< mid) break;
}
}
return 0;
}
int main(){
int n,k,i,j,a,b,sp,ans;
while(scanf("%d",&n) && n){
n--;
ans = 0;
for(i = 0; i <= n; i++)
scanf("%d",&num[i]);
if(n < 10){
printf("0\n");
continue;
}
for(i = 0; i < n; i++)
num[i] =num[i+1]-num[i]+90;//保证每个值>0
num[n] = 0;
da(num,n+1,200);
calheight(num,n);
int l = 0,r =n,mid;
while(l <= r){
mid =(l+r)/2;
if(check(mid,n)){
l = mid+1;
ans = mid;
}else r =mid - 1;
}
if(ans < 4)printf("0\n");
else printf("%d\n",ans+1);
}
return 0;
}
例4：至少重复k次的最长子串(可重叠)(poj3261)
给定一个字符串，求至少出现K次的最长重复子串，这k个子串可以重叠。所谓出现K次就是在任意K个地方出现了这个子串，并不要求这K个是连续的。那么根据罗大神的论
| 给定一个字符串，求至少出现K次的最长重复子串，51自学网，这k个子串可以重叠。 所谓出现K次就是在任意K个地方出现了这个子串，并不要求这K个是连续的。
二分答案，51自学网，51自学网，然后用height数组判定是否至少出现了K个这种子串
[cpp]
例7：长度不小于k的公共子串的个数(poj3415) 统计两个串中长度>=k的公共子串的数量。 【题解】 很容易想到将两个串连接起来，求后缀数组，在根据k分组 但统计是件麻烦事。。。 以下是我自己搞出来的方法： 搞一个单调的栈维护lcp，然后维护单调性。 统计的时候要注意，只有不同的串的lcp才能统计。 这个方法严格的来说不能算是O(n)，但是很接近。 结合代码分析： |
1. #include <iostream>
2. #include <cstring>
3. using namespace std;
4. const int maxn=200010;
5.
6. char s[maxn],s2[maxn];
7. int sa[maxn],rk[maxn],wa[maxn],wb[maxn],w[maxn],v[maxn],h[maxn];
8. int n,len,kk;
9.
10. int cmp(int* r,int a,int b,int l)
11. {
12. return r[a]==r[b] && r[a+l]==r[b+l];
13. }
14.
15. void da(char* s,int* sa,int n,int m)
16. {
17. int *x=wa,*y=wb,*t,i,j,p;
18. for (i=0;i<m;i++) w[i]=0;
19. for (i=0;i<n;i++) w[x[i]=s[i]]++;
20. for (i=1;i<m;i++) w[i]+=w[i-1];
21. for (i=n-1;i>=0;i--) sa[--w[x[i]]]=i;
22. for (j=1,p=1;p<n;j*=2,m=p)
23. {
24. for (p=0,i=n-j;i<n;i++) y[p++]=i;
25. for (i=0;i<n;i++) if (sa[i]>=j) y[p++]=sa[i]-j;
26. for (i=0;i<n;i++) v[i]=x[y[i]];
27. for (i=0;i<m;i++) w[i]=0;
28. for (i=0;i<n;i++) w[v[i]]++;
29. for (i=1;i<m;i++) w[i]+=w[i-1];
30. for (i=n-1;i>=0;i--) sa[--w[v[i]]]=y[i];
31. for (t=x,x=y,y=t,p=1,x[sa[0]]=0,i=1;i<n;i++)
32. x[sa[i]]=cmp(y,sa[i-1],sa[i],j)?p-1:p++;
33. }
34. }
35.
36. void calheight()
37. {
38. int i,j,k=0;
39. for (i=1;i<=n;i++) rk[sa[i]]=i;
40. for (i=0;i<n;h[rk[i++]]=k)
41. for (k?k--:k=0,j=sa[rk[i]-1];s[i+k]==s[j+k];k++);
42. }
43.
44. long long work()
45. {
46. int i,top=0,tmp;
47. long long ans=0,ss[3]={0,0,0};
48. for (i=1;i<=n;i++)
49. if (h[i]<kk)
50. {
51. top=ss[1]=ss[2]=ss[0]=0;
52. }
53. else
54. {
55. for (tmp=top;tmp && v[tmp]>h[i]-kk+1;tmp--)//维护单调性
56. {
57. ss[w[tmp]]+=h[i]-kk+1-v[tmp];
58. v[tmp]=h[i]-kk+1;
59. }
60. v[++top]=h[i]-kk+1;
61. if (sa[i-1]<len) w[top]=1;//注意是sa[i-1]，来判断前一个串是A串还是B串
62. if (sa[i-1]>len) w[top]=2;
63. ss[w[top]]+=h[i]-kk+1;
64. int t;
65. if (sa[i]<len) t=1;//自己是A串还是B串
66. if (sa[i]>len) t=2;
67. ans+=ss[3-t];//累加
68. }
69.
70. return ans;
71. }
72.
73. int main()
74. {
75. freopen("pin.txt","r",stdin);
76. freopen("pou.txt","w",stdout);
77. while (1)
78. {
79. cin >> kk;
80. if (kk==0) break;
81. cin >> s;
82. len=strlen(s);
83. cin >> s2;
84. strcat(s,"$");
85. strcat(s,s2);
86. n=strlen(s);
87. s[n]=0;
88. da(s,sa,n+1,128);
89. calheight();
90. cout << work() << endl;
91. }
92. return 0;
93. }
94.
例8：至少出现在k个串中的最长子串(poj3294)
1. /*
2. 倍增算法
3. */
4. #include<cstdio>
5. #include<cstring>
6. const int maxl=1010;
7. const int maxn=101;
8. bool vis[110];
9. int num[110];
10. #define max maxl*maxn
11. int w[max],wa[max],wb[max],wv[max];
12. int sa[max],height[max],rank[max];
13. int r[max],n,t,len,a[maxn];
14. char s[maxl];
15. int cmp(int *r,int a,int b,int l)
16. {
17. return r[a]==r[b]&&r[a+l]==r[b+l];
18. }
19. void da(int *r,int *sa,int n,int m)
20. {
21. int i,j,p,*x=wa,*y=wb,*t;
22. for (i=0; i<m; i++) w[i]=0;
23. for (i=0; i<n; i++) w[x[i]=r[i]]++;
24. for (i=1; i<m; i++) w[i]+=w[i-1];
25. for (i=n-1; i>=0; i--) sa[--w[x[i]]]=i;
26. for (p=1,j=1; p<n; m=p,j*=2)
27. {
28. for (p=0,i=n-j; i<n; i++) y[p++]=i;
29. for (i=0; i<n; i++) if (sa[i]>=j) y[p++]=sa[i]-j;
30. for (i=0; i<m; i++) w[i]=0;
31. for (i=0; i<n; i++) w[wv[i]=x[y[i]]]++;
32. for (i=1; i<m; i++) w[i]+=w[i-1];
33. for (i=n-1; i>=0; i--) sa[--w[wv[i]]]=y[i];
34. for (t=x,x=y,y=t,x[sa[0]]=0,p=1,i=1; i<n; i++)
35. x[sa[i]]=cmp(y,sa[i-1],sa[i],j)?p-1:p++;
36. }
37. return;
38. }
39. void cal(int *r,int *sa,int n)
40. {
41. int i,j,k=0;
42. for (i=1; i<=n; i++) rank[sa[i]]=i;
43. for (i=0; i<n; height[rank[i++]]=k)
44. for (k?k--:0,j=sa[rank[i]-1]; r[i+k]==r[j+k]; k++);
45. return;
46. }
47. int caa(int x)//找到子串的位置
48. {
49. for(int i=1,len=num[1]; i<=t; ++i)
50. {
51. if(x==len) return 0;
52. else if(x<len) return i;
53. len=len+num[i+1]+1;
54. }
55. }
56. bool check(int k)//注意绝对不要把自己加入进去区分两个串的字符算进去
57. {
58. int x,d=0,i;
59. memset(vis,0,sizeof(vis));
60. vis[0]=1;//放置自己加进去的特殊字符的
61. //printf("mid=%d n=%d t=%d/n",k,n,t);
62. x=caa(sa[0]);
63. if(vis[x]==0) d++;
64. for(i=1; i<=n; ++i)
65. {
66. if(height[i]<k)
67. {
68. d=0;
69. memset(vis,0,sizeof(vis));
70. vis[0]=1;
71. //printf("%d/n",i);
72. }
73. x=caa(sa[i]);
74. if(vis[x]==0) d++;
75. if(d>t/2) return 1;
76. //printf("x=%d/n",x);
77. vis[x]=1;
78. }
79. return 0;
80. }
81. void print(int k)
82. {
83. int x,d=0,i;
84. memset(vis,0,sizeof(vis));
85. vis[0]=1;//放置自己加进去的特殊字符的
86. //printf("mid=%d n=%d t=%d/n",k,n,t);
87. x=caa(sa[0]);//对第一个进行初始化
88. if(vis[x]==0) d++;
89. for(i=1; i<=n; ++i)
90. {
91. if(height[i]<k)
92. {
93. if(d>t/2)
94. {
95. for(int j=sa[i-1]; j<sa[i-1]+k; ++j) printf("%c",r[j]);
96. puts("");
97. }
98. d=0;
99. memset(vis,0,sizeof(vis));
100. vis[0]=1;
101. //printf("%d/n",i);
102. }
103. x=caa(sa[i]);
104. if(vis[x]==0) d++;
105. vis[x]=1;
106. }
107. }
108. int rmax;
109. void solve()
110. {
111. int left=1,right=rmax;
112. int ans=0;
113. while(left<=right)
114. {
115. int mid=(left+right)>>1;
116. if(check(mid))
117. {
118. ans=mid;
119. left=mid+1;
120. }
121. else
122. {
123. right=mid-1;
124. }
125. }
126. if(ans==0) printf("?/n");
127. else print(ans);
128. }
129. int main()
130. {
131. int cas=0;
132. while(scanf("%d",&t)!=EOF)
133. {
134. if(t==0) break;
135. memset(num,0,sizeof(num));
136. if(cas==0) cas=1;//输出空行
137. else printf("/n");
138. int len=0;
139. n=0;
140. rmax=-1;
141. for(int i=1; i<=t; ++i) //i<100,lower case letter
142. {
143. scanf("%s",s);
144. len=strlen(s);
145. num[i]=len;
146. if(len>rmax) rmax=len;
147. for(int j=0; j<len; ++j) r[n+j]=(int)s[j];
148. r[n+len]=i+ 200;
149. n=n+len+1;//用不可能出现的符号分割各个串
150. }
151. if(t==1)
152. {
153. printf("%s/n",s);
154. continue;
155. }
156.
157. da(r,sa,n+1,350);
158. //printf("2222222/n");
159. cal(r,sa,n);
160. //printf("1111111111111/n");
161. solve();
162. // puts("");
163. }
164. return 0;
165. }
585
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
