最长公共子序列(nlogn)

最长公共子序列（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对应序号为{4,0}，c对应序号为1，d对应为空集，生成的新序列为{2,4, 0, 1, 4, 0}，其最长严格递增子序列为{0,1,4}，对应的公共子序列为{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= A_j1, A_i2 = A_j2, ...., A_ik = A_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快至少一个数量级。以下为代码实现，可以参考：

/* ================================================================== *//* NlogN for LCS, using LIS *//* ================================================================== */struct lcs_sort {    int val;    int pos;};struct lcs_lis_pos {    int pos_val;    int pos_1;};/* first,  sort string2 */static int qsort_func(const void *s1, const void *s2){    struct lcs_sort *ls1 = (struct lcs_sort *)s1, *ls2 = (struct lcs_sort *)s2;    if (ls1->val > ls2->val)         return 1;    else if (ls1->val < ls2->val)        return -1;    else {        if (ls1->pos > ls2->pos)            return -1;        else            return 1;    }}static void sort_lcs_string(char *string2, int n2, struct lcs_sort **ls_pp){    int i = 0;    struct lcs_sort *ls = NULL;        ls = malloc(sizeof(struct lcs_sort) * n2);    for (i = 0;i < n2;i++) {        ls[i].val = string2[i];        ls[i].pos = i;    }    qsort(ls, n2, sizeof(struct lcs_sort), qsort_func);    *ls_pp = ls;}/* second, find match list for char in string1 and construct a new sequence */static int find_match_list(char ch, struct lcs_sort *ls, int l, int r, int *start){    int m = 0, len, len2;    int start1 = 0, start2 = 0;    if (l > r)        return 0;    else if (l == r) {        if (ch == ls[l].val) {            *start = r;            return 1;        }        return 0;    }    m = (l + r)>>1;    if (ls[m].val < ch)        len = find_match_list(ch, ls, m + 1, r, &start1);    else if (ls[m].val > ch)        len = find_match_list(ch, ls, l, m - 1, &start1);    else {        len = find_match_list(ch, ls, l, m - 1, &start1);        len2 = find_match_list(ch, ls, m + 1, r, &start2) + 1;        if (len == 0)             start1 = m;        len += len2;    }    *start = start1;    return len;}static int match_list_lcs(char *string1, int n1, struct lcs_sort *ls, int n2, struct lcs_lis_pos **lpos_pp){    int i = 0, j = 0, k = 0, len = 0, start = 0;    struct lcs_lis_pos *lpos = NULL;    lpos = malloc(sizeof(struct lcs_lis_pos) * n1 * n2);    for (i = 0;i < n1;i++) {        len = find_match_list(string1[i], ls, 0, n2 - 1, &start);        for (k = 0;k < len;k++) {            lpos[j].pos_val = ls[start + k].pos;            lpos[j++].pos_1 = i;        }    }        *lpos_pp = lpos;    return j;}/* third, find LIS for this new sequence */static inline int find_pos(int *B, int len, int value){    int left = 0, right = len - 1, middle = 0;    while (left <= right) {        /* must finding strictly increasing sub sequence */        middle = (left + right)>>1;        if (B[middle] < value)            left = middle + 1;        else            right = middle - 1;    }        return left;}static int lis_nlogn(struct lcs_lis_pos *p, int len, int *inc_seq){    int i = 0, pos = 0, curr_len = 0;    int *L = NULL, *prev = NULL, *M = NULL;    if (len <= 0)        return 0;    L = malloc(len * sizeof(int));    M = malloc(len * sizeof(int));        prev = malloc(len * sizeof(int));    L[0] = p[0].pos_val;    M[0] = 0;    prev[0] = -1;               /* the prev of the p[0] is NULL */     curr_len = 1;    /* Caculate prev and M */    for (i = 1;i < len;i++) {        pos = find_pos(L, curr_len, p[i].pos_val);        L[pos] = p[i].pos_val;        M[pos] = i;        if (pos > 0)            prev[i] = M[pos - 1];        else            prev[i] = -1;        if (pos + 1 > curr_len)            curr_len++;    }    /* Output increasing sequence */    pos = M[curr_len - 1];    for (i = curr_len - 1;i >= 0 && pos != -1;i--) {        inc_seq[i] = p[pos].pos_1;        pos = prev[pos];    }    return curr_len;}static int calc_lcs_using_lis(char *string1, int n1, char *string2, int n2){    int *inc_seq = NULL;    int list_len = 0, inc_len = 0, len = (n1 < n2) ? n1 : n2;    struct lcs_sort *ls = NULL;    struct lcs_lis_pos *lpos = NULL;    inc_seq = (int *)malloc(sizeof(int) * len);    /* first,  sort string2 */    sort_lcs_string(string2, n2, &ls);    /* second, find match list for char in string1 and construct a new sequence */    list_len = match_list_lcs(string1, n1, ls, n2, &lpos);    /* third, find LIS for this new sequence */    inc_len = lis_nlogn(lpos, list_len, inc_seq);    /* end, get LCS */#if 0    int i = 0;    printf("LCS is: ");    for (i = 0;i < inc_len;i++) {        inc_seq[i] = string1[inc_seq[i]];        printf("%c ", inc_seq[i]);    }    printf("\r\n");#endif    free(ls);    free(lpos);        return inc_len;}