Algorithms
This is my C++20 library Stormbreaker for competitive programming.
Table of Contents
- Standard Library Extensions
- Data Structures
- Graph Theory
- Mathematical Optimization
- Number Theory
- Numerical Analysis
- String Searching
- Computational Geometry
// Copyright 2018 Tyler Wang. All Rights Reserved.
// Author: Tyler Wang
#include <algorithm>
#include <array>
#include <bit>
#include <cassert>
#include <chrono>
#include <cmath>
#include <complex>
#include <cstdint>
#include <functional>
#include <initializer_list>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <limits>
#include <numeric>
#include <optional>
#include <queue>
#include <random>
#include <set>
#include <stack>
#include <string>
#include <string_view>
#include <tuple>
#include <type_traits>
#include <utility>
#include <vector>
#if __cplusplus < 202002L
namespace std {
template <typename C>
constexpr auto ssize(const C& c)
-> std::common_type_t<std::ptrdiff_t,
std::make_signed_t<decltype(c.size())>> { // C++20
using R = std::common_type_t<std::ptrdiff_t,
std::make_signed_t<decltype(c.size())>>;
return (R)c.size();
}
template <typename T, std::ptrdiff_t N>
constexpr std::ptrdiff_t ssize(const T (&)[N]) noexcept { // C++20
return N;
}
template <typename T>
constexpr std::enable_if_t<std::is_unsigned_v<T>, bool> has_single_bit(
T x) noexcept { // C++20
return x != 0 && (x & (x - 1)) == 0;
}
template <typename T>
constexpr std::enable_if_t<std::is_unsigned_v<T>, int> countl_zero(
T x) noexcept { // C++20
if constexpr (sizeof(T) <= sizeof(unsigned)) {
return __builtin_clz(x) - (std::numeric_limits<unsigned>::digits -
std::numeric_limits<T>::digits);
} else if constexpr (sizeof(T) <= sizeof(unsigned long)) {
return __builtin_clzl(x) - (std::numeric_limits<unsigned long>::digits -
std::numeric_limits<T>::digits);
} else {
static_assert(sizeof(T) <= sizeof(unsigned long long));
return __builtin_clzll(x) -
(std::numeric_limits<unsigned long long>::digits -
std::numeric_limits<T>::digits);
}
}
template <typename T>
constexpr std::enable_if_t<std::is_unsigned_v<T>, int> countr_zero(
T x) noexcept { // C++20
if constexpr (sizeof(T) <= sizeof(unsigned)) {
return __builtin_ctz(x);
} else if constexpr (sizeof(T) <= sizeof(unsigned long)) {
return __builtin_ctzl(x);
} else {
static_assert(sizeof(T) <= sizeof(unsigned long long));
return __builtin_ctzll(x);
}
}
template <typename T>
constexpr std::enable_if_t<std::is_unsigned_v<T>, int> popcount(
T x) noexcept { // C++20
if constexpr (sizeof(T) <= sizeof(unsigned)) {
return __builtin_popcount(x);
} else if constexpr (sizeof(T) <= sizeof(unsigned long long)) {
return __builtin_popcountl(x);
} else {
static_assert(sizeof(T) <= sizeof(unsigned long long));
return __builtin_popcountll(x);
}
}
template <typename T>
constexpr std::enable_if_t<std::is_unsigned_v<T>, T> bit_width(
T x) noexcept { // C++20
return (T)(std::numeric_limits<T>::digits - countl_zero(x));
}
template <typename T>
constexpr std::enable_if_t<std::is_unsigned_v<T>, T> bit_ceil(
T x) noexcept { // C++20
return (T)(x <= 1 ? 1 : (T)1 << bit_width((T)(x - 1)));
}
template <typename T>
constexpr std::enable_if_t<std::is_unsigned_v<T>, T> bit_floor(
T x) noexcept { // C++20
return (T)(x == 0 ? 0 : (T)1 << (bit_width(x) - 1));
}
} // namespace std
#endif
// Y combinator for recursive lambda functions.
// Reference: Derevenets, Yegor. "A Proposal to Add Y Combinator to the Standard
// Library." International Organization for Standardization (2016).
template <typename Fun>
class y_combinator_result {
public:
template <typename T>
explicit y_combinator_result(T&& fun) : fun_(std::forward<T>(fun)) {}
template <typename... Args>
decltype(auto) operator()(Args&&... args) {
return fun_(std::ref(*this), std::forward<Args>(args)...);
}
private:
Fun fun_;
};
template <typename Fun>
decltype(auto) y_combinator(Fun&& fun) {
return y_combinator_result<std::decay_t<Fun>>(std::forward<Fun>(fun));
}
template <typename T>
std::ostream& operator<<(std::ostream& os, const std::optional<T>& opt);
template <typename T1, typename T2>
std::ostream& operator<<(std::ostream& os, const std::pair<T1, T2>& p);
template <typename... Ts>
std::ostream& operator<<(std::ostream& os, const std::tuple<Ts...>& t);
template <typename Container,
std::enable_if_t<!std::is_convertible_v<Container, std::string_view>,
typename Container::const_iterator>* = nullptr>
std::ostream& operator<<(std::ostream& os, const Container& c);
template <typename ContainerAdaptor,
typename ContainerAdaptor::container_type* = nullptr>
std::ostream& operator<<(std::ostream& os, const ContainerAdaptor& a);
template <typename T>
std::ostream& operator<<(std::ostream& os, const std::optional<T>& opt) {
return opt ? os << *opt : os << "nullopt";
}
template <typename T1, typename T2>
std::ostream& operator<<(std::ostream& os, const std::pair<T1, T2>& p) {
return os << '{' << p.first << ", " << p.second << '}';
}
template <typename... Ts>
std::ostream& operator<<(std::ostream& os, const std::tuple<Ts...>& t) {
return std::apply(
[&os](const Ts&... args) -> std::ostream& {
os << '{';
int i = 0;
((os << (i++ == 0 ? "" : ", ") << args), ...);
return os << '}';
},
t);
}
template <typename Container,
std::enable_if_t<!std::is_convertible_v<Container, std::string_view>,
typename Container::const_iterator>*>
std::ostream& operator<<(std::ostream& os, const Container& c) {
os << '{';
for (auto it = c.begin(); it != c.end(); ++it) {
os << (it == c.begin() ? "" : ", ") << *it;
}
return os << '}';
}
template <typename ContainerAdaptor, typename ContainerAdaptor::container_type*>
std::ostream& operator<<(std::ostream& os, const ContainerAdaptor& a) {
class derived : public ContainerAdaptor {
public:
const auto& c() const { return ContainerAdaptor::c; }
};
return os << static_cast<const derived&>(a).c();
}
namespace detail {
template <typename... Ts>
void print(std::ostream& os, int line, std::string_view keys,
const Ts&... vals) {
if constexpr (constexpr int N = sizeof...(Ts); N == 0) {
os << " " << line << " |" << std::endl;
} else {
std::array<std::string_view, N> key_list;
int k = 0, first = 0, paren = 0;
for (int i = 0; i < (int)keys.size(); ++i) {
switch (keys[i]) {
case '(':
case '[':
case '{':
++paren;
break;
case ')':
case ']':
case '}':
--paren;
break;
case ',':
if (paren == 0) {
key_list[k++] = keys.substr(first, i - first);
first = i;
}
break;
}
}
assert(k == N - 1);
key_list[k] = keys.substr(first);
k = 0;
os << " " << line << " | ";
((os << key_list[k++] << " = " << vals), ...);
os << std::endl;
}
}
} // namespace detail
#if defined(TYLER) && defined(__clang__)
#define debug(...) \
detail::print(std::cerr, __LINE__, #__VA_ARGS__, ##__VA_ARGS__)
#elif defined(TYLER) && defined(__GNUC__)
#define debug(...) \
detail::print(std::cerr, __LINE__, #__VA_ARGS__ __VA_OPT__(, ) __VA_ARGS__)
#else
#define debug(...) ((void)0)
#endif
// Multidimensional array.
template <typename T, int Rank>
class mdarray {
static_assert(Rank > 0);
public:
using reference = typename std::vector<T>::reference;
using const_reference = typename std::vector<T>::const_reference;
mdarray(const std::array<int, Rank>& extents, const T& val = T()) {
std::copy(extents.begin(), extents.end(), extents_.begin());
strides_[Rank - 1] = 1;
for (int r = Rank - 2; r >= 0; --r) {
strides_[r] = strides_[r + 1] * extents_[r + 1];
}
data_.resize(strides_[0] * extents_[0], val);
}
static constexpr int rank() { return Rank; }
int size() const { return (int)data_.size(); }
int extent(int r) const { return extents_[r]; }
int stride(int r) const { return strides_[r]; }
template <typename... SizeType>
reference operator()(SizeType... indices) {
return data_[flatten(indices...)];
}
template <typename... SizeType>
const_reference operator()(SizeType... indices) const {
return data_[flatten(indices...)];
}
friend std::ostream& operator<<(std::ostream& os, const mdarray& a) {
for (int i = 0; i < a.size(); ++i) {
for (int r = a.rank() - 1; r >= 0; --r) {
if (i % (a.stride(r) * a.extent(r)) != 0) {
break;
}
os << '{';
}
os << a.data_[i];
for (int r = a.rank() - 1; r >= 0; --r) {
if ((i + 1) % (a.stride(r) * a.extent(r)) != 0) {
os << ", ";
break;
}
os << '}';
}
}
return os;
}
private:
std::array<int, Rank> extents_, strides_;
std::vector<T> data_;
template <typename... SizeType>
int flatten(SizeType... indices) const {
static_assert(sizeof...(SizeType) == Rank);
static_assert((std::is_integral_v<SizeType> && ...));
int pos = 0, r = 0;
((
#ifdef _GLIBCXX_DEBUG
assert(0 <= indices && indices < extents_[r]),
#endif
pos += indices * strides_[r++]),
...);
return pos;
}
};
#if defined(__GNUC__) && !defined(__clang__)
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
// Red-black order statistic tree.
// Supports find_by_order and order_of_key in O(log n).
template <typename Key, typename Compare = std::less<Key>,
typename Allocator = std::allocator<Key>>
using ordered_set =
__gnu_pbds::tree<Key, __gnu_pbds::null_type, Compare,
__gnu_pbds::rb_tree_tag,
__gnu_pbds::tree_order_statistics_node_update, Allocator>;
template <typename Key, typename T, typename Compare = std::less<Key>,
typename Allocator = std::allocator<std::pair<const Key, T>>>
using ordered_map =
__gnu_pbds::tree<Key, T, Compare, __gnu_pbds::rb_tree_tag,
__gnu_pbds::tree_order_statistics_node_update, Allocator>;
#endif
// Splay tree with order statistics node update.
// Amortized O(log n) lookup, insertion, and deletion.
// Reference: Sleator, Daniel Dominic, and Robert Endre Tarjan. "Self-adjusting
// binary search trees." Journal of the ACM (JACM) 32.3 (1985): 652-686.
template <typename Key, typename Compare = std::less<Key>>
class splay_set {
static_assert(std::is_same_v<
std::invoke_result_t<Compare, const Key&, const Key&>, bool>);
public:
using key_type = Key;
using value_type = Key;
using size_type = int;
using difference_type = int;
using key_compare = Compare;
using value_compare = Compare;
using reference = Key&;
using const_reference = const Key&;
using pointer = Key*;
using const_pointer = const Key*;
private:
struct node {
node *parent, *left, *right;
const Key key;
size_type size;
node(node* _parent, const Key& _key, int _size = 1)
: parent(_parent),
left(nullptr),
right(nullptr),
key(_key),
size(_size) {}
void update() {
size = 1 + (left == nullptr ? 0 : left->size) +
(right == nullptr ? 0 : right->size);
}
};
public:
struct iterator {
public:
using iterator_category = std::bidirectional_iterator_tag;
using value_type = Key;
using difference_type = int;
using pointer = const Key*;
using reference = const Key&;
iterator& operator++() {
if (ptr_->right == nullptr) {
while (ptr_->parent != nullptr && ptr_->parent->right == ptr_) {
ptr_ = ptr_->parent;
}
if ((ptr_ = ptr_->parent) != nullptr) {
tree_->splay(ptr_);
}
} else {
ptr_ = ptr_->right;
while (ptr_->left != nullptr) {
ptr_ = ptr_->left;
}
tree_->splay(ptr_);
}
return *this;
}
iterator operator++(int) {
iterator t = *this;
++*this;
return t;
}
iterator& operator--() {
if (ptr_ == nullptr) {
ptr_ = tree_->root_;
while (ptr_->right != nullptr) {
ptr_ = ptr_->right;
}
} else if (ptr_->left == nullptr) {
while (ptr_->parent->left == ptr_) {
ptr_ = ptr_->parent;
}
ptr_ = ptr_->parent;
} else {
ptr_ = ptr_->left;
while (ptr_->right != nullptr) {
ptr_ = ptr_->right;
}
}
tree_->splay(ptr_);
return *this;
}
iterator operator--(int) {
iterator t = *this;
--*this;
return t;
}
bool operator==(const iterator& other) const { return ptr_ == other.ptr_; }
bool operator!=(const iterator& other) const { return ptr_ != other.ptr_; }
reference operator*() const { return ptr_->key; }
private:
friend splay_set;
splay_set* tree_;
node* ptr_;
iterator(splay_set* tree, node* ptr) : tree_(tree), ptr_(ptr) {}
};
struct const_iterator {
public:
using iterator_category = std::bidirectional_iterator_tag;
using value_type = Key;
using difference_type = int;
using pointer = const Key*;
using reference = const Key&;
const_iterator(const iterator& it) : tree_(it.tree_), ptr_(it.ptr_) {}
const_iterator& operator++() {
if (ptr_->right == nullptr) {
while (ptr_->parent != nullptr && ptr_->parent->right == ptr_) {
ptr_ = ptr_->parent;
}
ptr_ = ptr_->parent;
} else {
ptr_ = ptr_->right;
while (ptr_->left != nullptr) {
ptr_ = ptr_->left;
}
}
return *this;
}
const_iterator operator++(int) {
const_iterator t = *this;
++*this;
return t;
}
const_iterator& operator--() {
if (ptr_ == nullptr) {
ptr_ = tree_->root_;
while (ptr_->right != nullptr) {
ptr_ = ptr_->right;
}
} else if (ptr_->left == nullptr) {
while (ptr_->parent->left == ptr_) {
ptr_ = ptr_->parent;
}
ptr_ = ptr_->parent;
} else {
ptr_ = ptr_->left;
while (ptr_->right != nullptr) {
ptr_ = ptr_->right;
}
}
return *this;
}
const_iterator operator--(int) {
const_iterator t = *this;
--*this;
return t;
}
bool operator==(const const_iterator& other) const {
return ptr_ == other.ptr_;
}
bool operator!=(const const_iterator& other) const {
return ptr_ != other.ptr_;
}
reference operator*() const { return ptr_->key; }
private:
friend splay_set;
const splay_set* tree_;
const node* ptr_;
const_iterator(const splay_set* tree, const node* ptr)
: tree_(tree), ptr_(ptr) {}
};
using reverse_iterator = std::reverse_iterator<iterator>;
using const_reverse_iterator = std::reverse_iterator<const_iterator>;
splay_set(const Compare& comp = Compare()) : comp_(comp), root_(nullptr) {}
splay_set(const splay_set& other) : comp_(other.comp_) {
if (other.root_ == nullptr) {
root_ = nullptr;
return;
}
root_ = new node(nullptr, other.root_->key, other.root_->size);
for (node *x = root_, *y = other.root_;;) {
if (x->left == nullptr && y->left != nullptr) {
y = y->left;
x = x->left = new node(x, y->key, y->size);
} else if (x->right == nullptr && y->right != nullptr) {
y = y->right;
x = x->right = new node(x, y->key, y->size);
} else {
if (x->parent == nullptr) {
return;
}
x = x->parent;
y = y->parent;
}
}
}
splay_set(splay_set&& other) : comp_(other.comp_), root_(other.root_) {
other.root_ = nullptr;
}
~splay_set() { clear(); }
splay_set& operator=(const splay_set& other) {
return *this = splay_set(other);
}
splay_set& operator=(splay_set&& other) {
clear();
comp_ = other.comp_;
root_ = other.root_;
other.root_ = nullptr;
return *this;
}
iterator begin() {
if (root_ == nullptr) {
return iterator(this, nullptr);
}
node* x = root_;
while (x->left != nullptr) {
x = x->left;
}
splay(x);
return iterator(this, x);
}
const_iterator begin() const {
if (root_ == nullptr) {
return const_iterator(this, nullptr);
}
node* x = root_;
while (x->left != nullptr) {
x = x->left;
}
return const_iterator(this, x);
}
const_iterator cbegin() const { return begin(); }
iterator end() { return iterator(this, nullptr); }
const_iterator end() const { return const_iterator(this, nullptr); }
const_iterator cend() const { return end(); }
reverse_iterator rbegin() { return reverse_iterator(end()); }
const_reverse_iterator rbegin() const {
return const_reverse_iterator(end());
}
const_reverse_iterator crbegin() const { return rbegin(); }
reverse_iterator rend() { return reverse_iterator(begin()); }
const_reverse_iterator rend() const {
return const_reverse_iterator(begin());
}
const_reverse_iterator crend() const { return rend(); }
bool empty() const { return root_ == nullptr; }
size_type size() const { return root_ == nullptr ? 0 : root_->size; }
void clear() {
if (root_ == nullptr) {
return;
}
for (node* x = root_;;) {
if (x->left != nullptr) {
x = x->left;
} else if (x->right != nullptr) {
x = x->right;
} else {
if (x->parent == nullptr) {
root_ = nullptr;
delete x;
return;
}
node* y = x;
x = x->parent;
if (x->left == y) {
x->left = nullptr;
} else {
x->right = nullptr;
}
delete y;
}
}
}
std::pair<iterator, bool> insert(const Key& key) {
if (root_ == nullptr) {
root_ = new node(nullptr, key);
return std::pair(iterator(this, root_), true);
}
for (node* x = root_;;) {
if (comp_(key, x->key)) {
if (x->left == nullptr) {
x = x->left = new node(x, key);
splay(x);
return std::pair(iterator(this, x), true);
}
x = x->left;
} else if (comp_(x->key, key)) {
if (x->right == nullptr) {
x = x->right = new node(x, key);
splay(x);
return std::pair(iterator(this, x), true);
}
x = x->right;
} else {
splay(x);
return std::pair(iterator(this, x), false);
}
}
}
iterator erase(iterator pos) {
node *x = pos.ptr_, *y, *z;
if (x->right == nullptr) {
if ((y = x->left) != nullptr) {
y->parent = x->parent;
}
z = x;
while (z->parent != nullptr && z->parent->right == z) {
z = z->parent;
}
z = z->parent;
} else {
y = x->right;
y->parent = nullptr;
while (y->left != nullptr) {
y = y->left;
}
splay(y);
if (x->left != nullptr) {
y->left = x->left;
y->left->parent = y;
y->update();
}
y->parent = x->parent;
z = y;
}
if (x->parent == nullptr) {
root_ = y;
} else {
if (x->parent->left == x) {
x->parent->left = y;
} else {
x->parent->right = y;
}
x->parent->update();
splay(x->parent);
}
delete x;
return iterator(this, z);
}
size_type erase(const Key& key) {
if (root_ == nullptr) {
return 0;
}
node* x = root_;
while (true) {
if (comp_(key, x->key)) {
if (x->left == nullptr) {
splay(x);
return 0;
}
x = x->left;
} else if (comp_(x->key, key)) {
if (x->right == nullptr) {
splay(x);
return 0;
}
x = x->right;
} else {
break;
}
}
node* y;
if (x->right == nullptr) {
if ((y = x->left) != nullptr) {
y->parent = x->parent;
}
} else {
y = x->right;
if (x->left != nullptr) {
y->parent = nullptr;
while (y->left != nullptr) {
y = y->left;
}
splay(y);
y->left = x->left;
y->left->parent = y;
y->update();
}
y->parent = x->parent;
}
if (x->parent == nullptr) {
root_ = y;
} else {
if (x->parent->left == x) {
x->parent->left = y;
} else {
x->parent->right = y;
}
x->parent->update();
splay(x->parent);
}
delete x;
return 1;
}
void swap(splay_set& other) {
std::swap(comp_, other.comp_);
std::swap(root_, other.root_);
}
size_type count(const Key& key) { return find(key) != end(); }
iterator find(const Key& key) {
if (root_ == nullptr) {
return iterator(this, nullptr);
}
for (node* x = root_;;) {
if (comp_(key, x->key)) {
if (x->left == nullptr) {
splay(x);
return iterator(this, nullptr);
}
x = x->left;
} else if (comp_(x->key, key)) {
if (x->right == nullptr) {
splay(x);
return iterator(this, nullptr);
}
x = x->right;
} else {
splay(x);
return iterator(this, x);
}
}
}
std::pair<iterator, iterator> equal_range(const Key& key) {
return std::pair(lower_bound(key), upper_bound(key));
}
iterator lower_bound(const Key& key) {
if (root_ == nullptr) {
return iterator(this, nullptr);
}
for (node *x = root_, *y = nullptr;;) {
if (comp_(key, x->key)) {
if (x->left == nullptr) {
splay(x);
return iterator(this, x);
}
y = x;
x = x->left;
} else if (comp_(x->key, key)) {
if (x->right == nullptr) {
splay(x);
return iterator(this, y);
}
x = x->right;
} else {
splay(x);
return iterator(this, x);
}
}
}
iterator upper_bound(const Key& key) {
if (root_ == nullptr) {
return iterator(this, nullptr);
}
for (node *x = root_, *y = nullptr;;) {
if (comp_(key, x->key)) {
if (x->left == nullptr) {
splay(x);
return iterator(this, x);
}
y = x;
x = x->left;
} else {
if (x->right == nullptr) {
splay(x);
return iterator(this, y);
}
x = x->right;
}
}
}
iterator find_by_order(size_type order) {
if (root_ == nullptr || order >= root_->size) {
return iterator(this, nullptr);
}
for (node* x = root_;;) {
if (size_type left_size = x->left == nullptr ? 0 : x->left->size;
order < left_size) {
x = x->left;
} else if (order > left_size) {
order -= left_size + 1;
x = x->right;
} else {
splay(x);
return iterator(this, x);
}
}
}
size_type order_of_key(const Key& key) {
if (root_ == nullptr) {
return 0;
}
int order = 0;
for (node* x = root_;;) {
if (comp_(key, x->key)) {
if (x->left == nullptr) {
splay(x);
return order;
}
x = x->left;
} else if (order += x->left == nullptr ? 0 : x->left->size;
comp_(x->key, key)) {
++order;
if (x->right == nullptr) {
splay(x);
return order;
}
x = x->right;
} else {
splay(x);
return order;
}
}
}
Compare key_comp() const { return comp_; }
Compare value_comp() const { return comp_; }
friend void swap(splay_set& lhs, splay_set& rhs) { lhs.swap(rhs); }
private:
Compare comp_;
node* root_;
void rotate_left(node* y) {
node* x = y->right;
if (y->parent != nullptr) {
if (y->parent->left == y) {
y->parent->left = x;
} else {
y->parent->right = x;
}
}
if ((y->right = x->left) != nullptr) {
y->right->parent = y;
}
x->left = y;
x->parent = y->parent;
y->parent = x;
y->update();
x->update();
}
void rotate_right(node* y) {
node* x = y->left;
if (y->parent != nullptr) {
if (y->parent->left == y) {
y->parent->left = x;
} else {
y->parent->right = x;
}
}
if ((y->left = x->right) != nullptr) {
y->left->parent = y;
}
x->right = y;
x->parent = y->parent;
y->parent = x;
y->update();
x->update();
}
void splay(node* x) {
while (x->parent != nullptr) {
if (x->parent->left == x) {
if (x->parent->parent == nullptr) {
rotate_right(x->parent);
} else if (x->parent->parent->left == x->parent) {
rotate_right(x->parent->parent);
rotate_right(x->parent);
} else {
rotate_right(x->parent);
rotate_left(x->parent);
}
} else {
if (x->parent->parent == nullptr) {
rotate_left(x->parent);
} else if (x->parent->parent->right == x->parent) {
rotate_left(x->parent->parent);
rotate_left(x->parent);
} else {
rotate_left(x->parent);
rotate_right(x->parent);
}
}
}
root_ = x;
}
};
// Disjoint-set data structure aka union-find data structure.
class disjoint_sets {
public:
disjoint_sets(int n) : count_(n), parent_or_size_(n, -1) {}
int count() const { return count_; }
int size(int x) { return -parent_or_size_[find(x)]; }
int find(int x) {
if (parent_or_size_[x] < 0) {
return x;
}
return parent_or_size_[x] = find(parent_or_size_[x]);
}
bool merge(int x, int y) {
if ((x = find(x)) == (y = find(y))) {
return false;
}
if (parent_or_size_[x] > parent_or_size_[y]) {
std::swap(x, y);
}
parent_or_size_[x] += parent_or_size_[y];
parent_or_size_[y] = x;
--count_;
return true;
}
std::vector<std::vector<int>> components() {
int n = (int)parent_or_size_.size();
std::vector<std::vector<int>> comps(count());
std::vector pos(n, comps.end());
auto it = comps.begin();
for (int i = 0; i < n; ++i) {
int p = find(i);
if (pos[p] == comps.end()) {
pos[p] = it++;
pos[p]->reserve(-parent_or_size_[p]);
}
pos[p]->push_back(i);
}
return comps;
}
friend std::ostream& operator<<(std::ostream& os, disjoint_sets ds) {
std::vector<std::vector<int>> comps = ds.components();
os << '{';
for (auto i = comps.begin(); i != comps.end(); ++i) {
os << (i == comps.begin() ? "" : ", ") << '{';
for (auto j = i->begin(); j != i->end(); ++j) {
os << (j == i->begin() ? "" : ", ") << *j;
}
os << '}';
}
return os << '}';
}
private:
int count_;
std::vector<int> parent_or_size_;
};
// Persistent disjoint-set data structure aka union-find data structure.
class persistent_disjoint_sets {
public:
persistent_disjoint_sets(int n) : count_(n), size_(n, 1) {
parent_.reserve(n);
for (int i = 0; i < n; ++i) {
parent_.push_back(i);
}
}
int count() const { return count_; }
int size(int x) const { return size_[find(x)]; }
int find(int x) const { return parent_[x] == x ? x : find(parent_[x]); }
bool merge(int x, int y) {
if ((x = find(x)) == (y = find(y))) {
return false;
}
if (size_[x] < size_[y]) {
std::swap(x, y);
}
edges_.emplace(x, y);
parent_[y] = x;
size_[x] += size_[y];
--count_;
return true;
}
bool rollback() {
if (edges_.empty()) {
return false;
}
auto [x, y] = edges_.top();
edges_.pop();
parent_[y] = y;
size_[x] -= size_[y];
++count_;
return true;
}
std::vector<std::vector<int>> components() const {
int n = (int)parent_.size();
std::vector<std::vector<int>> comps(count());
std::vector pos(n, comps.end());
auto it = comps.begin();
for (int i = 0; i < n; ++i) {
int p = find(i);
if (pos[p] == comps.end()) {
pos[p] = it++;
pos[p]->reserve(size_[p]);
}
pos[p]->push_back(i);
}
return comps;
}
friend std::ostream& operator<<(std::ostream& os,
const persistent_disjoint_sets& pds) {
std::vector<std::vector<int>> comps = pds.components();
os << '{';
for (auto i = comps.begin(); i != comps.end(); ++i) {
os << (i == comps.begin() ? "" : ", ") << '{';
for (auto j = i->begin(); j != i->end(); ++j) {
os << (j == i->begin() ? "" : ", ") << *j;
}
os << '}';
}
return os << '}';
}
private:
int count_;
std::vector<int> parent_, size_;
std::stack<std::pair<int, int>> edges_;
};
// Binary indexed tree aka Fenwick tree.
// O(n) construction, O(log n) update, and O(log n) query.
template <typename T>
class bi_tree {
public:
bi_tree(int n, const T& val = T()) : data_(n, val) {
for (int i = 0; i < size(); ++i) {
if (int j = i | (i + 1); j < size()) {
data_[j] += data_[i];
}
}
}
bi_tree(const std::vector<T>& data) : data_(data) {
for (int i = 0; i < size(); ++i) {
if (int j = i | (i + 1); j < size()) {
data_[j] += data_[i];
}
}
}
int size() const { return (int)data_.size(); }
void update(int pos, const T& delta) {
for (; pos < size(); pos |= pos + 1) {
data_[pos] += delta;
}
}
T query(int pos) const { return query(pos, pos + 1); }
// Sums the range [first, last).
T query(int first, int last) const {
assert(first < last);
T sum = T();
for (; last > first; last &= last - 1) {
sum += data_[last - 1];
}
for (; first > last; first &= first - 1) {
sum -= data_[first - 1];
}
return sum;
}
// Finds the first pos s.t. query(0, pos + 1) >= val.
int lower_bound(T val) const {
int low = -1;
for (int bit = std::bit_floor((unsigned)size()); bit != 0; bit >>= 1) {
if (int mid = low + bit; mid < size() && data_[mid] < val) {
val -= data_[low = mid];
}
}
return low + 1;
}
// Finds the first pos s.t. query(0, pos + 1) > val.
int upper_bound(T val) const {
int low = -1;
for (int bit = std::bit_floor((unsigned)size()); bit != 0; bit >>= 1) {
if (int mid = low + bit; mid < size() && data_[mid] <= val) {
val -= data_[low = mid];
}
}
return low + 1;
}
friend std::ostream& operator<<(std::ostream& os, const bi_tree& tree) {
os << '{';
for (int i = 0; i < tree.size(); ++i) {
os << (i == 0 ? "" : ", ") << tree.query(i);
}
return os << '}';
}
private:
std::vector<T> data_;
};
// Two-dimensional binary indexed tree.
template <typename T>
class bi_tree_2d {
public:
bi_tree_2d(int m, int n) : data_(m, std::vector<T>(n)) {}
int size1() const { return (int)data_.size(); }
int size2() const { return (int)data_[0].size(); }
void update(int x, int y, const T& delta) {
for (; x < size1(); x |= x + 1) {
for (int j = y; j < size2(); j |= j + 1) {
data_[x][j] += delta;
}
}
}
T query(int x, int y) const { return query(x, x + 1, y, y + 1); }
T query(int x1, int x2, int y1, int y2) const {
assert(x1 < x2 && y1 < y2);
T sum = T();
for (; x2 > x1; x2 &= x2 - 1) {
int first = y1, last = y2;
for (; last > first; last &= last - 1) {
sum += data_[x2 - 1][last - 1];
}
for (; first > last; first &= first - 1) {
sum -= data_[x2 - 1][first - 1];
}
}
for (; x1 > x2; x1 &= x1 - 1) {
int first = y1, last = y2;
for (; last > first; last &= last - 1) {
sum -= data_[x1 - 1][last - 1];
}
for (; first > last; first &= first - 1) {
sum += data_[x1 - 1][first - 1];
}
}
return sum;
}
friend std::ostream& operator<<(std::ostream& os, const bi_tree_2d& tree) {
os << '{';
for (int i = 0; i < tree.size1(); ++i) {
os << (i == 0 ? "" : ", ") << '{';
for (int j = 0; j < tree.size2(); ++j) {
os << (j == 0 ? "" : ", ") << tree.query(i, j);
}
os << '}';
}
return os << '}';
}
private:
std::vector<std::vector<T>> data_;
};
// Segment tree.
// O(n) construction, O(log n) update, and O(log n) query.
template <typename T, typename Merge>
class seg_tree {
static_assert(
std::is_same_v<std::invoke_result_t<Merge, const T&, const T&>, T>);
public:
seg_tree(int n, const T& val = T(), const Merge& merge = Merge())
: merge_(merge),
n_(n),
n_ceil_(std::bit_ceil((unsigned)n_)),
data_(n_ceil_ + n_, val) {
for (int first = n_ceil_, last = n_ceil_ + n_ - 1; first != 1;
first >>= 1, last >>= 1) {
data_[last >> 1] = data_[last];
for (int i = first; i < last; i += 2) {
data_[i >> 1] = merge_(data_[i], data_[i + 1]);
}
}
}
seg_tree(const std::vector<T>& data, const Merge& merge = Merge())
: merge_(merge),
n_((int)data.size()),
n_ceil_(std::bit_ceil((unsigned)n_)),
data_(n_ceil_ + n_) {
std::copy(data.begin(), data.end(), data_.begin() + n_ceil_);
for (int first = n_ceil_, last = n_ceil_ + n_ - 1; first != 1;
first >>= 1, last >>= 1) {
data_[last >> 1] = data_[last];
for (int i = first; i < last; i += 2) {
data_[i >> 1] = merge_(data_[i], data_[i + 1]);
}
}
}
int size() const { return n_; }
template <typename U, typename NodeUpdate>
void update(int pos, const U& val, NodeUpdate node_update) {
static_assert(
std::is_void_v<std::invoke_result_t<NodeUpdate, T&, const U&>>);
node_update(data_[pos += n_ceil_], val);
for (int last = n_ceil_ + n_ - 1; last != 1; pos >>= 1, last >>= 1) {
if ((pos | 1) <= last) {
data_[pos >> 1] = merge_(data_[pos & ~1], data_[pos | 1]);
} else {
data_[pos >> 1] = data_[pos];
}
}
}
T query(int pos) const { return data_[n_ceil_ + pos]; }
// Queries the range [first, last).
T query(int first, int last) const {
assert(first < last);
std::optional<T> left, right;
for (first += n_ceil_, last += n_ceil_ - 1; first <= last;
first >>= 1, last >>= 1) {
if ((first & 1) == 1) {
left = left ? merge_(*left, data_[first++]) : data_[first++];
}
if ((last & 1) == 0) {
right = right ? merge_(data_[last--], *right) : data_[last--];
}
}
return left ? (right ? merge_(*left, *right) : *left) : *right;
}
template <typename Contains>
int find(T val, Contains contains) const {
static_assert(
std::is_same_v<std::invoke_result_t<Contains, const T&, T&>, bool>);
if (!contains(data_[1], val)) {
return n_;
}
int pos = 1;
while (pos < n_ceil_) {
if (contains(data_[2 * pos], val)) {
pos = 2 * pos;
} else {
pos = 2 * pos + 1;
}
}
return pos - n_ceil_;
}
template <typename Contains>
int rfind(T val, Contains contains) const {
static_assert(
std::is_same_v<std::invoke_result_t<Contains, const T&, T&>, bool>);
if (!contains(data_[1], val)) {
return n_;
}
int pos = 1;
while (pos < n_ceil_) {
if (2 * pos + 1 < (int)data_.size() &&
contains(data_[2 * pos + 1], val)) {
pos = 2 * pos + 1;
} else {
pos = 2 * pos;
}
}
return pos - n_ceil_;
}
friend std::ostream& operator<<(std::ostream& os, const seg_tree& tree) {
os << '{';
for (int i = 0; i < tree.size(); ++i) {
os << (i == 0 ? "" : ", ") << tree.query(i);
}
return os << '}';
}
private:
Merge merge_;
const int n_, n_ceil_;
std::vector<T> data_;
};
// Segment tree with lazy propagation.
// O(n) construction, O(log n) update, and O(log n) query.
template <typename T, typename U, typename Merge, typename NodeUpdate,
typename LazyUpdate>
class lazy_seg_tree {
static_assert(
std::is_same_v<std::invoke_result_t<Merge, const T&, const T&>, T>);
static_assert(std::is_void_v<std::invoke_result_t<NodeUpdate, T&, const U&>>);
static_assert(std::is_void_v<std::invoke_result_t<LazyUpdate, U&, const U&>>);
public:
lazy_seg_tree(int n, const T& val = T(), const Merge& merge = Merge(),
const NodeUpdate& node_update = NodeUpdate(),
const LazyUpdate& lazy_update = LazyUpdate())
: merge_(merge),
node_update_(node_update),
lazy_update_(lazy_update),
n_(n),
lazy_(std::bit_ceil((unsigned)n_), std::nullopt),
data_(2 * lazy_.size(), val) {
init(1, 0, n_);
}
lazy_seg_tree(const std::vector<T>& data, const Merge& merge = Merge(),
const NodeUpdate& node_update = NodeUpdate(),
const LazyUpdate& lazy_update = LazyUpdate())
: merge_(merge),
node_update_(node_update),
lazy_update_(lazy_update),
n_((int)data.size()),
lazy_(std::bit_ceil((unsigned)n_), std::nullopt),
data_(2 * lazy_.size()) {
init(1, 0, n_, data);
}
int size() const { return n_; }
void update(int pos, const U& val) { update(pos, pos + 1, val); }
// Updates the range [first, last) with val.
void update(int first, int last, const U& val) {
assert(first < last);
update(1, 0, n_, first, last, val);
}
T query(int pos) { return query(pos, pos + 1); }
// Queries the range [first, last).
T query(int first, int last) {
assert(first < last);
return query(1, 0, n_, first, last);
}
template <typename Contains>
int find(T val, Contains contains) {
static_assert(
std::is_same_v<std::invoke_result_t<Contains, const T&, T&>, bool>);
if (!contains(data_[1], val)) {
return n_;
}
int node = 1, low = 0, high = n_ - 1;
while (low < high) {
push(node, low, high);
int mid = (low + high) >> 1;
if (contains(data_[2 * node], val)) {
node = 2 * node;
high = mid;
} else {
node = 2 * node + 1;
low = mid + 1;
}
}
return low;
}
template <typename Contains>
int rfind(T val, Contains contains) {
static_assert(
std::is_same_v<std::invoke_result_t<Contains, const T&, T&>, bool>);
if (!contains(data_[1], val)) {
return n_;
}
int node = 1, low = 0, high = n_ - 1;
while (low < high) {
push(node, low, high);
int mid = (low + high) >> 1;
if (contains(data_[2 * node + 1], val)) {
node = 2 * node + 1;
low = mid + 1;
} else {
node = 2 * node;
high = mid;
}
}
return low;
}
friend std::ostream& operator<<(std::ostream& os, lazy_seg_tree tree) {
os << '{';
for (int i = 0; i < tree.size(); ++i) {
os << (i == 0 ? "" : ", ") << tree.query(i);
}
return os << '}';
}
private:
Merge merge_;
NodeUpdate node_update_;
LazyUpdate lazy_update_;
const int n_;
std::vector<std::optional<U>> lazy_;
std::vector<T> data_;
void init(int node, int low, int high) {
if (high - low == 1) {
return;
}
int mid = (low + high) >> 1;
init(2 * node, low, mid);
init(2 * node + 1, mid, high);
data_[node] = merge_(data_[2 * node], data_[2 * node + 1]);
}
void init(int node, int low, int high, const std::vector<T>& data) {
if (high - low == 1) {
data_[node] = data[low];
return;
}
int mid = (low + high) >> 1;
init(2 * node, low, mid, data);
init(2 * node + 1, mid, high, data);
data_[node] = merge_(data_[2 * node], data_[2 * node + 1]);
}
void update_node(int node, int low, int high, const U& val) {
node_update_(data_[node], val);
if (high - low == 1) {
return;
}
if (!lazy_[node]) {
lazy_[node] = val;
return;
}
lazy_update_(*lazy_[node], val);
}
void push(int node, int low, int high) {
if (!lazy_[node]) {
return;
}
int mid = (low + high) >> 1;
update_node(2 * node, low, mid, *lazy_[node]);
update_node(2 * node + 1, mid, high, *lazy_[node]);
lazy_[node].reset();
}
void update(int node, int low, int high, int first, int last, const U& val) {
if (first <= low && high <= last) {
update_node(node, low, high, val);
return;
}
push(node, low, high);
int mid = (low + high) >> 1;
if (first < mid) {
update(2 * node, low, mid, first, last, val);
}
if (mid < last) {
update(2 * node + 1, mid, high, first, last, val);
}
data_[node] = merge_(data_[2 * node], data_[2 * node + 1]);
}
T query(int node, int low, int high, int first, int last) {
if (first <= low && high <= last) {
return data_[node];
}
push(node, low, high);
int mid = (low + high) >> 1;
if (last <= mid) {
return query(2 * node, low, mid, first, last);
}
if (mid <= first) {
return query(2 * node + 1, mid, high, first, last);
}
return merge_(query(2 * node, low, mid, first, last),
query(2 * node + 1, mid, high, first, last));
}
};
template <typename T, typename U, typename Merge, typename NodeUpdate,
typename LazyUpdate>
lazy_seg_tree<T, U, Merge, NodeUpdate, LazyUpdate> make_lazy_seg_tree(
int n, const T& val, const Merge& merge, const NodeUpdate& node_update,
const LazyUpdate& lazy_update) {
return lazy_seg_tree<T, U, Merge, NodeUpdate, LazyUpdate>(
n, val, merge, node_update, lazy_update);
}
template <typename T, typename U, typename Merge, typename NodeUpdate,
typename LazyUpdate>
lazy_seg_tree<T, U, Merge, NodeUpdate, LazyUpdate> make_lazy_seg_tree(
const std::vector<T>& data, const Merge& merge,
const NodeUpdate& node_update, const LazyUpdate& lazy_update) {
return lazy_seg_tree<T, U, Merge, NodeUpdate, LazyUpdate>(
data, merge, node_update, lazy_update);
}
// Heavy-light decomposition.
// O(n) construction, O(log^2 n) update, and O(log^2 n) query.
template <typename T, typename U, typename Merge, typename NodeUpdate,
typename LazyUpdate>
class heavy_light {
public:
heavy_light(const std::vector<std::vector<int>>& graph, int root,
const T& val = T(), const Merge& merge = Merge(),
const NodeUpdate& node_update = NodeUpdate(),
const LazyUpdate& lazy_update = LazyUpdate())
: merge_(merge),
tree_((int)graph.size(), val, merge, node_update, lazy_update),
parent_(graph.size(), -1),
depth_(graph.size(), 0),
size_(graph.size(), 1),
head_(graph.size(), root),
pos_(graph.size()) {
std::vector<int> heavy(graph.size(), -1);
dfs1(graph, root, heavy);
int index = 0;
dfs2(graph, root, heavy, index);
}
heavy_light(const std::vector<std::vector<int>>& graph, int root,
const std::vector<T>& data, const Merge& merge = Merge(),
const NodeUpdate& node_update = NodeUpdate(),
const LazyUpdate& lazy_update = LazyUpdate())
: merge_(merge),
tree_(data, merge, node_update, lazy_update),
parent_(graph.size(), -1),
depth_(graph.size(), 0),
size_(graph.size(), 1),
head_(graph.size(), root),
pos_(graph.size()) {
assert(graph.size() == data.size());
std::vector<int> heavy(graph.size(), -1);
dfs1(graph, root, heavy);
int index = 0;
dfs2(graph, root, heavy, index);
}
void update_node(int node, const U& val) { tree_.update(pos_[node], val); }
void update_subtree(int root, const U& val) {
tree_.update(pos_[root], pos_[root] + size_[root], val);
}
void update_path(int u, int v, const U& val) {
while (head_[u] != head_[v]) {
if (depth_[head_[u]] > depth_[head_[v]]) {
std::swap(u, v);
}
tree_.update(pos_[head_[v]], pos_[v] + 1, val);
v = parent_[head_[v]];
}
if (depth_[u] > depth_[v]) {
std::swap(u, v);
}
tree_.update(pos_[u], pos_[v] + 1, val);
}
T query_node(int node) { return tree_.query(pos_[node]); }
T query_subtree(int root) {
return tree_.query(pos_[root], pos_[root] + size_[root]);
}
T query_path(int u, int v) {
std::optional<T> val;
while (head_[u] != head_[v]) {
if (depth_[head_[u]] > depth_[head_[v]]) {
std::swap(u, v);
}
T t = tree_.query(pos_[head_[v]], pos_[v] + 1);
val = val ? merge_(*val, t) : t;
v = parent_[head_[v]];
}
if (depth_[u] > depth_[v]) {
std::swap(u, v);
}
T t = tree_.query(pos_[u], pos_[v] + 1);
return val ? merge_(*val, t) : t;
}
private:
Merge merge_;
lazy_seg_tree<T, U, Merge, NodeUpdate, LazyUpdate> tree_;
std::vector<int> parent_, depth_, size_, head_, pos_;
void dfs1(const std::vector<std::vector<int>>& graph, int u,
std::vector<int>& heavy) {
for (int v : graph[u]) {
if (v == parent_[u]) {
continue;
}
parent_[v] = u;
depth_[v] = depth_[u] + 1;
dfs1(graph, v, heavy);
size_[u] += size_[v];
if (heavy[u] == -1 || size_[v] > size_[heavy[u]]) {
heavy[u] = v;
}
}
}
void dfs2(const std::vector<std::vector<int>>& graph, int u,
const std::vector<int>& heavy, int& index) {
pos_[u] = index++;
if (heavy[u] != -1) {
head_[heavy[u]] = head_[u];
dfs2(graph, heavy[u], heavy, index);
}
for (int v : graph[u]) {
if (v == parent_[u] || v == heavy[u]) {
continue;
}
head_[v] = v;
dfs2(graph, v, heavy, index);
}
}
};
template <typename T, typename U, typename Merge, typename NodeUpdate,
typename LazyUpdate>
heavy_light<T, U, Merge, NodeUpdate, LazyUpdate> make_heavy_light(
const std::vector<std::vector<int>>& graph, int root, const T& val,
const Merge& merge, const NodeUpdate& node_update,
const LazyUpdate& lazy_update) {
return heavy_light<T, U, Merge, NodeUpdate, LazyUpdate>(
graph, root, val, merge, node_update, lazy_update);
}
template <typename T, typename U, typename Merge, typename NodeUpdate,
typename LazyUpdate>
heavy_light<T, U, Merge, NodeUpdate, LazyUpdate> make_heavy_light(
const std::vector<std::vector<int>>& graph, int root,
const std::vector<T>& data, const Merge& merge,
const NodeUpdate& node_update, const LazyUpdate& lazy_update) {
return heavy_light<T, U, Merge, NodeUpdate, LazyUpdate>(
graph, root, data, merge, node_update, lazy_update);
}
// Sparse table.
// O(n log n) construction, O(1) relaxed query, and O(log n) strict query.
template <typename T, typename Merge>
class sparse_table {
static_assert(
std::is_same_v<std::invoke_result_t<Merge, const T&, const T&>, T>);
public:
sparse_table(const std::vector<T>& data, const Merge& merge = Merge())
: merge_(merge), data_(std::bit_width(data.size())) {
int n = (int)data.size();
data_[0] = data;
for (int lg = 1; lg < (int)data_.size(); ++lg) {
data_[lg].reserve(n - (1 << lg) + 1);
for (int i = 0; i + (1 << lg) <= n; ++i) {
data_[lg].push_back(
merge_(data_[lg - 1][i], data_[lg - 1][i + (1 << (lg - 1))]));
}
}
}
int size() const { return (int)data_[0].size(); }
T query(int pos) const { return data_[0][pos]; }
// Queries the range [first, last). Some elements are merged twice.
T query(int first, int last) const {
assert(first < last);
int lg = std::bit_width((unsigned)(last - first)) - 1;
return merge_(data_[lg][first], data_[lg][last - (1 << lg)]);
}
// Queries the range [first, last). Each element is merged once.
T strict_query(int first, int last) const {
assert(first < last);
int lg = std::bit_width((unsigned)(last - first)) - 1;
T val = data_[lg][first];
while ((first += 1 << lg) < last) {
lg = std::bit_width((unsigned)(last - first)) - 1;
val = merge_(val, data_[lg][first]);
}
return val;
}
friend std::ostream& operator<<(std::ostream& os, const sparse_table& table) {
os << '{';
for (int i = 0; i < table.size(); ++i) {
os << (i == 0 ? "" : ", ") << table.query(i);
}
return os << '}';
}
private:
Merge merge_;
std::vector<std::vector<T>> data_;
};
// Lowest common ancestor with sparse table.
// O(n log n) preprocessing + O(1) query.
class tree_lca {
public:
tree_lca(const std::vector<std::vector<int>>& graph, int root)
: depth_(graph.size()),
first_(graph.size()),
table_(euler_tour(graph, root),
[this](int u, int v) { return depth_[u] < depth_[v] ? u : v; }) {
}
const std::vector<int>& depth() const { return depth_; }
int solve(int u, int v) const {
if ((u = first_[u]) > (v = first_[v])) {
std::swap(u, v);
}
return table_.query(u, v + 1);
}
int dist(int u, int v) const {
return depth_[u] + depth_[v] - 2 * depth_[solve(u, v)];
}
private:
std::vector<int> depth_, first_;
sparse_table<int, std::function<int(int, int)>> table_;
std::vector<int> euler_tour(const std::vector<std::vector<int>>& graph,
int root) {
std::vector<int> tour;
tour.reserve(2 * graph.size() - 1);
dfs(graph, root, -1, tour);
return tour;
}
void dfs(const std::vector<std::vector<int>>& graph, int u, int parent,
std::vector<int>& tour) {
first_[u] = (int)tour.size();
tour.push_back(u);
for (int v : graph[u]) {
if (v == parent) {
continue;
}
depth_[v] = depth_[u] + 1;
dfs(graph, v, u, tour);
tour.push_back(u);
}
}
};
// Kahn's topological sort algorithm. O(V + E).
std::vector<int> topo_sort(const std::vector<std::vector<int>>& graph) {
int n = (int)graph.size();
std::vector<int> degree(n, 0);
for (const std::vector<int>& edges : graph) {
for (int v : edges) {
++degree[v];
}
}
std::vector<int> sorted;
sorted.reserve(n);
for (int v = 0; v < n; ++v) {
if (degree[v] == 0) {
sorted.push_back(v);
}
}
for (auto u = sorted.begin(); u != sorted.end(); ++u) {
for (int v : graph[*u]) {
if (--degree[v] == 0) {
sorted.push_back(v);
}
}
}
if ((int)sorted.size() < n) { // cycle
return std::vector<int>();
}
return sorted;
}
// Tarjan's biconnected components algorithm. O(V + E).
class biconnected_components {
public:
biconnected_components(int n) : m_(0), graph_(n) {}
int add_edge(int u, int v) {
graph_[u].emplace_back(v, m_);
graph_[v].emplace_back(u, m_);
return m_++;
}
int solve() {
comp_.assign(m_, -1);
bridges_.clear();
int n = (int)graph_.size();
int index = 0;
std::vector<int> pos(n, -1);
std::stack<int> estack;
int count = 0;
for (int u = 0; u < n; ++u) {
if (pos[u] == -1) {
dfs(u, -1, index, pos, estack, count);
}
}
return count;
}
const std::vector<int>& component() const { return comp_; }
const std::vector<int>& bridges() const { return bridges_; }
private:
struct edge {
const int to, id;
constexpr edge(int _to, int _id) : to(_to), id(_id) {}
};
int m_;
std::vector<std::vector<edge>> graph_;
std::vector<int> comp_;
std::vector<int> bridges_;
int dfs(int u, int parent, int& index, std::vector<int>& pos,
std::stack<int>& estack, int& count) {
int u_low = pos[u] = index++;
for (const edge& e : graph_[u]) {
if (e.to == parent) {
continue;
}
if (pos[e.to] != -1) {
u_low = std::min(u_low, pos[e.to]);
if (pos[e.to] < pos[u]) {
estack.push(e.id);
}
continue;
}
int size = (int)estack.size();
int v_low = dfs(e.to, u, index, pos, estack, count);
u_low = std::min(u_low, v_low);
if (v_low > pos[u]) {
bridges_.push_back(e.id);
continue;
}
estack.push(e.id);
if (v_low == pos[u]) {
do {
comp_[estack.top()] = count;
estack.pop();
} while ((int)estack.size() > size);
++count;
}
}
return u_low;
}
};
namespace detail {
int art_dfs(const std::vector<std::vector<int>>& graph, int u, int parent,
int& index, std::vector<int>& pos, std::vector<int>& points) {
int u_low = pos[u] = index++, children = 0;
bool articulate = false;
for (int v : graph[u]) {
if (v == parent) {
continue;
}
if (pos[v] != -1) {
u_low = std::min(u_low, pos[v]);
continue;
}
int v_low = art_dfs(graph, v, u, index, pos, points);
u_low = std::min(u_low, v_low);
if (v_low >= pos[u]) {
articulate = true;
}
++children;
}
if (articulate || (parent == -1 && children > 1)) {
points.push_back(u);
}
return u_low;
}
} // namespace detail
// Tarjan's articulation points algorithm. O(V + E).
std::vector<int> articulation_points(
const std::vector<std::vector<int>>& graph) {
int n = (int)graph.size();
int index = 0;
std::vector<int> pos(n, -1), points;
for (int u = 0; u < n; ++u) {
if (pos[u] == -1) {
detail::art_dfs(graph, u, -1, index, pos, points);
}
}
return points;
}
namespace detail {
int scc_dfs(const std::vector<std::vector<int>>& graph, int u,
std::vector<int>& pos, std::stack<int>& vstack, int& count,
std::vector<int>& comp) {
int low = pos[u] = (int)vstack.size();
vstack.push(u);
for (int v : graph[u]) {
if (pos[v] == -1) {
low = std::min(low, scc_dfs(graph, v, pos, vstack, count, comp));
} else if (comp[v] == -1) {
low = std::min(low, pos[v]);
}
}
if (low == pos[u]) {
int v;
do {
v = vstack.top();
vstack.pop();
comp[v] = count;
} while (v != u);
++count;
}
return low;
}
} // namespace detail
// Tarjan's strongly connected components algorithm. O(V + E).
std::pair<int, std::vector<int>> strong_components(
const std::vector<std::vector<int>>& graph) {
int n = (int)graph.size();
std::vector<int> pos(n, -1);
std::stack<int> vstack;
int count = 0;
std::vector<int> comp(n, -1);
for (int u = 0; u < n; ++u) {
if (pos[u] == -1) {
detail::scc_dfs(graph, u, pos, vstack, count, comp);
}
}
return std::pair(count, comp);
}
// 2-satisfiability. O(V + E).
class two_sat {
public:
two_sat(int n) : graph_(2 * n), assignment_(n) {}
static constexpr int negate(int a) { return ~a; }
void add_true(int a) {
a = std::max(2 * a, 2 * ~a + 1);
graph_[a ^ 1].push_back(a);
}
void add_or(int a, int b) {
a = std::max(2 * a, 2 * ~a + 1);
b = std::max(2 * b, 2 * ~b + 1);
graph_[a ^ 1].push_back(b);
graph_[b ^ 1].push_back(a);
}
bool solve() {
std::vector<int> comp = strong_components(graph_).second;
for (int i = 0; i < (int)graph_.size(); i += 2) {
if (comp[i] == comp[i + 1]) {
return false;
}
assignment_[i >> 1] = comp[i] > comp[i + 1];
}
return true;
}
const std::vector<bool>& assignment() const { return assignment_; }
private:
std::vector<std::vector<int>> graph_;
std::vector<bool> assignment_;
};
// Eulerian path. O(V + E).
class euler_path {
public:
euler_path(int n) : m_(0), graph_(n) {}
int add_undirected(int u, int v) {
graph_[u].emplace_back(v, m_);
graph_[v].emplace_back(u, m_);
return m_++;
}
int add_directed(int u, int v) {
graph_[u].emplace_back(v, m_);
return m_++;
}
std::vector<int> trail(int src) const { return solve<false>(src); }
std::vector<int> cycle(int src) const { return solve<true>(src); }
private:
struct edge {
const int to, id;
constexpr edge(int _to, int _id) : to(_to), id(_id) {}
};
int m_;
std::vector<std::vector<edge>> graph_;
template <bool Cycle>
std::vector<int> solve(int src) const {
int n = (int)graph_.size();
std::vector<bool> used(m_, false);
std::vector<int> degree(n, 0), path;
if constexpr (!Cycle) {
++degree[src];
}
path.reserve(m_ + 1);
std::vector<std::vector<edge>::const_iterator> curr;
curr.reserve(n);
for (const std::vector<edge>& edges : graph_) {
curr.push_back(edges.begin());
}
std::stack<edge> estack;
estack.emplace(src, -1);
while (!estack.empty()) {
const edge& u = estack.top();
if (curr[u.to] == graph_[u.to].end()) {
path.push_back(u.id);
estack.pop();
continue;
}
const edge& v = *curr[u.to]++;
if (!used[v.id]) {
used[v.id] = true;
estack.push(v);
--degree[u.to];
++degree[v.to];
}
}
path.pop_back();
if ((int)path.size() < m_ || std::any_of(degree.begin(), degree.end(),
[](int d) { return d < 0; })) {
return std::vector<int>();
}
std::reverse(path.begin(), path.end());
return path;
}
};
// Bellman-Ford algorithm for single-source shortest paths. O(V E).
template <typename W>
std::pair<std::vector<W>, std::vector<int>> bellman_shortest_paths(
const std::vector<std::vector<std::pair<int, W>>>& graph, int src,
W inf = std::numeric_limits<W>::max()) {
int n = (int)graph.size();
std::vector<W> dist(n, inf);
dist[src] = 0;
std::vector<int> pred(n, -1), count(n, 0);
pred[src] = src;
count[src] = n - 1;
std::queue<int> q;
q.push(src);
std::vector<bool> active(n, false);
active[src] = true;
while (!q.empty()) {
int u = q.front();
q.pop();
active[u] = false;
for (const auto& [v, w] : graph[u]) {
if (pred[v] == -1 || dist[v] > dist[u] + w) {
dist[v] = dist[u] + w;
pred[v] = u;
if (!active[v]) {
if (++count[v] == n) { // negative-weight cycle
return {{}, {}};
}
q.push(v);
active[v] = true;
}
}
}
}
return std::pair(dist, pred);
}
// Dijkstra's algorithm for single-source shortest paths. O(E log V).
// Edge weights must be non-negative.
template <typename W>
std::pair<std::vector<W>, std::vector<int>> dijkstra_shortest_paths(
const std::vector<std::vector<std::pair<int, W>>>& graph, int src,
W inf = std::numeric_limits<W>::max()) {
int n = (int)graph.size();
std::vector<W> dist(n, inf);
dist[src] = 0;
std::vector<int> pred(n, -1);
pred[src] = src;
constexpr auto comp = [](const std::pair<int, W>& lhs,
const std::pair<int, W>& rhs) {
return lhs.second > rhs.second;
};
std::priority_queue<std::pair<int, W>, std::vector<std::pair<int, W>>,
decltype(comp)>
pq(comp);
pq.emplace(src, 0);
while (!pq.empty()) {
auto [u, d_u] = pq.top();
pq.pop();
if (d_u > dist[u]) {
continue;
}
for (const auto& [v, w] : graph[u]) {
assert(w >= 0);
if (pred[v] == -1 || dist[v] > dist[u] + w) {
pq.emplace(v, dist[v] = dist[u] + w);
pred[v] = u;
}
}
}
return std::pair(dist, pred);
}
// Johnson's algorithm for all-pairs shortest paths. O(V E log V).
template <typename W>
std::pair<std::vector<std::vector<W>>, std::vector<std::vector<int>>>
johnson_shortest_paths(std::vector<std::vector<std::pair<int, W>>> graph,
W inf = std::numeric_limits<W>::max()) {
int n = (int)graph.size();
std::vector<std::pair<int, W>>& graph_n = graph.emplace_back();
graph_n.reserve(n);
for (int i = 0; i < n; ++i) {
graph_n.emplace_back(i, 0);
}
auto [dist_n, pred_n] = bellman_shortest_paths(graph, n, inf);
if (dist_n.empty()) { // negative-weight cycle
return {{}, {}};
}
for (int u = 0; u < n; ++u) {
for (auto& [v, w] : graph[u]) {
w += dist_n[u] - dist_n[v];
}
}
graph.pop_back();
std::vector<std::vector<W>> dist(n);
std::vector<std::vector<int>> pred(n);
for (int u = 0; u < n; ++u) {
std::tie(dist[u], pred[u]) = dijkstra_shortest_paths(graph, u, inf);
for (int v = 0; v < n; ++v) {
if (pred[u][v] != -1) {
dist[u][v] += dist_n[v] - dist_n[u];
}
}
}
return std::pair(dist, pred);
}
// Floyd-Warshall algorithm for all-pairs shortest paths. O(V ^ 3).
template <typename W>
std::pair<std::vector<std::vector<W>>, std::vector<std::vector<int>>>
floyd_shortest_paths(const std::vector<std::vector<std::pair<int, W>>>& graph,
W inf = std::numeric_limits<W>::max()) {
int n = (int)graph.size();
std::vector<std::vector<W>> dist(n, std::vector<W>(n, inf));
std::vector<std::vector<int>> pred(n, std::vector<int>(n, -1));
for (int u = 0; u < n; ++u) {
for (const auto& [v, w] : graph[u]) {
dist[u][v] = w;
pred[u][v] = u;
}
}
for (int v = 0; v < n; ++v) {
dist[v][v] = 0;
pred[v][v] = v;
}
for (int k = 0; k < n; ++k) {
for (int i = 0; i < n; ++i) {
if (pred[i][k] != -1) {
for (int j = 0; j < n; ++j) {
if (pred[k][j] != -1 &&
(pred[i][j] == -1 || dist[i][j] > dist[i][k] + dist[k][j])) {
dist[i][j] = dist[i][k] + dist[k][j];
pred[i][j] = pred[k][j];
}
}
}
}
}
for (int v = 0; v < n; ++v) {
if (dist[v][v] < 0) { // negative-weight cycle
return {{}, {}};
}
}
return std::pair(dist, pred);
}
template <typename W>
struct edge {
int from, to;
W cost;
constexpr edge() : from(-1), to(-1), cost() {}
constexpr edge(int _from, int _to, W _cost)
: from(_from), to(_to), cost(_cost) {}
friend constexpr bool operator<(const edge& lhs, const edge& rhs) {
return lhs.cost < rhs.cost;
}
friend std::ostream& operator<<(std::ostream& os, const edge& e) {
return os << '{' << e.from << ", " << e.to << ", " << e.cost << '}';
}
};
// Kruskal's minimum spanning tree algorithm. O(E log V).
template <typename W>
std::pair<W, std::vector<int>> min_span_tree(
int n, const std::vector<edge<W>>& edges) {
int m = (int)edges.size();
std::vector<int> ids;
ids.reserve(m);
for (int i = 0; i < m; ++i) {
ids.push_back(i);
}
std::sort(ids.begin(), ids.end(),
[&edges](int lhs, int rhs) { return edges[lhs] < edges[rhs]; });
disjoint_sets forest(n);
W weight = 0;
std::vector<int> chosen;
chosen.reserve(n - 1);
for (int i : ids) {
if (forest.merge(edges[i].from, edges[i].to)) {
weight += edges[i].cost;
chosen.push_back(i);
}
}
return std::pair(weight, chosen);
}
// Hopcroft-Karp algorithm for maximum bipartite matching. O(E sqrt(V)).
class bipartite_match {
public:
bipartite_match(int m, int n)
: graph_(m), succ_(m, -1), pred_(n, -1), dist_(m), matches_(0) {}
void add_edge(int from, int to) { graph_[from].push_back(to); }
int solve() {
while (bfs()) {
for (int u = 0; u < (int)succ_.size(); ++u) {
if (succ_[u] == -1 && dfs(u)) {
++matches_;
}
}
}
return matches_;
}
const std::vector<int>& succ() const { return succ_; }
const std::vector<int>& pred() const { return pred_; }
std::vector<int> vertex_cover() {
std::vector<int> cover;
int m = (int)succ_.size();
for (int u = 0; u < m; ++u) {
if (dist_[u] == -1) {
cover.push_back(u);
}
}
int n = (int)pred_.size();
for (int v = 0; v < n; ++v) {
if (pred_[v] != -1 && dist_[pred_[v]] != -1) {
cover.push_back(m + v);
}
}
assert((int)cover.size() == matches_);
return cover;
}
private:
std::vector<std::vector<int>> graph_;
std::vector<int> succ_, pred_, dist_;
int matches_;
bool bfs() {
std::queue<int> q;
for (int u = 0; u < (int)succ_.size(); ++u) {
if (succ_[u] == -1) {
dist_[u] = 0;
q.push(u);
} else {
dist_[u] = -1;
}
}
while (!q.empty()) {
int u = q.front();
q.pop();
for (int v : graph_[u]) {
if (pred_[v] == -1) {
return true;
}
if (dist_[pred_[v]] == -1) {
dist_[pred_[v]] = dist_[u] + 1;
q.push(pred_[v]);
}
}
}
return false;
}
bool dfs(int u) {
for (int v : graph_[u]) {
if (pred_[v] == -1 ||
(dist_[pred_[v]] == dist_[u] + 1 && dfs(pred_[v]))) {
succ_[u] = v;
pred_[v] = u;
return true;
}
}
dist_[u] = -1;
return false;
}
};
// Highest-label push-relabel maximum flow algorithm. O(V ^ 2 * sqrt(E)).
template <typename Flow>
class max_flow {
public:
max_flow(int n) : graph_(n) {}
void add_edge(int from, int to, Flow cap, Flow rev_cap = 0) {
graph_[from].emplace_back(to, (int)graph_[to].size(), cap);
graph_[to].emplace_back(from, (int)graph_[from].size() - 1, rev_cap);
}
Flow solve(int source, int sink) {
for (std::vector<edge>& edges : graph_) {
for (edge& e : edges) {
e.flow = 0;
}
}
int n = (int)graph_.size();
std::vector<int> height(n, 0), count(n, 0);
height[source] = n;
count[0] = n - 1;
std::vector<Flow> excess(n, 0);
std::vector<typename std::vector<edge>::iterator> curr;
curr.reserve(n);
for (std::vector<edge>& edges : graph_) {
curr.push_back(edges.begin());
}
auto comp = [&height](int lhs, int rhs) {
return height[lhs] < height[rhs];
};
std::priority_queue<int, std::vector<int>, decltype(comp)> pq(comp);
for (edge& e : graph_[source]) {
if (e.cap > 0) {
if (e.to != sink && excess[e.to] == 0) {
pq.push(e.to);
}
e.flow = e.cap;
graph_[e.to][e.rev].flow = -e.cap;
excess[source] -= e.cap;
excess[e.to] += e.cap;
}
}
while (!pq.empty()) {
int u = pq.top();
pq.pop();
while (excess[u] > 0) {
if (curr[u] == graph_[u].end()) { // relabel
if (height[u] < n && --count[height[u]] == 0) { // gap heuristic
for (int& h : height) {
if (height[u] < h && h < n) {
--count[h];
h = n + 1;
}
}
}
height[u] = 2 * n;
for (auto e = graph_[u].begin(); e != graph_[u].end(); ++e) {
if (e->flow < e->cap && height[u] > height[e->to] + 1) {
height[u] = height[e->to] + 1;
curr[u] = e;
}
}
if (height[u] < n) {
++count[height[u]];
}
} else if (edge& e = *curr[u];
e.flow < e.cap && height[u] == height[e.to] + 1) { // push
if (e.to != sink && excess[e.to] == 0) {
pq.push(e.to);
}
Flow delta = std::min(excess[u], e.cap - e.flow);
e.flow += delta;
graph_[e.to][e.rev].flow -= delta;
excess[u] -= delta;
excess[e.to] += delta;
} else {
++curr[u];
}
}
}
return excess[sink];
}
private:
struct edge {
const int to, rev;
const Flow cap;
Flow flow;
constexpr edge(int _to, int _rev, Flow _cap)
: to(_to), rev(_rev), cap(_cap) {}
};
std::vector<std::vector<edge>> graph_;
};
// Successive shortest paths algorithm for minimum-cost flow.
template <typename Flow, typename Cost>
class min_cost_flow {
public:
min_cost_flow(int n) : m_(0), graph_(n) {}
void add_edge(int from, int to, Flow cap, Cost cost) {
graph_[from].emplace_back(to, (int)graph_[to].size(), cap, cost);
graph_[to].emplace_back(from, (int)graph_[from].size() - 1, 0, -cost);
m_ += 2;
}
std::pair<Flow, Cost> solve(
int source, int sink, Flow max_flow = std::numeric_limits<Flow>::max()) {
static constexpr Cost INF_COST = std::numeric_limits<Cost>::max();
std::vector<Cost> costs;
costs.reserve(m_);
for (std::vector<edge>& edges : graph_) {
for (edge& e : edges) {
e.flow = 0;
costs.push_back(e.cost);
}
}
int n = (int)graph_.size();
std::vector<Cost> dist(n, INF_COST);
dist[source] = 0;
std::vector<std::pair<int, edge*>> pred(n);
std::queue<int> q;
q.push(source);
std::vector<bool> active(n, false);
active[source] = true;
while (!q.empty()) { // Bellman-Ford shortest paths
int u = q.front();
q.pop();
active[u] = false;
for (edge& e : graph_[u]) {
if (e.flow < e.cap && dist[u] + e.cost < dist[e.to]) {
dist[e.to] = dist[u] + e.cost;
pred[e.to] = {u, &e};
if (!active[e.to]) {
q.push(e.to);
active[e.to] = true;
}
}
}
}
constexpr auto comp = [](const std::pair<int, Cost>& lhs,
const std::pair<int, Cost>& rhs) {
return lhs.second > rhs.second;
};
std::priority_queue<std::pair<int, Cost>, std::vector<std::pair<int, Cost>>,
decltype(comp)>
pq(comp);
Flow total_flow = 0;
Cost total_cost = 0, curr_cost = 0;
while (dist[sink] != INF_COST) {
Flow delta = max_flow - total_flow;
for (int v = sink; v != source;) {
const auto& [u, e] = pred[v];
delta = std::min(delta, e->cap - e->flow);
v = u;
}
for (int v = sink; v != source;) {
const auto& [u, e] = pred[v];
e->flow += delta;
graph_[v][e->rev].flow -= delta;
v = u;
}
total_cost += delta * (curr_cost += dist[sink]);
if ((total_flow += delta) == max_flow) {
break;
}
for (int u = 0; u < n; ++u) {
for (edge& e : graph_[u]) {
e.cost += dist[u] - dist[e.to];
}
}
std::fill(dist.begin(), dist.end(), INF_COST);
dist[source] = 0;
pq.emplace(source, 0);
while (!pq.empty()) { // Dijkstra's shortest paths
auto [u, d_u] = pq.top();
pq.pop();
if (d_u > dist[u]) {
continue;
}
for (edge& e : graph_[u]) {
if (e.flow < e.cap && dist[u] + e.cost < dist[e.to]) {
pq.emplace(e.to, dist[e.to] = dist[u] + e.cost);
pred[e.to] = {u, &e};
}
}
}
}
int k = 0;
for (std::vector<edge>& edges : graph_) {
for (edge& e : edges) {
e.cost = costs[k++];
}
}
return std::pair(total_flow, total_cost);
}
private:
struct edge {
const int to, rev;
const Flow cap;
Flow flow;
Cost cost;
constexpr edge(int _to, int _rev, Flow _cap, Cost _cost)
: to(_to), rev(_rev), cap(_cap), cost(_cost) {}
};
int m_;
std::vector<std::vector<edge>> graph_;
};
// Hungarian algorithm for the assignment problem. O(m ^ 2 * n).
template <typename Cost>
std::pair<Cost, std::vector<int>> min_cost_assign(
const std::vector<std::vector<Cost>>& costs) {
if (costs.empty() || costs[0].empty()) {
return std::pair(Cost(0), std::vector<int>(costs.size(), -1));
}
int m = (int)costs.size(), n = (int)costs[0].size();
bool transpose = m > n;
if (transpose) {
std::swap(m, n);
}
Cost total_cost = 0;
std::vector<bool> tight(n);
std::vector<int> col_match(n, -1), pred(n);
std::vector<Cost> row_pot(m, 0), col_pot(n, 0), slack(n);
for (int i = 0; i < m; ++i) {
std::fill(tight.begin(), tight.end(), false);
int y = -1;
for (int x = i; x != -1; x = col_match[y]) {
Cost delta = 0;
int next_y = -1;
for (int j = 0; j < n; ++j) {
if (tight[j]) {
continue;
}
Cost t =
(transpose ? costs[j][x] : costs[x][j]) - row_pot[x] - col_pot[j];
if (y == -1 || t < slack[j]) {
slack[j] = t;
pred[j] = y;
}
if (next_y == -1 || slack[j] < delta) {
delta = slack[j];
next_y = j;
}
}
y = next_y;
total_cost += delta;
row_pot[i] += delta;
for (int j = 0; j < n; ++j) {
if (tight[j]) {
row_pot[col_match[j]] += delta;
col_pot[j] -= delta;
} else {
slack[j] -= delta;
}
}
tight[y] = true;
}
while (pred[y] != -1) {
col_match[y] = col_match[pred[y]];
y = pred[y];
}
col_match[y] = i;
}
if (transpose) {
return std::pair(total_cost, col_match);
}
std::vector<int> row_match(m);
for (int j = 0; j < n; ++j) {
if (col_match[j] != -1) {
row_match[col_match[j]] = j;
}
}
return std::pair(total_cost, row_match);
}
// Simplex algorithm for linear programming.
// Maximizes c^T x subject to a x <= b and x >= 0.
template <typename T>
class linear_program {
static_assert(!std::is_integral_v<T>);
public:
linear_program(const std::vector<std::vector<T>>& a, const std::vector<T>& b,
const std::vector<T>& c,
T eps = std::sqrt(std::numeric_limits<T>::epsilon()))
: eps_(eps),
m_((int)b.size()),
n_((int)c.size()),
feasible_(true),
bounded_(true),
basic_(m_),
nonbasic_(n_ + 1),
tableau_(m_ + 2, std::vector<T>(n_ + 2, 0)) {
std::iota(basic_.begin(), basic_.end(), n_);
std::iota(nonbasic_.begin(), nonbasic_.begin() + n_, 0);
nonbasic_[n_] = -1;
for (int i = 0; i < m_; ++i) {
std::copy(a[i].begin(), a[i].end(), tableau_[i].begin());
tableau_[i][n_ + 1] = b[i];
}
int p = (int)(min_element(b.begin(), b.end()) - b.begin());
if (b[p] >= -eps_) {
for (int j = 0; j < n_; ++j) {
tableau_[m_][j] = -c[j];
}
simplex();
return;
}
for (int i = 0; i < m_; ++i) {
tableau_[i][n_] = -1;
}
tableau_[m_][n_] = 1;
for (int j = 0; j < n_; ++j) {
tableau_[m_ + 1][j] = -c[j];
}
pivot(p, n_);
simplex();
if (tableau_[m_][n_ + 1] < -eps_) {
feasible_ = false;
return;
}
int q;
if ((p = (int)(std::find(basic_.begin(), basic_.end(), -1) -
basic_.begin())) != m_) {
q = -1;
for (int j = 0; j <= n_; ++j) {
if (tableau_[p][j] < -eps_ &&
(q == -1 || nonbasic_[j] < nonbasic_[q])) {
q = j;
}
}
pivot(p, q);
} else {
q = (int)(std::find(nonbasic_.begin(), nonbasic_.end(), -1) -
nonbasic_.begin());
}
for (int i = 0; i < m_; ++i) {
tableau_[i][q] = 0;
}
std::fill(tableau_[m_].begin(), tableau_[m_].end(), 0);
tableau_[m_ + 1][q] = 0;
std::swap(tableau_[m_], tableau_[m_ + 1]);
simplex();
}
bool feasible() const { return feasible_; }
bool bounded() const { return bounded_; }
T solve() const { return tableau_[m_][n_ + 1]; }
std::vector<T> primal() const {
std::vector<T> x(n_, 0);
for (int i = 0; i < m_; ++i) {
if (basic_[i] < n_) {
x[basic_[i]] = tableau_[i][n_ + 1];
}
}
return x;
}
std::vector<T> dual() const {
std::vector<T> y(m_, 0);
for (int j = 0; j <= n_; ++j) {
if (nonbasic_[j] >= n_) {
y[nonbasic_[j] - n_] = tableau_[m_][j];
}
}
return y;
}
private:
const T eps_;
const int m_, n_;
bool feasible_, bounded_;
std::vector<int> basic_, nonbasic_;
std::vector<std::vector<T>> tableau_;
void pivot(int p, int q) {
int size1 = (int)tableau_.size(), size2 = (int)tableau_[0].size();
T inv_pq = 1 / tableau_[p][q];
for (int i = 0; i < size1; ++i) {
if (i != p) {
for (int j = 0; j < size2; ++j) {
if (j != q) {
tableau_[i][j] -= tableau_[i][q] * tableau_[p][j] * inv_pq;
}
}
}
}
for (int i = 0; i < size1; ++i) {
if (i != p) {
tableau_[i][q] *= -inv_pq;
}
}
for (int j = 0; j < size2; ++j) {
if (j != q) {
tableau_[p][j] *= inv_pq;
}
}
tableau_[p][q] = inv_pq;
std::swap(basic_[p], nonbasic_[q]);
}
void simplex() {
while (true) {
int q = -1;
for (int j = 0; j <= n_; ++j) {
if (tableau_[m_][j] < -eps_ &&
(q == -1 || nonbasic_[j] < nonbasic_[q])) {
q = j;
}
}
if (q == -1) {
return;
}
int p = -1;
T diff;
for (int i = 0; i < m_; ++i) {
if (tableau_[i][q] > eps_ &&
(p == -1 ||
(diff = tableau_[i][n_ + 1] / tableau_[i][q] -
tableau_[p][n_ + 1] / tableau_[p][q]) < -eps_ ||
(diff <= eps_ && basic_[i] < basic_[p]))) {
p = i;
}
}
if (p == -1) {
bounded_ = false;
return;
}
pivot(p, q);
}
}
};
constexpr int MODULUS = 1'000'000'007;
template <const int& Modulus = MODULUS>
class mod_int {
public:
constexpr mod_int() : val_(0) {}
template <typename T>
constexpr mod_int(T val) : val_((int)(val % Modulus)) {
static_assert(std::is_integral_v<T>);
if (val_ < 0) {
val_ += Modulus;
}
}
static constexpr int modulus() { return Modulus; }
template <typename T>
constexpr explicit operator T() const {
return (T)val_;
}
constexpr mod_int& operator++() {
if (++val_ == Modulus) {
val_ = 0;
}
return *this;
}
constexpr mod_int operator++(int) {
mod_int t = *this;
++*this;
return t;
}
constexpr mod_int& operator--() {
if (val_ == 0) {
val_ = Modulus;
}
--val_;
return *this;
}
constexpr mod_int operator--(int) {
mod_int t = *this;
--*this;
return t;
}
constexpr mod_int& operator+=(mod_int other) {
if ((val_ += other.val_ - Modulus) < 0) {
val_ += Modulus;
}
return *this;
}
constexpr mod_int& operator-=(mod_int other) {
if ((val_ -= other.val_) < 0) {
val_ += Modulus;
}
return *this;
}
constexpr mod_int& operator*=(mod_int other) {
val_ = (int)((long long)val_ * other.val_ % Modulus);
return *this;
}
constexpr mod_int& operator/=(mod_int other) { return *this *= inv(other); }
friend constexpr mod_int operator+(mod_int mi) { return mi; }
friend constexpr mod_int operator-(mod_int mi) {
if (mi.val_ != 0) {
mi.val_ = Modulus - mi.val_;
}
return mi;
}
friend constexpr mod_int operator+(mod_int lhs, mod_int rhs) {
return lhs += rhs;
}
friend constexpr mod_int operator-(mod_int lhs, mod_int rhs) {
return lhs -= rhs;
}
friend constexpr mod_int operator*(mod_int lhs, mod_int rhs) {
return lhs *= rhs;
}
friend constexpr mod_int operator/(mod_int lhs, mod_int rhs) {
return lhs /= rhs;
}
friend constexpr bool operator==(mod_int lhs, mod_int rhs) {
return lhs.val_ == rhs.val_;
}
friend constexpr bool operator!=(mod_int lhs, mod_int rhs) {
return lhs.val_ != rhs.val_;
}
friend constexpr bool operator<(mod_int lhs, mod_int rhs) {
return lhs.val_ < rhs.val_;
}
friend constexpr bool operator<=(mod_int lhs, mod_int rhs) {
return lhs.val_ <= rhs.val_;
}
friend constexpr bool operator>(mod_int lhs, mod_int rhs) {
return lhs.val_ > rhs.val_;
}
friend constexpr bool operator>=(mod_int lhs, mod_int rhs) {
return lhs.val_ >= rhs.val_;
}
friend std::ostream& operator<<(std::ostream& os, mod_int mi) {
return os << mi.val_;
}
friend std::istream& operator>>(std::istream& is, mod_int& mi) {
long long val;
is >> val;
mi = mod_int(val);
return is;
}
friend constexpr mod_int inv(mod_int mi) {
int a = mi.val_, b = Modulus;
int x = 0, xx = 1;
while (a != 0) {
int q = b / a;
int t = a;
a = b - q * a;
b = t;
t = xx;
xx = x - q * xx;
x = t;
}
assert(b == 1);
return mod_int(x);
}
template <typename T>
friend constexpr mod_int pow(mod_int base, T expo) {
static_assert(std::is_integral_v<T>);
if (expo < 0) {
base = inv(base);
expo = -expo;
}
mod_int power = 1;
for (; expo != 0; expo >>= 1) {
if (expo & 1) {
power *= base;
}
base *= base;
}
return power;
}
private:
int val_;
};
template <const int& Modulus = MODULUS>
class mod_comb {
public:
mod_comb(int n) : fact_{1}, inv_fact_{1} { resize(n); }
mod_int<Modulus> factorial(int n) {
if (n >= (int)fact_.size()) {
resize(n);
}
return fact_[n];
}
mod_int<Modulus> permute(int n, int k) {
if (k < 0 || k > n) {
return 0;
}
if (n >= (int)fact_.size()) {
resize(n);
}
return fact_[n] * inv_fact_[n - k];
}
mod_int<Modulus> choose(int n, int k) {
if (k < 0 || k > n) {
return 0;
}
if (n >= (int)fact_.size()) {
resize(n);
}
return fact_[n] * inv_fact_[n - k] * inv_fact_[k];
}
private:
std::vector<mod_int<Modulus>> fact_, inv_fact_;
void resize(int n) {
int m = (int)fact_.size();
fact_.reserve(n + 1);
for (int i = m; i <= n; ++i) {
fact_.push_back(fact_[i - 1] * i);
}
inv_fact_.resize(n + 1);
inv_fact_[n] = inv(fact_[n]);
for (int i = n; i > m; --i) {
inv_fact_[i - 1] = inv_fact_[i] * i;
}
}
};
// Extended Euclidean algorithm.
// a x + b y = gcd(a, b).
template <typename T>
constexpr T gcd(T a, T b, T& x, T& y) {
static_assert(std::is_integral_v<T>);
T xx = y = 1;
T yy = x = 0;
while (a != 0) {
T q = b / a;
T t = a;
a = b - q * a;
b = t;
t = xx;
xx = x - q * xx;
x = t;
t = yy;
yy = y - q * yy;
y = t;
}
if (b < 0) {
b = -b;
x = -x;
y = -y;
}
return b;
}
// Chinese remainder theorem.
// Returns {x, n} s.t. x = r[i] (mod m[i]) and n = lcm(m[i]).
// m[i] need not be pairwise coprime.
template <typename T>
std::optional<std::pair<T, T>> chinese_remainder(
const std::vector<std::pair<T, T>>& rem_mod) {
static_assert(std::is_integral_v<T>);
T x = 0, n = 1;
for (const auto& [r, m] : rem_mod) {
T inv_m, inv_n, g = gcd(m, n, inv_m, inv_n);
if ((r - x) % g != 0) {
return std::nullopt;
}
x += (r - x) % m * inv_n % m / g * n;
n *= m / g;
if (x < 0) {
x += n;
}
}
return std::pair(x, n);
}
// Miller-Rabin primality test.
constexpr bool is_prime(int n) {
if (n == 2) {
return true;
}
if (n <= 1 || (n & 1) == 0) {
return false;
}
int r = std::countr_zero((unsigned)(n - 1));
int d = (n - 1) >> r;
for (long long a : {2, 3, 5, 7}) {
if ((a == 3 && n < 2'047) || (a == 5 && n < 1'373'653) ||
(a == 7 && n < 25'326'001)) {
return true;
}
long long x = 1;
for (int b = d; b != 0; b >>= 1) {
if (b & 1) {
x = x * a % n;
}
a = a * a % n;
}
if (x == 1) {
continue;
}
for (int i = 1; x != n - 1; ++i) {
if (i == r || (x = x * x % n) == 1) {
return false;
}
}
}
return true;
}
#ifdef __SIZEOF_INT128__
// Miller-Rabin primality test.
constexpr bool is_prime(long long n) {
if (n == 2) {
return true;
}
if (n <= 1 || (n & 1) == 0) {
return false;
}
int r = std::countr_zero((unsigned long long)(n - 1));
long long d = (n - 1) >> r;
for (__int128 a : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
if ((a == 3 && n < 2'047LL) || (a == 5 && n < 1'373'653LL) ||
(a == 7 && n < 25'326'001LL) || (a == 11 && n < 3'215'031'751LL) ||
(a == 13 && n < 2'152'302'898'747LL) ||
(a == 17 && n < 3'474'749'660'383LL) ||
(a == 19 && n < 341'550'071'728'321LL) ||
(a == 29 && n < 3'825'123'056'546'413'051LL)) {
return true;
}
__int128 x = 1;
for (long long b = d; b != 0; b >>= 1) {
if (b & 1) {
x = x * a % n;
}
a = a * a % n;
}
if (x == 1) {
continue;
}
for (int i = 1; x != n - 1; ++i) {
if (i == r || (x = x * x % n) == 1) {
return false;
}
}
}
return true;
}
#endif
// Pollard's rho algorithm for integer factorization.
// n must be composite.
int factor(int n) {
std::mt19937 rng((std::uint_fast32_t)std::chrono::high_resolution_clock::now()
.time_since_epoch()
.count());
std::uniform_int_distribution<int> dist_c(1, n - 1), dist_x(0, n - 1);
int fact;
do {
int c = dist_c(rng), x = dist_x(rng), y = x;
int i = 1, k = 2;
do {
x = (int)(((long long)x * x + c) % n);
fact = std::gcd(x - y, n);
if (++i == k) {
y = x;
k <<= 1;
}
} while (fact == 1);
} while (fact == n);
return fact;
}
#ifdef __SIZEOF_INT128__
// Pollard's rho algorithm for integer factorization.
// n must be composite.
long long factor(long long n) {
std::mt19937_64 rng(
std::chrono::high_resolution_clock::now().time_since_epoch().count());
std::uniform_int_distribution<long long> dist_c(1, n - 1), dist_x(0, n - 1);
long long fact;
do {
long long c = dist_c(rng), x = dist_x(rng), y = x;
int i = 1, k = 2;
do {
x = (long long)(((__int128)x * x + c) % n);
fact = std::gcd(x - y, n);
if (++i == k) {
y = x;
k <<= 1;
}
} while (fact == 1);
} while (fact == n);
return fact;
}
#endif
// Sieve of Eratosthenes. O(n log log n) construction.
class sieve {
public:
sieve(int n) {
least_factor_.reserve(n + 1);
for (int i = 0; i <= n; ++i) {
least_factor_.push_back(i);
}
for (int i = 2; i * i <= n; ++i) {
if (least_factor_[i] == i) {
for (int j = i * i; j <= n; j += i) {
if (least_factor_[j] == j) {
least_factor_[j] = i;
}
}
}
}
for (int i = 2; i <= n; ++i) {
if (least_factor_[i] == i) {
primes_.push_back(i);
}
}
}
bool is_prime(int x) { return x >= 2 && least_factor_[x] == x; }
const std::vector<int>& primes() const { return primes_; }
template <typename T>
std::vector<std::pair<T, int>> factorize(T x) const {
static_assert(std::is_integral_v<T>);
std::vector<std::pair<T, int>> facts;
for (T prime : primes_) {
if (x < (int)least_factor_.size()) {
while (x > 1) {
prime = least_factor_[x];
int count = 0;
do {
x /= prime;
++count;
} while (least_factor_[x] == prime);
facts.emplace_back(prime, count);
}
return facts;
}
if (prime * prime > x) {
facts.emplace_back(x, 1);
return facts;
}
if (x % prime == 0) {
int count = 0;
do {
x /= prime;
++count;
} while (x % prime == 0);
facts.emplace_back(prime, count);
}
}
int k = (int)facts.size();
while (x > 1) {
T prime = x;
while (!::is_prime(prime)) {
prime = factor(prime);
}
int count = 0;
do {
x /= prime;
++count;
} while (x % prime == 0);
facts.emplace_back(prime, count);
}
std::sort(facts.begin() + k, facts.end());
return facts;
}
template <typename T>
std::vector<T> factors(T x) const {
std::vector<T> facts = {1};
for (const auto& [prime, count] : factorize(x)) {
int n = count * (int)facts.size();
for (int i = 0; i < n; ++i) {
facts.push_back(facts[i] * prime);
}
}
return facts;
}
private:
std::vector<int> primes_, least_factor_;
};
#ifdef __SIZEOF_INT128__
template <>
struct std::is_integral<__int128> : public std::true_type {};
namespace std {
constexpr __int128 abs(__int128 x) { return x < 0 ? -x : x; }
} // namespace std
template <>
constexpr __int128 std::gcd(__int128 a, __int128 b) {
a = std::abs(a);
b = std::abs(b);
while (a != 0) {
__int128 t = a;
a = b % a;
b = t;
}
return b;
}
template <>
constexpr __int128 std::lcm(__int128 a, __int128 b) {
if (a == 0 || b == 0) {
return 0;
}
return std::abs(a) / std::gcd(a, b) * std::abs(b);
}
std::ostream& operator<<(std::ostream& os, __int128 x) {
static constexpr long long BASE = 1'000'000'000'000'000'000LL;
char f = os.fill();
int w = (int)os.width();
if (__int128 mid = x / BASE; mid != 0) {
x = std::abs(x - mid * BASE);
if (int high = (int)(mid / BASE); high != 0) {
mid = std::abs(mid - (__int128)high * BASE);
os << high << std::setfill('0') << std::setw(18);
}
os << (long long)mid << std::setfill('0') << std::setw(18);
}
return os << (long long)x << std::setfill(f) << std::setw(w);
}
#endif
template <typename T>
class fraction {
static_assert(std::is_integral_v<T>);
public:
constexpr fraction(T num = 0) : num_(num), den_(1) {}
constexpr fraction(T num, T den) {
assert(den != 0);
T g = std::gcd(num, den);
if (den < 0) {
g = -g;
}
num_ = num / g;
den_ = den / g;
}
fraction(const std::vector<T>& cont) : num_(cont.back()), den_(1) {
for (int i = (int)cont.size() - 2; i >= 0; --i) {
std::swap(num_, den_);
num_ += cont[i] * den_;
}
}
constexpr T num() const { return num_; }
constexpr T den() const { return den_; }
template <typename U>
constexpr explicit operator U() const {
return (U)num_ / (U)den_;
}
constexpr fraction& operator+=(const fraction& other) {
return *this = *this + other;
}
constexpr fraction& operator-=(const fraction& other) {
return *this = *this - other;
}
constexpr fraction& operator*=(const fraction& other) {
T g1 = std::gcd(num_, other.den_), g2 = std::gcd(den_, other.num_);
num_ = (num_ / g1) * (other.num_ / g2);
den_ = (den_ / g2) * (other.den_ / g1);
return *this;
}
constexpr fraction& operator/=(const fraction& other) {
assert(other.num_ != 0);
T g1 = std::gcd(num_, other.num_), g2 = std::gcd(den_, other.den_);
if (other.num_ < 0) {
g1 = -g1;
}
num_ = (num_ / g1) * (other.den_ / g2);
den_ = (den_ / g2) * (other.num_ / g1);
return *this;
}
friend constexpr fraction operator+(fraction f) { return f; }
friend constexpr fraction operator-(fraction f) {
f.num_ = -f.num_;
return f;
}
friend constexpr fraction operator+(const fraction& lhs,
const fraction& rhs) {
T g = std::gcd(lhs.den_, rhs.den_);
return fraction(rhs.den_ / g * lhs.num_ + lhs.den_ / g * rhs.num_,
lhs.den_ / g * rhs.den_);
}
friend constexpr fraction operator-(const fraction& lhs,
const fraction& rhs) {
T g = std::gcd(lhs.den_, rhs.den_);
return fraction(rhs.den_ / g * lhs.num_ - lhs.den_ / g * rhs.num_,
lhs.den_ / g * rhs.den_);
}
friend constexpr fraction operator*(fraction lhs, const fraction& rhs) {
return lhs *= rhs;
}
friend constexpr fraction operator/(fraction lhs, const fraction& rhs) {
return lhs /= rhs;
}
friend constexpr bool operator==(const fraction& lhs, const fraction& rhs) {
return lhs.num_ == rhs.num_ && lhs.den_ == rhs.den_;
}
friend constexpr bool operator!=(const fraction& lhs, const fraction& rhs) {
return !(lhs == rhs);
}
friend constexpr bool operator<(const fraction& lhs, const fraction& rhs) {
return lhs.num_ * rhs.den_ < rhs.num_ * lhs.den_;
}
friend constexpr bool operator<=(const fraction& lhs, const fraction& rhs) {
return lhs.num_ * rhs.den_ <= rhs.num_ * lhs.den_;
}
friend constexpr bool operator>(const fraction& lhs, const fraction& rhs) {
return lhs.num_ * rhs.den_ > rhs.num_ * lhs.den_;
}
friend constexpr bool operator>=(const fraction& lhs, const fraction& rhs) {
return lhs.num_ * rhs.den_ >= rhs.num_ * lhs.den_;
}
friend std::ostream& operator<<(std::ostream& os, const fraction& f) {
return os << f.num_ << '/' << f.den_;
}
friend constexpr fraction inv(fraction f) {
assert(f.num_ != 0);
if (T t = f.num_; t < 0) {
f.num_ = -f.den_;
f.den_ = -t;
} else {
f.num_ = f.den_;
f.den_ = t;
}
return f;
}
friend std::vector<T> continued(fraction f) {
T q = f.num_ < 0 ? (f.num_ - f.den_ + 1) / f.den_ : f.num_ / f.den_;
f.num_ -= q * f.den_;
std::vector<T> cont = {q};
while (f.num_ != 0) {
std::swap(f.num_, f.den_);
q = f.num_ / f.den_;
f.num_ -= q * f.den_;
cont.push_back(q);
}
return cont;
}
private:
T num_, den_;
};
template <typename T>
class matrix {
public:
matrix(int m, int n, const T& val = T()) : m_(m), n_(n), data_(m * n, val) {}
matrix(std::initializer_list<std::initializer_list<T>> data)
: m_((int)data.size()), n_((int)data.begin()->size()) {
data_.reserve(m_ * n_);
for (const std::initializer_list<T>& row : data) {
assert((int)row.size() == n_);
std::copy(row.begin(), row.end(), std::back_inserter(data_));
}
}
static matrix identity(int n) {
matrix mat(n, n);
for (int i = 0; i < n; ++i) {
mat[i][i] = 1;
}
return mat;
}
int size1() const { return m_; }
int size2() const { return n_; }
T* operator[](int i) { return data_.data() + i * n_; }
const T* operator[](int i) const { return data_.data() + i * n_; }
matrix& operator+=(const matrix& other) {
assert(m_ == other.m_ && n_ == other.n_);
for (int i = 0; i < (int)data_.size(); ++i) {
data_[i] += other.data_[i];
}
return *this;
}
matrix& operator-=(const matrix& other) {
assert(m_ == other.m_ && n_ == other.n_);
for (int i = 0; i < (int)data_.size(); ++i) {
data_[i] -= other.data_[i];
}
return *this;
}
matrix& operator*=(const matrix& other) { return *this = *this * other; }
matrix& operator*=(const T& val) {
for (T& d : data_) {
d *= val;
}
return *this;
}
matrix& operator/=(const T& val) {
for (T& d : data_) {
d /= val;
}
return *this;
}
friend matrix operator+(matrix lhs, const matrix& rhs) { return lhs += rhs; }
friend matrix operator-(matrix lhs, const matrix& rhs) { return lhs -= rhs; }
friend matrix operator*(const matrix& lhs, const matrix& rhs) {
assert(lhs.n_ == rhs.m_);
matrix prod(lhs.m_, rhs.n_);
for (int i = 0; i < lhs.m_; ++i) {
for (int k = 0; k < lhs.n_; ++k) {
for (int j = 0; j < rhs.n_; ++j) {
prod[i][j] += lhs[i][k] * rhs[k][j];
}
}
}
return prod;
}
friend matrix operator*(matrix lhs, const T& val) { return lhs *= val; }
friend matrix operator*(const T& val, matrix rhs) { return rhs *= val; }
friend matrix operator/(matrix lhs, const T& val) { return lhs /= val; }
friend bool operator==(const matrix& lhs, const matrix& rhs) {
return lhs.m_ == rhs.m_ && lhs.data_ == rhs.data_;
}
friend bool operator!=(const matrix& lhs, const matrix& rhs) {
return !(lhs == rhs);
}
friend std::ostream& operator<<(std::ostream& os, const matrix& mat) {
os << '{';
for (int i = 0; i < mat.m_; ++i) {
os << (i == 0 ? "" : ", ") << '{';
for (int j = 0; j < mat.n_; ++j) {
os << (j == 0 ? "" : ", ") << mat[i][j];
}
os << '}';
}
return os << '}';
}
template <typename U>
friend matrix pow(matrix base, U expo) {
static_assert(std::is_integral_v<U>);
matrix power = identity(base.m_);
for (; expo > 0; expo >>= 1) {
if (expo & 1) {
power *= base;
}
base *= base;
}
return power;
}
private:
int m_, n_;
std::vector<T> data_;
};
// Berlekamp-Massey algorithm for linear recurrence.
template <const int& Modulus = MODULUS>
class linear_recur {
public:
linear_recur(const std::vector<mod_int<Modulus>>& seq)
: coeff_(0, 0), seq_(0, 0) {
int prev_i = -1;
mod_int<Modulus> prev_diff = 1;
std::vector<mod_int<Modulus>> prev_coeff = {1}, coeff = {1};
for (int i = 0; i < (int)seq.size(); ++i) {
mod_int<Modulus> diff;
for (int j = 0; j < (int)coeff.size(); ++j) {
diff += coeff[j] * seq[i - j];
}
if (diff == 0) {
continue;
}
int offset = i - prev_i;
mod_int<Modulus> q = diff / prev_diff;
if (prev_coeff.size() + offset <= coeff.size()) {
for (int j = 0; j < (int)prev_coeff.size(); ++j) {
coeff[j + offset] -= q * prev_coeff[j];
}
continue;
}
std::vector<mod_int<Modulus>> next_coeff(prev_coeff.size() + offset);
std::copy(coeff.begin(), coeff.end(), next_coeff.begin());
for (int j = 0; j < (int)prev_coeff.size(); ++j) {
next_coeff[j + offset] -= q * prev_coeff[j];
}
prev_i = i;
prev_diff = diff;
prev_coeff = std::move(coeff);
coeff = std::move(next_coeff);
}
int n = (int)coeff.size() - 1;
coeff_ = matrix<mod_int<Modulus>>(n, n);
for (int i = 0; i < n - 1; ++i) {
coeff_[i][i + 1] = 1;
}
for (int j = 0; j < n; ++j) {
coeff_[n - 1][j] = -coeff[n - j];
}
seq_ = matrix<mod_int<Modulus>>(n, 1);
std::copy_n(seq.begin(), n, seq_[0]);
}
mod_int<Modulus> solve(long long n) const {
int m = seq_.size1();
if (n < m) {
return seq_[(int)n][0];
}
return (pow(coeff_, n - m + 1) * seq_)[m - 1][0];
}
private:
matrix<mod_int<Modulus>> coeff_, seq_;
};
// Lagrange polynomial interpolation. O(n).
template <typename T>
T poly_interpolate(const std::vector<T>& vals, long long x) {
int n = (int)vals.size();
if (x < n) {
return vals[x];
}
T basis = 1;
for (int i = 1; i < n; ++i) {
basis *= T(x - i) / T(-i);
}
T y = vals[0] * basis;
for (int i = 1; i < n; ++i) {
basis *= (T(x - (i - 1)) / T(x - i)) * (T(i - n) / T(i));
y += vals[i] * basis;
}
return y;
}
// Five-point stencil for numerical differentiation.
template <typename T, typename Function>
constexpr auto derivative(Function f, T h) {
static_assert(std::is_floating_point_v<T>);
return [f, h](T x) {
return (-f(x + 2 * h) + 8 * f(x + h) - 8 * f(x - h) + f(x - 2 * h)) /
(12 * h);
};
}
// Newton's method for root finding.
template <typename T, typename Function, typename Derivative>
T find_root(Function f, Derivative df, T x, T min, T max, T eps, int max_iter) {
static_assert(std::is_floating_point_v<T>);
std::mt19937_64 thanos(
std::chrono::high_resolution_clock::now().time_since_epoch().count());
std::uniform_real_distribution<T> snap(min, max);
for (int resources = max_iter;;) {
T dx = f(x) / df(x);
if (!std::isfinite(dx) || resources-- == 0) { // needs correction
x = snap(thanos);
resources = max_iter;
continue;
}
x -= dx;
if (-eps <= dx && dx <= eps) {
return x;
}
}
}
namespace detail {
template <int ArgSign, typename T>
void fast_fourier_transform(std::vector<std::complex<T>>& a) {
static_assert(ArgSign == 1 || ArgSign == -1);
static_assert(std::is_floating_point_v<T>);
assert(std::has_single_bit(a.size()));
static constexpr T PI = (T)3.141592653589793238462643383279502884L;
int n = (int)a.size();
for (int i = 1, rev = 0; i < n; ++i) {
int bit = n;
do {
rev ^= (bit >>= 1);
} while ((rev & bit) == 0);
if (i < rev) {
std::swap(a[i], a[rev]);
}
}
std::vector<std::complex<T>> w(n >> 1);
w[0] = std::complex<T>(1.0, 0.0);
for (int k = 1; k < n; k <<= 1) {
T arg = ArgSign * PI / k;
std::complex<T> w_k(std::cos(arg), std::sin(arg));
for (int j = k - 2; j >= 0; j -= 2) {
w[j + 1] = w[j >> 1] * w_k;
w[j] = w[j >> 1];
}
for (int i = 0; i < n; i += k << 1) {
for (int j = 0; j < k; ++j) {
std::complex<T> t = w[j] * a[i + j + k];
a[i + j + k] = a[i + j] - t;
a[i + j] += t;
}
}
}
}
} // namespace detail
// Fast Fourier transform. O(n log n).
// a.size() must be a power of 2.
template <typename T>
void fft(std::vector<std::complex<T>>& a) {
detail::fast_fourier_transform<-1>(a);
}
// Inverse fast Fourier transform. O(n log n).
// a.size() must be a power of 2.
template <typename T>
void ifft(std::vector<std::complex<T>>& a) {
detail::fast_fourier_transform<1>(a);
T n = (T)a.size();
for (std::complex<T>& z : a) {
z /= n;
}
}
// Convolution. O(n log n).
template <typename T>
std::vector<T> convolve(const std::vector<T>& a, const std::vector<T>& b) {
static_assert(std::is_floating_point_v<T>);
int m = (int)a.size() + (int)b.size() - 1;
int n = std::bit_ceil((unsigned)m);
std::vector<std::complex<T>> c(n);
for (int i = 0; i < (int)a.size(); ++i) {
c[i].real(a[i]);
}
for (int i = 0; i < (int)b.size(); ++i) {
c[i].imag(b[i]);
}
fft(c);
c[0] = std::complex<T>(c[0].real() * c[0].imag(), 0.0);
for (int i = 1; i < n >> 1; ++i) {
std::complex<T> t = c[i] * c[i] - std::conj(c[n - i] * c[n - i]);
c[i] = std::complex<T>(0.25 * t.imag(), -0.25 * t.real());
c[n - i] = std::conj(c[i]);
}
c[n >> 1] = std::complex<T>(c[n >> 1].real() * c[n >> 1].imag(), 0.0);
ifft(c);
std::vector<T> conv;
conv.reserve(m);
for (int i = 0; i < m; ++i) {
conv.push_back(c[i].real());
}
return conv;
}
// Randomized polynomial rolling hash.
template <typename SequenceContainer>
class str_hasher {
static_assert(std::is_integral_v<typename SequenceContainer::value_type>);
public:
str_hasher(const SequenceContainer& s = SequenceContainer())
: muls_(rand2()) {
init(s);
}
str_hasher(const str_hasher& other, const SequenceContainer& s)
: muls_(other.muls_) {
init(s);
}
int size() const { return (int)pows_[0].size() - 1; }
void push_back(typename SequenceContainer::value_type c) {
for (int k = 0; k < 2; ++k) {
pows_[k].push_back(pows_[k].back() * muls_[k]);
prefs_[k].push_back(prefs_[k].back() * muls_[k] + c);
}
}
void pop_back() {
for (int k = 0; k < 2; ++k) {
pows_[k].pop_back();
prefs_[k].pop_back();
}
}
// Hashes the range [first, last).
long long hash(int first, int last) const {
assert(first < last);
std::array<mod_int<MOD>, 2> vals;
for (int k = 0; k < 2; ++k) {
vals[k] = prefs_[k][last] - prefs_[k][first] * pows_[k][last - first];
}
return ((long long)vals[0] << 32) + (long long)vals[1];
}
long long hash(const SequenceContainer& other) const {
std::array<mod_int<MOD>, 2> vals = {0, 0};
for (auto c : other) {
for (int k = 0; k < 2; ++k) {
vals[k] = vals[k] * muls_[k] + c;
}
}
return ((long long)vals[0] << 32) + (long long)vals[1];
}
private:
static constexpr int MOD = 2'147'483'647; // 2 ^ 31 - 1
const std::array<mod_int<MOD>, 2> muls_;
std::array<std::vector<mod_int<MOD>>, 2> pows_, prefs_;
static std::array<mod_int<MOD>, 2> rand2() {
std::mt19937 rng(
(std::uint_fast32_t)std::chrono::high_resolution_clock::now()
.time_since_epoch()
.count());
std::uniform_int_distribution<int> dist(1'000'000, 1'000'000'000);
return {dist(rng), dist(rng)};
}
void init(const SequenceContainer& s) {
int n = (int)s.size();
for (int k = 0; k < 2; ++k) {
pows_[k].reserve(n + 1);
pows_[k].emplace_back(1);
prefs_[k].reserve(n + 1);
prefs_[k].emplace_back(0);
}
for (auto c : s) {
push_back(c);
}
}
};
// Knuth-Morris-Pratt string-searching algorithm.
// O(m) preprocessing and O(n) matching.
template <typename SequenceContainer>
class str_searcher {
public:
str_searcher(const SequenceContainer& patt) : pattern_(patt) {
int m = (int)pattern_.size();
table_.reserve(m);
table_.push_back(0);
for (int i = 1, k = 0; i < m; ++i) {
while (k > 0 && pattern_[k] != pattern_[i]) {
k = table_[k - 1];
}
if (pattern_[k] == pattern_[i]) {
++k;
}
table_.push_back(k);
}
}
const SequenceContainer& pattern() const { return pattern_; }
const std::vector<int>& kmp_table() const { return table_; }
template <typename SequenceContainer::value_type Base = 'a', int Size = 26>
std::vector<std::array<int, Size>> automaton() const {
int m = (int)pattern_.size();
std::vector<std::array<int, Size>> transition(m + 1);
transition[0][pattern_[0] - Base] = 1;
for (int i = 1; i <= m; ++i) {
for (int j = 0; j < Size; ++j) {
if (i < m && pattern_[i] == Base + j) {
transition[i][j] = i + 1;
} else {
transition[i][j] = transition[table_[i - 1]][j];
}
}
}
return transition;
}
std::vector<int> search(const SequenceContainer& text) const {
int m = (int)pattern_.size(), n = (int)text.size();
std::vector<int> matches;
for (int i = 0, k = 0; i < n; ++i) {
while (k > 0 && pattern_[k] != text[i]) {
k = table_[k - 1];
}
if (pattern_[k] == text[i]) {
++k;
}
if (k == m) {
matches.push_back(i - m + 1);
k = table_[k - 1];
}
}
return matches;
}
private:
const SequenceContainer& pattern_;
std::vector<int> table_;
};
// Aho-Corasick string-searching algorithm.
template <typename SequenceContainer,
typename SequenceContainer::value_type Base = 'a', int Size = 26>
class multistr_searcher {
public:
multistr_searcher(const std::vector<SequenceContainer>& patterns)
: patterns_(patterns), trie_(1, node()) {
for (int i = 0; i < (int)patterns.size(); ++i) {
int v = 0;
for (auto c : patterns[i]) {
int j = c - Base;
if (trie_[v].child[j] == 0) {
trie_[v].child[j] = (int)trie_.size();
trie_.emplace_back();
}
v = trie_[v].child[j];
}
++trie_[v].count;
trie_[v].match = i;
}
std::vector<int> suffix(trie_.size(), 0);
std::queue<int> q;
for (int v : trie_[0].child) {
if (v != 0) {
q.push(v);
}
}
while (!q.empty()) {
int u = q.front();
q.pop();
for (int i = 0; i < Size; ++i) {
int& v = trie_[u].child[i];
int w = trie_[suffix[u]].child[i];
if (v == 0) {
v = w;
continue;
}
suffix[v] = w;
q.push(v);
trie_[v].count += trie_[w].count;
if (trie_[v].match == -1) {
trie_[v].match = trie_[w].match;
trie_[v].prev = trie_[w].prev;
} else if (trie_[w].match != -1) {
trie_[v].prev = w;
}
}
}
}
static constexpr int root() { return 0; }
int size() const { return (int)trie_.size(); }
const std::vector<SequenceContainer>& patterns() const { return patterns_; }
int next(int state, typename SequenceContainer::value_type c) const {
return trie_[state].child[c - Base];
}
int count(int state) const { return trie_[state].count; }
int match(int state) const { return trie_[state].match; }
std::vector<int> search(int state) const {
std::vector<int> matches;
if (match(state) != -1) {
do {
matches.push_back(match(state));
state = trie_[state].prev;
} while (state != -1);
}
return matches;
}
std::vector<int> search(const SequenceContainer& text) const {
std::vector<int> matches;
std::vector<bool> visited(trie_.size(), false);
int state = 0;
for (auto c : text) {
state = next(state, c);
if (match(state) != -1) {
for (int v = state; v != -1 && !visited[v]; v = trie_[v].prev) {
matches.push_back(match(v));
visited[v] = true;
}
}
}
return matches;
}
private:
struct node {
std::array<int, Size> child;
int count, match, prev;
constexpr node() : child(), count(0), match(-1), prev(-1) {}
};
const std::vector<SequenceContainer>& patterns_;
std::vector<node> trie_;
};
// Suffix array. O(n log n).
std::vector<int> suffix_array(std::string_view s) {
int n = (int)s.size();
std::vector<int> count(128, 0);
for (char c : s) {
++count[c];
}
std::partial_sum(count.begin(), count.end(), count.begin());
std::vector<int> suff(n), rank(n);
for (int i = 0; i < n; ++i) {
suff[--count[s[i]]] = i;
}
rank[suff[0]] = 0;
for (int i = 1; i < n; ++i) {
if (s[suff[i]] == s[suff[i - 1]]) {
rank[suff[i]] = rank[suff[i - 1]];
} else {
rank[suff[i]] = i;
}
}
std::vector<int> prev(n), pos(n);
for (int len = 1; len < n; len <<= 1) {
prev.swap(suff);
std::iota(pos.begin(), pos.end(), 0);
for (int i = n - len; i < n; ++i) {
suff[pos[rank[i]]++] = i;
}
for (int i = 0; i < n; ++i) {
if (int j = prev[i] - len; j >= 0) {
suff[pos[rank[j]]++] = j;
}
}
prev.swap(rank);
rank[suff[0]] = 0;
for (int i = 1; i < n; ++i) {
if (suff[i] + len < n && suff[i - 1] + len < n &&
prev[suff[i]] == prev[suff[i - 1]] &&
prev[suff[i] + len] == prev[suff[i - 1] + len]) {
rank[suff[i]] = rank[suff[i - 1]];
} else {
rank[suff[i]] = i;
}
}
}
return suff;
}
// Longest common prefix array. O(n).
std::vector<int> lcp_array(std::string_view s, const std::vector<int>& suff) {
int n = (int)s.size();
std::vector<int> rank(n);
for (int i = 0; i < n; ++i) {
rank[suff[i]] = i;
}
std::vector<int> lcp(n - 1);
for (int i = 0, k = 0; i < n; ++i) {
if (rank[i] == n - 1) {
k = 0;
continue;
}
int j = suff[rank[i] + 1];
while (i + k < n && j + k < n && s[i + k] == s[j + k]) {
++k;
}
lcp[rank[i]] = k;
if (k > 0) {
--k;
}
}
return lcp;
}
// Z-Function. O(n).
// Returns maximum lengths z s.t. s.substr(i, z[i]) == s.substr(0, z[i]).
template <typename SequenceContainer>
std::vector<int> z_fun(const SequenceContainer& s) {
int n = (int)s.size();
std::vector<int> z(n, 0);
int left = 0, right = 1;
for (int i = 1; i < n; ++i) {
if (i < right) {
z[i] = std::min(right - i, z[i - left]);
}
while (i + z[i] < n && s[i + z[i]] == s[z[i]]) {
++z[i];
}
if (i + z[i] > right) {
left = i;
right = i + z[i];
}
}
return z;
}
// Manacher's algorithm for the longest palindromic substring. O(n).
// Returns the half lengths of the palindromes for all 2n - 1 centers.
template <typename SequenceContainer>
std::vector<int> palindromic_substr(const SequenceContainer& s) {
int n = (int)s.size();
std::vector<int> p(2 * n - 1, 0);
int left = 0, right = 0;
for (int i = 0; i < 2 * n - 1; ++i) {
int l = (i + 1) >> 1, r = i >> 1;
if (l < right) {
p[i] = std::min(right - l, p[2 * (left + right) - i]);
}
while (0 <= l - p[i] - 1 && r + p[i] + 1 < n &&
s[l - p[i] - 1] == s[r + p[i] + 1]) {
++p[i];
}
if (r + p[i] > right) {
left = l - p[i];
right = r + p[i];
}
}
return p;
}
template <typename T, int Dim = 2>
struct point;
template <typename T>
struct point<T, 2> {
T x, y;
constexpr point(const T& _x = T(), const T& _y = T()) : x(_x), y(_y) {}
constexpr point& operator+=(const point& other) {
x += other.x;
y += other.y;
return *this;
}
constexpr point& operator-=(const point& other) {
x -= other.x;
y -= other.y;
return *this;
}
constexpr point& operator*=(const T& val) {
x *= val;
y *= val;
return *this;
}
constexpr point& operator/=(const T& val) {
x /= val;
y /= val;
return *this;
}
friend constexpr point operator+(const point& p) { return p; }
friend constexpr point operator-(const point& p) { return point(-p.x, -p.y); }
friend constexpr point operator+(const point& lhs, const point& rhs) {
return point(lhs.x + rhs.x, lhs.y + rhs.y);
}
friend constexpr point operator-(const point& lhs, const point& rhs) {
return point(lhs.x - rhs.x, lhs.y - rhs.y);
}
friend constexpr point operator*(const point& lhs, const T& val) {
return point(lhs.x * val, lhs.y * val);
}
friend constexpr point operator*(const T& val, const point& rhs) {
return point(val * rhs.x, val * rhs.y);
}
friend constexpr point operator/(const point& lhs, const T& val) {
return point(lhs.x / val, lhs.y / val);
}
friend constexpr bool operator==(const point& lhs, const point& rhs) {
return lhs.x == rhs.x && lhs.y == rhs.y;
}
friend constexpr bool operator!=(const point& lhs, const point& rhs) {
return !(lhs == rhs);
}
friend constexpr bool operator<(const point& lhs, const point& rhs) {
return lhs.x != rhs.x ? lhs.x < rhs.x : lhs.y < rhs.y;
}
friend std::ostream& operator<<(std::ostream& os, const point& p) {
return os << '(' << p.x << ',' << p.y << ')';
}
friend constexpr T norm(const point& p) { return p.x * p.x + p.y * p.y; }
friend constexpr T dot(const point& lhs, const point& rhs) {
return lhs.x * rhs.x + lhs.y * rhs.y;
}
friend constexpr T cross(const point& lhs, const point& rhs) {
return lhs.x * rhs.y - lhs.y * rhs.x;
}
};
template <typename T>
struct point<T, 3> {
T x, y, z;
constexpr point(const T& _x = T(), const T& _y = T(), const T& _z = T())
: x(_x), y(_y), z(_z) {}
constexpr point& operator+=(const point& other) {
x += other.x;
y += other.y;
z += other.z;
return *this;
}
constexpr point& operator-=(const point& other) {
x -= other.x;
y -= other.y;
z -= other.z;
return *this;
}
constexpr point& operator*=(const T& val) {
x *= val;
y *= val;
z *= val;
return *this;
}
constexpr point& operator/=(const T& val) {
x /= val;
y /= val;
z /= val;
return *this;
}
friend constexpr point operator+(const point& p) { return p; }
friend constexpr point operator-(const point& p) {
return point(-p.x, -p.y, -p.z);
}
friend constexpr point operator+(const point& lhs, const point& rhs) {
return point(lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z);
}
friend constexpr point operator-(const point& lhs, const point& rhs) {
return point(lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z);
}
friend constexpr point operator*(const point& lhs, const T& val) {
return point(lhs.x * val, lhs.y * val, lhs.z * val);
}
friend constexpr point operator*(const T& val, const point& rhs) {
return point(val * rhs.x, val * rhs.y, val * rhs.z);
}
friend constexpr point operator/(const point& lhs, const T& val) {
return point(lhs.x / val, lhs.y / val, lhs.z / val);
}
friend constexpr bool operator==(const point& lhs, const point& rhs) {
return lhs.x == rhs.x && lhs.y == rhs.y && lhs.z == rhs.z;
}
friend constexpr bool operator!=(const point& lhs, const point& rhs) {
return !(lhs == rhs);
}
friend constexpr bool operator<(const point& lhs, const point& rhs) {
if (lhs.x != rhs.x) {
return lhs.x < rhs.x;
}
if (lhs.y != rhs.y) {
return lhs.y < rhs.y;
}
return lhs.z < rhs.z;
}
friend std::ostream& operator<<(std::ostream& os, const point& p) {
return os << '(' << p.x << ',' << p.y << ',' << p.z << ')';
}
friend constexpr T norm(const point& p) {
return p.x * p.x + p.y * p.y + p.z * p.z;
}
friend constexpr T dot(const point& lhs, const point& rhs) {
return lhs.x * rhs.x + lhs.y * rhs.y + lhs.z * rhs.z;
}
friend constexpr point cross(const point& lhs, const point& rhs) {
return point(lhs.y * rhs.z - lhs.z * rhs.y, lhs.z * rhs.x - lhs.x * rhs.z,
lhs.x * rhs.y - lhs.y * rhs.x);
}
};
template <typename T>
int sign(T x, T eps = std::sqrt(std::numeric_limits<T>::epsilon())) {
return x > eps ? 1 : (x < -eps ? -1 : 0);
}
template <typename T>
struct line {
point<T, 2> first, second;
constexpr line() = default;
constexpr line(const point<T, 2>& _first, const point<T, 2>& _second)
: first(_first), second(_second) {}
friend std::ostream& operator<<(std::ostream& os, const line& l) {
return os << '{' << l.first << ", " << l.second << '}';
}
friend constexpr point<T, 2> proj(const line& l, const point<T, 2>& q) {
const auto& [p1, p2] = l;
point<T, 2> seg = p2 - p1;
return p1 + dot(q - p1, seg) * seg / norm(seg);
}
friend constexpr bool parallel(
const line& l1, const line& l2,
T eps = std::sqrt(std::numeric_limits<T>::epsilon())) {
const auto& [p1, p2] = l1;
const auto& [p3, p4] = l2;
T c = cross(p2 - p1, p4 - p3);
return -eps <= c && c <= eps;
}
friend constexpr bool intersects(
const line& l1, const line& l2,
T eps = std::sqrt(std::numeric_limits<T>::epsilon())) {
const auto& [p1, p2] = l1;
const auto& [p3, p4] = l2;
T c1 = cross(p3 - p1, p2 - p1);
T c2 = cross(p4 - p1, p2 - p1);
if (-eps <= c1 && c1 <= eps && -eps <= c2 && c2 <= eps) {
return std::min(p1.x, p2.x) <= std::max(p3.x, p4.x) &&
std::min(p3.x, p4.x) <= std::max(p1.x, p2.x);
} else {
return sign(c1, eps) != sign(c2, eps) &&
sign(cross(p1 - p3, p4 - p3), eps) !=
sign(cross(p2 - p3, p4 - p3), eps);
}
}
friend constexpr point<T, 2> intersection(const line& l1, const line& l2) {
const auto& [p1, p2] = l1;
const auto& [p3, p4] = l2;
T num1 = p1.x * p2.y - p1.y * p2.x;
T num2 = p3.x * p4.y - p3.y * p4.x;
T den = (p1.x - p2.x) * (p3.y - p4.y) - (p1.y - p2.y) * (p3.x - p4.x);
return point<T, 2>((num1 * (p3.x - p4.x) - num2 * (p1.x - p2.x)) / den,
((num1 * (p3.y - p4.y) - num2 * (p1.y - p2.y)) / den));
}
};
// Andrew's monotone chain convex hull algorithm. O(n log n).
template <typename T>
std::vector<point<T, 2>> lower_hull(std::vector<point<T, 2>> points) {
std::sort(points.begin(), points.end());
std::vector<point<T, 2>> hull;
for (const point<T, 2>& p : points) {
while ((int)hull.size() >= 2 &&
cross(hull.back() - hull.end()[-2], p - hull.back()) <= 0) {
hull.pop_back();
}
hull.push_back(p);
}
return hull;
}
// Andrew's monotone chain convex hull algorithm. O(n log n).
template <typename T>
std::vector<point<T, 2>> upper_hull(std::vector<point<T, 2>> points) {
std::sort(points.begin(), points.end());
std::vector<point<T, 2>> hull;
for (const point<T, 2>& p : points) {
while ((int)hull.size() >= 2 &&
cross(hull.back() - hull.end()[-2], p - hull.back()) >= 0) {
hull.pop_back();
}
hull.push_back(p);
}
return hull;
}
// Andrew's monotone chain convex hull algorithm. O(n log n).
// Returns convex hull vertices in counterclockwise order.
template <typename T>
std::vector<point<T, 2>> convex_hull(std::vector<point<T, 2>> points) {
std::sort(points.begin(), points.end());
int n = (int)points.size();
std::vector<point<T, 2>> hull;
for (int i = 0; i < n; ++i) { // lower hull
while ((int)hull.size() >= 2 &&
cross(hull.back() - hull.end()[-2], points[i] - hull.back()) <= 0) {
hull.pop_back();
}
hull.push_back(points[i]);
}
int k = (int)hull.size() + 1;
for (int i = n - 2; i >= 0; --i) { // upper hull
while ((int)hull.size() >= k &&
cross(hull.back() - hull.end()[-2], points[i] - hull.back()) <= 0) {
hull.pop_back();
}
hull.push_back(points[i]);
}
hull.pop_back();
return hull;
}
template <typename T>
class monotone_hull {
public:
using iterator = typename std::deque<point<T, 2>>::iterator;
using const_iterator = typename std::deque<point<T, 2>>::const_iterator;
monotone_hull() = default;
monotone_hull(std::vector<point<T, 2>> lines) {
std::sort(lines.begin(), lines.end());
for (const point<T, 2>& line : lines) {
push_back(line);
}
}
int size() const { return (int)lines_.size(); }
iterator begin() { return lines_.begin(); }
const_iterator begin() const { return lines_.begin(); }
iterator end() { return lines_.end(); }
const_iterator end() const { return lines_.end(); }
void push_back(const point<T, 2>& line) {
while ((int)lines_.size() >= 2 &&
cross(lines_.back() - lines_.end()[-2], line - lines_.back()) >= 0) {
lines_.pop_back();
}
lines_.push_back(line);
}
void emplace_back(const T& slope, const T& intercept) {
push_back(point<T, 2>(slope, intercept));
}
T linear_max(const T& x) {
point<T, 2> p(x, 1);
while ((int)lines_.size() >= 2 && dot(lines_[0], p) < dot(lines_[1], p)) {
lines_.pop_front();
}
return dot(lines_[0], p);
}
private:
std::deque<point<T, 2>> lines_;
};
template <typename T>
class dynamic_hull {
public:
struct line {
public:
constexpr operator const point<T, 2>&() const { return m_b_; }
constexpr const point<T, 2>& get() const { return m_b_; }
friend constexpr bool operator<(const line& lhs, const line& rhs) {
return lhs.m_b_ < rhs.m_b_;
}
friend constexpr bool operator<(const line& lhs, const T& x) {
return lhs.next_ != nullptr &&
dot(lhs.m_b_, point<T, 2>(x, 1)) <
dot(lhs.next_->m_b_, point<T, 2>(x, 1));
}
friend std::ostream& operator<<(std::ostream& os, const line& l) {
return os << l.m_b_;
}
private:
friend dynamic_hull;
point<T, 2> m_b_;
mutable const line* next_;
constexpr line(const point<T, 2>& m_b) : m_b_(m_b), next_(nullptr) {}
};
using iterator = typename std::set<line, std::less<void>>::iterator;
using const_iterator =
typename std::set<line, std::less<void>>::const_iterator;
dynamic_hull() = default;
dynamic_hull(const dynamic_hull& other) : lines_(other.lines_) {
for (auto it = lines_.begin(); std::next(it) != lines_.end(); ++it) {
it->next_ = &*std::next(it);
}
}
dynamic_hull(dynamic_hull&&) = default;
dynamic_hull& operator=(const dynamic_hull& other) {
lines_ = other.lines_;
for (auto it = lines_.begin(); std::next(it) != lines_.end(); ++it) {
it->next_ = &*std::next(it);
}
return *this;
}
dynamic_hull& operator=(dynamic_hull&&) = default;
int size() const { return (int)lines_.size(); }
iterator begin() { return lines_.begin(); }
const_iterator begin() const { return lines_.begin(); }
iterator end() { return lines_.end(); }
const_iterator end() const { return lines_.end(); }
void insert(const point<T, 2>& slope_intercept) {
auto [it, inserted] = lines_.insert(slope_intercept);
if (!inserted) {
return;
}
if (it != lines_.begin() && std::next(it) != lines_.end() &&
cross(it->m_b_ - std::prev(it)->m_b_, std::next(it)->m_b_ - it->m_b_) >=
0) {
lines_.erase(it);
return;
}
if (it != lines_.begin()) {
auto left = std::prev(it);
while (left != lines_.begin() && cross(left->m_b_ - std::prev(left)->m_b_,
it->m_b_ - left->m_b_) >= 0) {
left = std::prev(lines_.erase(left));
}
left->next_ = &*it;
}
if (auto right = std::next(it); right != lines_.end()) {
while (std::next(right) != lines_.end() &&
cross(right->m_b_ - it->m_b_,
std::next(right)->m_b_ - right->m_b_) >= 0) {
right = lines_.erase(right);
}
it->next_ = &*right;
}
}
void emplace(const T& slope, const T& intercept) {
insert(point<T, 2>(slope, intercept));
}
T linear_max(const T& x) const {
return dot(lines_.lower_bound(x)->m_b_, point<T, 2>(x, 1));
}
private:
std::set<line, std::less<void>> lines_;
};
int main() {
auto start_time = std::chrono::steady_clock::now();
std::mt19937 rng((std::uint_fast32_t)start_time.time_since_epoch().count());
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
std::cout.precision(24);
#ifdef TYLER
std::cerr << "Duration: "
<< std::chrono::duration_cast<std::chrono::milliseconds>(
std::chrono::steady_clock::now() - start_time)
.count()
<< "ms\n";
#endif
return 0;
}