m1une's library

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

View on GitHub

:warning: Algorithms All
(algo/all.hpp)

Overview

algo/all.hpp includes one-shot, domain-neutral algorithms. The public namespace is m1une::algo; subdirectories are only browsing categories.

Convex optimization helpers live under convex/.

If a header builds an object and then answers repeated queries, it belongs in ds/ instead, even when it is static. For example, cumulative sums live in ds/range_query/.

Included Headers

Header Contents
algo/sequence/all.hpp Sequence and array algorithms such as LIS, inversion count, non-adjacent selection, run-length encoding, and subset sum.
algo/search/all.hpp Search-over-answer and unimodal optimization helpers.
algo/offline/all.hpp Offline query processing such as Mo’s algorithm.
algo/enumeration/all.hpp Combinatorial traversal helpers such as Gray-code enumeration.
algo/dp/all.hpp Domain-neutral DP helpers such as knapsack routines.

Depends on

Code

#ifndef M1UNE_ALGO_ALL_HPP
#define M1UNE_ALGO_ALL_HPP 1

#include "dp/all.hpp"
#include "enumeration/all.hpp"
#include "offline/all.hpp"
#include "search/all.hpp"
#include "sequence/all.hpp"

#endif  // M1UNE_ALGO_ALL_HPP
#line 1 "algo/all.hpp"



#line 1 "algo/dp/all.hpp"



#line 1 "algo/dp/knapsack.hpp"



#include <algorithm>
#include <bit>
#include <cassert>
#include <cstddef>
#include <deque>
#include <limits>
#include <vector>

namespace m1une {
namespace algo {

namespace internal {

using SubsetSumWord = unsigned long long;

inline std::vector<SubsetSumWord> subset_sum_reachability_bits(
    const std::vector<int>& weights,
    int limit
) {
    assert(0 <= limit);
    using Word = SubsetSumWord;
    constexpr int word_bits = std::numeric_limits<Word>::digits;

    const std::size_t bit_count = std::size_t(limit) + 1;
    std::vector<Word> bits((bit_count + word_bits - 1) / word_bits, Word(0));
    bits[0] = Word(1);

    auto trim = [&]() {
        const int extra = int(bit_count % word_bits);
        if (extra != 0) {
            bits.back() &= (Word(1) << extra) - Word(1);
        }
    };

    for (int weight : weights) {
        assert(0 <= weight);
        if (weight == 0 || limit < weight) continue;

        const std::size_t word_shift = std::size_t(weight / word_bits);
        const int bit_shift = weight % word_bits;
        for (std::size_t i = bits.size() - word_shift; i-- > 0;) {
            const Word source = bits[i];
            if (source == Word(0)) continue;
            const std::size_t target = i + word_shift;
            bits[target] |= source << bit_shift;
            if (bit_shift != 0 && target + 1 < bits.size()) {
                bits[target + 1] |= source >> (word_bits - bit_shift);
            }
        }
        trim();
    }
    return bits;
}

}  // namespace internal

inline std::vector<char> subset_sum_reachable(const std::vector<int>& weights, int limit) {
    using Word = internal::SubsetSumWord;
    constexpr int word_bits = std::numeric_limits<Word>::digits;
    const std::vector<Word> bits =
        internal::subset_sum_reachability_bits(weights, limit);

    std::vector<char> reachable(std::size_t(limit) + 1, 0);
    for (int sum = 0; sum <= limit; ++sum) {
        reachable[sum] = char((bits[std::size_t(sum / word_bits)] >> (sum % word_bits)) & Word(1));
    }
    return reachable;
}

// Returns the maximum subset sum not exceeding limit.
inline int subset_sum_max_value(const std::vector<int>& weights, int limit) {
    using Word = internal::SubsetSumWord;
    constexpr int word_bits = std::numeric_limits<Word>::digits;
    const std::vector<Word> bits =
        internal::subset_sum_reachability_bits(weights, limit);

    for (std::size_t i = bits.size(); i-- > 0;) {
        if (bits[i] != Word(0)) {
            return int(i * word_bits + std::bit_width(bits[i]) - 1);
        }
    }
    return 0;
}

template <typename Value = long long>
std::vector<Value> zero_one_knapsack_max_value(
    const std::vector<int>& weights,
    const std::vector<Value>& values,
    int capacity,
    Value neg_inf = std::numeric_limits<Value>::lowest() / Value(4)
) {
    assert(weights.size() == values.size());
    assert(0 <= capacity);

    std::vector<Value> dp(std::size_t(capacity) + 1, neg_inf);
    dp[0] = Value{};
    for (std::size_t item = 0; item < weights.size(); ++item) {
        const int weight = weights[item];
        assert(0 <= weight);
        for (int current = capacity; weight <= current; --current) {
            if (dp[current - weight] == neg_inf) continue;
            dp[current] = std::max(dp[current], dp[current - weight] + values[item]);
        }
    }

    for (int current = 1; current <= capacity; ++current) {
        dp[current] = std::max(dp[current], dp[current - 1]);
    }
    return dp;
}

template <typename Value = long long>
std::vector<Value> bounded_knapsack_max_value(
    const std::vector<int>& weights,
    const std::vector<Value>& values,
    const std::vector<int>& counts,
    int capacity,
    Value neg_inf = std::numeric_limits<Value>::lowest() / Value(4)
) {
    assert(weights.size() == values.size());
    assert(weights.size() == counts.size());
    assert(0 <= capacity);

    std::vector<Value> dp(std::size_t(capacity) + 1, neg_inf);
    dp[0] = Value{};

    for (std::size_t item = 0; item < weights.size(); ++item) {
        const int weight = weights[item];
        const Value value = values[item];
        const int count = counts[item];
        assert(0 <= weight);
        assert(0 <= count);
        if (count == 0) continue;

        if (weight == 0) {
            if (Value{} < value) {
                const Value gain = value * Value(count);
                for (Value& current : dp) {
                    if (current != neg_inf) current += gain;
                }
            }
            continue;
        }

        std::vector<Value> next = dp;
        for (int residue = 0; residue < weight && residue <= capacity; ++residue) {
            std::deque<int> indices;
            std::deque<Value> bases;
            int k = 0;
            for (int current = residue; current <= capacity; current += weight, ++k) {
                if (dp[current] != neg_inf) {
                    const Value base = dp[current] - Value(k) * value;
                    while (!bases.empty() && bases.back() <= base) {
                        bases.pop_back();
                        indices.pop_back();
                    }
                    bases.push_back(base);
                    indices.push_back(k);
                }

                while (!indices.empty() && indices.front() < k - count) {
                    indices.pop_front();
                    bases.pop_front();
                }
                if (!bases.empty()) {
                    next[current] = std::max(next[current], bases.front() + Value(k) * value);
                }
            }
        }
        dp.swap(next);
    }

    for (int current = 1; current <= capacity; ++current) {
        dp[current] = std::max(dp[current], dp[current - 1]);
    }
    return dp;
}

template <typename Weight = long long>
std::vector<Weight> zero_one_knapsack_min_weight_for_value(
    const std::vector<Weight>& weights,
    const std::vector<int>& values,
    int value_limit,
    Weight inf = std::numeric_limits<Weight>::max() / Weight(4)
) {
    assert(weights.size() == values.size());
    assert(0 <= value_limit);

    std::vector<Weight> dp(std::size_t(value_limit) + 1, inf);
    dp[0] = Weight{};
    for (std::size_t item = 0; item < weights.size(); ++item) {
        assert(Weight{} <= weights[item]);
        assert(0 <= values[item]);
        for (int value = value_limit; values[item] <= value; --value) {
            if (dp[value - values[item]] == inf) continue;
            dp[value] = std::min(dp[value], dp[value - values[item]] + weights[item]);
        }
    }
    return dp;
}

}  // namespace algo
}  // namespace m1une


#line 5 "algo/dp/all.hpp"


#line 1 "algo/enumeration/all.hpp"



#line 1 "algo/enumeration/combination.hpp"



#line 5 "algo/enumeration/combination.hpp"
#include <concepts>
#include <cstdint>
#line 8 "algo/enumeration/combination.hpp"
#include <type_traits>

namespace m1une {
namespace algo {

namespace internal {

template <std::unsigned_integral UInt>
requires(!std::same_as<std::remove_cv_t<UInt>, bool>)
UInt combination_low_bits(int bit_count) {
    constexpr int digits = std::numeric_limits<UInt>::digits;
    assert(0 <= bit_count && bit_count <= digits);
    if (bit_count == digits) return ~UInt(0);
    return (UInt(1) << bit_count) - UInt(1);
}

}  // namespace internal

template <std::unsigned_integral UInt = std::uint64_t>
requires(!std::same_as<std::remove_cv_t<UInt>, bool>)
UInt first_combination_mask(int bit_count, int choose) {
    constexpr int digits = std::numeric_limits<UInt>::digits;
    assert(0 <= choose && choose <= bit_count && bit_count <= digits);
    if (choose == 0) return UInt(0);
    if (choose == bit_count) return internal::combination_low_bits<UInt>(bit_count);
    return (UInt(1) << choose) - UInt(1);
}

template <std::unsigned_integral UInt>
requires(!std::same_as<std::remove_cv_t<UInt>, bool>)
bool next_combination_mask(UInt& mask, int bit_count) {
    const UInt universe = internal::combination_low_bits<UInt>(bit_count);
    assert((mask & ~universe) == 0);
    if (mask == 0) return false;

    const UInt lowest = mask & (~mask + UInt(1));
    const UInt ripple = mask + lowest;
    if (ripple == 0 || (ripple & ~universe) != 0) return false;

    const UInt next = (((ripple ^ mask) >> 2) / lowest) | ripple;
    if ((next & ~universe) != 0) return false;
    mask = next;
    return true;
}

template <std::unsigned_integral UInt = std::uint64_t, class F>
requires(!std::same_as<std::remove_cv_t<UInt>, bool>)
void for_each_combination_mask(int bit_count, int choose, F f) {
    constexpr int digits = std::numeric_limits<UInt>::digits;
    assert(0 <= choose && choose <= bit_count && bit_count <= digits);
    UInt mask = first_combination_mask<UInt>(bit_count, choose);
    while (true) {
        f(mask);
        if (!next_combination_mask(mask, bit_count)) break;
    }
}

}  // namespace algo
}  // namespace m1une


#line 1 "algo/enumeration/gray_code.hpp"



#line 11 "algo/enumeration/gray_code.hpp"

namespace m1une {
namespace algo {

// Converts a binary value to its binary-reflected Gray code.
template <std::unsigned_integral UInt>
requires(!std::same_as<std::remove_cv_t<UInt>, bool>)
constexpr UInt gray_encode(UInt value) noexcept {
    return value ^ (value >> 1);
}

// Converts a binary-reflected Gray code to the corresponding binary value.
template <std::unsigned_integral UInt>
requires(!std::same_as<std::remove_cv_t<UInt>, bool>)
constexpr UInt gray_decode(UInt code) noexcept {
    for (int shift = 1; shift < std::numeric_limits<UInt>::digits;
         shift <<= 1) {
        code ^= code >> shift;
    }
    return code;
}

// Returns all bit_count-bit binary-reflected Gray codes in traversal order.
template <std::unsigned_integral UInt = std::uint64_t>
requires(!std::same_as<std::remove_cv_t<UInt>, bool>)
std::vector<UInt> gray_code_sequence(int bit_count) {
    constexpr int uint_digits = std::numeric_limits<UInt>::digits;
    constexpr int size_digits = std::numeric_limits<std::size_t>::digits;
    assert(0 <= bit_count);
    assert(bit_count <= uint_digits);
    assert(bit_count < size_digits);
    if (bit_count < 0 || uint_digits < bit_count || size_digits <= bit_count) {
        return {};
    }

    const std::size_t size = std::size_t(1) << bit_count;
    std::vector<UInt> result(size);
    for (std::size_t index = 0; index < size; ++index) {
        result[index] = gray_encode(static_cast<UInt>(index));
    }
    return result;
}

}  // namespace algo
}  // namespace m1une


#line 1 "algo/enumeration/permutation_lexicographical_order.hpp"



#line 8 "algo/enumeration/permutation_lexicographical_order.hpp"
#include <optional>
#line 10 "algo/enumeration/permutation_lexicographical_order.hpp"
#include <utility>
#line 12 "algo/enumeration/permutation_lexicographical_order.hpp"

namespace m1une {
namespace algo {

namespace internal {

struct PermutationOrderFenwick {
    std::vector<int> data;

    explicit PermutationOrderFenwick(int size) : data(size + 1) {}

    void add(int index, int value) {
        for (index++; index < int(data.size()); index += index & -index) {
            data[index] += value;
        }
    }

    int prefix_sum(int right) const {
        int result = 0;
        for (; 0 < right; right -= right & -right) result += data[right];
        return result;
    }

    int kth(int order) const {
        int index = 0;
        int accumulated = 0;
        int step = 1;
        while (step < int(data.size())) step <<= 1;
        for (; 0 < step; step >>= 1) {
            const int next = index + step;
            if (next < int(data.size()) &&
                accumulated + data[next] <= order) {
                index = next;
                accumulated += data[next];
            }
        }
        return index;
    }
};

}  // namespace internal

// Returns the zero-based lexicographical rank of a permutation of [0, n).
// Returns nullopt when the sequence is invalid or the rank does not fit in UInt.
template <
    std::unsigned_integral UInt = std::uint64_t,
    class Permutation
>
requires(!std::same_as<std::remove_cv_t<UInt>, bool>)
std::optional<UInt> checked_permutation_lexicographical_rank(
    const Permutation& permutation
) {
    const int size = int(permutation.size());
    internal::PermutationOrderFenwick fenwick(size);
    for (int value = 0; value < size; value++) fenwick.add(value, 1);

    UInt rank = 0;
    constexpr UInt limit = std::numeric_limits<UInt>::max();
    for (int index = 0; index < size; index++) {
        const auto& value_reference = permutation[index];
        using Value = std::remove_cvref_t<decltype(value_reference)>;
        static_assert(std::integral<Value>);
        static_assert(!std::same_as<Value, bool>);

        if (std::cmp_less(value_reference, 0) ||
            std::cmp_greater_equal(value_reference, size)) {
            return std::nullopt;
        }
        const int value = int(value_reference);
        if (fenwick.prefix_sum(value + 1) == fenwick.prefix_sum(value)) {
            return std::nullopt;
        }

        const std::uintmax_t smaller = fenwick.prefix_sum(value);
        const std::uintmax_t remaining = size - index;
        if (smaller > std::uintmax_t(limit) ||
            std::uintmax_t(rank) >
                (std::uintmax_t(limit) - smaller) / remaining) {
            return std::nullopt;
        }
        rank = UInt(std::uintmax_t(rank) * remaining + smaller);
        fenwick.add(value, -1);
    }
    return rank;
}

// Every value must occur exactly once in [0, n), and the rank must fit in UInt.
template <
    std::unsigned_integral UInt = std::uint64_t,
    class Permutation
>
requires(!std::same_as<std::remove_cv_t<UInt>, bool>)
UInt permutation_lexicographical_rank(const Permutation& permutation) {
    const std::optional<UInt> result =
        checked_permutation_lexicographical_rank<UInt>(permutation);
    assert(result.has_value());
    return result.value_or(UInt(0));
}

// Returns the permutation of [0, size) with the given zero-based rank.
// Returns nullopt when size is negative or rank is at least size factorial.
template <std::unsigned_integral UInt = std::uint64_t>
requires(!std::same_as<std::remove_cv_t<UInt>, bool>)
std::optional<std::vector<int>> checked_kth_lexicographical_permutation(
    int size,
    UInt rank
) {
    if (size < 0) return std::nullopt;

    std::vector<int> lehmer_code(size);
    UInt remaining_rank = rank;
    for (int base = 1; base <= size && remaining_rank != 0; base++) {
        lehmer_code[size - base] = int(remaining_rank % UInt(base));
        remaining_rank /= UInt(base);
    }
    if (remaining_rank != 0) return std::nullopt;

    internal::PermutationOrderFenwick fenwick(size);
    for (int value = 0; value < size; value++) fenwick.add(value, 1);

    std::vector<int> permutation(size);
    for (int index = 0; index < size; index++) {
        const int value = fenwick.kth(lehmer_code[index]);
        permutation[index] = value;
        fenwick.add(value, -1);
    }
    return permutation;
}

// Rank must be less than size factorial.
template <std::unsigned_integral UInt = std::uint64_t>
requires(!std::same_as<std::remove_cv_t<UInt>, bool>)
std::vector<int> kth_lexicographical_permutation(int size, UInt rank) {
    std::optional<std::vector<int>> result =
        checked_kth_lexicographical_permutation(size, rank);
    assert(result.has_value());
    return result.value_or(std::vector<int>());
}

}  // namespace algo
}  // namespace m1une


#line 1 "algo/enumeration/segtree_range.hpp"



#line 10 "algo/enumeration/segtree_range.hpp"

namespace m1une {
namespace algo {

// Splits [left, right) into maximal segment-tree ranges from left to right.
template <std::integral Int>
requires(!std::same_as<std::remove_cv_t<Int>, bool>)
std::vector<std::pair<Int, Int>> split_segtree_range(Int left, Int right) {
    if constexpr (std::signed_integral<Int>) assert(Int(0) <= left);
    assert(left <= right);
    if constexpr (std::signed_integral<Int>) {
        if (left < 0) return {};
    }
    if (right < left) return {};

    using UInt = std::make_unsigned_t<Int>;
    UInt position = static_cast<UInt>(left);
    const UInt end = static_cast<UInt>(right);
    std::vector<std::pair<Int, Int>> result;
    if (position == end) return result;
    result.reserve(2 * std::bit_width(end - position));

    while (position < end) {
        UInt length = std::bit_floor(end - position);
        if (position != 0) {
            const UInt alignment = position & (~position + UInt(1));
            if (alignment < length) length = alignment;
        }
        const UInt next = position + length;
        result.emplace_back(
            static_cast<Int>(position), static_cast<Int>(next)
        );
        position = next;
    }
    return result;
}

}  // namespace algo
}  // namespace m1une


#line 1 "algo/enumeration/submask.hpp"



#line 8 "algo/enumeration/submask.hpp"

namespace m1une {
namespace algo {

namespace internal {

template <std::unsigned_integral UInt>
requires(!std::same_as<std::remove_cv_t<UInt>, bool>)
UInt submask_low_bits(int bit_count) {
    constexpr int digits = std::numeric_limits<UInt>::digits;
    assert(0 <= bit_count && bit_count <= digits);
    if (bit_count == digits) return ~UInt(0);
    return (UInt(1) << bit_count) - UInt(1);
}

}  // namespace internal

template <std::unsigned_integral UInt, class F>
requires(!std::same_as<std::remove_cv_t<UInt>, bool>)
void for_each_submask(UInt mask, F f) {
    UInt submask = mask;
    while (true) {
        f(submask);
        if (submask == 0) break;
        submask = (submask - 1) & mask;
    }
}

template <std::unsigned_integral UInt, class F>
requires(!std::same_as<std::remove_cv_t<UInt>, bool>)
void for_each_nonzero_submask(UInt mask, F f) {
    for (UInt submask = mask; submask != 0; submask = (submask - 1) & mask) {
        f(submask);
    }
}

template <std::unsigned_integral UInt, class F>
requires(!std::same_as<std::remove_cv_t<UInt>, bool>)
void for_each_supermask(UInt mask, int bit_count, F f) {
    const UInt universe = internal::submask_low_bits<UInt>(bit_count);
    assert((mask & ~universe) == 0);
    const UInt free_bits = universe ^ mask;
    for_each_submask(free_bits, [&](UInt added_bits) {
        f(mask | added_bits);
    });
}

}  // namespace algo
}  // namespace m1une


#line 9 "algo/enumeration/all.hpp"


#line 1 "algo/offline/all.hpp"



#line 1 "algo/offline/cdq_divide_and_conquer.hpp"



#line 5 "algo/offline/cdq_divide_and_conquer.hpp"

namespace m1une {
namespace algo {

template <class SolveCross>
void cdq_divide_and_conquer(int left, int right, SolveCross solve_cross) {
    assert(left <= right);

    auto dfs = [&](auto& self, int l, int r) -> void {
        if (r - l <= 1) return;
        const int middle = l + (r - l) / 2;
        self(self, l, middle);
        self(self, middle, r);
        solve_cross(l, middle, r);
    };
    dfs(dfs, left, right);
}

template <class SolveCross>
void cdq_divide_and_conquer(int n, SolveCross solve_cross) {
    assert(0 <= n);
    cdq_divide_and_conquer(0, n, solve_cross);
}

}  // namespace algo
}  // namespace m1une


#line 1 "algo/offline/mo.hpp"



#line 6 "algo/offline/mo.hpp"
#include <cmath>
#include <numeric>
#line 9 "algo/offline/mo.hpp"

namespace m1une {
namespace algo {

// Offline Mo's algorithm for half-open array ranges.
struct Mo {
    struct Query {
        int left;
        int right;
        int id;
    };

   private:
    int _n;
    std::vector<Query> _queries;

   public:
    Mo() : _n(0) {}

    explicit Mo(int n) : _n(n) {
        assert(0 <= n);
    }

    int size() const {
        return _n;
    }

    int query_count() const {
        return int(_queries.size());
    }

    bool empty() const {
        return _queries.empty();
    }

    const std::vector<Query>& queries() const {
        return _queries;
    }

    void reserve(int query_capacity) {
        assert(0 <= query_capacity);
        _queries.reserve(query_capacity);
    }

    void clear() {
        _queries.clear();
    }

    // Adds [left, right) and returns its insertion-order ID.
    int add_query(int left, int right) {
        assert(0 <= left && left <= right && right <= _n);
        int id = query_count();
        _queries.push_back(Query{left, right, id});
        return id;
    }

    // Returns query IDs in Mo order. A non-positive block size selects one
    // automatically.
    std::vector<int> order(int block_size = 0) const {
        int query_size = query_count();
        std::vector<int> result(query_size);
        std::iota(result.begin(), result.end(), 0);
        if (query_size == 0) return result;

        if (block_size <= 0) {
            block_size = std::max(1, int(_n / std::sqrt(static_cast<double>(query_size))));
        }

        std::sort(result.begin(), result.end(), [&](int first, int second) {
            const Query& a = _queries[first];
            const Query& b = _queries[second];
            int first_block = a.left / block_size;
            int second_block = b.left / block_size;
            if (first_block != second_block) {
                return first_block < second_block;
            }
            if (first_block & 1) return a.right > b.right;
            return a.right < b.right;
        });
        return result;
    }

    // Maintains [left, right). Each movement callback receives the array index
    // being inserted or erased. `answer(query_id)` stores or reports a result.
    template <class AddLeft, class AddRight, class RemoveLeft, class RemoveRight, class Answer>
    void run(AddLeft add_left, AddRight add_right, RemoveLeft remove_left, RemoveRight remove_right, Answer answer,
             int block_size = 0) const {
        int left = 0;
        int right = 0;
        for (int query_index : order(block_size)) {
            const Query& query = _queries[query_index];
            while (query.left < left) add_left(--left);
            while (right < query.right) add_right(right++);
            while (left < query.left) remove_left(left++);
            while (query.right < right) remove_right(--right);
            answer(query.id);
        }
    }

    // Convenience overload for statistics whose update is independent of
    // which side moves.
    template <class Add, class Remove, class Answer>
    void run(Add add, Remove remove, Answer answer, int block_size = 0) const {
        run(add, add, remove, remove, answer, block_size);
    }
};

}  // namespace algo
}  // namespace m1une


#line 1 "algo/offline/parallel_binary_search.hpp"



#line 6 "algo/offline/parallel_binary_search.hpp"

namespace m1une {
namespace algo {

template <class Apply, class Check, class Reset>
std::vector<int> parallel_binary_search(
    int query_count,
    int event_count,
    Apply apply,
    Check check,
    Reset reset
) {
    assert(0 <= query_count);
    assert(0 <= event_count);

    std::vector<int> low(query_count, -1);
    std::vector<int> high(query_count, event_count + 1);
    std::vector<std::vector<int>> bucket(event_count + 1);

    while (true) {
        bool active = false;
        for (auto& queries : bucket) queries.clear();

        for (int query = 0; query < query_count; ++query) {
            if (high[query] - low[query] <= 1) continue;
            const int middle = low[query] + (high[query] - low[query]) / 2;
            bucket[middle].push_back(query);
            active = true;
        }
        if (!active) break;

        reset();
        int applied = 0;
        for (int middle = 0; middle <= event_count; ++middle) {
            while (applied < middle) {
                apply(applied);
                ++applied;
            }
            for (int query : bucket[middle]) {
                if (check(query)) {
                    high[query] = middle;
                } else {
                    low[query] = middle;
                }
            }
        }
    }

    return high;
}

}  // namespace algo
}  // namespace m1une


#line 7 "algo/offline/all.hpp"


#line 1 "algo/search/all.hpp"



#line 1 "algo/search/bisect.hpp"



#line 5 "algo/search/bisect.hpp"

namespace m1une {
namespace algo {

template <typename F>
long long first_true(long long ng, long long ok, F pred) {
    auto distance = [](long long a, long long b) {
        return a > b ? static_cast<__int128_t>(a) - b : static_cast<__int128_t>(b) - a;
    };
    while (distance(ng, ok) > 1) {
        long long mid = std::midpoint(ng, ok);
        if (pred(mid)) {
            ok = mid;
        } else {
            ng = mid;
        }
    }
    return ok;
}

template <typename F>
long long last_true(long long ok, long long ng, F pred) {
    auto distance = [](long long a, long long b) {
        return a > b ? static_cast<__int128_t>(a) - b : static_cast<__int128_t>(b) - a;
    };
    while (distance(ok, ng) > 1) {
        long long mid = std::midpoint(ok, ng);
        if (pred(mid)) {
            ok = mid;
        } else {
            ng = mid;
        }
    }
    return ok;
}

template <typename F>
double real_first_true(double ng, double ok, F pred, int iterations = 80) {
    for (int i = 0; i < iterations; ++i) {
        double mid = (ng + ok) / 2.0;
        if (pred(mid)) {
            ok = mid;
        } else {
            ng = mid;
        }
    }
    return ok;
}

}  // namespace algo
}  // namespace m1une


#line 1 "algo/search/golden_section_search.hpp"



#line 10 "algo/search/golden_section_search.hpp"

namespace m1une {
namespace algo {

namespace detail {

template <std::integral Int, class F, class Compare>
Int integer_golden_section_search(Int left, Int right, F f, Compare comp) {
    assert(left < right);

    using UInt = std::make_unsigned_t<Int>;
    using Uint128 = unsigned __int128;
    const Uint128 n = static_cast<Uint128>(static_cast<UInt>(right) - static_cast<UInt>(left));

    auto add_offset = [left](Uint128 offset) -> Int {
        if constexpr (std::signed_integral<Int>) {
            if (left < 0) {
                const Uint128 negative_count = static_cast<Uint128>(-(left + 1)) + 1;
                if (offset < negative_count) {
                    return static_cast<Int>(left + static_cast<Int>(offset));
                }
                return static_cast<Int>(offset - negative_count);
            }
        }
        return static_cast<Int>(left + static_cast<Int>(offset));
    };

    using Value = std::decay_t<decltype(f(left))>;
    struct Evaluated {
        Uint128 pos;
        const Value* value;
    };

    Uint128 fib0 = 1;
    Uint128 fib1 = 1;
    Uint128 fib2 = 2;
    int k = 2;
    while (fib2 < n) {
        fib0 = fib1;
        fib1 = fib2;
        fib2 = fib0 + fib1;
        ++k;
    }

    std::vector<std::pair<Uint128, Value>> cache;
    cache.reserve(static_cast<unsigned>(k) + 4);

    auto find_cached = [&](Uint128 pos) -> const Value* {
        for (const auto& [cached_pos, value] : cache) {
            if (cached_pos == pos) return &value;
        }
        return nullptr;
    };

    auto advance_fibonacci = [&]() {
        const Uint128 old0 = fib0;
        const Uint128 old1 = fib1;
        fib0 = old1 - old0;
        fib1 = old0;
        fib2 = old1;
        --k;
    };

    auto eval = [&](Uint128 pos) -> Evaluated {
        if (pos >= n) return Evaluated{pos, nullptr};
        if (const Value* value = find_cached(pos)) return Evaluated{pos, value};
        cache.emplace_back(pos, f(add_offset(pos)));
        return Evaluated{pos, &cache.back().second};
    };

    auto get_value = [&](Uint128 pos) -> const Value& {
        if (const Value* value = find_cached(pos)) return *value;
        cache.emplace_back(pos, f(add_offset(pos)));
        return cache.back().second;
    };

    auto scan = [&](Uint128 scan_left, Uint128 scan_right) -> Int {
        Int best = add_offset(scan_left);
        const Value* best_value = &get_value(scan_left);
        for (Uint128 pos = scan_left + 1; pos <= scan_right; ++pos) {
            Int x = add_offset(pos);
            const Value& value = get_value(pos);
            if (comp(value, *best_value)) {
                best = x;
                best_value = &value;
            }
        }
        return best;
    };

    if (n <= 3) return scan(0, n - 1);

    auto better = [&](const Evaluated& a, const Evaluated& b) -> bool {
        if ((a.value != nullptr) != (b.value != nullptr)) return a.value != nullptr;
        if (a.value == nullptr) return false;
        return comp(*a.value, *b.value);
    };

    Uint128 left_pos = 0;
    Uint128 right_pos = fib2 - 1;
    Uint128 x1 = left_pos + fib0 - 1;
    Uint128 x2 = left_pos + fib1 - 1;
    Evaluated y1 = eval(x1);
    Evaluated y2 = eval(x2);

    while (k > 2) {
        if (better(y2, y1)) {
            left_pos = x1 + 1;
            x1 = x2;
            y1 = y2;
            advance_fibonacci();
            if (k == 2) break;
            x2 = left_pos + fib1 - 1;
            y2 = eval(x2);
        } else {
            right_pos = x2;
            x2 = x1;
            y2 = y1;
            advance_fibonacci();
            if (k == 2) break;
            x1 = left_pos + fib0 - 1;
            y1 = eval(x1);
        }
    }

    const Uint128 last_valid = n - 1;
    if (right_pos > last_valid) right_pos = last_valid;
    assert(left_pos <= right_pos);
    return scan(left_pos, right_pos);
}

}  // namespace detail

template <std::integral Int, class F>
Int golden_section_search_argmin(Int left, Int right, F f) {
    return detail::integer_golden_section_search(left, right, f, [](const auto& a, const auto& b) { return a < b; });
}

template <std::integral Int, class F>
Int golden_section_search_argmax(Int left, Int right, F f) {
    return detail::integer_golden_section_search(left, right, f, [](const auto& a, const auto& b) { return b < a; });
}

template <class F>
double golden_section_search_argmin(double left, double right, F f, int iterations = 100) {
    assert(left <= right);
    assert(0 <= iterations);
    if (left == right || iterations == 0) return std::midpoint(left, right);

    constexpr double inv_phi = 0.6180339887498948482045868343656381177203;
    double x1 = right - (right - left) * inv_phi;
    double x2 = left + (right - left) * inv_phi;
    auto y1 = f(x1);
    auto y2 = f(x2);

    for (int i = 1; i < iterations; ++i) {
        if (y2 < y1) {
            left = x1;
            x1 = x2;
            y1 = std::move(y2);
            x2 = left + (right - left) * inv_phi;
            y2 = f(x2);
        } else {
            right = x2;
            x2 = x1;
            y2 = std::move(y1);
            x1 = right - (right - left) * inv_phi;
            y1 = f(x1);
        }
    }

    if (y2 < y1) {
        left = x1;
    } else {
        right = x2;
    }
    return std::midpoint(left, right);
}

template <class F>
double golden_section_search_argmax(double left, double right, F f, int iterations = 100) {
    assert(left <= right);
    assert(0 <= iterations);
    if (left == right || iterations == 0) return std::midpoint(left, right);

    constexpr double inv_phi = 0.6180339887498948482045868343656381177203;
    double x1 = right - (right - left) * inv_phi;
    double x2 = left + (right - left) * inv_phi;
    auto y1 = f(x1);
    auto y2 = f(x2);

    for (int i = 1; i < iterations; ++i) {
        if (y1 < y2) {
            left = x1;
            x1 = x2;
            y1 = std::move(y2);
            x2 = left + (right - left) * inv_phi;
            y2 = f(x2);
        } else {
            right = x2;
            x2 = x1;
            y2 = std::move(y1);
            x1 = right - (right - left) * inv_phi;
            y1 = f(x1);
        }
    }

    if (y1 < y2) {
        left = x1;
    } else {
        right = x2;
    }
    return std::midpoint(left, right);
}

}  // namespace algo
}  // namespace m1une


#line 1 "algo/search/ternary_search.hpp"



#line 6 "algo/search/ternary_search.hpp"

namespace m1une {
namespace algo {

template <std::integral Int, class F>
Int ternary_search_argmin(Int left, Int right, F f) {
    assert(left < right);
    while (right - left > 3) {
        const Int third = (right - left) / 3;
        const Int middle_left = left + third;
        const Int middle_right = right - third;
        if (f(middle_right) < f(middle_left)) {
            left = middle_left + 1;
        } else {
            right = middle_right;
        }
    }

    Int best = left;
    auto best_value = f(best);
    for (Int x = left + 1; x < right; ++x) {
        auto value = f(x);
        if (value < best_value) {
            best = x;
            best_value = value;
        }
    }
    return best;
}

template <std::integral Int, class F>
Int ternary_search_argmax(Int left, Int right, F f) {
    assert(left < right);
    while (right - left > 3) {
        const Int third = (right - left) / 3;
        const Int middle_left = left + third;
        const Int middle_right = right - third;
        if (f(middle_left) < f(middle_right)) {
            left = middle_left + 1;
        } else {
            right = middle_right;
        }
    }

    Int best = left;
    auto best_value = f(best);
    for (Int x = left + 1; x < right; ++x) {
        auto value = f(x);
        if (best_value < value) {
            best = x;
            best_value = value;
        }
    }
    return best;
}

template <class F>
double real_ternary_search_argmin(double left, double right, F f, int iterations = 100) {
    assert(left <= right);
    assert(0 <= iterations);
    for (int i = 0; i < iterations; ++i) {
        const double middle_left = (left * 2.0 + right) / 3.0;
        const double middle_right = (left + right * 2.0) / 3.0;
        if (f(middle_right) < f(middle_left)) {
            left = middle_left;
        } else {
            right = middle_right;
        }
    }
    return (left + right) / 2.0;
}

template <class F>
double real_ternary_search_argmax(double left, double right, F f, int iterations = 100) {
    assert(left <= right);
    assert(0 <= iterations);
    for (int i = 0; i < iterations; ++i) {
        const double middle_left = (left * 2.0 + right) / 3.0;
        const double middle_right = (left + right * 2.0) / 3.0;
        if (f(middle_left) < f(middle_right)) {
            left = middle_left;
        } else {
            right = middle_right;
        }
    }
    return (left + right) / 2.0;
}

}  // namespace algo
}  // namespace m1une


#line 7 "algo/search/all.hpp"


#line 1 "algo/sequence/all.hpp"



#line 1 "algo/sequence/inversion_count.hpp"



#line 5 "algo/sequence/inversion_count.hpp"

namespace m1une {
namespace algo {

// Returns the number of pairs (i, j) with i < j and a[i] > a[j].
// The vector is taken by value because merge sort rearranges it.
template <typename T>
long long inversion_count(std::vector<T> a) {
    const int n = int(a.size());
    std::vector<T> temp = a;

    auto merge_sort = [&](auto& self, int l, int r) -> long long {
        if (r - l <= 1) return 0;

        const int m = l + (r - l) / 2;
        long long inv = self(self, l, m) + self(self, m, r);

        int i = l;
        int j = m;
        int k = l;
        while (i < m && j < r) {
            if (!(a[j] < a[i])) {
                temp[k++] = a[i++];
            } else {
                temp[k++] = a[j++];
                inv += m - i;
            }
        }

        while (i < m) temp[k++] = a[i++];
        while (j < r) temp[k++] = a[j++];

        for (int p = l; p < r; ++p) {
            a[p] = temp[p];
        }

        return inv;
    };

    return merge_sort(merge_sort, 0, n);
}

}  // namespace algo
}  // namespace m1une


#line 1 "algo/sequence/lis.hpp"



#line 5 "algo/sequence/lis.hpp"
#include <iterator>
#line 7 "algo/sequence/lis.hpp"

namespace m1une {
namespace algo {

// Returns the zero-based indices of a longest increasing subsequence.
// If `strict` is false, equal adjacent values are also allowed.
template <typename T>
std::vector<int> lis(const std::vector<T>& a, bool strict = true) {
    const int n = int(a.size());
    std::vector<T> tails;
    std::vector<int> tail_positions;
    std::vector<int> predecessor(n, -1);
    tails.reserve(n);
    tail_positions.reserve(n);

    for (int i = 0; i < n; ++i) {
        auto it = strict ? std::lower_bound(tails.begin(), tails.end(), a[i])
                         : std::upper_bound(tails.begin(), tails.end(), a[i]);
        const int length = int(std::distance(tails.begin(), it));

        if (it == tails.end()) {
            tails.push_back(a[i]);
            tail_positions.push_back(i);
        } else {
            *it = a[i];
            tail_positions[length] = i;
        }

        if (length > 0) {
            predecessor[i] = tail_positions[length - 1];
        }
    }

    if (tail_positions.empty()) return {};

    std::vector<int> result;
    result.reserve(tail_positions.size());
    int current = tail_positions.back();
    while (current != -1) {
        result.push_back(current);
        current = predecessor[current];
    }
    std::reverse(result.begin(), result.end());
    return result;
}

}  // namespace algo
}  // namespace m1une


#line 1 "algo/sequence/merge_intervals.hpp"



#line 9 "algo/sequence/merge_intervals.hpp"

namespace m1une {
namespace algo {

// Returns the union of half-open intervals as sorted, disjoint intervals.
template <typename T>
std::vector<std::pair<T, T>> merge_intervals(
    std::vector<std::pair<T, T>> intervals
) {
    for (const auto& [left, right] : intervals) {
        if (right < left) assert(false);
    }

    std::sort(
        intervals.begin(),
        intervals.end(),
        [](const auto& lhs, const auto& rhs) {
            if (lhs.first < rhs.first) return true;
            if (rhs.first < lhs.first) return false;
            return lhs.second < rhs.second;
        }
    );

    std::size_t result_size = 0;
    for (std::size_t index = 0; index < intervals.size(); ++index) {
        auto& [left, right] = intervals[index];
        if (!(left < right)) continue;
        if (result_size == 0 || intervals[result_size - 1].second < left) {
            if (result_size != index) {
                intervals[result_size] = std::move(intervals[index]);
            }
            ++result_size;
        } else if (intervals[result_size - 1].second < right) {
            intervals[result_size - 1].second = std::move(right);
        }
    }
    intervals.erase(intervals.begin() + result_size, intervals.end());
    return intervals;
}

}  // namespace algo
}  // namespace m1une


#line 1 "algo/sequence/mex.hpp"



#line 9 "algo/sequence/mex.hpp"

namespace m1une {
namespace algo {

// Returns the smallest nonnegative integer absent from values.
template <class T>
int mex(const std::vector<T>& values) {
    static_assert(
        std::is_integral_v<T> && sizeof(T) <= sizeof(std::uintmax_t),
        "mex requires standard integral values"
    );
    assert(values.size() <= static_cast<std::size_t>(std::numeric_limits<int>::max()));
    const int n = int(values.size());
    std::vector<unsigned char> present(n, 0);
    for (T value : values) {
        if constexpr (std::is_signed_v<T>) {
            if (value < 0) continue;
        }
        if (static_cast<std::uintmax_t>(value) < static_cast<std::uintmax_t>(n)) {
            present[int(value)] = 1;
        }
    }
    int answer = 0;
    while (answer < n && present[answer]) ++answer;
    return answer;
}

}  // namespace algo
}  // namespace m1une


#line 1 "algo/sequence/non_adjacent_selection.hpp"



#include <functional>
#include <queue>
#line 7 "algo/sequence/non_adjacent_selection.hpp"

namespace m1une {
namespace algo {

namespace detail {

template <typename T>
struct NonAdjacentSelectionEntry {
    T value;
    int index;
};

template <typename T, typename Better>
struct NonAdjacentSelectionCompare {
    Better better;

    bool operator()(
        const NonAdjacentSelectionEntry<T>& lhs,
        const NonAdjacentSelectionEntry<T>& rhs
    ) const {
        if (better(lhs.value, rhs.value)) return false;
        if (better(rhs.value, lhs.value)) return true;
        return lhs.index > rhs.index;
    }
};

template <typename T, typename Better>
std::vector<T> non_adjacent_selection_sums(const std::vector<T>& values, Better better) {
    const int n = int(values.size());
    std::vector<T> weight = values;
    std::vector<int> left(n), right(n);
    std::vector<char> alive(n, true);
    for (int i = 0; i < n; ++i) {
        left[i] = i - 1;
        right[i] = (i + 1 == n ? -1 : i + 1);
    }

    using Entry = NonAdjacentSelectionEntry<T>;
    using Compare = NonAdjacentSelectionCompare<T, Better>;
    std::priority_queue<Entry, std::vector<Entry>, Compare> heap(Compare{better});
    for (int i = 0; i < n; ++i) heap.push(Entry{weight[i], i});

    std::vector<T> result;
    result.reserve((n + 1) / 2);
    T sum{};
    while (int(result.size()) < (n + 1) / 2) {
        while (!alive[heap.top().index]) heap.pop();
        const int current = heap.top().index;
        heap.pop();

        sum += weight[current];
        result.push_back(sum);

        const int l = left[current];
        const int r = right[current];
        if (l != -1 && r != -1) {
            weight[current] = weight[l] + weight[r] - weight[current];

            const int ll = left[l];
            const int rr = right[r];
            alive[l] = false;
            alive[r] = false;
            left[current] = ll;
            right[current] = rr;
            if (ll != -1) right[ll] = current;
            if (rr != -1) left[rr] = current;
            heap.push(Entry{weight[current], current});
        } else {
            const int ll = (l == -1 ? -1 : left[l]);
            const int rr = (r == -1 ? -1 : right[r]);
            alive[current] = false;
            if (l != -1) alive[l] = false;
            if (r != -1) alive[r] = false;
            if (ll != -1) right[ll] = rr;
            if (rr != -1) left[rr] = ll;
        }
    }
    return result;
}

}  // namespace detail

// Entry k - 1 is the maximum sum obtained by selecting exactly k values, with
// no two selected indices adjacent.
template <typename T>
std::vector<T> maximum_non_adjacent_selection_sums(const std::vector<T>& values) {
    return detail::non_adjacent_selection_sums(values, std::greater<T>{});
}

// Entry k - 1 is the minimum sum obtained by selecting exactly k values, with
// no two selected indices adjacent.
template <typename T>
std::vector<T> minimum_non_adjacent_selection_sums(const std::vector<T>& values) {
    return detail::non_adjacent_selection_sums(values, std::less<T>{});
}

}  // namespace algo
}  // namespace m1une


#line 1 "algo/sequence/number_of_subsequences.hpp"



#line 6 "algo/sequence/number_of_subsequences.hpp"

namespace m1une {
namespace algo {

// Returns the number of distinct nonempty subsequences.
template <class Mint, class T>
Mint number_of_distinct_subsequences(const std::vector<T>& values) {
    std::vector<T> compressed = values;
    std::sort(compressed.begin(), compressed.end());
    compressed.erase(
        std::unique(compressed.begin(), compressed.end()),
        compressed.end()
    );

    std::vector<Mint> previous_total(compressed.size(), Mint(0));
    Mint total = 1;
    for (const T& value : values) {
        int rank = int(
            std::lower_bound(
                compressed.begin(),
                compressed.end(),
                value
            ) - compressed.begin()
        );
        Mint old_total = total;
        total = total + total - previous_total[rank];
        previous_total[rank] = old_total;
    }
    return total - Mint(1);
}

template <class Mint, class T>
Mint number_of_subsequences(const std::vector<T>& values) {
    return number_of_distinct_subsequences<Mint>(values);
}

}  // namespace algo
}  // namespace m1une


#line 1 "algo/sequence/run_length_encoding.hpp"



#line 7 "algo/sequence/run_length_encoding.hpp"

namespace m1une {
namespace algo {

template <typename Container>
auto run_length_encoding(const Container& values) {
    using T = typename Container::value_type;
    std::vector<std::pair<T, long long>> result;

    auto it = std::begin(values);
    auto last = std::end(values);
    if (it == last) {
        return result;
    }

    T current = *it;
    long long count = 0;
    for (; it != last; ++it) {
        if (*it == current) {
            ++count;
        } else {
            result.emplace_back(current, count);
            current = *it;
            count = 1;
        }
    }
    result.emplace_back(current, count);
    return result;
}

}  // namespace algo
}  // namespace m1une


#line 1 "algo/sequence/subset_sum.hpp"



#line 8 "algo/sequence/subset_sum.hpp"

namespace m1une {
namespace algo {

namespace internal {

template <typename T>
std::vector<T> enumerate_sorted_subset_sums(
    const std::vector<T>& values,
    int left,
    int right
) {
    std::vector<T> sums(1, T{});
    std::vector<T> merged;

    for (int i = left; i < right; ++i) {
        const std::size_t size = sums.size();
        merged.clear();
        merged.reserve(size * 2);

        std::size_t without = 0;
        std::size_t with = 0;
        while (without < size && with < size) {
            const T with_current = sums[with] + values[i];
            if (with_current < sums[without]) {
                merged.push_back(with_current);
                ++with;
            } else {
                merged.push_back(sums[without]);
                ++without;
            }
        }
        while (without < size) {
            merged.push_back(sums[without]);
            ++without;
        }
        while (with < size) {
            merged.push_back(sums[with] + values[i]);
            ++with;
        }
        sums.swap(merged);
    }

    return sums;
}

}  // namespace internal

// Returns the sorted subset sums of values[0, n / 2) and values[n / 2, n).
template <typename T>
std::pair<std::vector<T>, std::vector<T>> enumerate_half_subset_sums(
    const std::vector<T>& values
) {
    const int n = int(values.size());
    const int middle = n / 2;
    return {
        internal::enumerate_sorted_subset_sums(values, 0, middle),
        internal::enumerate_sorted_subset_sums(values, middle, n)
    };
}

// Returns the maximum subset sum not exceeding limit.
template <typename T>
T maximum_subset_sum(const std::vector<T>& values, const T& limit) {
    assert(!(limit < T{}));
    auto [left_sums, right_sums] = enumerate_half_subset_sums(values);

    T answer{};
    std::size_t right_count = right_sums.size();
    for (const T& left : left_sums) {
        while (
            right_count > 0 &&
            limit < left + right_sums[right_count - 1]
        ) {
            --right_count;
        }
        if (right_count == 0) break;

        const T candidate = left + right_sums[right_count - 1];
        if (answer < candidate) answer = candidate;
    }
    return answer;
}

}  // namespace algo
}  // namespace m1une


#line 12 "algo/sequence/all.hpp"


#line 9 "algo/all.hpp"
Back to top page