:heavy_check_mark: verify/geometry/UNIT_halfplane_intersection.test.cpp

Depends on

Code

#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"

#include "template/template.hpp"
#include "geometry/halfplane-intersection.hpp"

Polygon halfplane_cut_naive(const Polygon& P, const Line& l) {
  Polygon ret;
  for (int i = 0; i < (int)P.size(); ++i) {
    Point a = P[i], b = P[(i + 1) % P.size()];
    int sa = sign(cross(l.b - l.a, a - l.a));
    int sb = sign(cross(l.b - l.a, b - l.a));
    if (sa >= 0) ret.push_back(a);
    if (sa * sb < 0) ret.push_back(cross_point_ll(Line{a, b}, l));
  }
  return ret;
}

int main() {
  {
    Lines ls{{{0, 0}, {2, 0}},
             {{2, 2}, {0, 2}},
             {{0, 2}, {0, 0}},
             {{2, 0}, {2, 2}},
             {{-1, -1}, {3, -1}}};
    Polygon P = halfplane_intersection(ls);
    assert(P.size() == 4);
    assert(abs(area(P) - 4) < EPS);
  }
  {
    Lines empty{{{1, 1}, {1, 0}}, {{0, 0}, {0, 1}}};
    Lines unbounded{{{0, 0}, {1, 0}}, {{0, 1}, {0, 0}}};
    assert(halfplane_intersection(empty).empty());
    assert(halfplane_intersection(unbounded).empty());
  }
  mt19937 rng(7654321);
  uniform_int_distribution<int> coord(-50, 50);
  for (int tc = 0; tc < 1000; ++tc) {
    vector<Point> ps(30);
    for (auto& p : ps) p = Point{(Real)coord(rng), (Real)coord(rng)};
    Polygon C = convex_hull(ps);
    if (C.size() < 3) continue;
    Lines ls;
    for (int i = 0; i < (int)C.size(); ++i)
      ls.emplace_back(C[i], C[(i + 1) % C.size()]);
    shuffle(ls.begin(), ls.end(), rng);
    Polygon P = halfplane_intersection(ls);
    assert(P.size() == C.size());
    assert(abs(area(P) - area(C)) < 1e-8);
  }
  uniform_int_distribution<int> small(-80, 80);
  for (int tc = 0; tc < 2000; ++tc) {
    Polygon want{{-100, -100}, {100, -100}, {100, 100}, {-100, 100}};
    Lines ls{{{-100, -100}, {100, -100}},
             {{100, -100}, {100, 100}},
             {{100, 100}, {-100, 100}},
             {{-100, 100}, {-100, -100}}};
    for (int i = 0; i < 12; ++i) {
      Point a{(Real)small(rng), (Real)small(rng)};
      Point b{(Real)small(rng), (Real)small(rng)};
      if (a == b) b.x += 1;
      ls.emplace_back(a, b);
      if (!want.empty()) want = halfplane_cut_naive(want, ls.back());
    }
    shuffle(ls.begin(), ls.end(), rng);
    Polygon got = halfplane_intersection(ls);
    Real want_area = want.empty() ? 0 : area(want);
    Real got_area = got.empty() ? 0 : area(got);
    assert(abs(want_area - got_area) < 1e-7);
  }
  int a, b;
  in(a, b);
  out(a + b);
}
#line 1 "verify/geometry/UNIT_halfplane_intersection.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"

#line 2 "template/template.hpp"
#include <bits/stdc++.h>
using namespace std;

#line 2 "template/macro.hpp"
#define rep(i, a, b) for (int i = (a); i < (int)(b); i++)
#define rrep(i, a, b) for (int i = (int)(b) - 1; i >= (a); i--)
#define ALL(v) (v).begin(), (v).end()
#define UNIQUE(v) sort(ALL(v)), (v).erase(unique(ALL(v)), (v).end())
#define SZ(v) (int)v.size()
#define MIN(v) *min_element(ALL(v))
#define MAX(v) *max_element(ALL(v))
#define LB(v, x) int(lower_bound(ALL(v), (x)) - (v).begin())
#define UB(v, x) int(upper_bound(ALL(v), (x)) - (v).begin())
#define YN(b) cout << ((b) ? "YES" : "NO") << "\n";
#define Yn(b) cout << ((b) ? "Yes" : "No") << "\n";
#define yn(b) cout << ((b) ? "yes" : "no") << "\n";
#line 6 "template/template.hpp"

#line 2 "template/util.hpp"
using uint = unsigned int;
using ll = long long int;
using ull = unsigned long long;
using i128 = __int128_t;
using u128 = __uint128_t;
template <class T>
using priority_queue_asc = priority_queue<T, vector<T>, greater<T>>;

template <class T, class S = T>
S SUM(const vector<T>& a) {
  return accumulate(ALL(a), S(0));
}
template <class T1, class T2>
inline bool chmin(T1& a, T2 b) {
  if (a > b) {
    a = b;
    return true;
  }
  return false;
}
template <class T1, class T2>
inline bool chmax(T1& a, T2 b) {
  if (a < b) {
    a = b;
    return true;
  }
  return false;
}
template <class T1, class T2>
inline bool chmin_opt(optional<T1>& a, T2 b) {
  if (!a || a > b) {
    a = b;
    return true;
  }
  return false;
}
template <class T1, class T2>
inline bool chmax_opt(optional<T1>& a, T2 b) {
  if (!a || a < b) {
    a = b;
    return true;
  }
  return false;
}

template <class T>
int popcnt(T x) {
  return __builtin_popcountll(x);
}
template <class T>
int topbit(T x) {
  return (x == 0 ? -1 : 63 - __builtin_clzll(x));
}
template <class T>
int lowbit(T x) {
  return (x == 0 ? -1 : __builtin_ctzll(x));
}
#line 8 "template/template.hpp"

#line 2 "template/inout.hpp"
struct Fast {
  Fast() {
    cin.tie(nullptr);
    ios_base::sync_with_stdio(false);
    cout << fixed << setprecision(15);
  }
} fast;

ostream& operator<<(ostream& os, __uint128_t x) {
  char buf[40];
  size_t k = 0;
  while (x > 0) buf[k++] = (char)(x % 10 + '0'), x /= 10;
  if (k == 0) buf[k++] = '0';
  while (k) os << buf[--k];
  return os;
}
ostream& operator<<(ostream& os, __int128_t x) {
  return x < 0 ? (os << '-' << (__uint128_t)(-x)) : (os << (__uint128_t)x);
}
template <class T, size_t N>
ostream& operator<<(ostream& os, const array<T, N>& a);
template <class T1, class T2>
istream& operator>>(istream& is, pair<T1, T2>& p) {
  return is >> p.first >> p.second;
}
template <class T1, class T2>
ostream& operator<<(ostream& os, const pair<T1, T2>& p) {
  return os << p.first << " " << p.second;
}
template <class T>
istream& operator>>(istream& is, vector<T>& a) {
  for (auto& v : a) is >> v;
  return is;
}
template <class T>
ostream& operator<<(ostream& os, const vector<T>& a) {
  for (auto it = a.begin(); it != a.end();) {
    os << *it;
    if (++it != a.end()) os << " ";
  }
  return os;
}
template <class T, size_t N>
ostream& operator<<(ostream& os, const array<T, N>& a) {
  for (auto it = a.begin(); it != a.end();) {
    os << *it;
    if (++it != a.end()) os << " ";
  }
  return os;
}
template <class T>
ostream& operator<<(ostream& os, const set<T>& st) {
  os << "{";
  for (auto it = st.begin(); it != st.end();) {
    os << *it;
    if (++it != st.end()) os << ",";
  }
  os << "}";
  return os;
}
template <class T1, class T2>
ostream& operator<<(ostream& os, const map<T1, T2>& mp) {
  os << "{";
  for (auto it = mp.begin(); it != mp.end();) {
    os << it->first << ":" << it->second;
    if (++it != mp.end()) os << ",";
  }
  os << "}";
  return os;
}

void in() {}
template <typename T, class... U>
void in(T& t, U&... u) {
  cin >> t;
  in(u...);
}
template <class... T>
void in_zip(int n, T&... t) {
  assert(n >= 0 && ((size(t) >= static_cast<size_t>(n)) && ...));
  for (int i = 0; i < n; i++) in(t[i]...);
}
void out() { cout << "\n"; }
template <typename T, class... U, char sep = ' '>
void out(const T& t, const U&... u) {
  cout << t;
  if (sizeof...(u)) cout << sep;
  out(u...);
}
template <class T, class U>
void out_opt(const optional<T>& opt, const U& fallback, ostream& os = cout) {
  if (opt.has_value())
    os << opt.value();
  else
    os << fallback;
  os << "\n";
}
template <class T, class U>
void out_opt(const vector<optional<T>>& vec, const U& fallback, ostream& os = cout) {
  for (auto it = vec.begin(); it != vec.end();) {
    if ((*it).has_value())
      os << (*it).value();
    else
      os << fallback;
    if (++it != vec.end()) os << " ";
  }
  os << "\n";
}

namespace IO {
template <class T, class... U>
T read(U&&... u) {
  T t = T(forward<U>(u)...);
  in(t);
  return t;
}
namespace Graph {
vector<vector<int>> unweighted(int n, int m, bool directed = false, int offset = 1) {
  vector<vector<int>> g(n);
  for (int i = 0; i < m; i++) {
    int u, v;
    cin >> u >> v;
    u -= offset, v -= offset;
    g[u].push_back(v);
    if (!directed) g[v].push_back(u);
  }
  return g;
}
template <class T>
vector<vector<pair<int, T>>> weighted(int n, int m, bool directed = false, int offset = 1) {
  vector<vector<pair<int, T>>> g(n);
  for (int i = 0; i < m; i++) {
    int u, v;
    T w;
    cin >> u >> v >> w;
    u -= offset, v -= offset;
    g[u].push_back({v, w});
    if (!directed) g[v].push_back({u, w});
  }
  return g;
}
}  // namespace Graph
namespace Tree {
vector<vector<int>> unweighted(int n, bool directed = false, int offset = 1) {
  return Graph::unweighted(n, n - 1, directed, offset);
}
template <class T>
vector<vector<pair<int, T>>> weighted(int n, bool directed = false, int offset = 1) {
  return Graph::weighted<T>(n, n - 1, directed, offset);
}
vector<vector<int>> rooted(int n, bool to_root = true, bool to_leaf = true, int offset = 1) {
  vector<vector<int>> g(n);
  for (int i = 1; i < n; i++) {
    int p;
    cin >> p;
    p -= offset;
    if (to_root) g[i].push_back(p);
    if (to_leaf) g[p].push_back(i);
  }
  return g;
}
}  // namespace Tree
}  // namespace IO
#line 10 "template/template.hpp"

#line 2 "template/debug.hpp"
#ifdef LOCAL
#define debug 1
#define show(...) _show(0, #__VA_ARGS__, __VA_ARGS__)
#else
#define debug 0
#define show(...) true
#endif
template <class T>
void _show(int, T) {
  cerr << '\n';
}
template <class T1, class T2, class... T3>
void _show(int i, const T1& a, const T2& b, const T3&... c) {
  for (; a[i] != ',' && a[i] != '\0'; i++) cerr << a[i];
  cerr << ":" << b << " ";
  _show(i + 1, a, c...);
}
#line 2 "geometry/halfplane-intersection.hpp"

#line 2 "geometry/line.hpp"

#line 2 "geometry/geometry-base.hpp"

#line 4 "geometry/geometry-base.hpp"

using Real = long double;
constexpr Real EPS = 1e-10;
constexpr Real PI = 3.141592653589793238462643383279L;
bool equals(Real x, Real y) { return fabs(x - y) < EPS; }
int sign(Real a) { return equals(a, 0) ? 0 : (a > 0 ? 1 : -1); }

template <class R>
struct PointBase {
  using P = PointBase;
  R x, y;
  PointBase() : x(0), y(0) {}
  PointBase(R _x, R _y) : x(_x), y(_y) {}
  template <typename T, typename U>
  PointBase(const pair<T, U>& p) : x(p.first), y(p.second) {}

  P operator+(const P& r) const { return P{x + r.x, y + r.y}; }
  P operator-(const P& r) const { return P{x - r.x, y - r.y}; }
  P operator-() const { return P{-x, -y}; }
  P operator*(R r) const { return P{x * r, y * r}; }
  P operator/(R r) const { return P{x / r, y / r}; }

  P& operator+=(const P& r) { return (*this) = (*this) + r; }
  P& operator-=(const P& r) { return (*this) = (*this) - r; }
  P& operator*=(R r) { return (*this) = (*this) * r; }
  P& operator/=(R r) { return (*this) = (*this) / r; }

  bool operator<(const P& r) const { return x != r.x ? x < r.x : y < r.y; }
  bool operator==(const P& r) const { return x == r.x and y == r.y; }
  bool operator!=(const P& r) const { return !((*this) == r); }

  P rotate(R rad) const {
    return {x * cos(rad) - y * sin(rad), x * sin(rad) + y * cos(rad)};
  }
  P rotate90() const { return {-y, x}; }

  R real() const { return x; }
  R imag() const { return y; }

  friend P operator*(R r, const P& p) { return p * r; }
  friend R real(const P& p) { return p.x; }
  friend R imag(const P& p) { return p.y; }
  friend R dot(const P& l, const P& r) { return l.x * r.x + l.y * r.y; }
  friend R cross(const P& l, const P& r) { return l.x * r.y - l.y * r.x; }
  friend R abs(const P& p) { return sqrt(p.x * p.x + p.y * p.y); }
  friend R norm(const P& p) { return p.x * p.x + p.y * p.y; }
  friend R arg(const P& p) { return atan2(p.y, p.x); }

  friend istream& operator>>(istream& is, P& p) {
    R a, b;
    is >> a >> b;
    p = P{a, b};
    return is;
  }
  friend ostream& operator<<(ostream& os, const P& p) {
    return os << p.x << " " << p.y;
  }
};
using Point = PointBase<Real>;
using Points = vector<Point>;

// relative position of c from a->b
int ccw(const Point& a, const Point& b, const Point& c) {
  Point x = b - a, y = c - a;
  if (cross(x, y) > EPS) return +1;        // counter-clockwise
  if (cross(x, y) < -EPS) return -1;       // clockwise
  if (dot(x, y) < -EPS) return +2;         // collinear in the order c-a-b
  if (norm(x) + EPS < norm(y)) return -2;  // collinear in the order a-b-c
  return 0;                                // collinear in the order a-c-b
}

/**
 * @brief 二次元幾何の基本要素
 * @docs docs/geometry/geometry-base.md
 */
#line 2 "geometry/polygon.hpp"

#line 4 "geometry/polygon.hpp"

// 反時計周り
using Polygon = vector<Point>;

// 多角形の内部に点があるか
// OUT : 0, ON : 1, IN : 2
// 凸でなくてもよい, 自己交差はもたない
// crossing number algorithm : x 軸正方向に伸ばした半直線が奇数回交われば内部
int contains_polygon(const Polygon& P, const Point& p) {
  bool in = false;
  for (int i = 0; i < (int)P.size(); i++) {
    Point a = P[i] - p, b = P[(i + 1) % P.size()] - p;
    if (a.y > b.y) swap(a, b);
    if (equals(cross(a, b), 0) && sign(dot(a, b)) <= 0) return 1;
    if (sign(a.y) <= 0 && 0 < sign(b.y) && sign(cross(a, b)) < 0)
      in = !in;
  }
  return in ? 2 : 0;
}

// 多角形の面積
Real area(const Polygon& P) {
  Real A = 0;
  for (int i = 0; i < (int)P.size(); i++) {
    A += cross(P[i], P[(i + 1) % P.size()]);
  }
  return A * 0.5;
}

bool is_convex(const Polygon& P) {
  int n = P.size();
  if (n < 3) return true;
  for (int i = 0; i < n; ++i) {
    if (ccw(P[i], P[(i + 1) % n], P[(i + 2) % n]) == -1) return false;
  }
  return true;
}

// 頂点集合から凸包を生成
// boundary : 周上の点も列挙する場合 true
template <bool boundary = false>
Polygon convex_hull(vector<Point> ps) {
  sort(begin(ps), end(ps));
  ps.erase(unique(begin(ps), end(ps)), end(ps));
  int n = ps.size(), k = 0;
  if (n <= 2) return ps;
  if constexpr (boundary) {
    bool collinear = true;
    for (int i = 2; i < n; ++i) {
      if (!equals(cross(ps[1] - ps[0], ps[i] - ps[0]), 0)) {
        collinear = false;
        break;
      }
    }
    if (collinear) return ps;
  }
  vector<Point> ch(2 * n);
  // 反時計周り
  const Real th = boundary ? -EPS : +EPS;
  for (int i = 0; i < n; ch[k++] = ps[i++]) {
    while (k >= 2 && cross(ch[k - 1] - ch[k - 2], ps[i] - ch[k - 1]) < th) --k;
  }
  for (int i = n - 2, t = k + 1; i >= 0; ch[k++] = ps[i--]) {
    while (k >= t && cross(ch[k - 1] - ch[k - 2], ps[i] - ch[k - 1]) < th) --k;
  }
  ch.resize(k - 1);
  return ch;
}

// 凸包の内部に点があるか
// OUT : 0, ON : 1, IN : 2
int contains_convex(const Polygon& C, const Point& p) {
  int N = C.size();
  if (N == 0) return 0;
  if (N == 1) return abs(C[0] - p) < EPS ? 1 : 0;
  if (N == 2) return ccw(C[0], C[1], p) == 0 ? 1 : 0;
  auto b1 = cross(C[1] - C[0], p - C[0]);
  auto b2 = cross(C[N - 1] - C[0], p - C[0]);
  if (b1 < -EPS or b2 > EPS) return 0;
  int L = 1, R = N - 1;
  while (L + 1 < R) {
    int M = (L + R) / 2;
    (cross(p - C[0], C[M] - C[0]) >= 0 ? R : L) = M;
  }
  auto v = cross(C[L] - p, C[R] - p);
  if (equals(v, 0)) {
    return 1;
  } else if (v > 0) {
    return equals(b1, 0) or equals(b2, 0) ? 1 : 2;
  } else {
    return 0;
  }
}

// 凸多角形の最遠点対を返す
// 返り値：頂点番号のペア
pair<int, int> convex_polygon_diameter(const Polygon& C) {
  int N = (int)C.size();
  assert(N > 0);
  if (N == 1) return {0, 0};
  if (N == 2) return {0, 1};
  int is = 0, js = 0;
  for (int i = 1; i < N; i++) {
    if (C[i].y > C[is].y) is = i;
    if (C[i].y < C[js].y) js = i;
  }
  Real maxdis = norm(C[is] - C[js]);

  int maxi, maxj, i, j;
  i = maxi = is;
  j = maxj = js;
  do {
    if (cross(C[(i + 1) % N] - C[i], C[(j + 1) % N] - C[j]) >= 0) {
      j = (j + 1) % N;
    } else {
      i = (i + 1) % N;
    }
    if (norm(C[i] - C[j]) > maxdis) {
      maxdis = norm(C[i] - C[j]);
      maxi = i;
      maxj = j;
    }
  } while (i != is || j != js);
  return minmax(maxi, maxj);
}

// min argmax_i dot(C[i], p)
int argmax_dot(const Polygon& C, const Point& p) {
  int n = C.size();
  assert(n > 0);
  assert(p.x != 0 || p.y != 0);
  if (n <= 2) {
    if (n == 1 || dot(C[0], p) >= dot(C[1], p)) return 0;
    return 1;
  }
  auto arg_half = [](const Point& a) {
    return a.y < 0 || (a.y == 0 && a.x < 0);
  };
  auto arg_less = [&](const Point& a, const Point& b) {
    int ha = arg_half(a), hb = arg_half(b);
    return ha != hb ? ha < hb : cross(a, b) > 0;
  };
  auto edge = [&](int i) { return C[(i + 1) % n] - C[i]; };
  int l = 0, r = n - 1;
  while (l < r) {
    int m = (l + r) / 2;
    if (arg_less(edge(r), edge(m)))
      l = m + 1;
    else
      r = m;
  }
  int s = l;
  Point q = p.rotate90();
  l = 0, r = n;
  while (l < r) {
    int m = (l + r) / 2;
    if (arg_less(edge((s + m) % n), q))
      l = m + 1;
    else
      r = m;
  }
  int i = (s + (l == n ? 0 : l)) % n;
  Point e = edge(i);
  if (!arg_less(e, q) && !arg_less(q, e)) {
    int low = l;
    l = low, r = n;
    while (l < r) {
      int m = (l + r) / 2;
      if (!arg_less(q, edge((s + m) % n)))
        l = m + 1;
      else
        r = m;
    }
    int count = l - low;
    return i + count >= n ? 0 : i;
  }
  return i;
}

struct ConvexPolygonCutResult {
  int first = 0, last = 0;
  optional<Point> front, back;

  bool empty() const { return first == last && !front && !back; }
  Polygon to_polygon(const Polygon& C) const {
    Polygon ret;
    auto push = [&](const Point& p) { if (ret.empty() || abs(ret.back() - p) >= EPS) ret.push_back(p); };
    if (front) push(*front);
    for (int i = first; i < last; ++i) push(C[i % C.size()]);
    if (back) push(*back);
    if (ret.size() >= 2 && abs(ret.front() - ret.back()) < EPS) ret.pop_back();
    return ret;
  }
};

ConvexPolygonCutResult convex_polygon_cut_info(const Polygon& C, const Point& a, const Point& b) {
  int n = C.size();
  assert(n > 0 && a != b);
  Point d = b - a, normal = d.rotate90();
  auto side = [&](int i) { return sign(cross(d, C[i % n] - a)); };
  int mx = argmax_dot(C, normal), mn = argmax_dot(C, -normal);
  if (side(mx) < 0) return {};
  if (side(mn) >= 0) return {0, n, nullopt, nullopt};

  int l = mn, r = mx;
  if (r < l) r += n;
  while (l + 1 < r) {
    int m = (l + r) / 2;
    (side(m) >= 0 ? r : l) = m;
  }
  int first = r;
  l = mx, r = mn;
  if (r <= l) r += n;
  while (l + 1 < r) {
    int m = (l + r) / 2;
    (side(m) < 0 ? r : l) = m;
  }
  int last = r;
  if (last < first) last += n;
  int count = last - first;
  first %= n;
  ConvexPolygonCutResult ret{first, first + count, nullopt, nullopt};
  auto cross_point = [&](const Point& p, const Point& q) {
    return p + (q - p) * (cross(d, a - p) / cross(d, q - p));
  };
  if (side(first) > 0) ret.front = cross_point(C[(first + n - 1) % n], C[first]);
  int back = (ret.last - 1) % n;
  if (side(back) > 0) ret.back = cross_point(C[back], C[(back + 1) % n]);
  return ret;
}

Polygon convex_polygon_cut(const Polygon& C, const Point& a, const Point& b) {
  return convex_polygon_cut_info(C, a, b).to_polygon(C);
}

/**
 * @brief 多角形
 * @docs docs/geometry/polygon.md
 */
#line 5 "geometry/line.hpp"

struct Line {
  using P = Point;
  P a, b;

  Line() = default;
  Line(const P& _a, const P& _b) : a(_a), b(_b) {}

  // Ax + By = C
  Line(Real A, Real B, Real C) {
    assert(!equals(A, 0) || !equals(B, 0));
    if (!equals(B, 0)) {
      a = P{0, C / B};
      b = P{1, (C - A) / B};
    } else {
      a = P{C / A, 0};
      b = P{C / A, 1};
    }
  }

  friend istream& operator>>(istream& is, Line& l) { return is >> l.a >> l.b; }
  friend ostream& operator<<(ostream& os, const Line& l) {
    return os << l.a << " to " << l.b;
  }
};

using Lines = vector<Line>;

bool is_parallel(const Line& l, const Line& m) {
  return equals(cross(l.b - l.a, m.b - m.a), 0);
}

bool is_orthogonal(const Line& l, const Line& m) {
  return equals(dot(l.b - l.a, m.b - m.a), 0);
}

bool is_intersect_ll(const Line& l, const Line& m) {
  return !is_parallel(l, m) || equals(cross(l.b - l.a, m.a - l.a), 0);
}

// Parallel distinct lines must not be given.
Point cross_point_ll(const Line& l, const Line& m) {
  Point u = l.b - l.a, v = m.b - m.a;
  Real d = cross(u, v);
  if (equals(d, 0)) {
    assert(equals(cross(u, m.a - l.a), 0));
    return m.a;
  }
  return l.a + u * (cross(m.a - l.a, v) / d);
}

Point projection(const Line& l, const Point& p) {
  Point d = l.b - l.a;
  return l.a + d * (dot(p - l.a, d) / norm(d));
}

Point reflection(const Line& l, const Point& p) {
  return projection(l, p) * 2 - p;
}

Real distance_lp(const Line& l, const Point& p) {
  return abs(p - projection(l, p));
}

// Returns the left side of the directed line l.a -> l.b.
Polygon convex_polygon_cut(const Polygon& P, const Line& l) {
  return convex_polygon_cut(P, l.a, l.b);
}

ConvexPolygonCutResult convex_polygon_cut_info(const Polygon& P,
                                               const Line& l) {
  return convex_polygon_cut_info(P, l.a, l.b);
}

/**
 * @brief 直線
 * @docs docs/geometry/line.md
 */
#line 4 "geometry/halfplane-intersection.hpp"

// Returns the bounded, positive-area intersection of left half-planes.
// Returns an empty polygon if the intersection is empty, degenerate, or unbounded.
Polygon halfplane_intersection(Lines ls) {
  auto direction = [](const Line& l) { return l.b - l.a; };
  auto arg_half = [](const Point& a) {
    return a.y < 0 || (a.y == 0 && a.x < 0);
  };
  auto arg_less = [&](const Line& l, const Line& r) {
    Point a = direction(l), b = direction(r);
    int ha = arg_half(a), hb = arg_half(b);
    return ha != hb ? ha < hb : cross(a, b) > 0;
  };
  for (const auto& l : ls) assert(l.a != l.b);
  sort(ls.begin(), ls.end(), arg_less);
  Lines lines;
  for (const auto& l : ls) {
    if (!lines.empty() && is_parallel(lines.back(), l) &&
        dot(direction(lines.back()), direction(l)) > 0) {
      if (sign(cross(direction(lines.back()), l.a - lines.back().a)) > 0)
        lines.back() = l;
    } else {
      lines.push_back(l);
    }
  }
  auto outside = [&](const Line& l, const Point& p) {
    return sign(cross(direction(l), p - l.a)) < 0;
  };
  deque<Line> dq;
  for (const auto& l : lines) {
    while (dq.size() >= 2 &&
           outside(l, cross_point_ll(dq[dq.size() - 2], dq.back())))
      dq.pop_back();
    while (dq.size() >= 2 && outside(l, cross_point_ll(dq[0], dq[1])))
      dq.pop_front();
    if (!dq.empty() && is_parallel(dq.back(), l)) return {};
    dq.push_back(l);
  }
  while (dq.size() >= 3 &&
         outside(dq.front(), cross_point_ll(dq[dq.size() - 2], dq.back())))
    dq.pop_back();
  while (dq.size() >= 3 &&
         outside(dq.back(), cross_point_ll(dq[0], dq[1])))
    dq.pop_front();
  if (dq.size() < 3) return {};
  int n = dq.size();
  for (int i = 0; i < n; ++i) {
    if (sign(cross(direction(dq[i]), direction(dq[(i + 1) % n]))) <= 0)
      return {};
  }
  Polygon ret;
  for (int i = 0; i < n; ++i)
    ret.push_back(cross_point_ll(dq[i], dq[(i + 1) % n]));
  return ret;
}

/**
 * @brief 半平面交差
 * @docs docs/geometry/halfplane-intersection.md
 */
#line 5 "verify/geometry/UNIT_halfplane_intersection.test.cpp"

Polygon halfplane_cut_naive(const Polygon& P, const Line& l) {
  Polygon ret;
  for (int i = 0; i < (int)P.size(); ++i) {
    Point a = P[i], b = P[(i + 1) % P.size()];
    int sa = sign(cross(l.b - l.a, a - l.a));
    int sb = sign(cross(l.b - l.a, b - l.a));
    if (sa >= 0) ret.push_back(a);
    if (sa * sb < 0) ret.push_back(cross_point_ll(Line{a, b}, l));
  }
  return ret;
}

int main() {
  {
    Lines ls{{{0, 0}, {2, 0}},
             {{2, 2}, {0, 2}},
             {{0, 2}, {0, 0}},
             {{2, 0}, {2, 2}},
             {{-1, -1}, {3, -1}}};
    Polygon P = halfplane_intersection(ls);
    assert(P.size() == 4);
    assert(abs(area(P) - 4) < EPS);
  }
  {
    Lines empty{{{1, 1}, {1, 0}}, {{0, 0}, {0, 1}}};
    Lines unbounded{{{0, 0}, {1, 0}}, {{0, 1}, {0, 0}}};
    assert(halfplane_intersection(empty).empty());
    assert(halfplane_intersection(unbounded).empty());
  }
  mt19937 rng(7654321);
  uniform_int_distribution<int> coord(-50, 50);
  for (int tc = 0; tc < 1000; ++tc) {
    vector<Point> ps(30);
    for (auto& p : ps) p = Point{(Real)coord(rng), (Real)coord(rng)};
    Polygon C = convex_hull(ps);
    if (C.size() < 3) continue;
    Lines ls;
    for (int i = 0; i < (int)C.size(); ++i)
      ls.emplace_back(C[i], C[(i + 1) % C.size()]);
    shuffle(ls.begin(), ls.end(), rng);
    Polygon P = halfplane_intersection(ls);
    assert(P.size() == C.size());
    assert(abs(area(P) - area(C)) < 1e-8);
  }
  uniform_int_distribution<int> small(-80, 80);
  for (int tc = 0; tc < 2000; ++tc) {
    Polygon want{{-100, -100}, {100, -100}, {100, 100}, {-100, 100}};
    Lines ls{{{-100, -100}, {100, -100}},
             {{100, -100}, {100, 100}},
             {{100, 100}, {-100, 100}},
             {{-100, 100}, {-100, -100}}};
    for (int i = 0; i < 12; ++i) {
      Point a{(Real)small(rng), (Real)small(rng)};
      Point b{(Real)small(rng), (Real)small(rng)};
      if (a == b) b.x += 1;
      ls.emplace_back(a, b);
      if (!want.empty()) want = halfplane_cut_naive(want, ls.back());
    }
    shuffle(ls.begin(), ls.end(), rng);
    Polygon got = halfplane_intersection(ls);
    Real want_area = want.empty() ? 0 : area(want);
    Real got_area = got.empty() ? 0 : area(got);
    assert(abs(want_area - got_area) < 1e-7);
  }
  int a, b;
  in(a, b);
  out(a + b);
}
Back to top page