最长回文子串求解

题目：给定一个字符串，求其的最大回文子串。例如：字符串：owwoshisbsiha，它的最大回文子串是：hisbsih。

求解方法：暴力枚举、动态规划、后缀数组、线性算法

方法一：暴力枚举

最简单的方法当然就是对字符串的每一个子串进行回文判断。一个字符串有O(n²)个子串，然后判断是否回文复杂度是O(n)，所以该算法的算法复杂度是O(n³)。

方法二：动态规划

动态规划之所以能改进算法是因为该算法保留之前计算过的情况，这样就后面的情况就能转化为在前面已有的结果上进行求解，这也是动态规划和递归的本质的区别。废话少说直接进入主题，我觉得解这个题目的一个很自然的想法就是从把字符串的每一字符当做回文子串的中心点向两边延伸来计算得到以该字符为中心的最长回文子串，如上面字符串以b为中心就会得到最长回文子串hisbsih，当然这里遇到回文子串长度是偶数（如：owwo)的还不能解决，稍后在说这个问题。假设现在有回文子串X了，那么sXs也是回文，如果X不是回文子串，则sXs也不能是回文。可以看出是以某点为中心子串长度从小到大顺序来构建最长回文子串。使用 dp[i][j] 记录字符串从位置i到位置j的回文子串长度，然后在用两个变量mx和str_begin分别记录最大回文子串的长度和回文子串的在字符串的起始位置。

下面来解决回文子串长度是偶数的解决方法：在字符与字符之间插入一个特殊字符如#，来间隔开字符，这样回文子串的长度都会变成奇数了。那就剩实现了，这个方法的关键是动态规划怎样进行两层循环，才能到达动态规划——转化为能使用前面结果的情况或者说是避免计算过的情况再次计算的效果。

第一层循环：回文子串的长度（二分之一）的从 1 到字符串长度的½

第二层循环：回文子串的中心点的移动

这样做，后面的情况才能在前面已有的结果上进行求解，当然还有初始化条件：dp[i][i]=1，其实i=0,1,…,n-1;其他值都初始化为0，初始化的意图是每一个中心点就是一个回文子串且长度为1。

╔

Cpp代码  

1. int mx;  
2. bool dp[SIZE][SIZE]={};  
3. int str_begin=0;  
4. void LPS_dp(char * str, int length)   // 略去测试X合法性  
5. {  
6.     maxlen = 1;  
7.    
8.     for(int i = 0; i < length; ++i) // 初始化  
9.     {  
10.         dp[i][i] = 1;    // 单字符为回文  
11.     }  
12.    
13.     for(int len = 1; len <= length; ++len) //长度  
14.     {  
15.         for(int begin = 0; begin < length+1-len; ++begin)  
16.         {  
17.             int end = begin + len; // 从长度为2开始，首尾  
18.    
19.             if((X[begin]==X[end]) && (dp[begin+1][end-1]==1))  
20.             {  
21.                 dp[begin][end] = 1;  
22.                 if(end - begin + 1 > mx)  
23.                 {  
24.                     mx = end - begin + 1;  
25.                     str_begin=begin;  
26.                 }  
27.             }  
28.         }  
29.     }  
30. }  

╝①  

动态规划的算法复杂度都是O(n²)；这样也是减少重复计算的效果。

方法三：后缀数组

后缀数组，顾名思义就是从字符串某一个位置开始到结尾，例如：字符串dsqiu的后缀数组是dsqiu,sqiu,qiu,iu,u。然后对后缀数组进行排序（可以只以首字母来排序，规则可以自定义），排序之后后缀数组变为：dsqiu,iu,qiu,sqiu,u，排序的目的是方便进行枚举比较。

╔

后缀数组的应用：

 

例1：最长公共前缀

    给定一个串，求任意两个后缀的最长公共前缀。

解：先根据rank确定这两个后缀的排名i和j(i<j)，在height数组i+1和j之间寻找最小值。(可以用rmq优化)

 例2：最长重复子串(不重叠)(poj1743)

解：二分长度，根据长度len分组，若某组里SA的最大值与最小值的差>=len，则说明存在长度为len的不重叠的重复子串。

 例3：最长重复子串(可重叠)

解：height数组里的最大值。这个问题等价于求两个后缀之间的最长公共前缀。

 例4：至少重复k次的最长子串(可重叠)(poj3261)

解：二分长度，根据长度len分组，若某组里的个数>=k，则说明存在长度为len的至少重复k次子串。

 例5：最长回文子串(ural1297)

     给定一个串，对于它的某个子串，正过来写和反过来写一样，称为回文子串。

解：枚举每一位，计算以这个位为中心的的最长回文子串(注意串长要分奇数和偶数考虑)。将整个字符串反转写在原字符串后面，中间用$分隔。这样把问题转化为求某两个后缀的最长公共前缀。

 例6：最长公共子串(poj2774)

   给定两个字符串s1和s2，求出s1和s2的最长公共子串。

解：将s2连接到s1后，中间用$分隔开。这样就转化为求两个后缀的最长公共前缀，注意不是height里的最大值，是要满足sa[i-1]和sa[i]不能同时属于s1或者s2。

 例7：长度不小于k的公共子串的个数(poj3415)

给定两个字符串s1和s2，求出s1和s2的长度不小于k的公共子串的个数(可以相同)。

解：将两个字符串连接，中间用$分隔开。扫描一遍，每遇到一个s2的后缀就统计与前面的s1的后缀能产生多少个长度不小于k的公共子串，这里s1的后缀需要用单调栈来维护。然后对s1也这样做一次。

 例8：至少出现在k个串中的最长子串(poj3294)

给定n个字符串，求至少出现在n个串中k个的最长子串。

解：将n个字符串连接起来，中间用$分隔开。二分长度，根据长度len分组，判断每组的后缀是否出现在不小于k个原串中。

求解后缀数组的算法主要有两种：倍增算法和DC3算法。

╝②

算法实现：

1.反转字符串，并连接到原字符串后面，以一个特殊字符串（‘#’）间隔；

2.得到后缀字符串数组，并对后缀字符串数组进行快速排序；

3.枚举后缀字符串数组求解最大公共前缀（最大公共前缀：字符串从开头相同的子串）。

下面附上网友的有关倍增算法和DC3算法的代码（没有测试）

倍增算法

╔

 

Cpp代码  

1. const int N = 20005;//串A的最大长度  
2. const int MAX = 1000100;//串A的最大值  
3. //int n,m,k;  
4. int SA[N], rank[N], height[N], key[N];  
5. int A[N], C[MAX], t1[N+1], t2[N+1];  
6.   
7. //倍增法求sa[]-----待排序的字符串放在r 数组中，r[]为整型数组, 从r[0]到r[n-1]，长度为n，且最大值小于m  
8. //约定除r[n-1]外所有的r[i]都大于0, r[n-1]=0  
9. //结果放在sa 数组中，从sa[0]到sa[n-1]  
10. //先对所有后缀的第一个字符进行排序(采用挖坑式的基数排序,即统计每个字符的个数,以便在扫描时总能将字符放入合适的位置),放入sa中  
11. void da(int n, int m)  
12. {  
13.     int i, j, l, s,*t;  
14.     int *X = t1, *Y = t2;  
15.     memset(C, 0, sizeof(C));  
16.     for (i=0;i<n;i++) C[X[i] = A[i]]++;  
17.     for (i=1;i<m;i++) C[i] += C[i-1];  
18.     for (i=n-1;i>=0;i--) SA[--C[X[i]]] = i;  
19.     for (l=1; l<n; l<<=1)  
20.     {  
21.         for (i=n-l,j=0;i<n;i++) Y[j++] = i;  
22.         for (i=0;i<n;i++) if (SA[i] >= l) Y[j++] = SA[i] - l;  
23.         for (i=0;i<n;i++) key[i] = X[Y[i]];  
24.         memset(C, 0, sizeof(C));  
25.         for (i=0;i<n;i++) C[key[i]]++;  
26.         for (i=1;i<m;i++) C[i] += C[i-1];  
27.         for (i=n-1;i>=0;i--) SA[--C[key[i]]] = Y[i];  
28.         t = X;  
29.         X = Y;  
30.         Y = t;  
31.         X[SA[0]] = j = 0;  
32.         for (i=1;i<n;i++)  
33.         {  
34.             if (Y[SA[i]] != Y[SA[i-1]] || Y[SA[i]+l] != Y[SA[i-1]+l])  
35.                 j++;  
36.             X[SA[i]] = j;  
37.         }  
38.         m = j + 1;  
39.         if (m==n) break;  
40.     }  
41.    
42.     for (i=0;i<n;i++)  
43.         rank[SA[i]] = i;  
44.     return;  
45. }  
46.   
47. //height[i]:suffix(sa[i-1])与suffix(sa[i])的最长公共前缀，即排名相邻的两个后缀的最长公共前缀  
48. void calheight(int n)  
49. {  
50.     int i,j,k=0;  
51.     for(i=0; i<n; i++)  
52.     {  
53.         if (k > 0)  
54.             --k;  
55.         if(rank[i] == 0)  
56.             height[0] = 0;  
57.         else  
58.         {  
59.             j = SA[rank[i] - 1];  
60.             while(A[i+k] == A[j+k])  
61.                 k++;  
62.             height[rank[i]] = k;  
63.         }  
64.     }  
65. }  
66. //串A[0]...A[n-1]  
67. da(n,1000001);  //m的最大值不超过1,000,000  
68. calheight(n);  

 ╝②

 

DC3算法

╔

 

C代码  

1. #include "stdio.h"  
2. #include "string.h"  
3. #define maxn 2004  
4.   
5. #define F(x) ((x)/3+((x)%3==1?0:tb))  
6. #define G(x) ((x)<tb?(x)*3+1:((x)-tb)*3+2)  
7. int wa[maxn],wb[maxn],wv[maxn],ws[maxn];  
8. int c0(int *r,int a,int b)  
9. {return r[a]==r[b]&&r[a+1]==r[b+1]&&r[a+2]==r[b+2];}  
10. int c12(int k,int *r,int a,int b)  
11. {if(k==2) return r[a]<r[b]||r[a]==r[b]&&c12(1,r,a+1,b+1);  
12.  else return r[a]<r[b]||r[a]==r[b]&&wv[a+1]<wv[b+1];}  
13. void sort(int *r,int *a,int *b,int n,int m)  
14. {  
15.      int i;  
16.      for(i=0;i<n;i++) wv[i]=r[a[i]];  
17.      for(i=0;i<m;i++) ws[i]=0;  
18.      for(i=0;i<n;i++) ws[wv[i]]++;  
19.      for(i=1;i<m;i++) ws[i]+=ws[i-1];  
20.      for(i=n-1;i>=0;i--) b[--ws[wv[i]]]=a[i];  
21.      return;  
22. }  
23. void dc3(int *r,int *sa,int n,int m)  
24. {  
25.      int i,j,*rn=r+n,*san=sa+n,ta=0,tb=(n+1)/3,tbc=0,p;  
26.      r[n]=r[n+1]=0;  
27.      for(i=0;i<n;i++) if(i%3!=0) wa[tbc++]=i;  
28.      sort(r+2,wa,wb,tbc,m);  
29.      sort(r+1,wb,wa,tbc,m);  
30.      sort(r,wa,wb,tbc,m);  
31.      for(p=1,rn[F(wb[0])]=0,i=1;i<tbc;i++)  
32.      rn[F(wb[i])]=c0(r,wb[i-1],wb[i])?p-1:p++;  
33.      if(p<tbc) dc3(rn,san,tbc,p);  
34.      else for(i=0;i<tbc;i++) san[rn[i]]=i;  
35.      for(i=0;i<tbc;i++) if(san[i]<tb) wb[ta++]=san[i]*3;  
36.      if(n%3==1) wb[ta++]=n-1;  
37.      sort(r,wb,wa,ta,m);  
38.      for(i=0;i<tbc;i++) wv[wb[i]=G(san[i])]=i;  
39.      for(i=0,j=0,p=0;i<ta && j<tbc;p++)  
40.      sa[p]=c12(wb[j]%3,r,wa[i],wb[j])?wa[i++]:wb[j++];  
41.      for(;i<ta;p++) sa[p]=wa[i++];  
42.      for(;j<tbc;p++) sa[p]=wb[j++];  
43.      return;  
44. }  
45. int rank[maxn],height[maxn];  
46. void calheight(int *r,int *sa,int n)  
47. {  
48.      int i,j,k=0;  
49.      for(i=1;i<=n;i++) rank[sa[i]]=i;  
50.      for(i=0;i<n;height[rank[i++]]=k)  
51.      for(k?k--:0,j=sa[rank[i]-1];r[i+k]==r[j+k];k++);  
52.      return;  
53. }  
54. int RMQ[maxn];  
55. int mm[maxn];  
56. int best[20][maxn];  
57. void initRMQ(int n)  
58. {  
59.      int i,j,a,b;  
60.      for(mm[0]=-1,i=1;i<=n;i++)  
61.      mm[i]=((i&(i-1))==0)?mm[i-1]+1:mm[i-1];  
62.      for(i=1;i<=n;i++) best[0][i]=i;  
63.      for(i=1;i<=mm[n];i++)  
64.      for(j=1;j<=n+1-(1<<i);j++)  
65.      {  
66.        a=best[i-1][j];  
67.        b=best[i-1][j+(1<<(i-1))];  
68.        if(RMQ[a]<RMQ[b]) best[i][j]=a;  
69.        else best[i][j]=b;  
70.      }  
71.      return;  
72. }  
73. int askRMQ(int a,int b)  
74. {  
75.     int t;  
76.     t=mm[b-a+1];b-=(1<<t)-1;  
77.     a=best[t][a];b=best[t][b];  
78.     return RMQ[a]<RMQ[b]?a:b;  
79. }  
80. int lcp(int a,int b)  
81. {  
82.     int t;  
83.     a=rank[a];b=rank[b];  
84.     if(a>b) {t=a;a=b;b=t;}  
85.     return(height[askRMQ(a+1,b)]);  
86. }  
87.   
88. char st[maxn];  
89. int r[maxn*3],sa[maxn*3];  
90. int main()  
91. {  
92.     int i,n,len,k,ans=0,w;  
93.     scanf("%s",st);  
94.     len=strlen(st);  
95.     for(i=0;i<len;i++) r[i]=st[i];  
96.     r[len]=1;  
97.     for(i=0;i<len;i++) r[i+len+1]=st[len-1-i];  
98.     n=len+len+1;  
99.     r[n]=0;  
100.     dc3(r,sa,n+1,128);  
101.     calheight(r,sa,n);  
102.     for(i=1;i<=n;i++) RMQ[i]=height[i];  
103.     initRMQ(n);  
104.     for(i=0;i<len;i++)  
105.     {  
106.       k=lcp(i,n-i);  
107.       if(k*2>ans) ans=k*2,w=i-k;  
108.       k=lcp(i,n-i-1);  
109.       if(k*2-1>ans) ans=k*2-1,w=i-k+1;  
110.     }  
111.     st[w+ans]=0;  
112.     printf("%s\n",st+w);  
113.     return 0;  
114. }  

 

 ╝③

拓展

后缀数组的算法复杂度是O(n㏒n)，主要是由排序引起的。那么，就会想到要是不经过排序的过程或者在构建后缀数组的过程就已经排好序，算法复杂度就会降到O(n)。这就得使用后缀树来完成构建后缀子串，这里后缀树就是子串中每一个字符都是后缀树的一个节点，如果两个前缀一样那么它们就拥有共同父亲节点。在构建后缀树的过程就记录从当前节点开始的最长公共前缀的长度，构建完成之后只要遍历一遍后缀树找到最长公共前缀，就是要找的最大回文子串的一半（如最长回文子串是abcdcba，得到的最长公共前缀是abcd）。这里说的比较简单，不过我觉得看到这里的理解应该都没问题吧，当然后缀树还有很多应用（如数据挖掘的FP-Growth Algorithm的FP tree）。

更多分析可以阅读参考④的内容。

 

方法四：线性算法

算法复杂度要尽可能小，一个优化方法就是避免之前情况的重复计算，正如前面动态规划对暴力枚举的改进——保留之前已经计算过的情况的结果，后面的情况转为在前面记录基础之上来算。无所不用其极，自然会想在动态规划有没有将已经计算结果充分利用殆尽。还是可以发现还有一个特征没有被利用，动态规划只保留前面情况的结果（利用仅此而已），其实真正的主角——回文还没有利用，就是只利用了回文子串的是不是的结果，但是没有利用回文子串对称的结果。

用一个数组 P[i] 来记录以字符S[i]为中心的最长回文子串向左/右扩张的长度，增加两个辅助变量（其实一个就够了，两个更清晰）id和mx，其中id表示最大回文子串中心的位置，mx则为id+P[id]，也就是最大回文子串的边界。

在计算以第 i 位置的字符为中心的最长回文子串时，所有在 i 前面的情况都计算过了，换而言之，p[j]都是已知，当j<i的时候，那么就充分利用这个性质吧。

看下图中，现在计算以第 i 位置的字符为中心的最长回文子串，前面已有的结果是最长回文子串的中心位置是 id 长度为 mx，j 是 i 关于 id 的对称点，如果mx>i，那么 i 就可以获得一个已知的回文子串的长度，即p[j]的值，就是在 j 的回文子串在 i 处同样会出现（对称嘛），这样就不用像动态规划那样每个位置begin的len都从1开始增大，而是直接从p[j]开始。

 

算法实现：

╔

 

Cpp代码  

1. /* O(n)解法 */  
2. #define MIN(a,b) ((a) < (b) ? (a) : (b))  
3. int mx;  
4. int maxid;        // 最长回文子串下标  
5. int LPS_rb[100];  // i为中心的回文子串右边界下标right border  
6. char str[100];    // 原字符串处理后的副本  
7.    
8. void LPS_linear(char * X, int xlen)  
9. {  
10.     mx= maxid = 0;  
11.    
12.     str[0] = '$';  // 将原串处理成所需的形式  
13.     char *p = str;  
14.     *(++p)++ = '#';  
15.     while((*p++ = *X++) != '\0')  
16.     {  
17.         *p++ = '#';  
18.     }  
19.    
20.     for(int i = 1; str[i]; ++i)  // 计算LPS_rb的值  
21.     {  
22.         if(maxlen > i)          // 初始化LPS[i]  
23.         {  
24.             LPS_rb[i] = MIN(LPS_rb[2*maxid-i],(mx));  
25.         }else  
26.         {  
27.             LPS_rb[i] = 1;  
28.         }  
29.         while(str[i-LPS_rb[i]] == str[i+LPS_rb[i]]) // 扩展  
30.         {  
31.             ++LPS_rb[i];  
32.         }  
33.         if(LPS_rb[i]-1-i > mx)  
34.         {  
35.             mx = LPS_rb[i]-1-i;  
36.             maxid = i;  
37.         }  
38.     }  
39. }  
40. 给出测试用例：  
41.   
42. void main()  
43. {  
44.     char X[30];  // 设串不超过30  
45.     /* test case 
46.      * aaaa 
47.      * abab 
48.      */  
49.     while(cin.getline(X,30))  
50.     {  
51.         /* 后缀数组方法 */  
52.         LPS_suffix(X,strlen(X));  
53.         printf("%d\n", maxlen);  
54.    
55.         /* O(n)方法 */  
56.         LPS_linear(X,strlen(X));  
57.         printf("%d\n", maxlen);  
58.     }  
59. }  

╝④

 

小结：

至此，将四种方法全面列举完毕，我觉得至少看出算法优化的一个范例，算法优化无非就是穷尽事物特征，无所不用其极，如果你有任何建议或者批评和补充，请留言指出，不胜感激，更多参考请移步互联网。

 

类似问题

最大公共子串，最大公共子序列，最长重复子串，最长不重复子串，最长递增序列，最大子数组和

参考：

 

①勇幸|Thinking： http://www.ahathinking.com/archives/132.html

②yzmduncan： http://yzmduncan.iteye.com/blog/979771

③小阮的菜田： http://www.cppblog.com/jericho/archive/2011/06/30/149862.aspx

④O興~： http://imlazy.ycool.com/post.2011818.html

⑤ felix021： http://www.felix021.com/blog/read.php?2040