static_range_lis_query.hpp¶
Answer LIS lengths on subarrays of a permutation. Seaweed doubling represents semi-local LCS against the sorted permutation as a subunit-Monge permutation. Unit-Monge distance products merge the two halves, and a wavelet matrix over the resulting critical points turns each interval LIS into one orthogonal counting query.
Verified by static_range_lis_query.
回答静态序列多个子区间的最长严格上升子序列长度，而不必为每个询问重新跑 LIS。
Implementation¶
#ifndef NOYA_STATIC_RANGE_LIS_QUERY_HPP
#define NOYA_STATIC_RANGE_LIS_QUERY_HPP 1
/// @complexity Time: O(n log^2 n) preprocessing and O(log n) per query.
/// Space: O(n log n).
#include <algorithm>
#include <cassert>
#include <climits>
#include <numeric>
#include <utility>
#include <vector>
namespace noya {
namespace static_range_lis_internal {
using uint = unsigned int;
using ll = long long;
static constexpr int word_bits = CHAR_BIT * sizeof(uint);
inline int popcount(uint value) {
#ifdef __GNUC__
return __builtin_popcount(value);
#else
static_assert(word_bits == 32);
value -= value >> 1 & 0x55555555;
value = (value & 0x33333333) + (value >> 2 & 0x33333333);
value = value + (value >> 4) & 0x0f0f0f0f;
return value * 0x01010101 >> 24 & 0x3f;
#endif
}
class bit_vector {
struct node {
uint bit = 0;
int sum = 0;
};
std::vector<node> data;
public:
explicit bit_vector(uint n) : data(n / word_bits + 1) {}
void set(uint index) {
data[index / word_bits].bit |= uint(1) << index % word_bits;
data[index / word_bits].sum++;
}
void build() {
for (int i = 1; i < int(data.size()); i++) {
data[i].sum += data[i - 1].sum;
}
}
int rank(uint index) const {
return data[index / word_bits].sum -
popcount(data[index / word_bits].bit &
(~uint(0) << (index % word_bits)));
}
int ones() const { return data.back().sum; }
};
class wavelet_matrix {
template <class Integer> static bool test(Integer value, int bit) {
return (value & (Integer(1) << bit)) != Integer(0);
}
std::vector<bit_vector> levels;
public:
template <class Integer>
wavelet_matrix(int bit_length, std::vector<Integer> values)
: levels(bit_length, bit_vector(values.size())) {
int n = int(values.size());
std::vector<Integer> zero;
zero.reserve(n);
for (int bit = bit_length - 1; bit >= 0; bit--) {
bit_vector &level = levels[bit];
auto one = values.begin();
for (int i = 0; i < n; i++) {
if (test(values[i], bit)) {
level.set(i);
*one++ = values[i];
} else {
zero.push_back(values[i]);
}
}
level.build();
std::copy(zero.begin(), zero.end(), one);
zero.clear();
}
}
int count_less_than(int left, int right, ll key) const {
int result = right - left;
for (int bit = int(levels.size()) - 1; bit >= 0; bit--) {
const bit_vector &level = levels[bit];
int rank_left = level.rank(left);
int rank_right = level.rank(right);
if (test(key, bit)) {
left = rank_left;
right = rank_right;
} else {
result -= rank_right - rank_left;
int ones = level.ones();
left += ones - rank_left;
right += ones - rank_right;
}
}
return result - (right - left);
}
};
using permutation = std::vector<int>;
using iterator = permutation::iterator;
static constexpr int none = -1;
inline permutation inverse(const permutation &p) {
permutation result(p.size(), none);
for (int i = 0; i < int(p.size()); i++) {
if (p[i] != none) {
result[p[i]] = i;
}
}
return result;
}
inline void unit_monge_distance_product(int n, iterator stack,
const iterator first,
const iterator second) {
if (n == 1) {
stack[0] = 0;
return;
}
const iterator result_row = stack;
stack += n;
const iterator result_column = stack;
stack += n;
const auto solve_part = [=](int length, const auto &belongs,
const auto &normalize) {
const iterator first_half = stack;
const iterator first_map = stack + length;
const iterator second_half = stack + 2 * length;
const iterator second_map = stack + 3 * length;
const auto split = [=](const iterator values, iterator half,
iterator map) {
for (int i = 0; i < n; i++) {
if (belongs(values[i])) {
*half++ = normalize(values[i]);
*map++ = i;
}
}
};
split(first, first_half, first_map);
split(second, second_half, second_map);
const iterator product = stack + 4 * length;
unit_monge_distance_product(length, product, first_half, second_half);
for (int i = 0; i < length; i++) {
int row = first_map[i];
int column = second_map[product[i]];
result_row[row] = column;
result_column[column] = row;
}
};
int middle = n / 2;
solve_part(middle, [middle](int value) { return value < middle; },
[](int value) { return value; });
solve_part(n - middle, [middle](int value) { return value >= middle; },
[middle](int value) { return value - middle; });
struct diagonal_iterator {
int delta = 0;
int column = 0;
};
int row = n;
const auto move_right = [&](diagonal_iterator &it) {
if (second[it.column] < middle) {
if (result_column[it.column] >= row) {
it.delta++;
}
} else if (result_column[it.column] < row) {
it.delta++;
}
it.column++;
};
const auto move_up = [&](diagonal_iterator &it) {
if (first[row] < middle) {
if (result_row[row] >= it.column) {
it.delta--;
}
} else if (result_row[row] < it.column) {
it.delta--;
}
};
diagonal_iterator negative, positive;
while (row != 0) {
while (positive.column != n) {
diagonal_iterator candidate = positive;
move_right(candidate);
if (candidate.delta != 0) {
break;
}
positive = candidate;
}
row--;
move_up(negative);
move_up(positive);
while (negative.delta != 0) {
move_right(negative);
}
if (negative.column > positive.column) {
result_row[row] = positive.column;
}
}
}
inline permutation subunit_monge_distance_product(permutation first,
permutation second) {
int n = int(first.size());
permutation inverse_first = inverse(first);
permutation inverse_second = inverse(second);
std::swap(second, inverse_second);
permutation first_map, second_map;
for (int i = n - 1; i >= 0; i--) {
if (first[i] != none) {
first_map.push_back(i);
first[n - int(first_map.size())] = first[i];
}
}
std::reverse(first_map.begin(), first_map.end());
{
int count = 0;
for (int i = 0; i < n; i++) {
if (inverse_first[i] == none) {
first[count++] = i;
}
}
}
for (int i = 0; i < n; i++) {
if (second[i] != none) {
second[second_map.size()] = second[i];
second_map.push_back(i);
}
}
{
int count = int(second_map.size());
for (int i = 0; i < n; i++) {
if (inverse_second[i] == none) {
second[count++] = i;
}
}
}
permutation workspace([](int length) {
int size = 0;
while (length > 1) {
size += 2 * length;
length = (length + 1) / 2;
size += 4 * length;
}
return size + 1;
}(n));
unit_monge_distance_product(n, workspace.begin(), first.begin(),
second.begin());
permutation result(n, none);
for (int i = 0; i < int(first_map.size()); i++) {
int column = workspace[n - int(first_map.size()) + i];
if (column < int(second_map.size())) {
result[first_map[i]] = second_map[column];
}
}
return result;
}
inline permutation seaweed_doubling(const permutation &p) {
int n = int(p.size());
if (n == 1) {
return {none};
}
int middle = n / 2;
permutation low, high, low_map, high_map;
for (int i = 0; i < n; i++) {
if (p[i] < middle) {
low.push_back(p[i]);
low_map.push_back(i);
} else {
high.push_back(p[i] - middle);
high_map.push_back(i);
}
}
low = seaweed_doubling(low);
high = seaweed_doubling(high);
permutation low_padded(n), high_padded(n);
std::iota(low_padded.begin(), low_padded.end(), 0);
std::iota(high_padded.begin(), high_padded.end(), 0);
for (int i = 0; i < middle; i++) {
low_padded[low_map[i]] = low[i] == none ? none : low_map[low[i]];
}
for (int i = 0; middle + i < n; i++) {
high_padded[high_map[i]] = high[i] == none ? none : high_map[high[i]];
}
return subunit_monge_distance_product(std::move(low_padded),
std::move(high_padded));
}
inline bool is_permutation(const permutation &p) {
std::vector<bool> used(p.size());
for (int value : p) {
if (value < 0 || value >= int(p.size()) || used[value]) {
return false;
}
used[value] = true;
}
return true;
}
inline wavelet_matrix build_wavelet_matrix(const permutation &p) {
assert(is_permutation(p));
int n = int(p.size());
permutation row;
if (n != 0) {
row = seaweed_doubling(p);
}
for (int &value : row) {
if (value == none) {
value = n;
}
}
int bit_length = 0;
for (int value = n; value > 0; value /= 2) {
bit_length++;
}
return wavelet_matrix(bit_length, std::move(row));
}
} // namespace static_range_lis_internal
/// @brief Answer LIS lengths on subarrays of a permutation.
/// Seaweed doubling represents semi-local LCS against the sorted permutation
/// as a subunit-Monge permutation. Unit-Monge distance products merge the two
/// halves, and a wavelet matrix over the resulting critical points turns each
/// interval LIS into one orthogonal counting query.
class static_range_lis_query {
int n;
static_range_lis_internal::wavelet_matrix matrix;
public:
static_range_lis_query() : static_range_lis_query(std::vector<int>{}) {}
explicit static_range_lis_query(const std::vector<int> &permutation)
: n(int(permutation.size())),
matrix(static_range_lis_internal::build_wavelet_matrix(permutation)) {}
int query(int left, int right) const {
assert(0 <= left && left <= right && right <= n);
return (right - left) - matrix.count_less_than(left, n, right);
}
};
} // namespace noya
#endif // NOYA_STATIC_RANGE_LIS_QUERY_HPP
#include <algorithm>
#include <cassert>
#include <climits>
#include <numeric>
#include <utility>
#include <vector>
/// @complexity Time: O(n log^2 n) preprocessing and O(log n) per query.
/// Space: O(n log n).
namespace noya {
namespace static_range_lis_internal {
using uint = unsigned int;
using ll = long long;
static constexpr int word_bits = CHAR_BIT * sizeof(uint);
inline int popcount(uint value) {
#ifdef __GNUC__
return __builtin_popcount(value);
#else
static_assert(word_bits == 32);
value -= value >> 1 & 0x55555555;
value = (value & 0x33333333) + (value >> 2 & 0x33333333);
value = value + (value >> 4) & 0x0f0f0f0f;
return value * 0x01010101 >> 24 & 0x3f;
#endif
}
class bit_vector {
struct node {
uint bit = 0;
int sum = 0;
};
std::vector<node> data;
public:
explicit bit_vector(uint n) : data(n / word_bits + 1) {}
void set(uint index) {
data[index / word_bits].bit |= uint(1) << index % word_bits;
data[index / word_bits].sum++;
}
void build() {
for (int i = 1; i < int(data.size()); i++) {
data[i].sum += data[i - 1].sum;
}
}
int rank(uint index) const {
return data[index / word_bits].sum -
popcount(data[index / word_bits].bit &
(~uint(0) << (index % word_bits)));
}
int ones() const { return data.back().sum; }
};
class wavelet_matrix {
template <class Integer> static bool test(Integer value, int bit) {
return (value & (Integer(1) << bit)) != Integer(0);
}
std::vector<bit_vector> levels;
public:
template <class Integer>
wavelet_matrix(int bit_length, std::vector<Integer> values)
: levels(bit_length, bit_vector(values.size())) {
int n = int(values.size());
std::vector<Integer> zero;
zero.reserve(n);
for (int bit = bit_length - 1; bit >= 0; bit--) {
bit_vector &level = levels[bit];
auto one = values.begin();
for (int i = 0; i < n; i++) {
if (test(values[i], bit)) {
level.set(i);
*one++ = values[i];
} else {
zero.push_back(values[i]);
}
}
level.build();
std::copy(zero.begin(), zero.end(), one);
zero.clear();
}
}
int count_less_than(int left, int right, ll key) const {
int result = right - left;
for (int bit = int(levels.size()) - 1; bit >= 0; bit--) {
const bit_vector &level = levels[bit];
int rank_left = level.rank(left);
int rank_right = level.rank(right);
if (test(key, bit)) {
left = rank_left;
right = rank_right;
} else {
result -= rank_right - rank_left;
int ones = level.ones();
left += ones - rank_left;
right += ones - rank_right;
}
}
return result - (right - left);
}
};
using permutation = std::vector<int>;
using iterator = permutation::iterator;
static constexpr int none = -1;
inline permutation inverse(const permutation &p) {
permutation result(p.size(), none);
for (int i = 0; i < int(p.size()); i++) {
if (p[i] != none) {
result[p[i]] = i;
}
}
return result;
}
inline void unit_monge_distance_product(int n, iterator stack,
const iterator first,
const iterator second) {
if (n == 1) {
stack[0] = 0;
return;
}
const iterator result_row = stack;
stack += n;
const iterator result_column = stack;
stack += n;
const auto solve_part = [=](int length, const auto &belongs,
const auto &normalize) {
const iterator first_half = stack;
const iterator first_map = stack + length;
const iterator second_half = stack + 2 * length;
const iterator second_map = stack + 3 * length;
const auto split = [=](const iterator values, iterator half,
iterator map) {
for (int i = 0; i < n; i++) {
if (belongs(values[i])) {
*half++ = normalize(values[i]);
*map++ = i;
}
}
};
split(first, first_half, first_map);
split(second, second_half, second_map);
const iterator product = stack + 4 * length;
unit_monge_distance_product(length, product, first_half, second_half);
for (int i = 0; i < length; i++) {
int row = first_map[i];
int column = second_map[product[i]];
result_row[row] = column;
result_column[column] = row;
}
};
int middle = n / 2;
solve_part(middle, [middle](int value) { return value < middle; },
[](int value) { return value; });
solve_part(n - middle, [middle](int value) { return value >= middle; },
[middle](int value) { return value - middle; });
struct diagonal_iterator {
int delta = 0;
int column = 0;
};
int row = n;
const auto move_right = [&](diagonal_iterator &it) {
if (second[it.column] < middle) {
if (result_column[it.column] >= row) {
it.delta++;
}
} else if (result_column[it.column] < row) {
it.delta++;
}
it.column++;
};
const auto move_up = [&](diagonal_iterator &it) {
if (first[row] < middle) {
if (result_row[row] >= it.column) {
it.delta--;
}
} else if (result_row[row] < it.column) {
it.delta--;
}
};
diagonal_iterator negative, positive;
while (row != 0) {
while (positive.column != n) {
diagonal_iterator candidate = positive;
move_right(candidate);
if (candidate.delta != 0) {
break;
}
positive = candidate;
}
row--;
move_up(negative);
move_up(positive);
while (negative.delta != 0) {
move_right(negative);
}
if (negative.column > positive.column) {
result_row[row] = positive.column;
}
}
}
inline permutation subunit_monge_distance_product(permutation first,
permutation second) {
int n = int(first.size());
permutation inverse_first = inverse(first);
permutation inverse_second = inverse(second);
std::swap(second, inverse_second);
permutation first_map, second_map;
for (int i = n - 1; i >= 0; i--) {
if (first[i] != none) {
first_map.push_back(i);
first[n - int(first_map.size())] = first[i];
}
}
std::reverse(first_map.begin(), first_map.end());
{
int count = 0;
for (int i = 0; i < n; i++) {
if (inverse_first[i] == none) {
first[count++] = i;
}
}
}
for (int i = 0; i < n; i++) {
if (second[i] != none) {
second[second_map.size()] = second[i];
second_map.push_back(i);
}
}
{
int count = int(second_map.size());
for (int i = 0; i < n; i++) {
if (inverse_second[i] == none) {
second[count++] = i;
}
}
}
permutation workspace([](int length) {
int size = 0;
while (length > 1) {
size += 2 * length;
length = (length + 1) / 2;
size += 4 * length;
}
return size + 1;
}(n));
unit_monge_distance_product(n, workspace.begin(), first.begin(),
second.begin());
permutation result(n, none);
for (int i = 0; i < int(first_map.size()); i++) {
int column = workspace[n - int(first_map.size()) + i];
if (column < int(second_map.size())) {
result[first_map[i]] = second_map[column];
}
}
return result;
}
inline permutation seaweed_doubling(const permutation &p) {
int n = int(p.size());
if (n == 1) {
return {none};
}
int middle = n / 2;
permutation low, high, low_map, high_map;
for (int i = 0; i < n; i++) {
if (p[i] < middle) {
low.push_back(p[i]);
low_map.push_back(i);
} else {
high.push_back(p[i] - middle);
high_map.push_back(i);
}
}
low = seaweed_doubling(low);
high = seaweed_doubling(high);
permutation low_padded(n), high_padded(n);
std::iota(low_padded.begin(), low_padded.end(), 0);
std::iota(high_padded.begin(), high_padded.end(), 0);
for (int i = 0; i < middle; i++) {
low_padded[low_map[i]] = low[i] == none ? none : low_map[low[i]];
}
for (int i = 0; middle + i < n; i++) {
high_padded[high_map[i]] = high[i] == none ? none : high_map[high[i]];
}
return subunit_monge_distance_product(std::move(low_padded),
std::move(high_padded));
}
inline bool is_permutation(const permutation &p) {
std::vector<bool> used(p.size());
for (int value : p) {
if (value < 0 || value >= int(p.size()) || used[value]) {
return false;
}
used[value] = true;
}
return true;
}
inline wavelet_matrix build_wavelet_matrix(const permutation &p) {
assert(is_permutation(p));
int n = int(p.size());
permutation row;
if (n != 0) {
row = seaweed_doubling(p);
}
for (int &value : row) {
if (value == none) {
value = n;
}
}
int bit_length = 0;
for (int value = n; value > 0; value /= 2) {
bit_length++;
}
return wavelet_matrix(bit_length, std::move(row));
}
} // namespace static_range_lis_internal
/// @brief Answer LIS lengths on subarrays of a permutation.
/// Seaweed doubling represents semi-local LCS against the sorted permutation
/// as a subunit-Monge permutation. Unit-Monge distance products merge the two
/// halves, and a wavelet matrix over the resulting critical points turns each
/// interval LIS into one orthogonal counting query.
class static_range_lis_query {
int n;
static_range_lis_internal::wavelet_matrix matrix;
public:
static_range_lis_query() : static_range_lis_query(std::vector<int>{}) {}
explicit static_range_lis_query(const std::vector<int> &permutation)
: n(int(permutation.size())),
matrix(static_range_lis_internal::build_wavelet_matrix(permutation)) {}
int query(int left, int right) const {
assert(0 <= left && left <= right && right <= n);
return (right - left) - matrix.count_less_than(left, n, right);
}
};
} // namespace noya