【算法】最长公共子序列(nlogn)

转载注明出处：http://blog.csdn.net/wdq347/article/details/9001005 (修正了一些错误，并自己重写了代码)

---

最长公共子序列（LCS）最常见的算法是时间复杂度为O(n^2)的动态规划(DP)算法，但在James W. Hunt和Thomas G. Szymansky 的论文"A Fast Algorithm for Computing Longest Common Subsequence"中，给出了O(nlogn)下限的一种算法。

 

定理：设序列A长度为n，{A(i)}，序列B长度为m，{B(i)}，考虑A中所有元素在B中的序号，即A某元素在B的序号为{P_k1,P_k2,..}，将这些序号按照降序排列，然后按照A中的顺序得到一个新序列，此新序列的最长严格递增子序列即对应为A、B的最长公共子序列。

 

举例来说，A={a,b,c,d,b}，B={b,c,a,b}，则a对应在B的序号为2，b对应序号为{3,0}，c对应序号为1，d对应为空集，生成的新序列为{2, 3, 0, 1, 3, 0}，其最长严格递增子序列为{0,1,3}，对应的公共子序列为{b, c, b}

 

原论文的证明过程较复杂，其实可以简单的通过一一对应来证明。即证明A、B的一个公共子序列和新序列的一个严格递增子序列一一对应。

(1) A、B的一个公共子序列对应新序列的一个严格递增子序列

假设A、B的某一个公共子序列长度为k，则其公共子序列在A和B中可以写为

{A_i1,A_i2, ..., A_ik}

{B_j1,B_j2, ..., B_jk}

 

如此有A_i1 = A_j1,A_i2 = A_j2, ...., A_ik = A_jk, 考虑元素B_j1在B中的序号P(B_j1)，则有

P(B_j1)< P(B_j2) < ... < P(B_jk)

注意此严格递增子序列属于新序列的一个子序列，因此得证

 

(2) 新序列的一个严格递增子序列对应A、B的一个公共子序列

设新序列的一个严格递增子序列{P₁,P₂, ..., P_k}，任意两个相同的P不可能属于A中同一个元素，因为A中某元素在B中的序号按照降序排列，但此序列为严格递增序列，矛盾。所以每个P均对应于A中不同位置的元素，设为{A_i1, A_i2, ..., A_ik}。

因为P是严格递增序列，则每个P也对应B中唯一的一个元素，假设为{B_j1,B_j2, ..., B_jk}，由P的定义可知A_i1= B_j1, A_i2 = B_j2, ...., A_ik = B_jk，因此得证。

 

实现上比较复杂，有以下几个步骤：

(1) 对序列B排序

(2) 计算A中每个元素在B中的序号，并构成新序列

(3) 使用LIS的方法计算最长严格递增子序列

(4) 获取最长公共子序列

性能分析：

(1) 排序复杂度为nlogn

(2) 获取一个元素在B中的序号的复杂度，最小为logn，最大为n，获取所有元素的复杂度为 nlogn === n*n

(3) LIS 复杂度为nlogn

因此总体复杂度在nlogn 到 n*n logn之间，但如果(2) 步骤中A中元素在B中的序号对数很少时，性能相当优越，在实际测试时，string 中均为小写字母，长度为10000的情况下，这种方法比普通的LCS快一倍以上；如果string 中的字符扩展成char，即0-255，则这种方法比普通的LCS快至少一个数量级。

---

以下是参考代码：

#include<cstdio>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<set>
#include<map>
#include<stack>
#include<vector>
#include<queue>
#include<string>
#include<sstream>
#define eps 1e-9
#define ALL(x) x.begin(),x.end()
#define INS(x) inserter(x,x.begin())
#define FOR(i,j,k) for(int i=j;i<=k;i++)
#define MAXN 1005
#define MAXM 40005
#define INF 0x3fffffff
using namespace std;
typedef long long LL;

struct node
{
    char c;
    int num;
} u[10005];

int i,j,k,n,m,x,y,T,ans,big,cas,num,len;
bool flag;

bool cmp(node a,node b)
{
    if (a.c==b.c) return a.num>b.num;
    return a.c<b.c;
}

vector <int> p;
char a[10005],b[10005],c[10005];
int lena,lenb,dp[10005];

int main()
{
    scanf("%s",a);//读入a串 
    scanf("%s",b);//读入b串 
    
    lena=strlen(a);
    lenb=strlen(b);
    
    for (i=0;i<lenb;i++)
    {
         u[i].c=b[i];
         u[i].num=i;
    }
    
    sort(u,u+lenb,cmp);//对b串排序 
    
    for (i=0;i<lenb;i++)//排序后存入字符串c中，便于使用lower_bound 
    {
        c[i]=u[i].c;
    }
    c[lenb]='\0';
     
     for (i=0;i<lena;i++)//计算A中每个元素在B中的序号 
     {
         k=lower_bound(c,c+lenb,a[i])-c;
         while (k<lenb && a[i]==c[k])
         {
             p.push_back(u[k].num);
            k++;
         }
     }
     
     n=p.size();
     
    memset(dp,0,sizeof(dp));//计算最长上升子序列 
    num=0;
    for (i=0;i<n;i++)
    {
        if (p[i]>dp[num])
        {
            dp[++num]=p[i];
            
        }else
        {
            k=lower_bound(dp+1,dp+1+num,p[i])-dp; 
            dp[k]=p[i];
            
        }
    }
     
     printf("%d\n",num);
     
    return 0;
}