longest_common_subsequence.hpp¶
求两个序列的一组最长公共子序列匹配下标；适合需要恢复方案而不只求长度的题目。
Complexity: Time: O(nm) for a certificate; bit-parallel length uses O(n log n + d ceil(n/64) + m(log d + ceil(n/64))) time and Space: O(nm) for a certificate or O(d ceil(n/64)) for bit-parallel length, where d is the number of distinct first-sequence values.
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(nm) for a certificate; bit-parallel length uses
/// O(n log n + d ceil(n/64) + m(log d + ceil(n/64))) time and
/// Space: O(nm) for a certificate or O(d ceil(n/64)) for bit-parallel length,
/// where d is the number of distinct first-sequence values.
#include <algorithm>
#include <bit>
#include <cstdint>
#include <utility>
#include <vector>
namespace noya {
/// @brief Return matched index pairs for one longest common subsequence in
/// O(nm) time and memory.
template <class First, class Second>
std::vector<std::pair<int, int>> longest_common_subsequence(const First &a,
const Second &b) {
int n = int(a.size());
int m = int(b.size());
std::vector dp(n + 1, std::vector<int>(m + 1));
for (int i = n - 1; i >= 0; i--) {
for (int j = m - 1; j >= 0; j--) {
if (a[i] == b[j]) {
dp[i][j] = dp[i + 1][j + 1] + 1;
} else {
dp[i][j] = std::max(dp[i + 1][j], dp[i][j + 1]);
}
}
}
std::vector<std::pair<int, int>> res;
int i = 0;
int j = 0;
while (i < n && j < m) {
if (a[i] == b[j] && dp[i][j] == dp[i + 1][j + 1] + 1) {
res.emplace_back(i++, j++);
} else if (dp[i + 1][j] >= dp[i][j + 1]) {
i++;
} else {
j++;
}
}
return res;
}
/// @brief Return only the LCS length using the bit-parallel DP-row-dif
/// update. Consecutive cells of a standard LCS DP row differ by zero or one,
/// so the complete row is represented by a bitset D of those differences.
/// For the next value, X marks equal positions in the first sequence. The
/// next row is `D = (D | X) & ~((D | X) - ((D << 1) | 1))`; the loop below
/// performs that subtraction word by word, carrying the shift and borrow
/// across 64-bit blocks. The popcount of the final dif row is the LCS
/// length. This is intended for long sequences when the quadratic certificate
/// table is too large.
template <class Sequence>
int longest_common_subsequence_length_bitset(const Sequence &a,
const Sequence &b) {
using value_type = typename Sequence::value_type;
const int n = int(a.size());
const int blk = (n + 63) / 64;
if (n == 0 || b.empty()) {
return 0;
}
std::vector<value_type> xs(a.begin(), a.end());
std::sort(xs.begin(), xs.end());
xs.erase(std::unique(xs.begin(), xs.end()), xs.end());
std::vector<std::vector<std::uint64_t>> mat(xs.size(),
std::vector<std::uint64_t>(blk));
for (int pos = 0; pos < n; pos++) {
int id = int(std::lower_bound(xs.begin(), xs.end(), a[pos]) - xs.begin());
mat[id][pos / 64] |= std::uint64_t(1) << (pos % 64);
}
std::vector<std::uint64_t> dif(blk);
for (const value_type &val : b) {
auto it = std::lower_bound(xs.begin(), xs.end(), val);
const std::vector<std::uint64_t> *eq =
it != xs.end() && *it == val ? &mat[it - xs.begin()] : nullptr;
std::uint64_t sc = 1;
std::uint64_t brw = 0;
for (int bl1 = 0; bl1 < blk; bl1++) {
std::uint64_t old = dif[bl1];
std::uint64_t sum = old | (eq ? (*eq)[bl1] : 0);
std::uint64_t sh = (old << 1) | sc;
sc = old >> 63;
std::uint64_t sub = sum - sh - brw;
std::uint64_t nb = sum < sh || (brw != 0 && sum == sh);
dif[bl1] = sum & ~sub;
brw = nb;
}
}
int res = 0;
for (std::uint64_t bl1 : dif) {
res += std::popcount(bl1);
}
return res;
}
} // namespace noya
#ifndef NOYA_LONGEST_COMMON_SUBSEQUENCE_HPP
#define NOYA_LONGEST_COMMON_SUBSEQUENCE_HPP 1
/// @complexity Time: O(nm) for a certificate; bit-parallel length uses
/// O(n log n + d ceil(n/64) + m(log d + ceil(n/64))) time and
/// Space: O(nm) for a certificate or O(d ceil(n/64)) for bit-parallel length,
/// where d is the number of distinct first-sequence values.
#include <algorithm>
#include <bit>
#include <cstdint>
#include <utility>
#include <vector>
namespace noya {
/// @brief Return matched index pairs for one longest common subsequence in
/// O(nm) time and memory.
template <class First, class Second>
std::vector<std::pair<int, int>> longest_common_subsequence(const First &a,
const Second &b) {
int n = int(a.size());
int m = int(b.size());
std::vector dp(n + 1, std::vector<int>(m + 1));
for (int i = n - 1; i >= 0; i--) {
for (int j = m - 1; j >= 0; j--) {
if (a[i] == b[j]) {
dp[i][j] = dp[i + 1][j + 1] + 1;
} else {
dp[i][j] = std::max(dp[i + 1][j], dp[i][j + 1]);
}
}
}
std::vector<std::pair<int, int>> res;
int i = 0;
int j = 0;
while (i < n && j < m) {
if (a[i] == b[j] && dp[i][j] == dp[i + 1][j + 1] + 1) {
res.emplace_back(i++, j++);
} else if (dp[i + 1][j] >= dp[i][j + 1]) {
i++;
} else {
j++;
}
}
return res;
}
/// @brief Return only the LCS length using the bit-parallel DP-row-dif
/// update. Consecutive cells of a standard LCS DP row differ by zero or one,
/// so the complete row is represented by a bitset D of those differences.
/// For the next value, X marks equal positions in the first sequence. The
/// next row is `D = (D | X) & ~((D | X) - ((D << 1) | 1))`; the loop below
/// performs that subtraction word by word, carrying the shift and borrow
/// across 64-bit blocks. The popcount of the final dif row is the LCS
/// length. This is intended for long sequences when the quadratic certificate
/// table is too large.
template <class Sequence>
int longest_common_subsequence_length_bitset(const Sequence &a,
const Sequence &b) {
using value_type = typename Sequence::value_type;
const int n = int(a.size());
const int blk = (n + 63) / 64;
if (n == 0 || b.empty()) {
return 0;
}
std::vector<value_type> xs(a.begin(), a.end());
std::sort(xs.begin(), xs.end());
xs.erase(std::unique(xs.begin(), xs.end()), xs.end());
std::vector<std::vector<std::uint64_t>> mat(xs.size(),
std::vector<std::uint64_t>(blk));
for (int pos = 0; pos < n; pos++) {
int id = int(std::lower_bound(xs.begin(), xs.end(), a[pos]) - xs.begin());
mat[id][pos / 64] |= std::uint64_t(1) << (pos % 64);
}
std::vector<std::uint64_t> dif(blk);
for (const value_type &val : b) {
auto it = std::lower_bound(xs.begin(), xs.end(), val);
const std::vector<std::uint64_t> *eq =
it != xs.end() && *it == val ? &mat[it - xs.begin()] : nullptr;
std::uint64_t sc = 1;
std::uint64_t brw = 0;
for (int bl1 = 0; bl1 < blk; bl1++) {
std::uint64_t old = dif[bl1];
std::uint64_t sum = old | (eq ? (*eq)[bl1] : 0);
std::uint64_t sh = (old << 1) | sc;
sc = old >> 63;
std::uint64_t sub = sum - sh - brw;
std::uint64_t nb = sum < sh || (brw != 0 && sum == sh);
dif[bl1] = sum & ~sub;
brw = nb;
}
}
int res = 0;
for (std::uint64_t bl1 : dif) {
res += std::popcount(bl1);
}
return res;
}
} // namespace noya
#endif // NOYA_LONGEST_COMMON_SUBSEQUENCE_HPP
#include <algorithm>
#include <bit>
#include <cstdint>
#include <utility>
#include <vector>
/// @complexity Time: O(nm) for a certificate; bit-parallel length uses
/// O(n log n + d ceil(n/64) + m(log d + ceil(n/64))) time and
/// Space: O(nm) for a certificate or O(d ceil(n/64)) for bit-parallel length,
/// where d is the number of distinct first-sequence values.
namespace noya {
/// @brief Return matched index pairs for one longest common subsequence in
/// O(nm) time and memory.
template <class First, class Second>
std::vector<std::pair<int, int>> longest_common_subsequence(const First &a,
const Second &b) {
int n = int(a.size());
int m = int(b.size());
std::vector dp(n + 1, std::vector<int>(m + 1));
for (int i = n - 1; i >= 0; i--) {
for (int j = m - 1; j >= 0; j--) {
if (a[i] == b[j]) {
dp[i][j] = dp[i + 1][j + 1] + 1;
} else {
dp[i][j] = std::max(dp[i + 1][j], dp[i][j + 1]);
}
}
}
std::vector<std::pair<int, int>> res;
int i = 0;
int j = 0;
while (i < n && j < m) {
if (a[i] == b[j] && dp[i][j] == dp[i + 1][j + 1] + 1) {
res.emplace_back(i++, j++);
} else if (dp[i + 1][j] >= dp[i][j + 1]) {
i++;
} else {
j++;
}
}
return res;
}
/// @brief Return only the LCS length using the bit-parallel DP-row-dif
/// update. Consecutive cells of a standard LCS DP row differ by zero or one,
/// so the complete row is represented by a bitset D of those differences.
/// For the next value, X marks equal positions in the first sequence. The
/// next row is `D = (D | X) & ~((D | X) - ((D << 1) | 1))`; the loop below
/// performs that subtraction word by word, carrying the shift and borrow
/// across 64-bit blocks. The popcount of the final dif row is the LCS
/// length. This is intended for long sequences when the quadratic certificate
/// table is too large.
template <class Sequence>
int longest_common_subsequence_length_bitset(const Sequence &a,
const Sequence &b) {
using value_type = typename Sequence::value_type;
const int n = int(a.size());
const int blk = (n + 63) / 64;
if (n == 0 || b.empty()) {
return 0;
}
std::vector<value_type> xs(a.begin(), a.end());
std::sort(xs.begin(), xs.end());
xs.erase(std::unique(xs.begin(), xs.end()), xs.end());
std::vector<std::vector<std::uint64_t>> mat(xs.size(),
std::vector<std::uint64_t>(blk));
for (int pos = 0; pos < n; pos++) {
int id = int(std::lower_bound(xs.begin(), xs.end(), a[pos]) - xs.begin());
mat[id][pos / 64] |= std::uint64_t(1) << (pos % 64);
}
std::vector<std::uint64_t> dif(blk);
for (const value_type &val : b) {
auto it = std::lower_bound(xs.begin(), xs.end(), val);
const std::vector<std::uint64_t> *eq =
it != xs.end() && *it == val ? &mat[it - xs.begin()] : nullptr;
std::uint64_t sc = 1;
std::uint64_t brw = 0;
for (int bl1 = 0; bl1 < blk; bl1++) {
std::uint64_t old = dif[bl1];
std::uint64_t sum = old | (eq ? (*eq)[bl1] : 0);
std::uint64_t sh = (old << 1) | sc;
sc = old >> 63;
std::uint64_t sub = sum - sh - brw;
std::uint64_t nb = sum < sh || (brw != 0 && sum == sh);
dif[bl1] = sum & ~sub;
brw = nb;
}
}
int res = 0;
for (std::uint64_t bl1 : dif) {
res += std::popcount(bl1);
}
return res;
}
} // namespace noya