基于位运算的最长公共子串算法

本文章来源

[常规动态规划算法]

设L[i , j]等于A[1..i] , B[1..j]的LCS.

则有L[i,j] = 1 + L[i-1 , j-1]                    如果 (A[i] = B[j])

              Max(L[i-1 , j] , L[i , j-1])   其他

复杂度为O(|A| * |B|)

[基于位运算的动态规划算法]

根据上面的动态规划算法，状态函数L具有如下性质：

L[i-1,j-1] ≤ L[i,j-1] , L[i-1,j] ≤ L[i,j]

  | L[i,j]-L[i-1,j-1] | ≤ 1

  对于L的每一行，相邻的两个元素的最多只相差1。这样一来，我们就可以用一个二进制的矩阵描述出L :    

#bits

 9    0 0 0 1 1 0 0 0 1 0 1 1 1 1 1 1 | T   string B

 9    0 1 0 0 1 0 0 0 1 0 1 1 1 1 1 1 | T 

 8    0 0 1 0 0 0 1 0 0 0 1 1 1 1 1 1 | C

 7    0 0 0 0 1 0 0 0 1 0 0 1 1 1 1 1 | T   - Row[11]

 7    1 0 0 0 0 0 0 1 0 0 0 1 1 1 1 1 | A   - Row[10]

 7    1 0 0 0 0 1 0 0 0 0 0 1 1 1 1 1 | G

 6    0 0 0 0 0 1 0 1 0 0 0 0 1 1 1 1 | A

 5    0 0 0 0 0 0 0 1 0 0 0 0 1 1 1 1 | A

 5    0 0 0 0 1 0 0 0 0 0 0 0 1 1 1 1 | T

 4    0 0 0 0 0 0 0 0 1 0 0 0 0 1 1 1 | T

 3    0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 1 | C

 3    1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 | G

 2    0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 | A

 1    0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 | T

     _________________________________.        matrix M_ij

String A:  G T C T T A C A T C C G T T C G 

 

这里，我们将串A从右往左写，串B从下往上写。Row[i]中的1的个数总是和Row[i-1]中的1的个数一样多或者恰好多一个。串A和串B的LCS即为最上面一行Row[|B|]中1的个数。

字符比较串表

这里我们定义一组称为字符比较串的二进制串。分别是字符集中的每一个字符与串A的比较结果(相同为1，不同为0)。

A： G T C T T A C A T C C G T T C G

‘A’-string: 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0

‘C’-string: 0 0 1 0 0 0 1 0 0 1 1 0 0 0 1 0

‘G’-string: 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1

‘T’-string: 0 1 0 1 1 0 0 0 1 0 0 0 1 1 0 0

预先计算这个字符比较串表，时间复杂度为O(|S|*|A|)。对于一个确定的字符集，时间复杂度为O(|A|)。对于一个不确定的字符集，最坏情况为O(|A|*|B|)。如果字符集较小，|S|<<|B|，预处理的时间复杂度可以被忽略。

矩阵M

为了计算Row[i]，我们需要用到字符比较串表中的B[i]-string。以Row[10]到Row[11]为例，我们来研究如何计算出Row[i]。

下面是Row[10] , 以及B[10]-string(‘T’-string). 按照Row[10]中的1分隔：

Row[10]: 1 0 0 0 0 0 0   1 0 0 0   1   1   1   1   1

T-string : 0 1 0 1 1 0 0   0 1 0 0   0   1   1   0   0

每一段都是从Row[i-1]的一个1的位开始往右延伸，直到下一个位置是1或者串结束。如果Row[i-1]的最左边的位置上是0，那么最左边的一段从B[i]-string的最左边的1的位置开始延伸直到下一个1。Row[i]的构成方式很简单：就是对于每一段，都是选择Row[i-1]或者B[i]-string最右边的1所在位置为1，其他的为0。如果这一段Row[i-1]和B[i]-string都是0，那么Row[i]这一段也为0。

‘T’-string Or Row[10]:

  1 1 0 1 1*0 0    1 1*0 0   1*  1* 1*  1*  1*

Row[11]: 0 0 0 0 1 0 0    0 1 0 0   1   1   1   1   1

*表示了进行或操作之后每一段最右边的1

附带一提，你可以假定在每个串的最左边（位置|A|+1）存在一个1，这样可以方便处理最左边一段全为0的情况。不过对于本算法并没有这个必要。

在Row[i-1]中的一个1的位置，代表了A中的一个最短前缀与B[1..i-1]的LCS达到了该长度。引进B[i]，最好的方法当然是在前面的最短基础上加入一个最短的与B[i]的匹配，也就是从那个1所在的位置往左找，找到的第一个能与B[i]匹配的(如果找得到的话)。

设 X = Row[i-1] Or B[i]-string

对于当前的例子，X = Row[10] Or ‘T’-string

    X:  1 1 0 1 1 0 0    1 1 0 0   1   1   1   1   1

把Row[i-1]往左移1位，并且设最低位为1，用X去减这个串：

X:         1 1 0 1 1 0 0    1 1 0 0   1   1   1   1   1

         - 0 0 0 0 0 0 1    0 0 0 1   1   1   1   1   1

           -------------    -------   -   -   -   -   -

           1 1 0 1 0 1 1    1 0 1 1   0   0   0   0   0

这个操作所表达的意思为对于每一段，把最右边的一个1变成0，把这个1右边的所有0变成1，设最低位为1的目的是为了可以同样处理最后一段。

再跟X 进行Xor(异或)操作，得到：

0 0 0 0 1 1 1    0 1 1 1   1   1   1   1   1

这步操作之后的结果就是，对于每一段，从最右边的1开始到段尾变成1，其他变成0。因为Xor为1当且仅当两个位上的数不同，而上一次操作我们修改过的位是从最右边的1直到段尾。

接下来就很显然了，只要再将上面的结果与X进行And操作即可。得到：

0 0 0 0 1 0 0    0 1 0 0   1   1   1   1   1

这样就得到了Row[i] , 即当前例子中的Row[11]。

综上，

Row[i] = X And ( (X – ((Row[i-1] << 1) + 1)) Xor X),   

其中X = Row[i-1] Or B[i]-string

算法到此结束。计算出Row[|B|]之后，数一下其中有多少个1的位，便是答案了。