m1une's library

This documentation is automatically generated by online-judge-tools/verification-helper

View on GitHub

:heavy_check_mark: Static Range LIS Query
(ds/range_query/range_lis_query.hpp)

Overview

RangeLisQuery<T> preprocesses a static array and returns the length of its strict longest increasing subsequence in any subarray [left, right). Equal values cannot both appear in the increasing subsequence.

The implementation converts the array to an equivalent permutation and builds its semi-local LIS (seaweed) matrix. A wavelet matrix then counts the seaweeds crossing each query rectangle.

Requirements

T must support < as a strict weak ordering. The array is static after construction. Duplicate values are supported.

Complexity

Operation Time Memory
Construction $O(N\log^2 N)$ $O(N\log N)$
Query $O(\log N)$ $O(1)$

Methods

Method Complexity Description
RangeLisQuery() $O(1)$ Constructs an empty structure.
explicit RangeLisQuery(const std::vector<T>& values) $O(N\log^2 N)$ Builds the structure for values.
int query(int left, int right) const $O(\log N)$ Returns the strict LIS length in [left, right).
int lis_length(int left, int right) const $O(\log N)$ Alias of query.
int size() const $O(1)$ Returns the array size.
bool empty() const $O(1)$ Returns whether the array is empty.

Example

#include "ds/range_query/range_lis_query.hpp"

#include <iostream>
#include <vector>

int main() {
    std::vector<int> values = {3, 1, 2, 2, 4};
    m1une::ds::RangeLisQuery<int> lis(values);

    std::cout << lis.query(0, 5) << "\n"; // 3: 1, 2, 4
    std::cout << lis.query(1, 4) << "\n"; // 2: 1, 2
}

Depends on

Verified with

Code

#ifndef M1UNE_DS_RANGE_QUERY_RANGE_LIS_QUERY_HPP
#define M1UNE_DS_RANGE_QUERY_RANGE_LIS_QUERY_HPP 1

#include "../wavelet_matrix/wavelet_matrix.hpp"

#include <algorithm>
#include <cassert>
#include <numeric>
#include <utility>
#include <vector>

namespace m1une {
namespace ds {

namespace range_lis_query_internal {

constexpr int none = -1;
using Permutation = std::vector<int>;
using Iterator = Permutation::iterator;

inline Permutation inverse(const Permutation& permutation) {
    int n = int(permutation.size());
    Permutation result(n, none);
    for (int i = 0; i < n; i++) {
        if (permutation[i] != none) result[permutation[i]] = i;
    }
    return result;
}

// Distance multiplication of two unit-Monge matrices, represented by their
// permutations. `workspace` must have the size used by subunit_monge_product.
inline void unit_monge_product(
    int n,
    Iterator workspace,
    Iterator first,
    Iterator second
) {
    if (n == 1) {
        workspace[0] = 0;
        return;
    }

    Iterator result_row = workspace;
    workspace += n;
    Iterator result_column = workspace;
    workspace += n;

    auto map_half = [=](int length, const auto& belongs, const auto& map) {
        Iterator first_half = workspace;
        Iterator first_position = workspace + length;
        Iterator second_half = workspace + 2 * length;
        Iterator second_position = workspace + 3 * length;

        auto split = [=](Iterator source, Iterator half, Iterator position) {
            for (int i = 0; i < n; i++) {
                if (belongs(source[i])) {
                    *half++ = map(source[i]);
                    *position++ = i;
                }
            }
        };
        split(first, first_half, first_position);
        split(second, second_half, second_position);

        Iterator product = workspace + 4 * length;
        unit_monge_product(
            length,
            product,
            first_half,
            second_half
        );
        for (int i = 0; i < length; i++) {
            int row = first_position[i];
            int column = second_position[product[i]];
            result_row[row] = column;
            result_column[column] = row;
        }
    };

    int middle = n / 2;
    map_half(
        middle,
        [middle](int value) { return value < middle; },
        [](int value) { return value; }
    );
    map_half(
        n - middle,
        [middle](int value) { return value >= middle; },
        [middle](int value) { return value - middle; }
    );

    struct DiagonalIterator {
        int delta = 0;
        int column = 0;
    };

    int row = n;
    auto move_right = [&](DiagonalIterator& iterator) {
        if (second[iterator.column] < middle) {
            if (result_column[iterator.column] >= row) iterator.delta++;
        } else {
            if (result_column[iterator.column] < row) iterator.delta++;
        }
        iterator.column++;
    };
    auto move_up = [&](DiagonalIterator& iterator) {
        if (first[row] < middle) {
            if (result_row[row] >= iterator.column) iterator.delta--;
        } else {
            if (result_row[row] < iterator.column) iterator.delta--;
        }
    };

    DiagonalIterator negative;
    DiagonalIterator positive;
    while (row != 0) {
        while (positive.column != n) {
            DiagonalIterator next = positive;
            move_right(next);
            if (next.delta != 0) break;
            positive = next;
        }
        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 int workspace_size(int n) {
    int result = 1;
    while (n > 1) {
        result += 2 * n;
        n = (n + 1) / 2;
        result += 4 * n;
    }
    return result;
}

inline Permutation subunit_monge_product(
    Permutation first,
    Permutation second
) {
    int n = int(first.size());
    Permutation first_inverse = inverse(first);
    Permutation second_inverse = inverse(second);
    std::swap(second, second_inverse);

    Permutation first_map;
    Permutation 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 missing = 0;
    for (int i = 0; i < n; i++) {
        if (first_inverse[i] == none) first[missing++] = i;
    }
    for (int i = 0; i < n; i++) {
        if (second[i] != none) {
            second[int(second_map.size())] = second[i];
            second_map.push_back(i);
        }
    }
    missing = int(second_map.size());
    for (int i = 0; i < n; i++) {
        if (second_inverse[i] == none) second[missing++] = i;
    }

    Permutation workspace(workspace_size(n));
    unit_monge_product(
        n,
        workspace.begin(),
        first.begin(),
        second.begin()
    );

    Permutation result(n, none);
    int first_count = int(first_map.size());
    int second_count = int(second_map.size());
    for (int i = 0; i < first_count; i++) {
        int mapped = workspace[n - first_count + i];
        if (mapped < second_count) {
            result[first_map[i]] = second_map[mapped];
        }
    }
    return result;
}

inline Permutation seaweed_doubling(const Permutation& permutation) {
    int n = int(permutation.size());
    if (n == 1) return Permutation(1, none);

    int middle = n / 2;
    Permutation lower;
    Permutation upper;
    Permutation lower_map;
    Permutation upper_map;
    for (int i = 0; i < n; i++) {
        int value = permutation[i];
        if (value < middle) {
            lower.push_back(value);
            lower_map.push_back(i);
        } else {
            upper.push_back(value - middle);
            upper_map.push_back(i);
        }
    }

    lower = seaweed_doubling(lower);
    upper = seaweed_doubling(upper);
    Permutation lower_padded(n);
    Permutation upper_padded(n);
    std::iota(lower_padded.begin(), lower_padded.end(), 0);
    std::iota(upper_padded.begin(), upper_padded.end(), 0);

    for (int i = 0; i < middle; i++) {
        lower_padded[lower_map[i]] =
            lower[i] == none ? none : lower_map[lower[i]];
    }
    for (int i = 0; middle + i < n; i++) {
        upper_padded[upper_map[i]] =
            upper[i] == none ? none : upper_map[upper[i]];
    }
    return subunit_monge_product(
        std::move(lower_padded),
        std::move(upper_padded)
    );
}

template <class T>
Permutation make_permutation(const std::vector<T>& values) {
    int n = int(values.size());
    Permutation order(n);
    std::iota(order.begin(), order.end(), 0);
    std::sort(order.begin(), order.end(), [&](int first, int second) {
        if (values[first] < values[second]) return true;
        if (values[second] < values[first]) return false;
        return first > second;
    });

    Permutation permutation(n);
    for (int rank = 0; rank < n; rank++) {
        permutation[order[rank]] = rank;
    }
    return permutation;
}

}  // namespace range_lis_query_internal

// Static strict-LIS lengths for arbitrary subarrays.
template <class T>
struct RangeLisQuery {
   private:
    int _n;
    WaveletMatrix<int> _seaweed;

   public:
    RangeLisQuery() : _n(0), _seaweed() {}

    explicit RangeLisQuery(const std::vector<T>& values)
        : _n(int(values.size())), _seaweed() {
        if (_n == 0) return;
        std::vector<int> row = range_lis_query_internal::seaweed_doubling(
            range_lis_query_internal::make_permutation(values)
        );
        for (int& value : row) {
            if (value == range_lis_query_internal::none) value = _n;
        }
        _seaweed = WaveletMatrix<int>(row);
    }

    int size() const {
        return _n;
    }

    bool empty() const {
        return _n == 0;
    }

    // Returns the strict LIS length of values[left, right).
    int query(int left, int right) const {
        assert(0 <= left && left <= right && right <= _n);
        if (left == right) return 0;
        return right - left - _seaweed.range_freq(left, _n, right);
    }

    int lis_length(int left, int right) const {
        return query(left, right);
    }
};

}  // namespace ds
}  // namespace m1une

#endif  // M1UNE_DS_RANGE_QUERY_RANGE_LIS_QUERY_HPP
#line 1 "ds/range_query/range_lis_query.hpp"



#line 1 "ds/wavelet_matrix/wavelet_matrix.hpp"



#include <algorithm>
#include <bit>
#include <cassert>
#include <concepts>
#include <cstdint>
#include <limits>
#include <optional>
#include <type_traits>
#include <utility>
#include <vector>

#if defined(__AVX2__) || defined(__BMI2__)
#include <immintrin.h>
#endif

namespace m1une {
namespace ds {

// A static wavelet matrix for integral values.
template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
struct WaveletMatrix {
    using value_type = T;
    using unsigned_type = std::make_unsigned_t<T>;

   private:
    static constexpr int value_bit_width =
        std::numeric_limits<unsigned_type>::digits;
    static constexpr unsigned_type sign_mask = [] {
        if constexpr (std::signed_integral<T>) {
            return unsigned_type(1) << (value_bit_width - 1);
        } else {
            return unsigned_type(0);
        }
    }();

    struct BitVector {
        std::vector<std::uint64_t> bits;
        std::vector<int> prefix;

        BitVector() = default;

        explicit BitVector(int n)
            : bits(((std::size_t(n) + 63) >> 6) + 1, 0),
              prefix(bits.size(), 0) {}

        void build() {
            for (std::size_t i = 0; i + 1 < bits.size(); i++) {
                prefix[i + 1] = prefix[i] + std::popcount(bits[i]);
            }
        }

        bool get(int p) const {
            return (bits[std::size_t(p) >> 6] >> (p & 63)) & 1;
        }

        int rank1(int r) const {
            std::size_t word = std::size_t(r) >> 6;
            int offset = r & 63;
            int result = prefix[word];
#if defined(__BMI2__)
            result += std::popcount(
                _bzhi_u64(bits[word], static_cast<unsigned int>(offset))
            );
#else
            if (offset != 0) {
                result += std::popcount(
                    bits[word] & ((std::uint64_t(1) << offset) - 1)
                );
            }
#endif
            return result;
        }
    };

    int _n;
    int _log;
    unsigned_type _key_prefix;
    unsigned_type _min_key;
    unsigned_type _max_key;
    std::vector<BitVector> _matrix;
    std::vector<int> _zero_count;

    static unsigned_type encode(T value) {
        unsigned_type bits;
        if constexpr (std::signed_integral<T>) {
            bits = std::bit_cast<unsigned_type>(value);
        } else {
            bits = value;
        }
        return bits ^ sign_mask;
    }

    static T decode(unsigned_type key) {
        unsigned_type bits = key ^ sign_mask;
        if constexpr (std::signed_integral<T>) {
            return std::bit_cast<T>(bits);
        } else {
            return bits;
        }
    }

    bool bit(unsigned_type value, int level) const {
        return (value >> (_log - 1 - level)) & unsigned_type(1);
    }

    static std::uint64_t extract_bits(
        const unsigned_type* values,
        int count,
        int shift
    ) {
        std::uint64_t result = 0;
        int i = 0;
#if defined(__AVX2__)
        if constexpr (sizeof(unsigned_type) == 8) {
            __m128i left = _mm_cvtsi32_si128(63 - shift);
            for (; i + 4 <= count; i += 4) {
                __m256i data = _mm256_loadu_si256(
                    reinterpret_cast<const __m256i*>(values + i)
                );
                data = _mm256_sll_epi64(data, left);
                int mask = _mm256_movemask_pd(_mm256_castsi256_pd(data));
                result |= std::uint64_t(mask) << i;
            }
        } else if constexpr (sizeof(unsigned_type) == 4) {
            __m128i left = _mm_cvtsi32_si128(31 - shift);
            for (; i + 8 <= count; i += 8) {
                __m256i data = _mm256_loadu_si256(
                    reinterpret_cast<const __m256i*>(values + i)
                );
                data = _mm256_sll_epi32(data, left);
                int mask = _mm256_movemask_ps(_mm256_castsi256_ps(data));
                result |= std::uint64_t(mask) << i;
            }
        }
#endif
        for (; i < count; i++) {
            result |= std::uint64_t((values[i] >> shift) & unsigned_type(1))
                      << i;
        }
        return result;
    }

    int count_less_encoded(int l, int r, unsigned_type upper) const {
        if (_n == 0 || upper <= _min_key) return 0;
        if (upper > _max_key) return r - l;

        int result = 0;
        for (int level = 0; level < _log; level++) {
            int l1 = _matrix[level].rank1(l);
            int r1 = _matrix[level].rank1(r);
            if (bit(upper, level)) {
                result += (r - l) - (r1 - l1);
                l = _zero_count[level] + l1;
                r = _zero_count[level] + r1;
            } else {
                l -= l1;
                r -= r1;
            }
        }
        return result;
    }

   public:
    WaveletMatrix()
        : _n(0),
          _log(0),
          _key_prefix(0),
          _min_key(0),
          _max_key(0) {}

    explicit WaveletMatrix(const std::vector<T>& values)
        : _n(int(values.size())),
          _log(0),
          _key_prefix(0),
          _min_key(0),
          _max_key(0) {
        std::vector<unsigned_type> current(_n);
        std::vector<unsigned_type> next(_n);
        for (int i = 0; i < _n; i++) current[i] = encode(values[i]);
        if (_n == 0) return;

        _min_key = current[0];
        _max_key = current[0];
        for (unsigned_type key : current) {
            if (key < _min_key) _min_key = key;
            if (_max_key < key) _max_key = key;
        }
        _log = int(std::bit_width(unsigned_type(_min_key ^ _max_key)));
        if (_log != value_bit_width) {
            _key_prefix = unsigned_type((_min_key >> _log) << _log);
        }
        _zero_count.assign(_log, 0);

        _matrix.reserve(_log);
        for (int level = 0; level < _log; level++) {
            _matrix.emplace_back(_n);
            BitVector& bit_vector = _matrix.back();
            int shift = _log - 1 - level;
            int zeros = 0;
            for (int base = 0; base < _n; base += 64) {
                int count = std::min(64, _n - base);
                std::uint64_t word = extract_bits(
                    current.data() + base,
                    count,
                    shift
                );
                bit_vector.bits[std::size_t(base) >> 6] = word;
                zeros += count - std::popcount(word);
            }
            bit_vector.build();

            _zero_count[level] = zeros;
            int zero_pos = 0;
            int one_pos = zeros;
            for (int base = 0; base < _n; base += 64) {
                int count = std::min(64, _n - base);
                std::uint64_t ones = bit_vector.bits[std::size_t(base) >> 6];
                std::uint64_t valid = count == 64
                                          ? ~std::uint64_t(0)
                                          : (std::uint64_t(1) << count) - 1;
                std::uint64_t zeroes = (~ones) & valid;
                while (zeroes != 0) {
                    int offset = std::countr_zero(zeroes);
                    next[zero_pos++] = current[base + offset];
                    zeroes &= zeroes - 1;
                }
                while (ones != 0) {
                    int offset = std::countr_zero(ones);
                    next[one_pos++] = current[base + offset];
                    ones &= ones - 1;
                }
            }
            current.swap(next);
        }
    }

    int size() const {
        return _n;
    }

    bool empty() const {
        return _n == 0;
    }

    T access(int p) const {
        assert(0 <= p && p < _n);
        unsigned_type key = _key_prefix;
        for (int level = 0; level < _log; level++) {
            int ones_before = _matrix[level].rank1(p);
            bool one = _matrix[level].get(p);
            if (one) {
                key |= unsigned_type(1) << (_log - 1 - level);
                p = _zero_count[level] + ones_before;
            } else {
                p -= ones_before;
            }
        }
        return decode(key);
    }

    T operator[](int p) const {
        return access(p);
    }

    int rank(T value, int r) const {
        assert(0 <= r && r <= _n);
        return rank(value, 0, r);
    }

    int rank(T value, int l, int r) const {
        assert(0 <= l && l <= r && r <= _n);
        unsigned_type key = encode(value);
        if (_n == 0 || key < _min_key || _max_key < key) return 0;
        for (int level = 0; level < _log; level++) {
            int l1 = _matrix[level].rank1(l);
            int r1 = _matrix[level].rank1(r);
            if (bit(key, level)) {
                l = _zero_count[level] + l1;
                r = _zero_count[level] + r1;
            } else {
                l -= l1;
                r -= r1;
            }
        }
        return r - l;
    }

    T kth_smallest(int l, int r, int k) const {
        assert(0 <= l && l <= r && r <= _n);
        assert(0 <= k && k < r - l);
        unsigned_type key = _key_prefix;
        for (int level = 0; level < _log; level++) {
            int l1 = _matrix[level].rank1(l);
            int r1 = _matrix[level].rank1(r);
            int l0 = l - l1;
            int r0 = r - r1;
            int zeros = r0 - l0;
            if (k < zeros) {
                l = l0;
                r = r0;
            } else {
                k -= zeros;
                key |= unsigned_type(1) << (_log - 1 - level);
                l = _zero_count[level] + l1;
                r = _zero_count[level] + r1;
            }
        }
        return decode(key);
    }

    T kth_largest(int l, int r, int k) const {
        assert(0 <= l && l <= r && r <= _n);
        assert(0 <= k && k < r - l);
        return kth_smallest(l, r, r - l - 1 - k);
    }

    int range_freq(int l, int r, T upper) const {
        assert(0 <= l && l <= r && r <= _n);
        return count_less_encoded(l, r, encode(upper));
    }

    int range_freq(int l, int r, T lower, T upper) const {
        assert(0 <= l && l <= r && r <= _n);
        if (upper <= lower) return 0;
        return range_freq(l, r, upper) - range_freq(l, r, lower);
    }

    std::optional<T> prev_value(int l, int r, T upper) const {
        assert(0 <= l && l <= r && r <= _n);
        int count = range_freq(l, r, upper);
        if (count == 0) return std::nullopt;
        return kth_smallest(l, r, count - 1);
    }

    std::optional<T> next_value(int l, int r, T lower) const {
        assert(0 <= l && l <= r && r <= _n);
        int count = range_freq(l, r, lower);
        if (count == r - l) return std::nullopt;
        return kth_smallest(l, r, count);
    }
};

}  // namespace ds
}  // namespace m1une


#line 5 "ds/range_query/range_lis_query.hpp"

#line 8 "ds/range_query/range_lis_query.hpp"
#include <numeric>
#line 11 "ds/range_query/range_lis_query.hpp"

namespace m1une {
namespace ds {

namespace range_lis_query_internal {

constexpr int none = -1;
using Permutation = std::vector<int>;
using Iterator = Permutation::iterator;

inline Permutation inverse(const Permutation& permutation) {
    int n = int(permutation.size());
    Permutation result(n, none);
    for (int i = 0; i < n; i++) {
        if (permutation[i] != none) result[permutation[i]] = i;
    }
    return result;
}

// Distance multiplication of two unit-Monge matrices, represented by their
// permutations. `workspace` must have the size used by subunit_monge_product.
inline void unit_monge_product(
    int n,
    Iterator workspace,
    Iterator first,
    Iterator second
) {
    if (n == 1) {
        workspace[0] = 0;
        return;
    }

    Iterator result_row = workspace;
    workspace += n;
    Iterator result_column = workspace;
    workspace += n;

    auto map_half = [=](int length, const auto& belongs, const auto& map) {
        Iterator first_half = workspace;
        Iterator first_position = workspace + length;
        Iterator second_half = workspace + 2 * length;
        Iterator second_position = workspace + 3 * length;

        auto split = [=](Iterator source, Iterator half, Iterator position) {
            for (int i = 0; i < n; i++) {
                if (belongs(source[i])) {
                    *half++ = map(source[i]);
                    *position++ = i;
                }
            }
        };
        split(first, first_half, first_position);
        split(second, second_half, second_position);

        Iterator product = workspace + 4 * length;
        unit_monge_product(
            length,
            product,
            first_half,
            second_half
        );
        for (int i = 0; i < length; i++) {
            int row = first_position[i];
            int column = second_position[product[i]];
            result_row[row] = column;
            result_column[column] = row;
        }
    };

    int middle = n / 2;
    map_half(
        middle,
        [middle](int value) { return value < middle; },
        [](int value) { return value; }
    );
    map_half(
        n - middle,
        [middle](int value) { return value >= middle; },
        [middle](int value) { return value - middle; }
    );

    struct DiagonalIterator {
        int delta = 0;
        int column = 0;
    };

    int row = n;
    auto move_right = [&](DiagonalIterator& iterator) {
        if (second[iterator.column] < middle) {
            if (result_column[iterator.column] >= row) iterator.delta++;
        } else {
            if (result_column[iterator.column] < row) iterator.delta++;
        }
        iterator.column++;
    };
    auto move_up = [&](DiagonalIterator& iterator) {
        if (first[row] < middle) {
            if (result_row[row] >= iterator.column) iterator.delta--;
        } else {
            if (result_row[row] < iterator.column) iterator.delta--;
        }
    };

    DiagonalIterator negative;
    DiagonalIterator positive;
    while (row != 0) {
        while (positive.column != n) {
            DiagonalIterator next = positive;
            move_right(next);
            if (next.delta != 0) break;
            positive = next;
        }
        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 int workspace_size(int n) {
    int result = 1;
    while (n > 1) {
        result += 2 * n;
        n = (n + 1) / 2;
        result += 4 * n;
    }
    return result;
}

inline Permutation subunit_monge_product(
    Permutation first,
    Permutation second
) {
    int n = int(first.size());
    Permutation first_inverse = inverse(first);
    Permutation second_inverse = inverse(second);
    std::swap(second, second_inverse);

    Permutation first_map;
    Permutation 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 missing = 0;
    for (int i = 0; i < n; i++) {
        if (first_inverse[i] == none) first[missing++] = i;
    }
    for (int i = 0; i < n; i++) {
        if (second[i] != none) {
            second[int(second_map.size())] = second[i];
            second_map.push_back(i);
        }
    }
    missing = int(second_map.size());
    for (int i = 0; i < n; i++) {
        if (second_inverse[i] == none) second[missing++] = i;
    }

    Permutation workspace(workspace_size(n));
    unit_monge_product(
        n,
        workspace.begin(),
        first.begin(),
        second.begin()
    );

    Permutation result(n, none);
    int first_count = int(first_map.size());
    int second_count = int(second_map.size());
    for (int i = 0; i < first_count; i++) {
        int mapped = workspace[n - first_count + i];
        if (mapped < second_count) {
            result[first_map[i]] = second_map[mapped];
        }
    }
    return result;
}

inline Permutation seaweed_doubling(const Permutation& permutation) {
    int n = int(permutation.size());
    if (n == 1) return Permutation(1, none);

    int middle = n / 2;
    Permutation lower;
    Permutation upper;
    Permutation lower_map;
    Permutation upper_map;
    for (int i = 0; i < n; i++) {
        int value = permutation[i];
        if (value < middle) {
            lower.push_back(value);
            lower_map.push_back(i);
        } else {
            upper.push_back(value - middle);
            upper_map.push_back(i);
        }
    }

    lower = seaweed_doubling(lower);
    upper = seaweed_doubling(upper);
    Permutation lower_padded(n);
    Permutation upper_padded(n);
    std::iota(lower_padded.begin(), lower_padded.end(), 0);
    std::iota(upper_padded.begin(), upper_padded.end(), 0);

    for (int i = 0; i < middle; i++) {
        lower_padded[lower_map[i]] =
            lower[i] == none ? none : lower_map[lower[i]];
    }
    for (int i = 0; middle + i < n; i++) {
        upper_padded[upper_map[i]] =
            upper[i] == none ? none : upper_map[upper[i]];
    }
    return subunit_monge_product(
        std::move(lower_padded),
        std::move(upper_padded)
    );
}

template <class T>
Permutation make_permutation(const std::vector<T>& values) {
    int n = int(values.size());
    Permutation order(n);
    std::iota(order.begin(), order.end(), 0);
    std::sort(order.begin(), order.end(), [&](int first, int second) {
        if (values[first] < values[second]) return true;
        if (values[second] < values[first]) return false;
        return first > second;
    });

    Permutation permutation(n);
    for (int rank = 0; rank < n; rank++) {
        permutation[order[rank]] = rank;
    }
    return permutation;
}

}  // namespace range_lis_query_internal

// Static strict-LIS lengths for arbitrary subarrays.
template <class T>
struct RangeLisQuery {
   private:
    int _n;
    WaveletMatrix<int> _seaweed;

   public:
    RangeLisQuery() : _n(0), _seaweed() {}

    explicit RangeLisQuery(const std::vector<T>& values)
        : _n(int(values.size())), _seaweed() {
        if (_n == 0) return;
        std::vector<int> row = range_lis_query_internal::seaweed_doubling(
            range_lis_query_internal::make_permutation(values)
        );
        for (int& value : row) {
            if (value == range_lis_query_internal::none) value = _n;
        }
        _seaweed = WaveletMatrix<int>(row);
    }

    int size() const {
        return _n;
    }

    bool empty() const {
        return _n == 0;
    }

    // Returns the strict LIS length of values[left, right).
    int query(int left, int right) const {
        assert(0 <= left && left <= right && right <= _n);
        if (left == right) return 0;
        return right - left - _seaweed.range_freq(left, _n, right);
    }

    int lis_length(int left, int right) const {
        return query(left, right);
    }
};

}  // namespace ds
}  // namespace m1une
Back to top page