LCS のワードサイズ高速化

概要

[1] で紹介されている LCS のワードサイズ高速化について、自分の理解を記します。特に新しい主張はないです。

内容は [2] を大きく参考にしています。

dp 配列とその差分

$\mathrm{dp} _ {i, j}$ を二つの文字列 $S _ {1 .. i}, T _ {1 .. j}$ の LCS の長さ、と定義します。次が成り立ちます。

\begin{align*} \mathrm{dp}_{i, j} = \begin{cases} 0 & (i = 0 \ \mathrm{or}\ j = 0) \\ \mathrm{dp}_{i-1, j-1} + 1 & (S_i = T_j) \\ \max( \mathrm{dp}_{i-1, j}, \mathrm{dp}_{i, j-1} ) & (\mathrm{Otherwise}) \end{cases} \end{align*}

$j$ についての差分 $D _ {i, j} := \mathrm{dp} _ {i, j} - \mathrm{dp} _ {i, j-1}$ を考えます。 $D _ {i, j} \in \lbrace 0, 1 \rbrace$ なので $D$ は bit 列として表現できます。

遷移の観察

$D _ {i-1}$ から $D _ {i}$ への遷移を考えます。

便宜上、 $D _ {i-1}$ の末尾に $1$ を、 $D _ {i}$ の末尾に $0$ を追加しています。

また、 $D _ {i-1, j} = 1$ なる添字列 $j _ 0 \lt j _ 1 \lt j _ 2 \lt \dots \lt j _ k (= |T|+1)$ について、 $\lbrack 1, j _ 0 \rbrack, \lbrack j _ 0 + 1, j _ 1 \rbrack, \lbrack j _ 1 + 1, j _ 2 \rbrack, \dots$ に分割して、二重線で囲んでいます。

各区間について、 $1$ の位置の変化に注目しましょう。同じ場所にあるか、左側に移動していることがわかります。移動しているのは

$D _ {i-1, j} = 1$ の場所
$S _ {i} = T _ {j}$ の場所

のうち、最も左にある場所です。

自然な解釈

この $1$ の移動は $D$ の意味を考えるとより自然かもしれません。 $D _ {i, j ^ { * } } = 1$ なる添字 $j ^ { * }$ を一つ取ります。このとき $\sum _ {j \leq j ^ { * } } D _ {i, j}$ は $\mathrm{dp} _ {i, j ^ *}$ (つまり $S _ {1 .. i}$ と $T _ {1 .. j ^ { * } }$ の LCS の長さ) に等しいです。

これは逆に LCS の長さが初めて $x$ になるような添字 $j$ の列を $D$ が管理していると考えることができます。文字 $\alpha$ が追加されると、各添字は (前の添字を追い越さないように) 左側に進んでいきます。番兵として $D$ の末尾に $1$ を追加しましたが、これが左側に移動することが、全体の LCS の長さが $1$ だけ長くなることに対応しています。

ビット演算

各区間における最も左の $1$ の場所に $1$ を立てたいです。

ここから先は整数に対する演算を考えるので、左右を反転して考えることにします。

D : 現在の配列。各区間の左端には 1 が立っている。
m : 新しく追加する文字と一致するところに 1 が立っている
各区間においてD|mの最も右の 1 の場所に 1 があるようにしたい

2 種類の計算式を説明します。

`D = (D | m) ^ ((D | m) & ((D | m) - (D << 1 | 1)))` [2]

右隣の区間の左端には 1 が立っているので、(D << 1 | 1)で各区間の右端に 1 がある状態になります（最も右の区間だけ 1 を追加する必要がある）。
(D | m) - x で (D | m) の最下位ビットが 0 に、それ以下のビットが 1 になったものが手に入ります。
(D | m) & x で (D | m) の最下位ビットだけが 0 になります。
(D | m) ^ x で (D | m) の最下位ビットだけが 1 になります。

                       D                       --- 1000000
                           m                   --- 0110100
                       D | m                   --- 1110100
                                 D << 1 | 1    --- 0000001
                      (D | m) - (D << 1 | 1)   --- 1110011
           (D | m) & ((D | m) - (D << 1 | 1))  --- 1110000
(D | m) ^ ((D | m) & ((D | m) - (D << 1 | 1))) --- 0000100

`D = (D | m) & (D - (m & ~D))` [1]

[1] では [2] で紹介されたビット演算を式変形によって簡略化しました。この演算を直接理解することもできます。

一つの区間で考えます。m が 0 のときは D はそのままです。

m が 0 でないとき、

m & ~D は m の中で D に含まれない 1 が立ちます。
D - (m & ~D) は、m を右から見ていって、最下位ビットから〜次に 1 が立つ場所の手前まで 1 が連続します。
(D | m) & (D - (m & ~D)) は m の最下位ビットだけが残ります。

 D                       --- 100000
     m                   --- 111010
 D | m                   --- 111010
                m & ~D   --- 011010
           D - (m & ~D)  --- 000110
(D | m) & (D - (m & ~D)) --- 000010

終わりに

実装例

#include <bits/stdc++.h>
using namespace std;

struct simple_bitset {
    int n;
    vector<unsigned long long> data;
    simple_bitset(int n_) : n((n_ + 64 - 1) / 64), data(n, 0) {} 

    using bs = simple_bitset;
    void set_one(const int i) {
        data[i / 64] |= 1ULL << (i % 64);
    }
    bs& operator&=(const bs& r) {
        for(int i = 0; i < n; i++) data[i] &= r.data[i];
        return *this;
    }
    bs& operator|=(const bs& r) {
        for(int i = 0; i < n; i++) data[i] |= r.data[i];
        return *this;
    }
    bs& operator-=(const bs& r) {
        for(int i = 0, c = 0, nc = 0; i < n; i++) {
            nc = data[i] < r.data[i] || (data[i] == r.data[i] && c);
            data[i] -= r.data[i] + c;
            c = nc;
        }
        return *this;
    }
    bs operator&(const bs& r) const { return bs(*this) &= r; }
    bs operator|(const bs& r) const { return bs(*this) |= r; }
    bs operator-(const bs& r) const { return bs(*this) -= r; }
    bs operator~() const {
        bs x = *this;
        for(int i = 0; i < n; i++) x.data[i] = ~x.data[i];
        return x;
    }
    int count() const {
        int cnt = 0;
        for(int i = 0; i < n; i++) cnt += __builtin_popcountll(data[i]);
        return cnt;
    }
};

int lcs(const string& s, const string& t) {
    const int n = s.size(), m = t.size();
    vector M(26, simple_bitset(m));
    for(int j = 0; j < m; j++) M[t[j] - 'a'].set_one(j);

    simple_bitset D(m);
    for(int i = 0; i < n; i++) {
        const simple_bitset& Ms = M[s[i] - 'a'];
        D = (D | Ms) & (D - (Ms & ~D));
    }
    return D.count();
}

int main() {
    string s, t; cin >> s >> t;
    cout << lcs(s, t) << endl;
}