convex_hull_trick.hpp

SECTIONOptimization INCLUDEnoya/convex_hull_trick.hpp

维护斜率单调加入的直线，并在查询横坐标也单调时求最优值。

Complexity: Time: O(log n) per insertion or monotone query. Space: O(n).

跳到代码 · GitHub ↗

Implementation

ReferenceSourceExpanded

当前头文件，省略 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.
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 ((__int128)(q1 - q2) * (p1 - p) < (__int128)(q - q1) * (p2 - p1)) {
        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;
    }
  }
};

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>::min(), amn = std::numeric_limits<T>::max();
  T bmx = std::numeric_limits<T>::min(), 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();
  }
  void add(T a, T b, int i = -1) {
    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>::min(), -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) {
    auto [f, i] = get_max(-x, -y);
    return {-f, i};
  }
};
} // namespace noya

#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.
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 ((__int128)(q1 - q2) * (p1 - p) < (__int128)(q - q1) * (p2 - p1)) {
        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;
    }
  }
};

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>::min(), amn = std::numeric_limits<T>::max();
  T bmx = std::numeric_limits<T>::min(), 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();
  }
  void add(T a, T b, int i = -1) {
    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>::min(), -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) {
    auto [f, i] = get_max(-x, -y);
    return {-f, i};
  }
};
} // namespace noya

#endif // NOYA_CONVEX_HULL_TRICK_HPP

#include <cassert>
#include <limits>
#include <map>
#include <set>
#include <type_traits>
#include <utility>
#include <vector>

/// @complexity Time: O(log n) per insertion or monotone query.
/// Space: O(n).

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.
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 ((__int128)(q1 - q2) * (p1 - p) < (__int128)(q - q1) * (p2 - p1)) {
        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;
    }
  }
};

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>::min(), amn = std::numeric_limits<T>::max();
  T bmx = std::numeric_limits<T>::min(), 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();
  }
  void add(T a, T b, int i = -1) {
    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>::min(), -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) {
    auto [f, i] = get_max(-x, -y);
    return {-f, i};
  }
};
} // namespace noya