位运算求最长公共子序列

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

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

简单动态规划

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

那么显然有：

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

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

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

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

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

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

\[\begin{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,k} \]

根据上述 \(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|\times |A|)\) 的空间。然而我们发现 \(p\) 中只有 \(0/1\)，所以我们可以用类似于 \(M\) 进行压位优化，那就只要 \(O\left(\frac{|\Sigma|\times |A|}{w}\right)\)，一般来说还是一个可承受的量级。

\(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}\ \texttt{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}\ \texttt{<<}\ 1) + 1))\ \texttt{xor}\ X)\ \texttt{and}\ X、 \]

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

复杂度

每次转移需要 \(O\left(\frac {|A|} w\right)\)，总时间复杂度为 \(O\left(\frac{|A|\times |B|}{w}\right)\)

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

参考代码

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

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\n", f.count());
  return 0;
}