LCS を時間 O(|S| |T| / w)、空間 O(|S| + |T|) で復元までやる

最長共通部分列 (Longest Common Subsequence; LCS) を時間計算量 $\Theta(|S| |T| / w)$、空間計算量 $\Theta(|S| + |T|)$ で復元までします（$w$ はワードサイズ）。

韓国語の記事はいくつか見つかるんですが日本語の記事は見つからないので韓国高度典型なのかも（なまじ音が近いせいか韓国語 → 日本語の機械翻訳精度が低いことが多くしんどい……）。

• 諸注意
  • この記事では「復元なし」「復元あり」とはそれぞれ以下のような問題設定を指します。
    • 復元なし：最長共通部分列の長さを求める
    • 復元あり：最長共通部分列（のうちどれかひとつ）を求める
  • 説明のためよく知られるアルゴリズムから順に触れていくので、適宜必要なところだけを読んでください。
  • 完全に理解しているわけではないので厳密な証明などは求めないでください。あくまで紹介記事と捉えていただけると嬉しいです。
  • 間違い（誤字等含む）を見つけたら気軽に教えてください。読まれていることが実感できて嬉しいまである。

---

---

時間 $\Theta(|S| |T|)$、空間 $\Theta(|S| |T|)$、復元なし

ジャッジ

• https://www.acmicpc.net/problem/9251
  • $|S|, |T| \le 1000$、復元なし

概要

詳しく知りたい方は以下の記事などを参考にしてください。

• 問題 8：最長共通部分列 (LCS) 問題 | アルゴ式

コード

int lcs(const string &s, const string &t) {
  int n = s.size(), m = t.size();
  vector<vector<int>> dp(n + 1, vector<int>(m + 1));
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < m; ++j) {
      dp[i + 1][j + 1] = max({ dp[i + 1][j + 1], dp[i + 1][j], dp[i][j + 1] });
      if (s[i] == t[j]) {
        dp[i + 1][j + 1] = max(dp[i + 1][j + 1], dp[i][j] + 1);
      }
    }
  }
  return dp[n][m];
}

---

時間 $\Theta(|S| |T|)$、空間 $\Theta(|S| |T|)$、復元あり

ジャッジ

• https://www.acmicpc.net/problem/9252
  • $|S|, |T| \le 1000$、復元あり
• https://atcoder.jp/contests/dp/tasks/dp_f
  • $|S|, |T| \le 3000$、復元あり

概要

復元なしの方で作ったテーブルを逆順に辿るか、テーブル作成時にどこから遷移したかをメモすることなどで求まります。

コード

string lcs(const string &s, const string &t) {
  int n = s.size(), m = t.size();
  vector<vector<int>> dp(n + 1, vector<int>(m + 1));
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < m; ++j) {
      dp[i + 1][j + 1] = max({ dp[i + 1][j + 1], dp[i + 1][j], dp[i][j + 1] });
      if (s[i] == t[j]) {
        dp[i + 1][j + 1] = max(dp[i + 1][j + 1], dp[i][j] + 1);
      }
    }
  }
  string res = "";
  int i = n, j = m;
  while (i > 0 && j > 0) {
    if (s[i - 1] == t[j - 1]) {
      res += s[i - 1];
      --i;
      --j;
    }
    else if (dp[i][j] == dp[i - 1][j]) {
      --i;
    }
    else {
      --j;
    }
  }
  reverse(res.begin(), res.end());
  return res;
}

---

時間 $\Theta(|S| |T|)$、空間 $\Theta(|S| + |T|)$、復元なし

ジャッジ

• なし
  • 時間 $\Theta(|S| |T|)$、復元なしで空間 $\Theta(|S| + |T|)$ と $\Theta(|S| |T|)$ を区別する問題は見つかりませんでした。あれば教えてください。

概要

テーブルを作成するときに一行前の結果しか使っていないので全部持つ必要はありません。
後々の都合上、s に対するループと t に対するループの順番を逆にしています。

コード

int lcs(const string &s, const string &t) {
  int n = s.size(), m = t.size();
  vector<int> dp(n + 1);
  for (int i = 0; i < m; ++i) {
    vector<int> nx = dp;
    for (int j = 0; j < n; ++j) {
      nx[j + 1] = max(nx[j + 1], nx[j]);
      if (t[i] == s[j]) {
        nx[j + 1] = max(nx[j + 1], dp[j] + 1);
      }
    }
    dp = nx;
  }
  return dp[n];
}

---

時間 $\Theta(|S| |T|)$、空間 $\Theta(|S| + |T|)$、復元あり

ジャッジ

• https://www.acmicpc.net/problem/18438
  • $|S|, |T| \le 7000$、復元あり、ML: 4 MB

概要

Hirschberg's algorithm というものが存在します。
以下がかなりわかりやすいです。

• Hirschberg's Algorithmについて - Learning Algorithms

お気持ちとしては分割統治で、$T$ を $T_1, T_2$ と前後半分に切り、$\mathrm{dp_1}[i] := S[:i]$ と $T_1$ の LCS の長さ、$\mathrm{dp_2}[i] := \mathrm{rev}(S)[:i]$ と $\mathrm{rev}(T_2)$ の LCS の長さ、という 2 つを求めることで通過すべき位置が求まり、サイズが半分の問題に帰着できます。

コード

vector<int> _lcs(const string &s, const string &t) {
  int n = s.size(), m = t.size();
  vector<int> dp(n + 1);
  for (int i = 0; i < m; ++i) {
    vector<int> nx = dp;
    for (int j = 0; j < n; ++j) {
      nx[j + 1] = max(nx[j + 1], nx[j]);
      if (t[i] == s[j]) {
        nx[j + 1] = max(nx[j + 1], dp[j] + 1);
      }
    }
    dp = nx;
  }
  return dp;
}

void lcs(string s, string t, string &result) {
  if (s.size() <= 1 || t.size() <= 1) {
    for (auto &i : s) {
      for (auto &j : t) {
        if (i == j) {
          result += i;
          return;
        }
      }
    }
    return;
  }
  int n = s.size();
  int l1 = t.size() / 2;
  int l2 = (int)t.size() - l1;
  auto lef = _lcs(s, t.substr(0, l1));
  reverse(s.begin(), s.end());
  reverse(t.begin(), t.end());
  auto rig = _lcs(s, t.substr(0, l2));
  reverse(s.begin(), s.end());
  reverse(t.begin(), t.end());
  int max_val = -1, max_idx = -1;
  for (int i = 0; i < n + 1; ++i) {
    if (max_val < lef[i] + rig[n - i]) {
      max_val = lef[i] + rig[n - i];
      max_idx = i;
    }
  }
  if (max_val <= 0) {
    return;
  }
  lcs(s.substr(0, max_idx), t.substr(0, l1), result);
  lcs(s.substr(max_idx), t.substr(l1), result);
}

---

時間 $\Theta(|S| |T| / w)$、空間 $\Theta(|S| + |T|)$、復元なし

ジャッジ

• https://www.acmicpc.net/problem/18439
  • $|S|, |T| \le 50000$、復元なし、ML: 8 MB

参考文献

• https://koosaga.com/245
• https://m.blog.naver.com/jinhan814/222545486610

概要

$\mathrm{dp}(S, T)[i] := S[:i]$ と $T$ の LCS の長さ、というテーブルを高速に求めたいです。
ここで、$\mathrm{dp}$ の階差を考えると、$\mathrm{dp}[i] \le \mathrm{dp}[i+1] \le \mathrm{dp}[i] + 1$ であることからこれは各要素が $0$ または $1$ の数列になっています。
階差から $\mathrm{dp}$ の復元は容易なので、階差がビット演算で効率的に求まれば良いです。
より具体的には、$B[j] := \mathrm{dp}(S, T[:j])$ の階差数列、としたとき、$B[j+1]$ が $B[j]$ から効率的に求まれば良いです。

……自分が自信を持って解説できるのはここまでで（短すぎる）、

結論だけ言うと、これは以下のように更新できます。

B[j+1] = (B[j] | m) ^ ((B[j] | m) & ((B[j] | m) - (B[j] << 1 | 1)))
（m は S[i] == T[j] を満たす i にビットが立ったビット列）

お気持ちをなんとなくは理解しているんですが説明できるほどではないです。 誰か完全理解したら正しいお気持ちを教えてください。 （上の参考文献にはもう少しまともな説明が書かれています）

ビット列同士の減算を行う必要があるので通常の bitset ではできず自作する必要はありますが、これで効率的に更新することができます。

……というのが韓国で知られている典型っぽいんですが、上のビット演算は少し式変形できて、（自分の実装だと）さらに少しだけ速くなります。

B[j+1] = (B[j] | m) ^ ((B[j] | m) & ((B[j] | m) - (B[j] << 1 | 1)))
       = (B[j] | m) ^ ((B[j] | m) & ((B[j] | m) - (B[j] << 1) - 1))
       = (B[j] | m) ^ ((B[j] | m) & ~((B[j] << 1) - (B[j] | m)))
       = (B[j] | m) ^ ((B[j] | m) & ~(B[j] + B[j] - (B[j] | m)))
       = (B[j] | m) ^ ((B[j] | m) & ~(B[j] - (m & ~B[j])))
       = (B[j] | m) ^ ((B[j] | m) & ((B[j] | m) ^ (B[j] - (m & ~B[j]))))
       = (B[j] | m) & ((B[j] | m) ^ (B[j] | m) ^ (B[j] - (m & ~B[j])))
       = (B[j] | m) & (B[j] - (m & ~B[j]))

コード

int lcs(const string &s, const string &t) {
  int n = s.size();
  int w = (n + 63) >> 6;
  vector<vector<uint64_t>> m(26, vector<uint64_t>(w));
  for (int i = 0; i < n; ++i) {
    m[s[i] - 'A'][i >> 6] |= 1ULL << (i & 63);
  }
  vector<uint64_t> b(w);
  for (auto &c : t) {
    const auto &mc = m[c - 'A'];
    // b = (b - (mc & ~b)) & (mc | b)
    for (int i = 0, borrow = 0; i < w; ++i) {
      uint64_t x = mc[i] & ~b[i];
      int nx = b[i] < x || (b[i] == x && borrow);
      b[i] = (b[i] - x - borrow) & (mc[i] | b[i]);
      borrow = nx;
    }
  }
  int res = 0;
  for (auto &i : b) {
    res += __builtin_popcountll(i);
  }
  return res;
}

---

時間 $\Theta(|S| |T| / w)$、空間 $\Theta(|S| + |T|)$、復元あり

ジャッジ

• https://www.acmicpc.net/problem/18440
  • $|S|, |T| \le 50000$、復元あり、ML: 8 MB

概要

bitset + Hirschberg's algorithm でできます。

コード

vector<uint64_t> _lcs(const string &s, const string &t) {
  int n = s.size();
  int w = (n + 63) >> 6;
  vector<vector<uint64_t>> m(26, vector<uint64_t>(w));
  for (int i = 0; i < n; ++i) {
    m[s[i] - 'A'][i >> 6] |= 1ULL << (i & 63);
  }
  vector<uint64_t> b(w);
  for (auto &c : t) {
    const auto &mc = m[c - 'A'];
    // b = (b - (mc & ~b)) & (mc | b)
    for (int i = 0, borrow = 0; i < w; ++i) {
      uint64_t x = mc[i] & ~b[i];
      int nx = b[i] < x || (b[i] == x && borrow);
      b[i] = (b[i] - x - borrow) & (mc[i] | b[i]);
      borrow = nx;
    }
  }
  return b;
}

void lcs(string s, string t, string &result) {
  if (s.size() <= 1 || t.size() <= 1) {
    for (auto &i : s) {
      for (auto &j : t) {
        if (i == j) {
          result += i;
          return;
        }
      }
    }
    return;
  }
  int n = s.size();
  int l1 = t.size() / 2;
  int l2 = (int)t.size() - l1;
  auto lef = _lcs(s, t.substr(0, l1));
  reverse(s.begin(), s.end());
  reverse(t.begin(), t.end());
  auto rig = _lcs(s, t.substr(0, l2));
  reverse(s.begin(), s.end());
  reverse(t.begin(), t.end());
  int lc = 0, rc = 0;
  for (auto &i : rig) {
    rc += __builtin_popcountll(i);
  }
  int max_val = lc + rc;
  int max_idx = 0;
  for (int i = 0; i < n; ++i) {
    lc += (lef[i >> 6] >> (i & 63)) & 1;
    rc -= (rig[(n - i - 1) >> 6] >> ((n - i - 1) & 63)) & 1;
    if (max_val < lc + rc) {
      max_val = lc + rc;
      max_idx = i + 1;
    }
  }
  if (max_val == 0) {
    return;
  }
  lcs(s.substr(0, max_idx), t.substr(0, l1), result);
  lcs(s.substr(max_idx), t.substr(l1), result);
}

あとがき

LCS とかいうド典型問題なのに全然知らなかった。

ネタバレにもなるのでここでは挙げてないですが、賢くやって $\Theta(|S| |T|)$ が想定解な問題を $\Theta(|S|^2 |T| / w)$ とかでわりと余裕もって殴れたりしてすごいです（$|S|, |T| \le 3000$ くらいまでは行けることがある雰囲気）。

最近は solved.ac を埋めるのにハマっています。よければぜひ。