convex_hull_trick.hpp¶
维护斜率单调加入的直线，并在查询横坐标也单调时求最优值。
Complexity: Time: O(log n) per insertion or monotone query. Space: O(n).
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(log n) per insertion or monotone query.
/// Space: O(n).
#include <cassert>
#include <limits>
#include <map>
#include <set>
#include <type_traits>
#include <utility>
#include <vector>
namespace noya {
/// @brief Monotone convex-hull trick for decreasing slopes and queries.
/// Obsolete middle lines are removed by cross multiplication. Because query
/// coordinates are also decreasing, the current optimum moves monotonically
/// along the same stack and can be popped permanently. Integral coefficients
/// may use the full range of types up to 64 bits; query evaluations must fit T.
template <class T> struct linear_cht_min {
std::vector<std::pair<T, T>> stk;
/// @brief Add line y = px + q. p must be decreasing.
void add(T p, T q) {
if (stk.size() >= 1) {
assert(stk.back().first >= p);
}
if (stk.size() >= 1 && stk.back().first == p) {
if (stk.back().second <= q) {
return;
} else {
stk.pop_back();
}
}
while ((int)stk.size() >= 2) {
const auto [p2, q2] = stk.end()[-2];
const auto [p1, q1] = stk.end()[-1];
if (keep_middle(p2, q2, p1, q1, p, q)) {
break;
}
stk.pop_back();
}
stk.emplace_back(p, q);
}
/// @brief Query minimum at x. x must be decreasing.
T get(T x) {
while ((int)stk.size() >= 2) {
const auto [p2, q2] = stk.end()[-2];
const auto [p1, q1] = stk.end()[-1];
if (p2 * x + q2 > p1 * x + q1) {
break;
} else {
stk.pop_back();
}
}
if (stk.empty()) {
return std::numeric_limits<T>::max();
} else {
auto [k, b] = stk.back();
return k * x + b;
}
}
private:
static bool keep_middle(T p2, T q2, T p1, T q1, T p, T q) {
if constexpr (std::numeric_limits<T>::is_integer) {
static_assert(sizeof(T) <= 8, "linear_cht_min supports up to 64-bit integers");
__int128 a = __int128(q1) - __int128(q2);
__int128 b = __int128(q) - __int128(q1);
// Slopes are strictly decreasing here. Each difference has at most 64
// magnitude bits, so its unsigned product fits even beyond signed 128.
if ((a < 0) != (b < 0)) {
return a < 0;
}
using U = unsigned __int128;
U left = U(a < 0 ? -a : a) * U(__int128(p1) - __int128(p));
U right = U(b < 0 ? -b : b) * U(__int128(p2) - __int128(p1));
return a < 0 ? left > right : left < right;
} else {
using W = long double;
return (W(q1) - W(q2)) * (W(p1) - W(p)) <
(W(q) - W(q1)) * (W(p2) - W(p1));
}
}
};
template <typename T> struct Line {
mutable T k, m, p;
bool operator<(const Line &o) const { return k < o.k; }
bool operator<(T x) const { return p < x; }
};
template <typename T> T lc_inf() { return std::numeric_limits<T>::max(); }
template <> inline long double lc_inf<long double>() { return 1 / .0; }
template <typename T> T lc_div(T a, T b) {
return a / b - ((a ^ b) < 0 and a % b);
}
template <> inline long double lc_div(long double a, long double b) { return a / b; }
template <> inline double lc_div(double a, double b) { return a / b; }
/// @brief Dynamic line container with logarithmic insertion and queries.
/// Lines are slope-ordered in a multiset and store the first coordinate where
/// they dominate the next line. Insertion repairs neighboring intersections;
/// heterogeneous lower_bound then locates the active interval for a query.
template <typename T, bool MIN = true>
struct line_container : std::multiset<Line<T>, std::less<>> {
using super = std::multiset<Line<T>, std::less<>>;
using super::begin, super::end, super::insert, super::erase;
using super::empty, super::lower_bound;
T inf = lc_inf<T>();
bool insect(typename super::iterator x, typename super::iterator y) {
if (y == end())
return x->p = inf, false;
if (x->k == y->k)
x->p = (x->m > y->m ? inf : -inf);
else
x->p = lc_div(y->m - x->m, x->k - y->k);
return x->p >= y->p;
}
void add(T k, T m) {
if (MIN) {
k = -k, m = -m;
}
auto z = insert({k, m, 0}), y = z++, x = y;
while (insect(y, z))
z = erase(z);
if (x != begin() and insect(--x, y))
insect(x, y = erase(y));
while ((y = x) != begin() and (--x)->p >= y->p)
insect(x, erase(y));
}
T query(T x) {
assert(!empty());
auto l = *lower_bound(x);
T v = (l.k * x + l.m);
return (MIN ? -v : v);
}
};
template <typename T> using cht_min = line_container<T, true>;
template <typename T> using cht_max = line_container<T, false>;
template <typename T> struct cht_xy {
static_assert(std::is_same_v<T, long long> || std::is_floating_point_v<T>);
using ld = long double;
cht_min<ld> cm;
cht_max<ld> cma;
T amx = std::numeric_limits<T>::lowest(), amn = std::numeric_limits<T>::max();
T bmx = std::numeric_limits<T>::lowest(), bmn = std::numeric_limits<T>::max();
int ix = -1, iy = -1;
int jx = -1, jy = -1;
bool emp = true;
std::map<std::pair<T, T>, int> MP;
void clear() {
emp = true;
cm.clear();
cma.clear();
MP.clear();
amx = bmx = std::numeric_limits<T>::lowest();
amn = bmn = std::numeric_limits<T>::max();
ix = iy = jx = jy = -1;
}
void add(T a, T b, int i = -1) {
if (emp) {
amx = amn = a;
bmx = bmn = b;
ix = iy = jx = jy = i;
}
emp = false;
cm.add(b, a);
cma.add(b, a);
std::pair<T, T> p = {a, b};
MP[p] = i;
if (amx < a)
amx = a, ix = i;
if (amn > a)
amn = a, iy = i;
if (bmx < b)
bmx = b, jx = i;
if (bmn > b)
bmn = b, jy = i;
}
std::pair<T, int> get_max(T x, T y) {
if (cm.empty())
return {std::numeric_limits<T>::lowest(), -1};
if (x == 0) {
if (bmx * y > bmn * y) {
return {bmx * y, jx};
}
return {bmn * y, jy};
}
ld z = ld(y) / x;
if (x > 0) {
auto l = cma.lower_bound(z);
T a = l->m, b = l->k;
std::pair<T, T> p = {a, b};
int idx = MP[p];
return {a * x + b * y, idx};
}
auto l = cm.lower_bound(z);
T a = -(l->m), b = -(l->k);
std::pair<T, T> p = {a, b};
int idx = MP[p];
return {a * x + b * y, idx};
}
std::pair<T, int> get_min(T x, T y) {
if (cm.empty())
return {std::numeric_limits<T>::max(), -1};
auto [f, i] = get_max(-x, -y);
return {-f, i};
}
};
} // namespace noya
// NOYA_CONVEX_HULL_TRICK_HPP
#ifndef NOYA_CONVEX_HULL_TRICK_HPP
#define NOYA_CONVEX_HULL_TRICK_HPP 1
/// @complexity Time: O(log n) per insertion or monotone query.
/// Space: O(n).
#include <cassert>
#include <limits>
#include <map>
#include <set>
#include <type_traits>
#include <utility>
#include <vector>
namespace noya {
/// @brief Monotone convex-hull trick for decreasing slopes and queries.
/// Obsolete middle lines are removed by cross multiplication. Because query
/// coordinates are also decreasing, the current optimum moves monotonically
/// along the same stack and can be popped permanently. Integral coefficients
/// may use the full range of types up to 64 bits; query evaluations must fit T.
template <class T> struct linear_cht_min {
std::vector<std::pair<T, T>> stk;
/// @brief Add line y = px + q. p must be decreasing.
void add(T p, T q) {
if (stk.size() >= 1) {
assert(stk.back().first >= p);
}
if (stk.size() >= 1 && stk.back().first == p) {
if (stk.back().second <= q) {
return;
} else {
stk.pop_back();
}
}
while ((int)stk.size() >= 2) {
const auto [p2, q2] = stk.end()[-2];
const auto [p1, q1] = stk.end()[-1];
if (keep_middle(p2, q2, p1, q1, p, q)) {
break;
}
stk.pop_back();
}
stk.emplace_back(p, q);
}
/// @brief Query minimum at x. x must be decreasing.
T get(T x) {
while ((int)stk.size() >= 2) {
const auto [p2, q2] = stk.end()[-2];
const auto [p1, q1] = stk.end()[-1];
if (p2 * x + q2 > p1 * x + q1) {
break;
} else {
stk.pop_back();
}
}
if (stk.empty()) {
return std::numeric_limits<T>::max();
} else {
auto [k, b] = stk.back();
return k * x + b;
}
}
private:
static bool keep_middle(T p2, T q2, T p1, T q1, T p, T q) {
if constexpr (std::numeric_limits<T>::is_integer) {
static_assert(sizeof(T) <= 8, "linear_cht_min supports up to 64-bit integers");
__int128 a = __int128(q1) - __int128(q2);
__int128 b = __int128(q) - __int128(q1);
// Slopes are strictly decreasing here. Each difference has at most 64
// magnitude bits, so its unsigned product fits even beyond signed 128.
if ((a < 0) != (b < 0)) {
return a < 0;
}
using U = unsigned __int128;
U left = U(a < 0 ? -a : a) * U(__int128(p1) - __int128(p));
U right = U(b < 0 ? -b : b) * U(__int128(p2) - __int128(p1));
return a < 0 ? left > right : left < right;
} else {
using W = long double;
return (W(q1) - W(q2)) * (W(p1) - W(p)) <
(W(q) - W(q1)) * (W(p2) - W(p1));
}
}
};
template <typename T> struct Line {
mutable T k, m, p;
bool operator<(const Line &o) const { return k < o.k; }
bool operator<(T x) const { return p < x; }
};
template <typename T> T lc_inf() { return std::numeric_limits<T>::max(); }
template <> inline long double lc_inf<long double>() { return 1 / .0; }
template <typename T> T lc_div(T a, T b) {
return a / b - ((a ^ b) < 0 and a % b);
}
template <> inline long double lc_div(long double a, long double b) { return a / b; }
template <> inline double lc_div(double a, double b) { return a / b; }
/// @brief Dynamic line container with logarithmic insertion and queries.
/// Lines are slope-ordered in a multiset and store the first coordinate where
/// they dominate the next line. Insertion repairs neighboring intersections;
/// heterogeneous lower_bound then locates the active interval for a query.
template <typename T, bool MIN = true>
struct line_container : std::multiset<Line<T>, std::less<>> {
using super = std::multiset<Line<T>, std::less<>>;
using super::begin, super::end, super::insert, super::erase;
using super::empty, super::lower_bound;
T inf = lc_inf<T>();
bool insect(typename super::iterator x, typename super::iterator y) {
if (y == end())
return x->p = inf, false;
if (x->k == y->k)
x->p = (x->m > y->m ? inf : -inf);
else
x->p = lc_div(y->m - x->m, x->k - y->k);
return x->p >= y->p;
}
void add(T k, T m) {
if (MIN) {
k = -k, m = -m;
}
auto z = insert({k, m, 0}), y = z++, x = y;
while (insect(y, z))
z = erase(z);
if (x != begin() and insect(--x, y))
insect(x, y = erase(y));
while ((y = x) != begin() and (--x)->p >= y->p)
insect(x, erase(y));
}
T query(T x) {
assert(!empty());
auto l = *lower_bound(x);
T v = (l.k * x + l.m);
return (MIN ? -v : v);
}
};
template <typename T> using cht_min = line_container<T, true>;
template <typename T> using cht_max = line_container<T, false>;
template <typename T> struct cht_xy {
static_assert(std::is_same_v<T, long long> || std::is_floating_point_v<T>);
using ld = long double;
cht_min<ld> cm;
cht_max<ld> cma;
T amx = std::numeric_limits<T>::lowest(), amn = std::numeric_limits<T>::max();
T bmx = std::numeric_limits<T>::lowest(), bmn = std::numeric_limits<T>::max();
int ix = -1, iy = -1;
int jx = -1, jy = -1;
bool emp = true;
std::map<std::pair<T, T>, int> MP;
void clear() {
emp = true;
cm.clear();
cma.clear();
MP.clear();
amx = bmx = std::numeric_limits<T>::lowest();
amn = bmn = std::numeric_limits<T>::max();
ix = iy = jx = jy = -1;
}
void add(T a, T b, int i = -1) {
if (emp) {
amx = amn = a;
bmx = bmn = b;
ix = iy = jx = jy = i;
}
emp = false;
cm.add(b, a);
cma.add(b, a);
std::pair<T, T> p = {a, b};
MP[p] = i;
if (amx < a)
amx = a, ix = i;
if (amn > a)
amn = a, iy = i;
if (bmx < b)
bmx = b, jx = i;
if (bmn > b)
bmn = b, jy = i;
}
std::pair<T, int> get_max(T x, T y) {
if (cm.empty())
return {std::numeric_limits<T>::lowest(), -1};
if (x == 0) {
if (bmx * y > bmn * y) {
return {bmx * y, jx};
}
return {bmn * y, jy};
}
ld z = ld(y) / x;
if (x > 0) {
auto l = cma.lower_bound(z);
T a = l->m, b = l->k;
std::pair<T, T> p = {a, b};
int idx = MP[p];
return {a * x + b * y, idx};
}
auto l = cm.lower_bound(z);
T a = -(l->m), b = -(l->k);
std::pair<T, T> p = {a, b};
int idx = MP[p];
return {a * x + b * y, idx};
}
std::pair<T, int> get_min(T x, T y) {
if (cm.empty())
return {std::numeric_limits<T>::max(), -1};
auto [f, i] = get_max(-x, -y);
return {-f, i};
}
};
} // namespace noya
#endif // NOYA_CONVEX_HULL_TRICK_HPP
/// @complexity Time: O(log n) per insertion or monotone query.
/// Space: O(n).
#include <cassert>
#include <limits>
#include <map>
#include <set>
#include <type_traits>
#include <utility>
#include <vector>
namespace noya {
/// @brief Monotone convex-hull trick for decreasing slopes and queries.
/// Obsolete middle lines are removed by cross multiplication. Because query
/// coordinates are also decreasing, the current optimum moves monotonically
/// along the same stack and can be popped permanently. Integral coefficients
/// may use the full range of types up to 64 bits; query evaluations must fit T.
template <class T> struct linear_cht_min {
std::vector<std::pair<T, T>> stk;
/// @brief Add line y = px + q. p must be decreasing.
void add(T p, T q) {
if (stk.size() >= 1) {
assert(stk.back().first >= p);
}
if (stk.size() >= 1 && stk.back().first == p) {
if (stk.back().second <= q) {
return;
} else {
stk.pop_back();
}
}
while ((int)stk.size() >= 2) {
const auto [p2, q2] = stk.end()[-2];
const auto [p1, q1] = stk.end()[-1];
if (keep_middle(p2, q2, p1, q1, p, q)) {
break;
}
stk.pop_back();
}
stk.emplace_back(p, q);
}
/// @brief Query minimum at x. x must be decreasing.
T get(T x) {
while ((int)stk.size() >= 2) {
const auto [p2, q2] = stk.end()[-2];
const auto [p1, q1] = stk.end()[-1];
if (p2 * x + q2 > p1 * x + q1) {
break;
} else {
stk.pop_back();
}
}
if (stk.empty()) {
return std::numeric_limits<T>::max();
} else {
auto [k, b] = stk.back();
return k * x + b;
}
}
private:
static bool keep_middle(T p2, T q2, T p1, T q1, T p, T q) {
if constexpr (std::numeric_limits<T>::is_integer) {
static_assert(sizeof(T) <= 8, "linear_cht_min supports up to 64-bit integers");
__int128 a = __int128(q1) - __int128(q2);
__int128 b = __int128(q) - __int128(q1);
// Slopes are strictly decreasing here. Each difference has at most 64
// magnitude bits, so its unsigned product fits even beyond signed 128.
if ((a < 0) != (b < 0)) {
return a < 0;
}
using U = unsigned __int128;
U left = U(a < 0 ? -a : a) * U(__int128(p1) - __int128(p));
U right = U(b < 0 ? -b : b) * U(__int128(p2) - __int128(p1));
return a < 0 ? left > right : left < right;
} else {
using W = long double;
return (W(q1) - W(q2)) * (W(p1) - W(p)) <
(W(q) - W(q1)) * (W(p2) - W(p1));
}
}
};
template <typename T> struct Line {
mutable T k, m, p;
bool operator<(const Line &o) const { return k < o.k; }
bool operator<(T x) const { return p < x; }
};
template <typename T> T lc_inf() { return std::numeric_limits<T>::max(); }
template <> inline long double lc_inf<long double>() { return 1 / .0; }
template <typename T> T lc_div(T a, T b) {
return a / b - ((a ^ b) < 0 and a % b);
}
template <> inline long double lc_div(long double a, long double b) { return a / b; }
template <> inline double lc_div(double a, double b) { return a / b; }
/// @brief Dynamic line container with logarithmic insertion and queries.
/// Lines are slope-ordered in a multiset and store the first coordinate where
/// they dominate the next line. Insertion repairs neighboring intersections;
/// heterogeneous lower_bound then locates the active interval for a query.
template <typename T, bool MIN = true>
struct line_container : std::multiset<Line<T>, std::less<>> {
using super = std::multiset<Line<T>, std::less<>>;
using super::begin, super::end, super::insert, super::erase;
using super::empty, super::lower_bound;
T inf = lc_inf<T>();
bool insect(typename super::iterator x, typename super::iterator y) {
if (y == end())
return x->p = inf, false;
if (x->k == y->k)
x->p = (x->m > y->m ? inf : -inf);
else
x->p = lc_div(y->m - x->m, x->k - y->k);
return x->p >= y->p;
}
void add(T k, T m) {
if (MIN) {
k = -k, m = -m;
}
auto z = insert({k, m, 0}), y = z++, x = y;
while (insect(y, z))
z = erase(z);
if (x != begin() and insect(--x, y))
insect(x, y = erase(y));
while ((y = x) != begin() and (--x)->p >= y->p)
insect(x, erase(y));
}
T query(T x) {
assert(!empty());
auto l = *lower_bound(x);
T v = (l.k * x + l.m);
return (MIN ? -v : v);
}
};
template <typename T> using cht_min = line_container<T, true>;
template <typename T> using cht_max = line_container<T, false>;
template <typename T> struct cht_xy {
static_assert(std::is_same_v<T, long long> || std::is_floating_point_v<T>);
using ld = long double;
cht_min<ld> cm;
cht_max<ld> cma;
T amx = std::numeric_limits<T>::lowest(), amn = std::numeric_limits<T>::max();
T bmx = std::numeric_limits<T>::lowest(), bmn = std::numeric_limits<T>::max();
int ix = -1, iy = -1;
int jx = -1, jy = -1;
bool emp = true;
std::map<std::pair<T, T>, int> MP;
void clear() {
emp = true;
cm.clear();
cma.clear();
MP.clear();
amx = bmx = std::numeric_limits<T>::lowest();
amn = bmn = std::numeric_limits<T>::max();
ix = iy = jx = jy = -1;
}
void add(T a, T b, int i = -1) {
if (emp) {
amx = amn = a;
bmx = bmn = b;
ix = iy = jx = jy = i;
}
emp = false;
cm.add(b, a);
cma.add(b, a);
std::pair<T, T> p = {a, b};
MP[p] = i;
if (amx < a)
amx = a, ix = i;
if (amn > a)
amn = a, iy = i;
if (bmx < b)
bmx = b, jx = i;
if (bmn > b)
bmn = b, jy = i;
}
std::pair<T, int> get_max(T x, T y) {
if (cm.empty())
return {std::numeric_limits<T>::lowest(), -1};
if (x == 0) {
if (bmx * y > bmn * y) {
return {bmx * y, jx};
}
return {bmn * y, jy};
}
ld z = ld(y) / x;
if (x > 0) {
auto l = cma.lower_bound(z);
T a = l->m, b = l->k;
std::pair<T, T> p = {a, b};
int idx = MP[p];
return {a * x + b * y, idx};
}
auto l = cm.lower_bound(z);
T a = -(l->m), b = -(l->k);
std::pair<T, T> p = {a, b};
int idx = MP[p];
return {a * x + b * y, idx};
}
std::pair<T, int> get_min(T x, T y) {
if (cm.empty())
return {std::numeric_limits<T>::max(), -1};
auto [f, i] = get_max(-x, -y);
return {-f, i};
}
};
} // namespace noya
// NOYA_CONVEX_HULL_TRICK_HPP