# Tối ưu thuật toán tìm xâu con chung dài nhất bằng bitset

**Người viết:**
> **Nguyễn Tường Duy - VNU University of Engineering and Technology**

## Bài toán

Xâu $S$ là xâu con của $T$ khi ta có thể xóa đi một số kí tự của $T$ và không thay đổi thứ tự của các kí tự còn lại để được xâu $S$.

Cho hai xâu $A$ độ dài $n$ và xâu $B$ độ dài $m$ **chỉ gồm những chữ cái in thường**, tìm xâu con chung dài nhất của $A$ và $B$.

## Tìm độ dài xâu con chung dài nhất

### 1. Công thức quy hoạch động phổ biến

Gọi $dp(i, j)$ là độ dài xâu con chung dài nhất của tiền tố độ dài $i$ của xâu $A$ và tiền tố độ dài $j$ của xâu $B$, ta có công thức quy hoạch động sau:

* Nếu $i = 0$ hoặc $j = 0$ thì $dp(i, j) = 0$.
* Nếu $A_i = B_j$ thì $dp(i, j) = dp(i - 1, j - 1) + 1$.
* Nếu $A_i \neq B_j$ thì $dp(i, j) = max(dp(i - 1, j), dp(i, j - 1))$.

Như vậy ta có thể giải bài toán trên trong thời gian $O(n \cdot m)$ và bộ nhớ $O(n \cdot m)$.

### 2. Ma trận bit

Xét ma trận $M$ với $M(i, j) = dp(i, j) - dp(i, j - 1)$, khi đó $dp(i, j) = \sum_{k = 1}^{j} M(i, k)$.

Vì $0 \leq M(i, j) \leq 1$ nên $M(i, j)$ là $1$ bit. Từ đó ta có ý tưởng lưu $W$ bit vào một số nguyên để tăng tốc độ tính toán.

### 3. Quy hoạch động trên ma trận bit

Ta có cách quy hoạch động như sau:

M[0] <- 00...0
Duyệt i từ 1 đến |A|:
    M[i] <- M[i - 1]
    j <- 1
    Lặp lại nếu j <= |B|:
        Nếu B[j] = A[i]:
            Xét hàng M[i], nếu không có bit 1 nào bên phải j thì:
                M[i][j] <- 1
                Thoát khỏi vòng lặp
            Nếu không thì:
                k <- vị trí bit 1 gần nhất bên phải j (k >= j)
                M[i][k] <- 0
                M[i][j] <- 1
                j <- k + 1
        Nếu không thì:
            j <- j + 1

### 4. Tối ưu bằng các phép toán trên bit

Theo như cách quy hoạch động ở trên, khi xét vị trí $j$ mà $B_j = A_i$, gọi $k$ là vị trí bit $1$ gần nhất bên phải $j$, ta cần chuyển bit ở vị trí $k$ thành $0$, chuyển bit ở vị trí $j$ thành $1$ và với các vị trí $l$ ($j < l < k$) mà $B_l = A_i$ thì ta giữ nguyên bit ở vị trí đó bằng $0$. Để thực hiện được việc này thì ta dùng phép trừ.

Gọi $dp$ là hàng thứ $i - 1$, $T$ là bitset gồm các vị trí mà $B_j = A_i$. Ban đầu ta có công thức tính hàng $i$ là $dp - T$.

Nhưng nếu có vị trí mà bit của $dp$ và $T$ đều là $1$ thì sẽ xóa thừa bit $1$ nên ta phải bỏ các bit chung của $dp$ và $T$. Như vậy ta có công thức $dp - (dp \oplus (dp \lor T))$.

Nhưng khi đó ta lại tính thừa bit $1$. Để giải quyết vấn đề này thì ta dùng phép $AND$ là xong: $(dp - (dp \oplus (dp \lor T))) \land (dp \lor T)$. Ta tối ưu tiếp công thức trên bằng cách bỏ những phép toán tính đi tính lại nhiều lần. Đặt $x = dp \lor T$, công thức là $(dp - (dp \oplus x)) \land x$.

Như vậy ta đã tối ưu được thời gian xuống $O(\frac{n \cdot m}{W})$.

### 5 . Tổng kết

Ta được thuật toán độ phức tạp thời gian là $O(\frac{n \cdot m}{W})$, độ phức tạp bộ nhớ là $O(\frac{n}{W})$.

::: spoiler Cài đặt

struct custom_bitset {
    int n;
    vector<long long> a;
    custom_bitset(int m): a(n = m / 63 + 1) {}
    void set(int i) {
        a[i / 63] |= 1LL << (i % 63);
    }
    int count() const {
        int res = 0;
        for (int i = 0; i < n; ++i) res += __builtin_popcountll(a[i]);
        return res;
    }
    custom_bitset operator & (const custom_bitset &o) const {
        custom_bitset res = *this;
        for (int i = 0; i < n; ++i) res.a[i] &= o.a[i];
        return res;
    }
    custom_bitset operator | (const custom_bitset &o) const {
        custom_bitset res = *this;
        for (int i = 0; i < n; ++i) res.a[i] |= o.a[i];
        return res;
    }
    custom_bitset operator ^ (const custom_bitset &o) const {
        custom_bitset res = *this;
        for (int i = 0; i < n; ++i) res.a[i] ^= o.a[i];
        return res;
    }
    custom_bitset operator - (const custom_bitset &o) const {
        custom_bitset res = *this;
        for (int i = 0, sub = 0; i < n; ++i) {
            res.a[i] -= o.a[i] + sub;
            sub = (res.a[i] < 0);
            res.a[i] &= LLONG_MAX;
        }
        return res;
    }
};

int lcs_length(const string &a, const string &b) {
    int n = a.size(), m = b.size();
    custom_bitset dp(m);
    vector<custom_bitset> T(26, custom_bitset(m));
    for (int i = 0; i < m; ++i) T[b[i] - 'a'].set(i);
    for (int i = 0; i < n; ++i) {
        custom_bitset x = dp | T[a[i] - 'a'];
        dp = (dp - (dp ^ x)) & x;
    }
    return dp.count();
}

:::

## Tìm xâu con chung dài nhất

### 1. Thuật toán ngây thơ

Nếu có đủ bộ nhớ để lưu ma trận $M$ thì ta có thể tìm xâu con chung dài nhất trong thời gian $O(\frac{n \cdot m}{W})$ và bộ nhớ $O(\frac{n \cdot m}{W})$.

Nhưng với $n$ và $m$ lớn (khoảng $10^5$) thì ta cần tới trên $1GB$ để lưu ma trận $M$. Giới hạn bộ nhớ trong CP thường là $1GB$ nên ta sẽ bị vượt quá giới hạn bộ nhớ.

### 2. Thuật toán tối ưu bộ nhớ

Ta sẽ chia bài toán thành các bài toán con. Ta chia xâu $A$ ra thành hai nửa bằng nhau, gọi là $A_1$, $A_2$. Ký hiệu $LCS(A, B)$ là độ dài xâu con chung dài nhất của $A$ và $B$, ta cần chia xâu $B$ ra thành hai phần $B_1$, $B_2$ sao cho $LCS(A1, B1) + LCS(A2, B2) = LCS(A, B)$.

Ký hiệu $A'$ là mảng đảo ngược của $A$, ta sẽ tìm độ dài xâu con chung dài nhất tất cả các tiền tố của $A_1$ và $B$, tất cả các tiền tố của $A_2'$ và $B'$, gọi là $dp1$, $dp2$. Khi đó vị trí để chia xâu $B$ ra thành hai phần là vị trí $i$ mà $dp1_i + dp2'_{i + 1}$ lớn nhất. Lưu ý các trường hợp nhỏ sau: xâu $B$ rỗng, xâu $A$ có độ dài là $1$.

Rõ ràng độ phức tạp bộ nhớ của thuật toán trên là $O(n + m)$, ta sẽ phân tích độ phức tạp thời gian. Phần tốn nhiều thời gian nhất của hàm đệ quy là tìm độ dài xâu con chung dài nhất, độ phức tạp là O($\frac{n \cdot m}{W}$). Xem xét cây đệ quy, ta thấy số lượng phép tính ở các đỉnh độ sâu $1$ là $\frac{n \cdot m}{W}$, độ sâu $2$ là $\frac{\frac{n}{2} \cdot m_1}{W} + \frac{\frac{n}{2} \cdot m_2}{W} = \frac{\frac{n}{2} \cdot (m_1 + m_2)}{W} = \frac{\frac{n}{2} \cdot m}{W}$, độ sâu $3$ là $\frac{\frac{n}{4} \cdot m}{W}$, $\dots$ Tổng số phép tính là $\frac{(n +\frac{n}{2} + \frac{n}{4} + \dots) \cdot m}{W} \approx \frac{n \cdot m}{W}$. Như vậy độ phức tạp thời gian là $O(\frac{n \cdot m}{W})$, giữ nguyên so với thuật toán ngây thơ.

### 3. Tổng kết

Ta được thuật toán độ phức tạp thời gian là $O(\frac{n \cdot m}{W})$, độ phức tạp bộ nhớ là $O(n + m)$.

::: spoiler Cài đặt

struct custom_bitset {
    int n;
    vector<long long> a;
    custom_bitset(int m): a(n = m / 63 + 1) {}
    void set(int i) {
        a[i / 63] |= 1LL << (i % 63);
    }
    bool test(int i) const {
        return (a[i / 63] >> (i % 63)) & 1;
    }
    int count() const {
        int res = 0;
        for (int i = 0; i < n; ++i) res += __builtin_popcountll(a[i]);
        return res;
    }
    custom_bitset operator & (const custom_bitset &o) const {
        custom_bitset res = *this;
        for (int i = 0; i < n; ++i) res.a[i] &= o.a[i];
        return res;
    }
    custom_bitset operator | (const custom_bitset &o) const {
        custom_bitset res = *this;
        for (int i = 0; i < n; ++i) res.a[i] |= o.a[i];
        return res;
    }
    custom_bitset operator ^ (const custom_bitset &o) const {
        custom_bitset res = *this;
        for (int i = 0; i < n; ++i) res.a[i] ^= o.a[i];
        return res;
    }
    custom_bitset operator - (const custom_bitset &o) const {
        custom_bitset res = *this;
        for (int i = 0, sub = 0; i < n; ++i) {
            res.a[i] -= o.a[i] + sub;
            sub = (res.a[i] < 0);
            res.a[i] &= LLONG_MAX;
        }
        return res;
    }
};

vector<int> get_dp_state(string a, string b, bool rev) {
    if (rev) {
        reverse(a.begin(), a.end());
        reverse(b.begin(), b.end());
    }
    int n = a.size(), m = b.size();
    custom_bitset dp(m);
    vector<custom_bitset> T(26, custom_bitset(m));
    for (int i = 0; i < m; ++i) T[b[i] - 'a'].set(i);
    for (int i = 0; i < n; ++i) {
        custom_bitset x = dp | T[a[i] - 'a'];
        dp = (dp - (dp ^ x)) & x;
    }
    vector<int> dp_state(m + 2, 0);
    for (int i = 0; i < m; ++i) dp_state[i + 1] = dp_state[i] + dp.test(i);
    if (rev) reverse(dp_state.begin(), dp_state.end());
    return dp_state;
}

string lcs(const string &a, const string &b) {
    int n = a.size(), m = b.size();
    if (!m) return "";
    if (n == 1) {
        for (auto &i: b) if (i == a[0]) return string(1, i);
        return "";
    }
    string a1 = a.substr(0, n / 2), a2 = a.substr(n / 2);
    vector<int> dp1 = get_dp_state(a1, b, 0), dp2 = get_dp_state(a2, b, 1);
    int ma = -1, j;
    for (int i = 0; i <= m; ++i) {
        if (ma < dp1[i] + dp2[i + 1]) {
            ma = dp1[i] + dp2[i + 1];
            j = i;
        }
    }
    return lcs(a1, b.substr(0, j)) + lcs(a2, b.substr(j));
}

:::

## Thuật toán có ứng dụng gì không?

Bạn không nên hỏi điều đó, đáng ra bạn nên nói "Thuật toán trên thật tuyệt!" hoặc "Blog rất thú vị!".

> **Thuật toán không có ứng dụng nào :sob:**

## Tham khảo

**A bit-string longest common subsequence
algorithm (Lloyd ALLISON, Trevor I. DIX)**

**A linear space algorithm for computing maximal common subsequences (D.S. Hirschberg)**