序列比对（24）最长公共子序列
本文介绍如何求解两个字符串的最长公共子序列。
最长公共子序列问题
前文《序列比对（23）最长公共子字符串》介绍了如何求解两个字符串的最长公共子字符串，本文将介绍如何求解两个字符串的最长公共子序列。二者听起来很像，所以我们首先得说明一下子字符串和子序列的区别。
与最长公共子字符串问题类似，最长公共子序列问题也是一种序列比对问题，可以用动态规划解决，只是在迭代时允许插入和缺失，而不允许错配而已。如果是匹配，得分为1，否则得分为0。其迭代公式如下：
效果如下：
动态规划求解最长公共子序列的代码
具体代码如下：
1#include <stdio.h>
2#include <stdlib.h>
3#include <string.h>
4#define MAXSEQ 1000
5#define max(a, b) ((a) > (b) ? (a) : (b))
6#define min(a, b) ((a) < (b) ? (a) : (b))
7
8struct Unit {
9 int W1; // 是否往上回溯一格
10 int W2; // 是否往左上回溯一格
11 int W3; // 是否往左回溯一格
12 int M; // 得分矩阵第(i, j)这个单元的分值，即序列s(1,...,i)与序列r(1,...,j)比对的最高得分
13};
14typedef struct Unit *pUnit;
15
16void strUpper(char *s);
17void printAlign(pUnit** a, const int i, const int j, char* s, char* r, int* is, int* ir, int n);
18void align(char *s, char *r);
19
20int main() {
21 char s[MAXSEQ];
22 char r[MAXSEQ];
23 printf("The 1st seq: ");
24 scanf("%s", s);
25 printf("The 2nd seq: ");
26 scanf("%s", r);
27 align(s, r);
28 return 0;
29}
30
31void strUpper(char *s) {
32 while (*s != '\0') {
33 if (*s >= 'a' && *s <= 'z') {
34 *s -= 32;
35 }
36 s++;
37 }
38}
39
40void printAlign(pUnit** a, const int i, const int j, char* s, char* r, int* is, int* ir, int n) {
41 int k;
42 pUnit p = a[i][j];
43 if (p->M == 0) { // 到值为0的矩阵单元就结束
44 printf("index of common seq in seq1: ");
45 for (k = n - 1; k >= 0; k--)
46 printf("%d ", is[k]);
47 printf("\n");
48 printf("index of common seq in seq2: ");
49 for (k = n - 1; k >= 0; k--)
50 printf("%d ", ir[k]);
51 printf("\n");
52 for (k = n - 1; k >= 0; k--)
53 printf("%c", s[is[k] - 1]);
54 printf("\n");
55 for (k = n - 1; k >= 0; k--)
56 printf("%c", r[ir[k] - 1]);
57 printf("\n\n");
58 return;
59 }
60 if (p->W1) // 向上回溯一格
61 printAlign(a, i - 1, j, s, r, is, ir, n);
62 if (p->W2) { // 向左上回溯一格
63 is[n] = i;
64 ir[n] = j;
65 printAlign(a, i - 1, j - 1, s, r, is, ir, n + 1);
66 }
67 if (p->W3) // 向左回溯一格
68 printAlign(a, i, j - 1, s, r, is, ir, n);
69}
70
71void align(char *s, char *r) {
72 int i, j;
73 int m = strlen(s);
74 int n = strlen(r);
75 int ml = min(m, n);
76 int m1, m2, m3, maxm;
77 pUnit **aUnit;
78 int* sidx;
79 int* ridx;
80 // 初始化
81 if ((aUnit = (pUnit **) malloc(sizeof(pUnit*) * (m + 1))) == NULL) {
82 fputs("Error: Out of space!\n", stderr);
83 exit(1);
84 }
85 for (i = 0; i <= m; i++) {
86 if ((aUnit[i] = (pUnit *) malloc(sizeof(pUnit) * (n + 1))) == NULL) {
87 fputs("Error: Out of space!\n", stderr);
88 exit(1);
89 }
90 for (j = 0; j <= n; j++) {
91 if ((aUnit[i][j] = (pUnit) malloc(sizeof(struct Unit))) == NULL) {
92 fputs("Error: Out of space!\n", stderr);
93 exit(1);
94 }
95 aUnit[i][j]->W1 = 0;
96 aUnit[i][j]->W2 = 0;
97 aUnit[i][j]->W3 = 0;
98 }
99 }
100 for (i = 0; i <= m; i++)
101 aUnit[i][0]->M = 0;
102 for (j = 1; j <= n; j++)
103 aUnit[0][j]->M = 0;
104 // 将字符串都变成大写
105 strUpper(s);
106 strUpper(r);
107 // 动态规划算法计算得分矩阵每个单元的分值
108 for (i = 1; i <= m; i++) {
109 for (j = 1; j <= n; j++) {
110 m1 = aUnit[i - 1][j]->M;
111 m2 = s[i - 1] == r[j - 1] ? aUnit[i - 1][j - 1]->M + 1 : -1;
112 m3 = aUnit[i][j - 1]->M;
113 maxm = max(max(m1, m2), m3);
114 aUnit[i][j]->M = maxm;
115 if (m1 == maxm) aUnit[i][j]->W1 = 1;
116 if (m2 == maxm) aUnit[i][j]->W2 = 1;
117 if (m3 == maxm) aUnit[i][j]->W3 = 1;
118 }
119 }
120/*
121 // 打印得分矩阵
122 for (i = 0; i <= m; i++) {
123 for (j = 0; j <= n; j++)
124 printf("%d ", aUnit[i][j]->M);
125 printf("\n");
126 }
127*/
128 printf("max score: %d\n", aUnit[m][n]->M);
129 // 打印最优比对结果，如果有多个，全部打印
130 // 递归法
131 if (aUnit[m][n]->M == 0) {
132 fputs("No common seq found.\n", stdout);
133 } else {
134 if ((sidx = (int*) malloc(sizeof(int) * ml)) == NULL || \
135 (ridx = (int*) malloc(sizeof(int) * ml)) == NULL) {
136 fputs("Error: Out of space!\n", stderr);
137 exit(1);
138 }
139 printAlign(aUnit, m, n, s, r, sidx, ridx, 0);
140 // 释放内存
141 free(sidx);
142 free(ridx);
143 }
144 for (i = 0; i <= m; i++) {
145 for (j = 0; j <= n; j++)
146 free(aUnit[i][j]);
147 free(aUnit[i]);
148 }
149 free(aUnit);
150}