位运算求最长公共子序列

最长公共子序列（LCS）问题

你有两个字符串 (A,B)，字符集为 (Sigma)，求 (A, B) 的最长公共子序列。

简单动态规划

首先有一个广为人知的 dp：(f_{i,j}) 为 (A) 的长度为 (j) 的前缀与 (B) 长度为 (i) 的前缀的 LCS。（注意 (i) 和 (j) 分别对于那个串）

那么显然有：

[f_{i,j} = egin{cases} f_{i-1, j-1} + 1 & (A_j = B_i) \ max(f_{i, j-1}, f_{i-1, j}) & (A_j e B_i) end{cases} ]

然而这是 (O(n^2)) 的，在略大的数据下就很容易 TLE。

还有一个 (O(nlog n)) 的算法，但只是针对排列的情况。

然后我们介绍一个基于 位运算 的优化方法。

这怎么就能位运算了呢？看着就不怎么 01。

但是有一个极其重要的性质：

[egin{cases} f_{i,j} ge f_{i-1,j} \ f_{i,j} ge f_{i,j-1} \ |f_{i,j} - f_{i,j-1}| le 1 end{cases} ]

即 (f) 的同一行内是 单调不减 并且 相邻两个相差不超过一。

矩阵 (M)

我们定义矩阵 (M) 为 (f) 数组每行分别 差分 的结果，即：

[f_{i,j} = sum_{k=1}^j M_{i,j} ]

根据上述 (f) 的性质，不难发现 (M) 是个 01矩阵。那么可以直接 压位（类似 std::bitset）。

然后考虑直接转移 (M_i) 整行，最后 (sum_{j}M_{|B|,j}) 就是答案。这就是优化的基本思想。

字符比较串 (p)

我们定义 (p(c)) 为字符 (c) 在字符串 (A) 中出现的所有位置的集合，(p(c)_i=1) 表示 (A_i=c)。这是我们转移的工具。

要预处理 (p) 我们需要 (O(|Sigma| imes |A|)) 的空间。然而我们发现 (p) 中只有 (0/1)，所以我们可以用类似于 (M) 进行压位优化，那就只要 (Oleft(frac{|Sigma| imes |A|}{w} ight))，一般来说还是一个可承受的量级。

(M) 的实际意义

上面只提到 (M) 是个差分数组，现在来考虑它的实际意义是什么，以便推出它的转移方式。

考虑一个 (M_{i,j}) 什么时候会是 (1)。观察原转移方程，发现 (f_{i,j-1}) 方向必然不会使 (f_{i,j}) 加一，唯一两个方向就是 (f_{i-1,j-1}) 或 (f_{i-1,j})。

如果是从 (f_{i-1,j-1}+1) 而来，那么说明这个位置 (A_j) 发生了配对，从而答案 (+1)；

如果是 (f_{i-1,j})，仔细思考一下还是一样的，在下面总有一个位置会和上面一条相同。

总而言之就是 (A_j) 被计入答案 了，但注意这不意味着 (M_i) 中所有的 (1) 都对应一个被选中的 (A_j)。

正确的理解是 (M_{i,j}) 如果为 (1)，设 (k) 为当前位到第一位之间 (1) 的个数，那就说明当前一个 LCS 长度为 (k) 的方案，最后的一位为 (j)。事实我们也是只需要考虑当前 LCS 的最后一位，添加时答案只要保证在当前方案的最后一位之后即可。

转移方式

对于一整行 (M_{i-1})，我们对其分段，每段有前面一个极长 (0) 段，由一个单独的 (1) 结尾，最后一整段 (0) 单独成段。

然后用当前 (B) 的字符 (p(B_i)) 与之比对（注意这里是倒着的）：

M[i - 1]: [1  0  0  0  0  0  0][1  0  0  0][1][1][1][1][1]
p[B[i]] : [0  1  0  1  1  0  0  0  1  0  0  0  1  1  0  0]
           ^                                            ^
           |                                            |
         j = |A|                                      j = 1

然后将两者做 按位或 操作，再对于每个段按位或的结果取 段中的最后一个 (1)，得：

M[i] :    [0  0  0  0  1  0  0  0  1  0  0  1  1  1  1  1]

这个过程相当于 (M_{i-1}) 借助 (p(B_i)) 将这些 (1) 尽量向字符串的开头移，以便为之后的匹配留足更大的机会。

至于其中的意义可以结合上面理解，大概就是对于每个长度的方案，都在不超过下一个长度的前提下前移。具体细节我也说不清楚

转移实现

上面的转移过于复杂，很难用我们熟知的位运算进行优化，于是尝试将它翻译成位运算。

我也不知道原论文作者怎么想到的，这里就说只一下做法吧。

我们记 (X = M_{i-1} exttt{OR} p(B_i))，然后我们需要取其中最后一位：

X :       [1  1  0  1  1  0  0  1  1  0  0  1  1  1  1  1]

然后将 (M_{i-1}) 右移一位，头部补上 (1)，并用 (X) 数值减 这个 01 串，得：

          [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]

这么做旨在将每段的末尾 (0) 段，然后将原来最右边的 (1) 变成 (0)。

然后和 (X) 进行 异或 操作：

          [0  0  0  0  1  1  1][0  1  1  1][1][1][1][1][1]

这样就使最开始的最右边的 (1) 到段尾变成 (1)，其余变成 (0)。

最后只要保留第一个 (1)，那么就刚好是 按位与 (X) 的结果。

于是得到：

[M_i = ((X-((M_{i-1} exttt{<<} 1) + 1)) exttt{xor} X) exttt{and} X、 ]

那么在实现时，只要手写一个 bitset，支持按位与、或、异或、数值相减、位移即可。

复杂度

每次转移需要 (Oleft(frac {|A|} w ight))，总时间复杂度为 (Oleft(frac{|A| imes |B|}{w} ight))

空间瓶颈为 (p) 集合，为 (Oleft(frac{|A| imes |Sigma|}{w} ight))，如果字符集 (Sigma) 不确定可以离散化，空间为 (Oleft( frac{|A|^2}{w} ight))。

参考代码

下面的代码实现 并不是倒着的（为了减法方便），于是位移什么的看着就有点诡异。

LOJ 提交入口

/*
 * Author : _Wallace_
 * Source : https://www.cnblogs.com/-Wallace-/
 * Problem : LOJ #6564. 最长公共子序列
 * Standard : GNU C++ 03
 * Optimal : -Ofast
 */
#include <algorithm>
#include <cstddef>
#include <cstdio>
#include <cstring>

typedef unsigned long long ULL;

const int N = 7e4 + 5;
int n, m, u;

struct bitset {
  ULL t[N / 64 + 5];

  bitset() {
    memset(t, 0, sizeof(t));
  }
  bitset(const bitset &rhs) {
    memcpy(t, rhs.t, sizeof(t));
  }

  bitset& set(int p) {
    t[p >> 6] |= 1llu << (p & 63);
    return *this;
  }
  bitset& shift() {
    ULL last = 0llu;
    for (int i = 0; i < u; i++) {
      ULL cur = t[i] >> 63;
      (t[i] <<= 1) |= last, last = cur;
    }
    return *this;
  }
  int count() {
    int ret = 0;
    for (int i = 0; i < u; i++)
      ret += __builtin_popcountll(t[i]);
    return ret;
  }

  bitset& operator = (const bitset &rhs) {
    memcpy(t, rhs.t, sizeof(t));
    return *this;
  }
  bitset& operator &= (const bitset &rhs) {
    for (int i = 0; i < u; i++) t[i] &= rhs.t[i];
    return *this;
  }
  bitset& operator |= (const bitset &rhs) {
    for (int i = 0; i < u; i++) t[i] |= rhs.t[i];
    return *this;
  }
  bitset& operator ^= (const bitset &rhs) {
    for (int i = 0; i < u; i++) t[i] ^= rhs.t[i];
    return *this;
  }

  friend bitset operator - (const bitset &lhs, const bitset &rhs) {
    ULL last = 0llu; bitset ret;
    for (int i = 0; i < u; i++){
      ULL cur = (lhs.t[i] < rhs.t[i] + last);
      ret.t[i] = lhs.t[i] - rhs.t[i] - last;
      last = cur;
    }
    return ret;
  }
} p[N], f, g;

signed main() {
  scanf("%d%d", &n, &m), u = n / 64 + 1;
  for (int i = 1, c; i <= n; i++)
    scanf("%d", &c), p[c].set(i);
  for (int i = 1, c; i <= m; i++) {
    scanf("%d", &c), (g = f) |= p[c];
    f.shift(), f.set(0);
    ((f = g - f) ^= g) &= g;
  }
  printf("%d
", f.count());
  return 0;
}

后记

• 原文地址：https://www.cnblogs.com/-Wallace-/p/14052964.html
• 本文作者：@-Wallace-
• 转载请附上出处。