convex_hull_trick.hpp¶
支持单调 CHT、任意顺序的动态直线容器与二维点积最值；整数交点精确比较，点积查询可返回最优点编号。
Complexity: Time: amortized \(O(1)\) per monotone insertion/query; amortized \(O(\log n)\) per dynamic insertion and \(O(\log n)\) per dynamic/XY query. Space: \(O(n)\).
Source: maspy (line container and XY queries)
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: amortized O(1) per monotone insertion/query;
/// amortized O(log n) per dynamic insertion and O(log n) per dynamic/XY query.
/// Space: O(n).
#include <cassert>
#include <cstdint>
#include <limits>
#include <set>
#include <type_traits>
#include <utility>
#include <vector>
namespace noya {
namespace cht_internal {
// Affine evaluations need twice the coefficient width. Unsigned products
// also fit in unsigned 128 bits after adding one 64-bit intercept.
template <class T>
using wide = std::conditional_t<std::numeric_limits<T>::is_integer,
std::conditional_t<std::is_signed_v<T>, __int128, unsigned __int128>, T>;
// Intersections have at most 64 magnitude bits in numerator/denominator.
// A zero denominator represents a signed infinity.
struct ratio {
std::uint64_t num = 0, den = 1;
bool neg = false;
ratio(__int128 n = 0, __int128 d = 1)
: num(n < 0 ? -n : n), den(d < 0 ? -d : d),
neg(n != 0 && ((n < 0) != (d < 0))) {}
ratio operator-() const {
ratio r = *this;
if (num) r.neg = !neg;
return r;
}
bool operator<(const ratio &r) const {
if (!den) return neg && (r.den || !r.neg);
if (!r.den) return !r.neg;
if (neg != r.neg) return neg;
using U = unsigned __int128;
U a = U(num) * r.den, b = U(r.num) * den;
return neg ? a > b : a < b;
}
};
template <class T, bool Fractional = true>
using point = std::conditional_t<std::numeric_limits<T>::is_integer,
std::conditional_t<Fractional, ratio, __int128>, T>;
} // namespace cht_internal
/// @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 results 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 (eval(p2, q2, x) > eval(p1, q1, x)) {
break;
} else {
stk.pop_back();
}
}
if (stk.empty()) {
return std::numeric_limits<T>::max();
} else {
auto [k, b] = stk.back();
return T(eval(k, b, x));
}
}
private:
static cht_internal::wide<T> eval(T k, T b, T x) {
return cht_internal::wide<T>(k) * x + b;
}
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, bool Min = false, bool Fractional = false> struct Line {
T k, m;
mutable cht_internal::point<T, Fractional> p;
int id = -1;
bool operator<(const Line &o) const { return Min ? k > o.k : k < o.k; }
bool operator<(const cht_internal::point<T, Fractional> &x) const { return p < x; }
};
template <typename T> T lc_inf() {
if constexpr (std::is_floating_point_v<T>) return std::numeric_limits<T>::infinity();
else return std::numeric_limits<T>::max();
}
template <typename T> T lc_div(T a, T b) {
if constexpr (std::is_floating_point_v<T>) return a / b;
else return a / b - ((a ^ b) < 0 && a % b);
}
/// @brief Dynamic line container with logarithmic insertion and queries.
/// Integer queries use floored intersections; Fractional retains exact ratios
/// for cht_xy, where rounding could discard a line optimal between integers.
/// Insertions use descending slopes for minima and ascending slopes for maxima.
/// Integral coefficients and coordinates may use up to 64 bits; results must fit T.
template <typename T, bool MIN = true, bool Fractional = false>
struct line_container : std::multiset<Line<T, MIN, Fractional>, std::less<>> {
static_assert(!std::numeric_limits<T>::is_integer || sizeof(T) <= 8);
using super = std::multiset<Line<T, MIN, Fractional>, std::less<>>;
using point_type = cht_internal::point<T, Fractional>;
using super::begin, super::end, super::insert, super::erase;
using super::empty, super::lower_bound;
point_type inf = [] {
if constexpr (std::numeric_limits<T>::is_integer && Fractional) return point_type(1, 0);
else return lc_inf<point_type>();
}();
bool insect(typename super::iterator x, typename super::iterator y) {
if (y == end()) return x->p = inf, false;
if (x->k == y->k) {
bool better = MIN ? x->m < y->m : x->m > y->m;
x->p = better ? inf : -inf;
} else if constexpr (std::numeric_limits<T>::is_integer) {
__int128 a = __int128(y->m) - x->m, b = __int128(x->k) - y->k;
if constexpr (Fractional) x->p = point_type(a, b);
else {
// Differences have at most 64 magnitude bits, so a 64-bit division suffices.
std::uint64_t u = a < 0 ? -a : a, v = b < 0 ? -b : b;
__int128 q = u / v;
x->p = (a < 0) != (b < 0) ? -q - (u % v != 0) : q;
}
} else {
x->p = (y->m - x->m) / (x->k - y->k);
}
return !(x->p < y->p);
}
void add(T k, T m, int id = -1) {
auto z = insert({k, m, {}, id}), y = z++, x = y;
while (insect(y, z)) z = erase(z);
if (x != begin() && insect(--x, y)) insect(x, y = erase(y));
while ((y = x) != begin() && !((--x)->p < y->p)) insect(x, erase(y));
}
T query(T x) const {
assert(!empty());
auto l = lower_bound(point_type(x));
return T(cht_internal::wide<T>(l->k) * x + l->m);
}
};
template <typename T> using cht_min = line_container<T, true>;
template <typename T> using cht_max = line_container<T, false>;
/// @brief Maintain points (a,b), querying min/max a*x+b*y and an attaining id.
/// Integer directions and intersections stay exact; returned values must fit T.
template <typename T> struct cht_xy {
static_assert(std::is_same_v<T, long long> || std::is_floating_point_v<T>);
using C = std::conditional_t<std::is_floating_point_v<T>, long double, T>;
using point_type = cht_internal::point<C>;
using W = cht_internal::wide<T>;
line_container<C, true, true> cm;
line_container<C, false, true> cma;
T bmx = std::numeric_limits<T>::lowest(), bmn = std::numeric_limits<T>::max();
int jx = -1, jy = -1;
bool emp = true;
void clear() {
emp = true;
cm.clear();
cma.clear();
bmx = std::numeric_limits<T>::lowest();
bmn = std::numeric_limits<T>::max();
jx = jy = -1;
}
void add(T a, T b, int i = -1) {
if (emp) {
bmx = bmn = b;
jx = jy = i;
}
emp = false;
cm.add(b, a, i);
cma.add(b, a, 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) const { return get<false>(x, y); }
std::pair<T, int> get_min(T x, T y) const { return get<true>(x, y); }
private:
template <bool Min> std::pair<T, int> get(T x, T y) const {
if (cm.empty())
return {Min ? std::numeric_limits<T>::max() : std::numeric_limits<T>::lowest(), -1};
if (x == 0) {
bool upper = Min ? y < 0 : y > 0;
return upper ? std::pair<T, int>{T(W(bmx) * y), jx}
: std::pair<T, int>{T(W(bmn) * y), jy};
}
point_type z;
if constexpr (std::numeric_limits<T>::is_integer) z = point_type(y, x);
else z = C(y) / x;
auto eval = [&](auto it) {
T a = it->m, b = it->k;
return std::pair<T, int>{T(W(a) * x + W(b) * y), it->id};
};
if ((x > 0) == Min) return eval(cm.lower_bound(z));
return eval(cma.lower_bound(z));
}
};
} // namespace noya
// NOYA_CONVEX_HULL_TRICK_HPP
#ifndef NOYA_CONVEX_HULL_TRICK_HPP
#define NOYA_CONVEX_HULL_TRICK_HPP 1
/// @complexity Time: amortized O(1) per monotone insertion/query;
/// amortized O(log n) per dynamic insertion and O(log n) per dynamic/XY query.
/// Space: O(n).
#include <cassert>
#include <cstdint>
#include <limits>
#include <set>
#include <type_traits>
#include <utility>
#include <vector>
namespace noya {
namespace cht_internal {
// Affine evaluations need twice the coefficient width. Unsigned products
// also fit in unsigned 128 bits after adding one 64-bit intercept.
template <class T>
using wide = std::conditional_t<std::numeric_limits<T>::is_integer,
std::conditional_t<std::is_signed_v<T>, __int128, unsigned __int128>, T>;
// Intersections have at most 64 magnitude bits in numerator/denominator.
// A zero denominator represents a signed infinity.
struct ratio {
std::uint64_t num = 0, den = 1;
bool neg = false;
ratio(__int128 n = 0, __int128 d = 1)
: num(n < 0 ? -n : n), den(d < 0 ? -d : d),
neg(n != 0 && ((n < 0) != (d < 0))) {}
ratio operator-() const {
ratio r = *this;
if (num) r.neg = !neg;
return r;
}
bool operator<(const ratio &r) const {
if (!den) return neg && (r.den || !r.neg);
if (!r.den) return !r.neg;
if (neg != r.neg) return neg;
using U = unsigned __int128;
U a = U(num) * r.den, b = U(r.num) * den;
return neg ? a > b : a < b;
}
};
template <class T, bool Fractional = true>
using point = std::conditional_t<std::numeric_limits<T>::is_integer,
std::conditional_t<Fractional, ratio, __int128>, T>;
} // namespace cht_internal
/// @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 results 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 (eval(p2, q2, x) > eval(p1, q1, x)) {
break;
} else {
stk.pop_back();
}
}
if (stk.empty()) {
return std::numeric_limits<T>::max();
} else {
auto [k, b] = stk.back();
return T(eval(k, b, x));
}
}
private:
static cht_internal::wide<T> eval(T k, T b, T x) {
return cht_internal::wide<T>(k) * x + b;
}
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, bool Min = false, bool Fractional = false> struct Line {
T k, m;
mutable cht_internal::point<T, Fractional> p;
int id = -1;
bool operator<(const Line &o) const { return Min ? k > o.k : k < o.k; }
bool operator<(const cht_internal::point<T, Fractional> &x) const { return p < x; }
};
template <typename T> T lc_inf() {
if constexpr (std::is_floating_point_v<T>) return std::numeric_limits<T>::infinity();
else return std::numeric_limits<T>::max();
}
template <typename T> T lc_div(T a, T b) {
if constexpr (std::is_floating_point_v<T>) return a / b;
else return a / b - ((a ^ b) < 0 && a % b);
}
/// @brief Dynamic line container with logarithmic insertion and queries.
/// Integer queries use floored intersections; Fractional retains exact ratios
/// for cht_xy, where rounding could discard a line optimal between integers.
/// Insertions use descending slopes for minima and ascending slopes for maxima.
/// Integral coefficients and coordinates may use up to 64 bits; results must fit T.
template <typename T, bool MIN = true, bool Fractional = false>
struct line_container : std::multiset<Line<T, MIN, Fractional>, std::less<>> {
static_assert(!std::numeric_limits<T>::is_integer || sizeof(T) <= 8);
using super = std::multiset<Line<T, MIN, Fractional>, std::less<>>;
using point_type = cht_internal::point<T, Fractional>;
using super::begin, super::end, super::insert, super::erase;
using super::empty, super::lower_bound;
point_type inf = [] {
if constexpr (std::numeric_limits<T>::is_integer && Fractional) return point_type(1, 0);
else return lc_inf<point_type>();
}();
bool insect(typename super::iterator x, typename super::iterator y) {
if (y == end()) return x->p = inf, false;
if (x->k == y->k) {
bool better = MIN ? x->m < y->m : x->m > y->m;
x->p = better ? inf : -inf;
} else if constexpr (std::numeric_limits<T>::is_integer) {
__int128 a = __int128(y->m) - x->m, b = __int128(x->k) - y->k;
if constexpr (Fractional) x->p = point_type(a, b);
else {
// Differences have at most 64 magnitude bits, so a 64-bit division suffices.
std::uint64_t u = a < 0 ? -a : a, v = b < 0 ? -b : b;
__int128 q = u / v;
x->p = (a < 0) != (b < 0) ? -q - (u % v != 0) : q;
}
} else {
x->p = (y->m - x->m) / (x->k - y->k);
}
return !(x->p < y->p);
}
void add(T k, T m, int id = -1) {
auto z = insert({k, m, {}, id}), y = z++, x = y;
while (insect(y, z)) z = erase(z);
if (x != begin() && insect(--x, y)) insect(x, y = erase(y));
while ((y = x) != begin() && !((--x)->p < y->p)) insect(x, erase(y));
}
T query(T x) const {
assert(!empty());
auto l = lower_bound(point_type(x));
return T(cht_internal::wide<T>(l->k) * x + l->m);
}
};
template <typename T> using cht_min = line_container<T, true>;
template <typename T> using cht_max = line_container<T, false>;
/// @brief Maintain points (a,b), querying min/max a*x+b*y and an attaining id.
/// Integer directions and intersections stay exact; returned values must fit T.
template <typename T> struct cht_xy {
static_assert(std::is_same_v<T, long long> || std::is_floating_point_v<T>);
using C = std::conditional_t<std::is_floating_point_v<T>, long double, T>;
using point_type = cht_internal::point<C>;
using W = cht_internal::wide<T>;
line_container<C, true, true> cm;
line_container<C, false, true> cma;
T bmx = std::numeric_limits<T>::lowest(), bmn = std::numeric_limits<T>::max();
int jx = -1, jy = -1;
bool emp = true;
void clear() {
emp = true;
cm.clear();
cma.clear();
bmx = std::numeric_limits<T>::lowest();
bmn = std::numeric_limits<T>::max();
jx = jy = -1;
}
void add(T a, T b, int i = -1) {
if (emp) {
bmx = bmn = b;
jx = jy = i;
}
emp = false;
cm.add(b, a, i);
cma.add(b, a, 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) const { return get<false>(x, y); }
std::pair<T, int> get_min(T x, T y) const { return get<true>(x, y); }
private:
template <bool Min> std::pair<T, int> get(T x, T y) const {
if (cm.empty())
return {Min ? std::numeric_limits<T>::max() : std::numeric_limits<T>::lowest(), -1};
if (x == 0) {
bool upper = Min ? y < 0 : y > 0;
return upper ? std::pair<T, int>{T(W(bmx) * y), jx}
: std::pair<T, int>{T(W(bmn) * y), jy};
}
point_type z;
if constexpr (std::numeric_limits<T>::is_integer) z = point_type(y, x);
else z = C(y) / x;
auto eval = [&](auto it) {
T a = it->m, b = it->k;
return std::pair<T, int>{T(W(a) * x + W(b) * y), it->id};
};
if ((x > 0) == Min) return eval(cm.lower_bound(z));
return eval(cma.lower_bound(z));
}
};
} // namespace noya
#endif // NOYA_CONVEX_HULL_TRICK_HPP
/// @complexity Time: amortized O(1) per monotone insertion/query;
/// amortized O(log n) per dynamic insertion and O(log n) per dynamic/XY query.
/// Space: O(n).
#include <cassert>
#include <cstdint>
#include <limits>
#include <set>
#include <type_traits>
#include <utility>
#include <vector>
namespace noya {
namespace cht_internal {
// Affine evaluations need twice the coefficient width. Unsigned products
// also fit in unsigned 128 bits after adding one 64-bit intercept.
template <class T>
using wide = std::conditional_t<std::numeric_limits<T>::is_integer,
std::conditional_t<std::is_signed_v<T>, __int128, unsigned __int128>, T>;
// Intersections have at most 64 magnitude bits in numerator/denominator.
// A zero denominator represents a signed infinity.
struct ratio {
std::uint64_t num = 0, den = 1;
bool neg = false;
ratio(__int128 n = 0, __int128 d = 1)
: num(n < 0 ? -n : n), den(d < 0 ? -d : d),
neg(n != 0 && ((n < 0) != (d < 0))) {}
ratio operator-() const {
ratio r = *this;
if (num) r.neg = !neg;
return r;
}
bool operator<(const ratio &r) const {
if (!den) return neg && (r.den || !r.neg);
if (!r.den) return !r.neg;
if (neg != r.neg) return neg;
using U = unsigned __int128;
U a = U(num) * r.den, b = U(r.num) * den;
return neg ? a > b : a < b;
}
};
template <class T, bool Fractional = true>
using point = std::conditional_t<std::numeric_limits<T>::is_integer,
std::conditional_t<Fractional, ratio, __int128>, T>;
} // namespace cht_internal
/// @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 results 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 (eval(p2, q2, x) > eval(p1, q1, x)) {
break;
} else {
stk.pop_back();
}
}
if (stk.empty()) {
return std::numeric_limits<T>::max();
} else {
auto [k, b] = stk.back();
return T(eval(k, b, x));
}
}
private:
static cht_internal::wide<T> eval(T k, T b, T x) {
return cht_internal::wide<T>(k) * x + b;
}
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, bool Min = false, bool Fractional = false> struct Line {
T k, m;
mutable cht_internal::point<T, Fractional> p;
int id = -1;
bool operator<(const Line &o) const { return Min ? k > o.k : k < o.k; }
bool operator<(const cht_internal::point<T, Fractional> &x) const { return p < x; }
};
template <typename T> T lc_inf() {
if constexpr (std::is_floating_point_v<T>) return std::numeric_limits<T>::infinity();
else return std::numeric_limits<T>::max();
}
template <typename T> T lc_div(T a, T b) {
if constexpr (std::is_floating_point_v<T>) return a / b;
else return a / b - ((a ^ b) < 0 && a % b);
}
/// @brief Dynamic line container with logarithmic insertion and queries.
/// Integer queries use floored intersections; Fractional retains exact ratios
/// for cht_xy, where rounding could discard a line optimal between integers.
/// Insertions use descending slopes for minima and ascending slopes for maxima.
/// Integral coefficients and coordinates may use up to 64 bits; results must fit T.
template <typename T, bool MIN = true, bool Fractional = false>
struct line_container : std::multiset<Line<T, MIN, Fractional>, std::less<>> {
static_assert(!std::numeric_limits<T>::is_integer || sizeof(T) <= 8);
using super = std::multiset<Line<T, MIN, Fractional>, std::less<>>;
using point_type = cht_internal::point<T, Fractional>;
using super::begin, super::end, super::insert, super::erase;
using super::empty, super::lower_bound;
point_type inf = [] {
if constexpr (std::numeric_limits<T>::is_integer && Fractional) return point_type(1, 0);
else return lc_inf<point_type>();
}();
bool insect(typename super::iterator x, typename super::iterator y) {
if (y == end()) return x->p = inf, false;
if (x->k == y->k) {
bool better = MIN ? x->m < y->m : x->m > y->m;
x->p = better ? inf : -inf;
} else if constexpr (std::numeric_limits<T>::is_integer) {
__int128 a = __int128(y->m) - x->m, b = __int128(x->k) - y->k;
if constexpr (Fractional) x->p = point_type(a, b);
else {
// Differences have at most 64 magnitude bits, so a 64-bit division suffices.
std::uint64_t u = a < 0 ? -a : a, v = b < 0 ? -b : b;
__int128 q = u / v;
x->p = (a < 0) != (b < 0) ? -q - (u % v != 0) : q;
}
} else {
x->p = (y->m - x->m) / (x->k - y->k);
}
return !(x->p < y->p);
}
void add(T k, T m, int id = -1) {
auto z = insert({k, m, {}, id}), y = z++, x = y;
while (insect(y, z)) z = erase(z);
if (x != begin() && insect(--x, y)) insect(x, y = erase(y));
while ((y = x) != begin() && !((--x)->p < y->p)) insect(x, erase(y));
}
T query(T x) const {
assert(!empty());
auto l = lower_bound(point_type(x));
return T(cht_internal::wide<T>(l->k) * x + l->m);
}
};
template <typename T> using cht_min = line_container<T, true>;
template <typename T> using cht_max = line_container<T, false>;
/// @brief Maintain points (a,b), querying min/max a*x+b*y and an attaining id.
/// Integer directions and intersections stay exact; returned values must fit T.
template <typename T> struct cht_xy {
static_assert(std::is_same_v<T, long long> || std::is_floating_point_v<T>);
using C = std::conditional_t<std::is_floating_point_v<T>, long double, T>;
using point_type = cht_internal::point<C>;
using W = cht_internal::wide<T>;
line_container<C, true, true> cm;
line_container<C, false, true> cma;
T bmx = std::numeric_limits<T>::lowest(), bmn = std::numeric_limits<T>::max();
int jx = -1, jy = -1;
bool emp = true;
void clear() {
emp = true;
cm.clear();
cma.clear();
bmx = std::numeric_limits<T>::lowest();
bmn = std::numeric_limits<T>::max();
jx = jy = -1;
}
void add(T a, T b, int i = -1) {
if (emp) {
bmx = bmn = b;
jx = jy = i;
}
emp = false;
cm.add(b, a, i);
cma.add(b, a, 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) const { return get<false>(x, y); }
std::pair<T, int> get_min(T x, T y) const { return get<true>(x, y); }
private:
template <bool Min> std::pair<T, int> get(T x, T y) const {
if (cm.empty())
return {Min ? std::numeric_limits<T>::max() : std::numeric_limits<T>::lowest(), -1};
if (x == 0) {
bool upper = Min ? y < 0 : y > 0;
return upper ? std::pair<T, int>{T(W(bmx) * y), jx}
: std::pair<T, int>{T(W(bmn) * y), jy};
}
point_type z;
if constexpr (std::numeric_limits<T>::is_integer) z = point_type(y, x);
else z = C(y) / x;
auto eval = [&](auto it) {
T a = it->m, b = it->k;
return std::pair<T, int>{T(W(a) * x + W(b) * y), it->id};
};
if ((x > 0) == Min) return eval(cm.lower_bound(z));
return eval(cma.lower_bound(z));
}
};
} // namespace noya
// NOYA_CONVEX_HULL_TRICK_HPP