geometry/HalfplaneSet.h

• View this file on GitHub
• Last update: 2025-11-21 16:12:02+07:00

Depends on

• geometry/Line.h
• geometry/LineIntersection.h
• geometry/Point.h

Code

#include "Line.h"
#include "LineIntersection.h"

struct HalfplaneSet : multiset<Line> {
  HalfplaneSet() {
    insert({+1, 0, INF});
    insert({0, +1, INF});
    insert({-1, 0, INF});
    insert({0, -1, INF});
  };

  auto adv(auto it, int z) {  // z = {-1, +1}
    return (z == -1 ? --(it == begin() ? end() : it)
                    : (++it == end() ? begin() : it));
  }
  bool chk(auto it) {
    Line l = *it, pl = *adv(it, -1), nl = *adv(it, +1);
    auto [x, y, d] = LineIntersection(pl, nl);
    T4 sat = l.a * x + l.b * y - (T4) l.c * d;
    if (d < 0 && sat < 0) return clear(), 0;  // unsat
    if ((d > 0 && sat <= 0) || (d == 0 && sat < 0)) return erase(it), 1;
    return 0;
  }
  void Cut(Line l) {  // add ax + by <= c
    if (empty()) return;
    auto it = insert(l);
    if (chk(it)) return;
    for (int z : {-1, +1})
      while (size() && chk(adv(it, z)));
  }
  ld Maximize(T a, T b) {  // max ax + by
    if (empty()) return -1 / 0.;
    auto it = lower_bound({a, b});
    if (it == end()) it = begin();
    auto [x, y, d] = LineIntersection(*adv(it, -1), *it);
    return (1.0 * a * x + 1.0 * b * y) / d;
  }
  ld Area() {  // half-plane intersection area
    ld total = 0.;
    for (auto it = begin(); it != end(); ++it) {
      auto [x1, y1, d1] = LineIntersection(*adv(it, -1), *it);
      auto [x2, y2, d2] = LineIntersection(*it, *adv(it, +1));
      total += (1.0 * x1 * y2 - 1.0 * x2 * y1) / d1 / d2;
    }
    return total * 0.5;
  }
};

#line 1 "geometry/Line.h"
using T = int; 
using T2 = long long;
using T4 = __int128_t;
const T2 INF = 4e18;

struct Line { T a, b; T2 c; };

bool half(Line m) { return m.a < 0 || m.a == 0 && m.b < 0; };
void normalize(Line& m) {
  T2 g =gcd((T2)gcd(abs(m.a), abs(m.b)), abs(m.c));
  if (half(m)) g *= -1;
  m.a /= g, m.b /= g, m.c /= g;
}
// Sorts halfplanes in clockwise order. 
// To sort lines, normalize first (gcd logic not needed).
bool operator<(Line m, Line n) {
  return make_pair(half(m), (T2)m.b * n.a) < 
         make_pair(half(n), (T2)m.a * n.b);
}
Line LineFromPoints(T x1, T y1, T x2, T y2) {
  T a = y1 - y2, b = x2 - x1;
  T2 c = (T2)a * x1 + (T2)b * y1;
  return {a, b, c}; // halfplane points to the left of vec.
}
#line 2 "geometry/Point.h"

template <class T>
int sgn(T x) { return (x > 0) - (x < 0); }
template <class T>
struct Point {
  typedef Point P;
  T x, y;
  explicit Point(T x = 0, T y = 0) : x(x), y(y) {}
  bool operator<(P p) const { return tie(x, y) < tie(p.x, p.y); }
  bool operator==(P p) const { return tie(x, y) == tie(p.x, p.y); }
  P operator+(P p) const { return P(x + p.x, y + p.y); }
  P operator-(P p) const { return P(x - p.x, y - p.y); }
  P operator*(T d) const { return P(x * d, y * d); }
  P operator/(T d) const { return P(x / d, y / d); }
  T dot(P p) const { return x * p.x + y * p.y; }
  T cross(P p) const { return x * p.y - y * p.x; }
  T cross(P a, P b) const { return (a - *this).cross(b - *this); }
  T dist2() const { return x * x + y * y; }
  T dist() const { return sqrt(dist2()); }
  // angle to x-axis in interval [-pi, pi]
  T angle() const { return atan2l(y, x); }
  P unit() const { return *this / dist(); }  // makes dist()=1
  P perp() const { return P(-y, x); }        // rotates +90 degrees
  P normal() const { return perp().unit(); }
  // returns point rotated 'a' radians ccw around the origin
  P rotate(ld a) const {
    return P(x * cos(a) - y * sin(a), x * sin(a) + y * cos(a));
  }
  friend ostream& operator<<(ostream& os, P p) {
    return os << "(" << p.x << "," << p.y << ")";
  }
};
#line 1 "geometry/Line.h"
using T = int; 
using T2 = long long;
using T4 = __int128_t;
const T2 INF = 4e18;

struct Line { T a, b; T2 c; };

bool half(Line m) { return m.a < 0 || m.a == 0 && m.b < 0; };
void normalize(Line& m) {
  T2 g =gcd((T2)gcd(abs(m.a), abs(m.b)), abs(m.c));
  if (half(m)) g *= -1;
  m.a /= g, m.b /= g, m.c /= g;
}
// Sorts halfplanes in clockwise order. 
// To sort lines, normalize first (gcd logic not needed).
bool operator<(Line m, Line n) {
  return make_pair(half(m), (T2)m.b * n.a) < 
         make_pair(half(n), (T2)m.a * n.b);
}
Line LineFromPoints(T x1, T y1, T x2, T y2) {
  T a = y1 - y2, b = x2 - x1;
  T2 c = (T2)a * x1 + (T2)b * y1;
  return {a, b, c}; // halfplane points to the left of vec.
}
#line 3 "geometry/LineIntersection.h"

template <class P>
pair<int, P> lineInter(P s1, P e1, P s2, P e2) {
  auto d = (e1 - s1).cross(e2 - s2);
  if (d == 0)  // if parallel
    return {-(s1.cross(e1, s2) == 0), P(0, 0)};
  auto p = s2.cross(e1, e2), q = s2.cross(e2, s1);
  return {1, (s1 * p + e1 * q) / d};
}

tuple<T4, T4, T2> LineIntersection(Line m, Line n) {
  T2 d = (T2)m.a * n.b - (T2)m.b * n.a; // assert(d);
  T4 x = (T4)m.c * n.b - (T4)m.b * n.c; 
  T4 y = (T4)m.a * n.c - (T4)m.c * n.a;
  return {x, y, d}; // (x/d, y/d) is intersection. 
}
#line 3 "geometry/HalfplaneSet.h"

struct HalfplaneSet : multiset<Line> {
  HalfplaneSet() {
    insert({+1, 0, INF});
    insert({0, +1, INF});
    insert({-1, 0, INF});
    insert({0, -1, INF});
  };

  auto adv(auto it, int z) {  // z = {-1, +1}
    return (z == -1 ? --(it == begin() ? end() : it)
                    : (++it == end() ? begin() : it));
  }
  bool chk(auto it) {
    Line l = *it, pl = *adv(it, -1), nl = *adv(it, +1);
    auto [x, y, d] = LineIntersection(pl, nl);
    T4 sat = l.a * x + l.b * y - (T4) l.c * d;
    if (d < 0 && sat < 0) return clear(), 0;  // unsat
    if ((d > 0 && sat <= 0) || (d == 0 && sat < 0)) return erase(it), 1;
    return 0;
  }
  void Cut(Line l) {  // add ax + by <= c
    if (empty()) return;
    auto it = insert(l);
    if (chk(it)) return;
    for (int z : {-1, +1})
      while (size() && chk(adv(it, z)));
  }
  ld Maximize(T a, T b) {  // max ax + by
    if (empty()) return -1 / 0.;
    auto it = lower_bound({a, b});
    if (it == end()) it = begin();
    auto [x, y, d] = LineIntersection(*adv(it, -1), *it);
    return (1.0 * a * x + 1.0 * b * y) / d;
  }
  ld Area() {  // half-plane intersection area
    ld total = 0.;
    for (auto it = begin(); it != end(); ++it) {
      auto [x1, y1, d1] = LineIntersection(*adv(it, -1), *it);
      auto [x2, y2, d2] = LineIntersection(*it, *adv(it, +1));
      total += (1.0 * x1 * y2 - 1.0 * x2 * y1) / d1 / d2;
    }
    return total * 0.5;
  }
};

Back to top page