longest_common_subsequence.hpp¶
Return matched index pairs for one longest common subsequence in O(nm) time and memory.
求两个序列的一组最长公共子序列匹配下标；适合需要恢复方案而不只求长度的题目。
Implementation¶
#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 &first, const Second &second) {
int n = int(first.size());
int m = int(second.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 (first[i] == second[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>> result;
int i = 0;
int j = 0;
while (i < n && j < m) {
if (first[i] == second[j] && dp[i][j] == dp[i + 1][j + 1] + 1) {
result.emplace_back(i++, j++);
} else if (dp[i + 1][j] >= dp[i][j + 1]) {
i++;
} else {
j++;
}
}
return result;
}
/// @brief Return only the LCS length using the bit-parallel DP-row-difference
/// 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 difference 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 &first,
const Sequence &second) {
using value_type = typename Sequence::value_type;
const int n = int(first.size());
const int blocks = (n + 63) / 64;
if (n == 0 || second.empty()) {
return 0;
}
std::vector<value_type> coordinates(first.begin(), first.end());
std::sort(coordinates.begin(), coordinates.end());
coordinates.erase(std::unique(coordinates.begin(), coordinates.end()),
coordinates.end());
std::vector<std::vector<std::uint64_t>> matches(
coordinates.size(), std::vector<std::uint64_t>(blocks));
for (int position = 0; position < n; position++) {
int id = int(std::lower_bound(coordinates.begin(), coordinates.end(),
first[position]) -
coordinates.begin());
matches[id][position / 64] |= std::uint64_t(1) << (position % 64);
}
std::vector<std::uint64_t> difference(blocks);
for (const value_type &value : second) {
auto iterator = std::lower_bound(coordinates.begin(), coordinates.end(),
value);
const std::vector<std::uint64_t> *equal =
iterator != coordinates.end() && *iterator == value
? &matches[iterator - coordinates.begin()]
: nullptr;
std::uint64_t shift_carry = 1;
std::uint64_t borrow = 0;
for (int block = 0; block < blocks; block++) {
std::uint64_t old = difference[block];
std::uint64_t combined = old | (equal ? (*equal)[block] : 0);
std::uint64_t shifted = (old << 1) | shift_carry;
shift_carry = old >> 63;
std::uint64_t subtracted = combined - shifted - borrow;
std::uint64_t next_borrow =
combined < shifted || (borrow != 0 && combined == shifted);
difference[block] = combined & ~subtracted;
borrow = next_borrow;
}
}
int result = 0;
for (std::uint64_t block : difference) {
result += std::popcount(block);
}
return result;
}
} // 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 &first, const Second &second) {
int n = int(first.size());
int m = int(second.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 (first[i] == second[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>> result;
int i = 0;
int j = 0;
while (i < n && j < m) {
if (first[i] == second[j] && dp[i][j] == dp[i + 1][j + 1] + 1) {
result.emplace_back(i++, j++);
} else if (dp[i + 1][j] >= dp[i][j + 1]) {
i++;
} else {
j++;
}
}
return result;
}
/// @brief Return only the LCS length using the bit-parallel DP-row-difference
/// 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 difference 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 &first,
const Sequence &second) {
using value_type = typename Sequence::value_type;
const int n = int(first.size());
const int blocks = (n + 63) / 64;
if (n == 0 || second.empty()) {
return 0;
}
std::vector<value_type> coordinates(first.begin(), first.end());
std::sort(coordinates.begin(), coordinates.end());
coordinates.erase(std::unique(coordinates.begin(), coordinates.end()),
coordinates.end());
std::vector<std::vector<std::uint64_t>> matches(
coordinates.size(), std::vector<std::uint64_t>(blocks));
for (int position = 0; position < n; position++) {
int id = int(std::lower_bound(coordinates.begin(), coordinates.end(),
first[position]) -
coordinates.begin());
matches[id][position / 64] |= std::uint64_t(1) << (position % 64);
}
std::vector<std::uint64_t> difference(blocks);
for (const value_type &value : second) {
auto iterator = std::lower_bound(coordinates.begin(), coordinates.end(),
value);
const std::vector<std::uint64_t> *equal =
iterator != coordinates.end() && *iterator == value
? &matches[iterator - coordinates.begin()]
: nullptr;
std::uint64_t shift_carry = 1;
std::uint64_t borrow = 0;
for (int block = 0; block < blocks; block++) {
std::uint64_t old = difference[block];
std::uint64_t combined = old | (equal ? (*equal)[block] : 0);
std::uint64_t shifted = (old << 1) | shift_carry;
shift_carry = old >> 63;
std::uint64_t subtracted = combined - shifted - borrow;
std::uint64_t next_borrow =
combined < shifted || (borrow != 0 && combined == shifted);
difference[block] = combined & ~subtracted;
borrow = next_borrow;
}
}
int result = 0;
for (std::uint64_t block : difference) {
result += std::popcount(block);
}
return result;
}
} // namespace noya