text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
#include "Angle.h"
#include "Line.h"
using namespace std;

// Computes the intersection of half-planes defined by the left side of
//   a set of lines (including the line itself)
// Points on the intersection are given in ccw order and may be identical
// Angle::pivot is set to (0, 0)
// Function Arguments:
//   lines: a vector of lines representing the half-planes defined by the
//     left side
// Return Value: a pair containing a return code and a vector of points with
//     the return codes defined as follows:
//   code == 0: the intersection is finite (possibly empty),
//     with the vector containing the points on the intersection;
//   code == 1: the intersection is infinite and bounded by a series of
//     line segments and rays (possibly a single infinite line), with the
//     vector P[0], P[1], ..., P[k] where the first ray is represented by
//     the point P[1] in the direction towards P[0], and the second ray is
//     represented by the point P[k - 1] in the direction towards P[k];
//     the remaining points P[1], ..., P[k - 1] represent the intersections
//     between the line segments and rays; if there are only 2 points, then
//     the intersection is bounded by the left side of the line from
//     P[0] to P[1]
//   code == 2: the intersection is infinite and bounded by two lines
//     (possibly coincident lines in opposite directions), with the vector
//     containing exactly 4 points (P[0], P[1], P[2], P[3]), where P[0] is a
//     point on the first bounding line, P[1] is another point on the first
//     line in the positive direction (based on the direction vector of the
//     line) from P[0], P[2] is a point on the second line, and P[3] is another
//     point on the second line in its positive direction from P[2]
//   code == 3: the intersection is the entire plane, with an empty vector
// In practice, has a moderate constant
// Time Complexity: O(N log N)
// Memory Complexity: O(N)
// Tested:
//   https://open.kattis.com/problems/bigbrother
//   https://open.kattis.com/problems/marshlandrescues
//   https://dmoj.ca/problem/ccoprep3p3
pair<int, vector<pt>> halfPlaneIntersection(vector<Line> lines) {
  Angle::setPivot(pt(0, 0));
  sort(lines.begin(), lines.end(), [&] (Line a, Line b) {
    Angle angA(a.v), angB(b.v);
    return angA == angB ? a.onLeft(b.proj(pt(0, 0))) < 0 : angA < angB;
  });
  lines.erase(unique(lines.begin(), lines.end(), [&] (Line a, Line b) {
    return Angle(a.v) == Angle(b.v);
  }), lines.end());
  int N = lines.size(); if (N == 0) return make_pair(3, vector<pt>());
  if (N == 1) {
    pt p = lines[0].proj(pt(0, 0));
    return make_pair(1, vector<pt>{p, p + lines[0].v});
  }
  int code = 0; for (int i = 0; code == 0 && i < N; i++) {
    Angle diff = Angle(lines[i].v) - Angle(lines[i == 0 ? N - 1 : i - 1].v);
    if (diff < Angle(pt(1, 0))) {
      rotate(lines.begin(), lines.begin() + i, lines.end()); code = 1;
      if (N == 2 && diff == Angle(pt(-1, 0))) {
        pt p = lines[0].proj(pt(0, 0));
        if (lines[1].onLeft(p) < 0) return make_pair(0, vector<pt>());
        pt q = lines[1].proj(pt(0, 0));
        return make_pair(2, vector<pt>{p, p + lines[0].v, q, q + lines[1].v});
      }
    }
  }
  vector<Line> q(N + 1, lines[0]); vector<pt> ret(N); pt inter;
  int front = 0, back = 0; for (int i = 1; i <= N - code; i++) {
    if (i == N) lines.push_back(q[front]);
    while (front < back && lines[i].onLeft(ret[back - 1]) < 0) back--;
    while (i != N && front < back && lines[i].onLeft(ret[front]) < 0) front++;
    if (lineLineIntersection(lines[i], q[back], inter) != 1) continue;
    ret[back++] = inter; q[back] = lines[i];
  }
  if (code == 0 && back - front < 3) return make_pair(code, vector<pt>());
  vector<pt> P(ret.begin() + front, ret.begin() + back); if (code == 1) {
    P.insert(P.begin(), P[0] - q[front].v); P.push_back(P.back() + q[back].v);
  }
  return make_pair(code, P);
}


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Supports adding lines in the form f(x) = mx + b, or line segments
//   in the form of f(x) = mx + b over l <= x <= r, and finding
//   the maximum value of f(x) at an integral point x where MN <= x <= MX
// Template Arguments:
//   IndexType: the type of x for queries (and l and r for addLineSegment)
//   T: the type of the slope (m) and intercept (b) of the line, as well as
//     the type of the function argument (x)
//   Cmp: the comparator to compare two f(x) values,
//       convention is same as std::priority_queue in STL
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
// Constructor Arguments:
//   MN: the minimum bound (inclusive) for the value of x, allowing queries
//     where MN <= x
//   MX: the maximum bound (inclusive) for the value of x, allowing queries
//     where x <= MX
//   cmp: an instance of the Cmp struct
//   INF: a value for positive infinity, must be negatable
// Functions:
//   addLine(m, b): adds a line in the form f(x) = mx + b to the set of lines
//   addLineSegment(m, b, l, r): adds a line segment in the form f(x) = mx + b
//     where l <= x <= r, to the set of lines
//   getMax(x): finds the maximum value of f(x) (based on the comparator)
//     for all inserted lines
//   clear(): removes all lines from the seg
// In practice, has a moderate constant, performance compared to
//   IncrementalConvexHullTrick (which uses multiset) and
//   IncrementalConvexHullTrickSqrtBuffer can vary, slower than LiChaoTree
// Time Complexity:
//   constructor: O(1)
//   addLine, getMax: O(log(MX - MN)) for the range [MX, MN]
//   addLineSegment: O(log(MX - MN) ^ 2) for the range [MX, MN]
//   clear: O(Q log(MX - MN) + U (log(MX - MN))^2) for the range
//            [MX, MN], Q addLine queries, and U addLineSegment queries
// Memory Complexity: O(Q log(MX - MN) + U (log(MX - MN))^2) for the range
//                      [MX, MN], Q addLine queries, and U addLineSegment
//                      queries
// Tested:
//   https://judge.yosupo.jp/problem/line_add_get_min
//   https://judge.yosupo.jp/problem/segment_add_get_min
//   https://open.kattis.com/problems/longestlife
//   https://www.spoj.com/problems/CHTPRAC/
template <class IndexType, class T, class Cmp = less<T>>
struct SparseLiChaoTree {
  static_assert(is_integral<IndexType>::value, "IndexType must be integeral");
  using Line = pair<T, T>;
  struct Node {
    Line line; int l, r; Node(T m, T b) : line(m, b), l(-1), r(-1) {}
  };
  IndexType MN, MX; Cmp cmp; T INF; int root; vector<Node> TR;
  T eval(Line l, IndexType x) const { return l.first * x + l.second; }
  bool majorize(Line a, Line b, IndexType l, IndexType r) {
    return !cmp(eval(a, l), eval(b, l)) && !cmp(eval(a, r), eval(b, r));
  }
  SparseLiChaoTree(IndexType MN, IndexType MX, Cmp cmp = Cmp(),
                   T inf = numeric_limits<T>::max())
      : MN(MN), MX(MX), cmp(cmp), INF(min(inf, -inf, cmp)), root(-1) {}
  int addLine(int k, IndexType tl, IndexType tr, Line line) {
    if (k == -1) { k = TR.size(); TR.emplace_back(T(), INF); }
    if (majorize(line, TR[k].line, tl, tr)) swap(line, TR[k].line);
    if (majorize(TR[k].line, line, tl, tr)) return k;
    if (cmp(eval(TR[k].line, tl), eval(line, tl))) swap(line, TR[k].line);
    IndexType m = tl + (tr - tl) / 2;
    if (!cmp(eval(line, m), eval(TR[k].line, m))) {
      swap(line, TR[k].line);
      int nl = addLine(TR[k].l, tl, m, line); TR[k].l = nl;
    } else { int nr = addLine(TR[k].r, m + 1, tr, line); TR[k].r = nr; }
    return k;
  }
  int addLineSegment(int k, IndexType tl, IndexType tr,
                      IndexType l, IndexType r, Line line) {
    if (r < tl || tr < l) return k;
    if (l <= tl && tr <= r) return addLine(k, tl, tr, line);
    if (k == -1) { k = TR.size(); TR.emplace_back(T(), INF); }
    IndexType m = tl + (tr - tl) / 2;
    int nl = addLineSegment(TR[k].l, tl, m, l, r, line); TR[k].l = nl;
    int nr = addLineSegment(TR[k].r, m + 1, tr, l, r, line); TR[k].r = nr;
    return k;
  }
  T getMax(int k, IndexType tl, IndexType tr, IndexType x) const {
    if (k == -1) return INF;
    T ret = eval(TR[k].line, x); IndexType m = tl + (tr - tl) / 2;
    if (x <= m) return max(ret, getMax(TR[k].l, tl, m, x), cmp);
    else return max(ret, getMax(TR[k].r, m + 1, tr, x), cmp);
  }
  void addLine(T m, T b) { root = addLine(root, MN, MX, Line(m, b)); }
  void addLineSegment(T m, T b, IndexType l, IndexType r) {
    root = addLineSegment(root, MN, MX, l, r, Line(m, b));
  }
  T getMax(IndexType x) const { return getMax(root, MN, MX, x); }
  void clear() { root = -1; TR.clear(); }
};

repo_name: wesley-a-leung/Resources

mime_type: text/x-c++


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
#include "ConvexHull.h"
using namespace std;

// Computes the maximmum area from 3 distinct points out of a set of N points
// Constructor Arguments:
//   P: the vector of points
// Fields:
//   PA: the first point in the maximum area triangle
//   PB: the second point in the maximum area triangle
//   PC: the third point in the maximum area triangle
//   maxArea2: twice the area of the maximum area triangle
//   hull: the points in the convex hull of P
// In practice, has a small constant
// Time Complexity:
//   constructor: O(N^2)
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/dmpg18g5
struct MaxTriangleArea {
  pt PA, PB, PC; T maxArea2; vector<pt> hull;
  MaxTriangleArea(const vector<pt> &P) : maxArea2(0), hull(convexHull(P)) {
    int H = hull.size(), a = 0, b = 1, c = 2; if (H < 3) return;
    maxArea2 = area2(PA = hull[a], PB = hull[b], PC = hull[c]);
    for (;; b = (a + 1) % H, c = (b + 1) % H) {
      for (; c != a; b = (b + 1) % H) {
        T A = area2(hull[a], hull[b], hull[c]), B = A;
        while ((B = area2(hull[a], hull[b], hull[(c + 1) % H])) >= A) {
          c = (c + 1) % H; A = B;
        }
        if (maxArea2 < A) {
          maxArea2 = A; PA = hull[a]; PB = hull[b]; PC = hull[c];
        }
      }
      if ((a = (a + 1) % H) == 0) break;
    }
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
#include "Angle.h"
#include "Line.h"
#include "Circle.h"
using namespace std;

// Functions for 2D polygons, represented by a vector of the N points in
//   the polygon

// Returns i mod n if i is in the range [0, 2n)
int mod(int i, int n) { return i < n ? i : i - n; }

// Determines twice the signed area of a simple polygon
// Function Arguments:
//   poly: the points of the simple polygon
// Return Value: twice the signed area of the polygon, positive if
//   ccw, negative if cw
// Time Complexity: O(N)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/polygonarea
//   https://open.kattis.com/problems/crane
T getArea2(const vector<pt> &poly) {
  T ret = 0; int n = poly.size();
  for (int i = 0; i < n; i++) ret += cross(poly[i], poly[mod(i + 1, n)]);
  return ret;
}

// Determines centroid of a simple polygon
// Function Arguments:
//   poly: the points of the simple polygon
// Return Value: the centroid of the polygon
// Time Complexity: O(N)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/crane
pt getCentroid(const vector<pt> &poly) {
  T A2 = 0; pt cen(0, 0); int n = poly.size(); for (int i = 0; i < n; i++) {
    T a = cross(poly[i], poly[mod(i + 1, n)]); A2 += a;
    cen += a * (poly[i] + poly[mod(i + 1, n)]);
  }
  return cen / A2 / T(3);
}

// Determines the orientation of a convex polygon
// Function Arguments:
//   poly: the points of the convex polygon
// Return Value: 1 if ccw, -1 if cw, 0 if a point or a line
// Time Complexity: O(1)
// Memory Complexity: O(1)
int isCcwConvexPolygon(const vector<pt> &poly) {
  return ccw(poly.back(), poly[0], poly[mod(1, poly.size())]);
}

// Determines the orientation of a simple polygon
// Function Arguments:
//   poly: the points of the simple polygon
// Return Value: 1 if ccw, -1 if cw, 0 if a point or a line
// Time Complexity: O(N)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/abstractart
int isCcwPolygon(const vector<pt> &poly) {
  int n = poly.size();
  int i = min_element(poly.begin(), poly.end()) - poly.begin();
  return ccw(poly[mod(i + n - 1, n)], poly[i], poly[mod(i + 1, n)]);
}

// Determines whether a point is inside a convex polygon
// Function Arguments:
//   poly: the points of the convex polygon in ccw order or cw order
//   p: the point to check
// Return Value: 1 if inside the polygon, 0 if on the edge, -1 if outside
// Time Complexity: O(log N)
// Memory Complexity: O(1)
// Tested:
//   https://codeforces.com/contest/166/problem/B
int isInConvexPolygon(const vector<pt> &poly, pt p) {
  int n = poly.size(), a = 1, b = n - 1;
  if (n < 3) return onSeg(p, poly[0], poly.back()) ? 0 : -1;
  if (ccw(poly[0], poly[a], poly[b]) > 0) swap(a, b);
  int o1 = ccw(poly[0], p, poly[a]), o2 = ccw(poly[0], p, poly[b]);
  if (o1 < 0 || o2 > 0) return -1;
  if (o1 == 0 || o2 == 0) return 0;
  while (abs(a - b) > 1) {
    int c = (a + b) / 2; (ccw(poly[0], p, poly[c]) < 0 ? b : a) = c;
  }
  return ccw(poly[a], p, poly[b]);
}

// Determines whether a point is inside a simple polygon
// Function Arguments:
//   poly: the points of the simple polygon in ccw order or cw order
//   p: the point to check
// Return Value: 1 if inside the polygon, 0 if on the edge, -1 if outside
// Time Complexity: O(N)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/pointinpolygon
int isInPolygon(const vector<pt> &poly, pt p) {
  int n = poly.size(), windingNumber = 0; for (int i = 0; i < n; i++) {
    pt a = poly[i], b = poly[mod(i + 1, n)]; if (lt(b.y, a.y)) swap(a, b);
    if (onSeg(p, a, b)) return 0;
    windingNumber ^= (!lt(p.y, a.y) && lt(p.y, b.y) && ccw(p, a, b) > 0);
  }
  return windingNumber == 0 ? -1 : 1;
}

// Finds an extreme vertex of a convex polygon (a vertex that is the furthest
//   point in that direction, selecting the rightmost vertex if there are
//   multiple)
// Function Arguments:
//   poly: the points of the convex polygon in ccw order
//   dir: the direction
// Return Value: the index of an extreme vertex
// Time Complexity: O(log N)
// Memory Complexity: O(1)
// Tested:
//   https://codeforces.com/contest/799/problem/G
//   https://www.acmicpc.net/problem/4225
int extremeVertex(const vector<pt> &poly, pt dir) {
  int n = poly.size(), lo = 0, hi = n; pt pp = perp(dir);
  auto cmp = [&] (int i, int j) {
    return sgn(cross(pp, poly[mod(i, n)] - poly[mod(j, n)]));
  };
  auto extr = [&] (int i) {
    return cmp(i + 1, i) >= 0 && cmp(i, i + n - 1) < 0;
  };
  if (extr(0)) return 0;
  while (lo + 1 < hi) {
    int m = lo + (hi - lo) / 2; if (extr(m)) return m;
    int ls = cmp(lo + 1, lo), ms = cmp(m + 1, m);
    (ls < ms || (ls == ms && ls == cmp(lo, m)) ? hi : lo) = m;
  }
  return lo;
}

// Finds the intersection of a convex polygon and a line
// Function Arguments:
//   poly: the points of the convex polygon in ccw order
//   l: the line
// Return Value: (-1, -1) if no collision
//               (i, -1) if touching corner i
//               (i, i) if along side (i, i + 1)
//               (i, j) if crossing sides (i, i + 1) and (j, j + 1)
//                 crossing corner i is treated as crossing side (i, i + 1)
//   first index in pair is guaranteed to NOT be on the left side of the line
// Time Complexity: O(log N)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
//   https://codeforces.com/contest/799/problem/G
pair<int, int> convexPolygonLineIntersection(const vector<pt> &poly, Line l) {
  int n = poly.size();
  if (n == 1) return make_pair(l.onLeft(poly[0]) == 0 ? 0 : -1, -1);
  if (n == 2) {
    int o0 = l.onLeft(poly[0]), o1 = l.onLeft(poly[1]);
    if (o0 == 0 && o1 == 0) return make_pair(0, 0);
    if (o0 == 0) return make_pair(0, -1);
    if (o1 == 0) return make_pair(1, -1);
    return o0 == o1 ? make_pair(-1, -1) : make_pair(0, 1);
  }
  int endA = extremeVertex(poly, -perp(l.v));
  int endB = extremeVertex(poly, perp(l.v));
  auto cmpL = [&] (int i) { return l.onLeft(poly[i]); };
  pair<int, int> ret(-1, -1); if (cmpL(endA) > 0 || cmpL(endB) < 0) return ret;
  for (int i = 0; i < 2; i++) {
    int lo = endB, hi = endA; while (mod(lo + 1, n) != hi) {
      int m = mod((lo + hi + (lo < hi ? 0 : n)) / 2, n);
      (cmpL(m) == cmpL(endB) ? lo : hi) = m;
    }
    (i ? ret.second : ret.first) = mod(lo + !cmpL(hi), n); swap(endA, endB);
  }
  if (ret.first == ret.second) return make_pair(ret.first, -1);
  if (!cmpL(ret.first) && !cmpL(ret.second)) {
    switch ((ret.first - ret.second + n + 1) % n) {
      case 0: return make_pair(ret.first, ret.first);
      case 2: return make_pair(ret.second, ret.second);
    }
  }
  return ret;
}

// Finds a single tangent of a convex polygon and a point strictly outside
//   the polygon
// Function Arguments:
//   poly: the points of the convex polygon in ccw order
//   p: the point strictly outside the polygon
//   left: whether the left or right tangent line is found if p is
//     considered to be below the polygon
// Return Value: the tangent index, with the closest point to p being
//   selected if there are multiple indices
// Time Complexity: O(log N)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
//   https://dmoj.ca/problem/coci19c2p5
int convexPolygonPointSingleTangent(const vector<pt> &poly, pt p, bool left) {
  int n = poly.size(), o = ccw(p, poly[0], poly.back());
  bool farSide = o ? o < 0 : lt(distSq(p, poly.back()), distSq(p, poly[0]));
  int lo = farSide != left, hi = lo + n - 2; while (lo <= hi) {
    int mid = lo + (hi - lo) / 2;
    if (ccw(p, poly[0], poly[mid]) == (left ? -1 : 1)) {
      if (farSide == left) hi = mid - 1;
      else lo = mid + 1;
    } else {
      if ((ccw(poly[mid], poly[mod(mid + 1, n)], p) < 0) == left) hi = mid - 1;
      else lo = mid + 1;
    }
  }
  return mod(lo, n);
}

// Finds the tangent of a convex polygon and a point strictly outside
//   the polygon
// Function Arguments:
//   poly: the points of the convex polygon in ccw order
//   p: the point strictly outside the polygon
// Return Value: a pair containing the tangent indices, with the first index
//   being the left tangent point and the second index being the right tangent
//   if p is considered to be below the polygon; all points strictly between
//   the tangent indices are strictly within the tangent lines, while all other
//   points are on or outside the tangent lines
// Time Complexity: O(log N)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
//   https://dmoj.ca/problem/coci19c2p5
pair<int, int> convexPolygonPointTangent(const vector<pt> &poly, pt p) {
  return make_pair(convexPolygonPointSingleTangent(poly, p, true),
                   convexPolygonPointSingleTangent(poly, p, false));
}

// Finds the tangent of a convex polygon and a circle strictly and completely
//   outside the polygon
// Function Arguments:
//   poly: the points of the convex polygon in ccw order
//   c: the circle strictly and completely outside the polygon
//   inner: whether to find the inner or outer tangents
// Return Value: a pair containing the tangent indices, with the first index
//   being the left tangent point and the second index being the right tangent
//   if c is considered to be below the polygon; all points strictly between
//   the tangent indices are strictly within the tangent lines, while all other
//   points are on or outside the tangent lines (same index means all points
//   are inside or on the tangent lines)
// Time Complexity: O(log N)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
pair<int, int> convexPolygonCircleTangent(const vector<pt> &poly,
                                          Circle c, bool inner) {
  int n = poly.size(), a = 0, b = 0; vector<pair<pt, pt>> t;
  for (int h = 0; h < 2; h++) {
    assert(circleCircleTangent(Circle(poly[0], 0), c, inner, t) == 1);
    pt q = t[h].second; int o = ccw(q, poly[0], poly.back());
    bool farSide = o ? o < 0 : lt(distSq(q, poly.back()), distSq(q, poly[0]));
    int lo = farSide == h, hi = lo + n - 2; while (lo <= hi) {
      int mid = lo + (hi - lo) / 2; t.clear();
      assert(circleCircleTangent(Circle(poly[mid], 0), c, inner, t) == 1);
      q = t[h].second; if (ccw(q, poly[0], poly[mid]) == (h ? 1 : -1)) {
        if (farSide != h) hi = mid - 1;
        else lo = mid + 1;
      } else {
        if ((ccw(poly[mid], poly[mod(mid + 1, n)], q) < 0) != h) hi = mid - 1;
        else lo = mid + 1;
      }
    }
    (h ? b : a) = mod(lo, n); t.clear();
  }
  return make_pair(a, b);
}

// Finds the tangents of two convex polygons that do not intersect
// Function Arguments:
//   poly1: the first convex polygon
//   poly2: the second convex polygon that does not intersect with the first
//   inner: whether to find the inner or outer tangents
// Return Value: a vector of pairs containing the tangent indices, with the
//   first element in each pair being the index in the first polygon, and the
//   second element being the index in the second polygon; first point in
//   each pair is the left tangent point of the first polygon if the second
//   polygon is considered to be below the first polygon
// Time Complexity: O(log N log M)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
vector<pair<int, int>> convexPolygonConvexPolygonTangent(
    const vector<pt> &poly1, const vector<pt> &poly2, bool inner) {
  int n = poly1.size(), a = 0, b = 0, c = 0, d = 0; vector<pair<int, int>> ret;
  for (int h = 0; h < 2; h++) {
    pt q = poly2[convexPolygonPointSingleTangent(poly2, poly1[0], inner ^ h)];
    int o = ccw(q, poly1[0], poly1.back());
    bool farSide = o ? o < 0
                     : lt(distSq(q, poly1.back()), distSq(q, poly1[0]));
    int lo = farSide == h, hi = lo + n - 2; while (lo <= hi) {
      int mid = lo + (hi - lo) / 2;
      q = poly2[convexPolygonPointSingleTangent(poly2, poly1[mid], inner ^ h)];
      if (ccw(q, poly1[0], poly1[mid]) == (h ? 1 : -1)) {
        if (farSide != h) hi = mid - 1;
        else lo = mid + 1;
      } else {
        if ((ccw(poly1[mid], poly1[mod(mid + 1, n)], q) < 0) != h)
          hi = mid - 1;
        else lo = mid + 1;
      }
    }
    (h ? b : a) = lo = mod(lo, n);
    (h ? d : c) = convexPolygonPointSingleTangent(poly2, poly1[lo], inner ^ h);
  }
  ret.emplace_back(a, c); ret.emplace_back(b, d); return ret;
}

// Finds the closest point on the edge of the polygon to a point strictly
//   outside the polygon
// Function Arguments:
//   poly: the points of the convex polygon in ccw order
//   p: the point strictly outside the polygon
// Return Value: the closest point on the edge of the polygon to the point p
// Time Complexity: O(log N)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
pt closestPointOnConvexPolygon(const vector<pt> &poly, pt p) {
  pair<int, int> tangent = convexPolygonPointTangent(poly, p);
  int n = poly.size(), len = tangent.second - tangent.first;
  if (len < 0) len += n;
  if (len == 0) return poly[tangent.first];
  int lo = 0, hi = len - 2; while (lo <= hi) {
    int mid = lo + (hi - lo) / 2, i = mod(tangent.first + mid, n);
    if (ptSegDist(p, poly[i], poly[mod(i + 1, n)])
        < ptSegDist(p, poly[mod(i + 1, n)], poly[mod(i + 2, n)]))
      hi = mid - 1;
    else lo = mid + 1;
  }
  int i = mod(tangent.first + lo, n);
  return closestPtOnSeg(p, poly[i], poly[mod(i + 1, n)]);
}

// Determines the intersection of a simple polygon and a half-plane defined
//   by the left side of a line (including the line itself)
// Function Arguments:
//   poly: the points of the simple polygon in ccw order
//   l: the line with the half-plane defined by the left side
// Return Value: the polygon defined by the intersection of the simple polygon
//   and the half-plane, assuming the result is a single simple polyon
// Time Complexity: O(N)
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/utso15p6
//   https://open.kattis.com/problems/canyon
vector<pt> polygonHalfPlaneIntersection(const vector<pt> &poly, Line l) {
  int n = poly.size(); vector<pt> ret; for (int i = 0; i < n; i++) {
    int j = mod(i + 1, n), o1 = l.onLeft(poly[i]), o2 = l.onLeft(poly[j]);
    if (o1 >= 0) ret.push_back(poly[i]);
    if (o1 && o2 && o1 != o2) {
      pt p; if (lineLineIntersection(l, Line(poly[i], poly[j]), p) == 1)
        ret.push_back(p);
    }
  }
  return ret;
}

// Computes the area of union of multiple polygons
// Function Arguments:
//   polys: a vector of the polygons represented by a vector of points given in
//     ccw order
// Return Value: the area of the union of all the polygons
// Time Complexity: O(N^2 log N) for N total points
// Memory Complexity: O(N) for N total points
// Tested:
//   https://open.kattis.com/problems/abstractart
T polygonUnion(const vector<vector<pt>> &polys) {
  auto rat = [&] (pt p, pt q) { return sgn(q.x) ? p.x / q.x : p.y / q.y; };
  T ret = 0; for (int i = 0; i < int(polys.size()); i++)
    for (int v = 0; v < int(polys[i].size()); v++) {
      pt a = polys[i][v], b = polys[i][mod(v + 1, polys[i].size())];
      vector<pair<T, int>> segs{make_pair(0, 0), make_pair(1, 0)};
      for (int j = 0; j < int(polys.size()); j++) if (i != j)
        for (int w = 0; w < int(polys[j].size()); w++) {
          pt c = polys[j][w], d = polys[j][mod(w + 1, polys[j].size())];
          int sc = ccw(a, b, c), sd = ccw(a, b, d); if (sc != sd) {
            if (min(sc, sd) < 0) {
              T sa = area2(c, d, a), sb = area2(c, d, b);
              segs.emplace_back(sa / (sa - sb), sgn(sc - sd));
            }
          } else if (j < i && !sc && !sd && sgn(dot(b - a, d - c)) > 0) {
            segs.emplace_back(rat(c - a, b - a), 1);
            segs.emplace_back(rat(d - a, b - a), -1);
          }
        }
      sort(segs.begin(), segs.end()); T sm = 0;
      for (auto &&s : segs) s.first = min(max(s.first, T(0)), T(1));
      for (int j = 1, cnt = segs[0].second; j < int(segs.size()); j++) {
        if (!cnt) sm += segs[j].first - segs[j - 1].first;
        cnt += segs[j].second;
      }
      ret += cross(a, b) * sm;
    }
  return ret / 2;
}

// Determines the area of the intersection of a simple polygon and a circle
// Function Arguments:
//   poly: the points of the simple polygon in ccw order
//   c: the circle
// Return Value: the area of the intersection of the simple polygon and
//   the circle
// Time Complexity: O(N)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/pizzacutting
//   https://open.kattis.com/problems/birthdaycake
T polygonCircleIntersectionArea(const vector<pt> &poly, Circle c) {
  T r2 = c.r * c.r / 2;
  auto tri = [&] (pt p, pt q) {
    pt d = q - p; T a = dot(d, p) / norm(d);
    T b = (norm(p) - c.r * c.r) / norm(d), det = a * a - b;
    if (!lt(0, det)) return ang(q, pt(0, 0), p) * r2;
    T s = max(T(0), -a - sqrt(det)), t = min(T(1), -a + sqrt(det));
    if (lt(t, 0) || !lt(s, 1)) return ang(q, pt(0, 0), p) * r2;
    pt u = p + d * s, v = p + d * t;
    return ang(u, pt(0, 0), p) * r2 + cross(u, v) / 2
        + ang(q, pt(0, 0), v) * r2;
  };
  T ret = 0; for (int n = poly.size(), i = 0; i < n; i++)
    ret += tri(poly[i] - c.o, poly[mod(i + 1, n)] - c.o);
  return ret;
}

// Computes the area of union of multiple polygons and multiple circles
// Assertion failure likely means there is precision error
// Function Arguments:
//   polys: a vector of the polygons represented by a vector of points given in
//     ccw order
//   circles: a vector of the circles
// Return Value: the area of the union of all the polygons and all the circles
// Time Complexity: O((N + M)^2 log (N + M)) for N total points and M circles
// Memory Complexity: O((N + M)^2) for N total points and M circles
// Tested:
//   https://open.kattis.com/problems/abstractart
//   https://www.spoj.com/problems/CIRU/
//   https://dmoj.ca/problem/noi05p6
T polygonCircleUnionArea(const vector<vector<pt>> &polys,
                         const vector<Circle> &circles) {
  int n = polys.size(), m = circles.size(); T ret = 0;
  auto rat = [&] (pt p, pt q) { return sgn(q.x) ? p.x / q.x : p.y / q.y; };
  for (int i = 0; i < n; i++) for (int v = 0; v < int(polys[i].size()); v++) {
    pt a = polys[i][v], b = polys[i][mod(v + 1, polys[i].size())];
    if (a == b) continue;
    vector<pair<T, int>> segs{make_pair(0, 0), make_pair(1, 0)};
    for (int j = 0; j < n; j++) if (i != j)
      for (int w = 0; w < int(polys[j].size()); w++) {
        pt c = polys[j][w], d = polys[j][mod(w + 1, polys[j].size())];
        int sc = ccw(a, b, c), sd = ccw(a, b, d); if (sc != sd) {
          if (min(sc, sd) < 0) {
            T sa = area2(c, d, a), sb = area2(c, d, b);
            segs.emplace_back(sa / (sa - sb), sgn(sc - sd));
          }
        } else if (j < i && !sc && !sd && sgn(dot(b - a, d - c)) > 0) {
          segs.emplace_back(rat(c - a, b - a), 1);
          segs.emplace_back(rat(d - a, b - a), -1);
        }
      }
    Line l(a, b); for (int j = 0; j < m; j++) {
      vector<pt> p = circleLineIntersection(circles[j], l);
      if (int(p.size()) == 2) {
        segs.emplace_back(rat(p[0] - a, b - a), 1);
        segs.emplace_back(rat(p[1] - a, b - a), -1);
      }
    }
    sort(segs.begin(), segs.end()); T sm = 0;
    for (auto &&s : segs) s.first = min(max(s.first, T(0)), T(1));
    for (int j = 1, cnt = segs[0].second; j < int(segs.size()); j++) {
      if (!cnt) sm += segs[j].first - segs[j - 1].first;
      cnt += segs[j].second;
    }
    ret += cross(a, b) * sm / 2;
  }
  for (int i = 0; i < m; i++) {
    vector<pair<Angle, int>> segs; Angle::setPivot(circles[i].o);
    segs.emplace_back(Angle(circles[i].o - pt(circles[i].r, 0)), 0);
    segs.emplace_back(Angle(circles[i].o + pt(circles[i].r, 0)), 0);
    segs.emplace_back(Angle(circles[i].o), 0); bool covered = false;
    for (int j = 0; j < m; j++) if (i != j) {
      int o = circles[j].contains(circles[i]);
      if (o > 0 || (o == 0 && (lt(circles[i].r, circles[j].r) || j < i))) {
        covered = true; break;
      }
      vector<pt> p; circleCircleIntersection(circles[i], circles[j], p);
      if (int(p.size()) == 2) {
        Angle a(p[0]), b(p[1]);
        segs.emplace_back(a, 1); segs.emplace_back(b, -1); if (a >= b) {
          segs.emplace_back(Angle(circles[i].o - pt(circles[i].r, 0)), 1);
          segs.emplace_back(Angle(circles[i].o), -1);
        }
      }
    }
    for (int j = 0; j < n && !covered; j++) {
      vector<pair<Angle, int>> tmp; bool hasInter = false;
      for (int w = 0; w < int(polys[j].size()); w++) {
        pt c = polys[j][w], d = polys[j][mod(w + 1, polys[j].size())];
        if (c == d) continue;
        int sc = circles[i].contains(c), sd = circles[i].contains(d);
        vector<pt> p = circleSegIntersection(circles[i], c, d); if (sc != sd) {
          if (min(sc, sd) > 0 && int(p.size()) == 1) {
            hasInter = true; tmp.emplace_back(Angle(p[0]), sgn(sc - sd));
          } else if (min(sc, sd) < 0 && int(p.size()) == 2) {
            hasInter = true; tmp.emplace_back(Angle(p[sc == 0]), sgn(sc - sd));
          }
        } else if (int(p.size()) == 2) {
          hasInter = true; tmp.emplace_back(Angle(p[1]), 1);
          tmp.emplace_back(Angle(p[0]), -1);
        }
      }
      auto cmp = [&] (const pair<Angle, int> &a, const pair<Angle, int> &b) {
        return make_pair(a.first, -a.second) < make_pair(b.first, -b.second);
      };
      sort(tmp.begin(), tmp.end(), cmp);
      if (!tmp.empty() && tmp[0].second < 0) {
        tmp.emplace_back(Angle(circles[i].o - pt(circles[i].r, 0)), 1);
        tmp.emplace_back(Angle(circles[i].o), -1);
      }
      segs.insert(segs.end(), tmp.begin(), tmp.end());
      if (!hasInter && isInPolygon(polys[j], circles[i].o) >= 0)
        for (auto &&p : polys[j])
          if (circles[i].contains(p) < 0) { covered = true; break; }
    }
    if (covered) continue;
    sort(segs.begin(), segs.end()); for (auto &&s : segs)
      if (s.first.p == circles[i].o) s.first.p -= pt(circles[i].r, 0);
    for (int j = 1, cnt = segs[0].second; j < int(segs.size()); j++) {
      if (!cnt) {
        pt a = segs[j - 1].first.p, b = segs[j].first.p; if (a != b) {
          ret += cross(a, b) / 2;
          ret += circleHalfPlaneIntersectionArea(circles[i], Line(b, a));
        }
      }
      cnt += segs[j].second;
    }
  }
  return ret;
}


text:
#pragma once
#include <bits/stdc++.h>
#include "../../search/BinarySearch.h"
using namespace std;

// Supports adding lines in the form f(x) = mx + b and finding
//   the maximum value of f(x) at any given x; this version allows for
//   updates and queries in arbitrary order
// Template Arguments:
//   T: the type of the slope (m) and intercept (b) of the line, as well as
//     the type of the function argument (x), must be able to store m * b
//   Cmp: the comparator to compare two f(x) values,
//       convention is same as std::priority_queue in STL
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
// Constructor Arguments:
//   cmp: an instance of the Cmp struct
//   SCALE: the value to scale sqrt by
// Functions:
//   addLine(m, b): adds a line in the form f(x) = mx + b to the set of lines
//   getMax(x): finds the maximum value of f(x) (based on the comparator)
//     for all inserted lines
//   size(): returns the number of lines in the convex hull
//   reserve(N): reserves space for N lines in the convex hull
// In practice, has a very small constant, performance compared to
//   IncrementalConvexHullTrick (which uses multiset) and SparseLiChao
//   can vary, slower than LiChao
// Time Complexity:
//   constructor: O(1)
//   addLine: O(1) amortized
//   getMax: O(sqrt(N) + log(N)) amortized for N lines in the convex hull
//   size: O(1)
//   reserve: O(N)
// Memory Complexity: O(N) for N lines in the convex hull
// Tested:
//   https://judge.yosupo.jp/problem/line_add_get_min
//   https://open.kattis.com/problems/longestlife
//   https://www.spoj.com/problems/CHTPRAC/
template <class T, class Cmp = less<T>>
struct IncrementalConvexHullTrickSqrtBuffer {
  struct Line {
    T m, b; Line(T m, T b) : m(m), b(b) {}
    T eval(T x) const { return m * x + b; }
  };
  vector<Line> large, small; Cmp cmp; double SCALE;
  IncrementalConvexHullTrickSqrtBuffer(Cmp cmp = Cmp(), double SCALE = 4)
      : cmp(cmp), SCALE(SCALE) {}
  bool ccw(Line a, Line b, Line c) {
    return (b.m - a.m) * (c.b - a.b) <= (b.b - a.b) * (c.m - a.m);
  }
  bool slope(Line a, Line b) {
    return !cmp(a.m, b.m) && !cmp(b.m, a.m) && !cmp(a.b, b.b);
  }
  void rebuildHull() {
    int back = 0; for (auto &&line : large) {
      while (back >= 1 && slope(line, large[back - 1])) back--;
      while (back >= 2 && ccw(line, large[back - 1], large[back - 2])) back--;
      large[back++] = line;
    }
    large.erase(large.begin() + back, large.end());
  }
  int size() const { return large.size() + small.size(); }
  void rebuild() {
    if (int(small.size()) > SCALE * sqrt(size())) {
      auto lcmp = [&] (Line a, Line b) { return cmp(a.m, b.m); };
      int lSz = large.size(); sort(small.begin(), small.end(), lcmp);
      large.insert(large.end(), small.begin(), small.end()); small.clear();
      inplace_merge(large.begin(), large.begin() + lSz, large.end(), lcmp);
      rebuildHull();
    }
  }
  void addLine(T m, T b) { small.emplace_back(m, b); }
  T getMax(T x) {
    rebuild(); int ind = bsearch<FIRST>(0, int(large.size()) - 1, [&] (int i) {
      return cmp(large[i + 1].eval(x), large[i].eval(x));
    });
    T mx = (large.empty() ? small[0] : large[ind]).eval(x);
    for (auto &&line : small) mx = max(mx, line.eval(x), cmp);
    return mx;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
#include "Angle.h"
#include "Line.h"
#include "Circle.h"
using namespace std;

// A ray represented by a point and line with a direction
// Fields:
//   p: the point
//   l: a line passing through p, representing the direction
//   nxt: a pointer to the next ray in the circular set
struct Ray {
  pt p; Line l; mutable const Ray *nxt;
  Ray(pt p, Line l = Line()) : p(p), l(l), nxt(nullptr) {}
  virtual bool cmp(const Ray &o) const {
    Angle::setPivot(pt(0, 0)); return Angle(o.l.v) < Angle(l.v);
  }
  bool operator < (const Ray &r) const { return r.cmp(*this); }
};

// Helper struct for isIn
struct IsInCmp : public Ray {
  pt q, r; IsInCmp(pt p, pt q, pt r) : Ray(p), q(q), r(r) {}
  bool cmp(const Ray &o) const override {
    return q != o.p && ccw(p, o.nxt->p, r) > 0;
  }
};

// Helper struct for pointTangents
struct PointTangentCmp : public Ray {
  pt q, r; bool left, farSide;
  PointTangentCmp(pt p, pt q, pt r, bool left, bool farSide)
      : Ray(p), q(q), r(r), left(left), farSide(farSide) {}
  bool cmp(const Ray &o) const override {
    if (farSide != left && o.p == p) return true;
    if (farSide == left && o.p == q) return false;
    if (ccw(r, p, o.p) == (left ? -1 : 1)) return farSide != left;
    else return (ccw(o.p, o.nxt->p, r) < 0) != left;
  }
};

// Helper struct for circleTangents
struct CircleTangentCmp : public Ray {
  pt q; Circle c; bool inner, h, farSide;
  CircleTangentCmp(pt p, pt q, Circle c, bool inner, bool h, bool farSide)
      : Ray(p), q(q), c(c), inner(inner), h(h), farSide(farSide) {}
  bool cmp(const Ray &o) const override {
    if (farSide == h && o.p == p) return true;
    if (farSide != h && o.p == q) return false;
    vector<pair<pt, pt>> t;
    assert(circleCircleTangent(Circle(o.p, 0), c, inner, t) == 1);
    pt q = t[h].second;
    if (ccw(q, p, o.p) == (h ? 1 : -1)) return farSide == h;
    else return (ccw(o.p, o.nxt->p, q) < 0) == h;
  }
};

// Helper struct for closestPt
struct ClosestPointCmp : public Ray {
  Angle lo, hi;
  ClosestPointCmp(pt p, Angle lo, Angle hi) : Ray(p), lo(lo), hi(hi) {}
  bool cmp(const Ray &o) const override {
    Angle::setPivot(pt(0, 0)); if (Angle(o.l.v) < lo) return true;
    if (hi < Angle(o.l.v)) return false;
    return ptSegDist(p, o.p, o.nxt->p)
        >= ptSegDist(p, o.nxt->p, o.nxt->nxt->p);
  }
};

// Helper struct for lineIntersection
struct LineIntersectionCmp : public Ray {
  Angle lo, hi;
  LineIntersectionCmp(Line l, Angle lo, Angle hi)
      : Ray(pt(0, 0), l), lo(lo), hi(hi) {}
  bool cmp(const Ray &o) const override {
    Angle::setPivot(pt(0, 0)); if (Angle(o.l.v) < lo) return true;
    if (hi < Angle(o.l.v)) return false;
    return l.onLeft(o.nxt->p) < 0;
  }
};

// Maintains the convex hull of points where points can be added at any time
// Points are stored as a circular ordered set of Rays, with an associated
//   Line in the direction towards the next point
// Points are strictly convex
// Fields:
//   a2: twice the area of the convex hull
// Functions:
//   clear(): clears the points in the convex hull
//   isIn(p): returns 1 if inside the polygon, 0 if on the edge, -1 if outside
//   singlePointTangent(p, left): returns an iterator pointing to the left or
//     right tangent, with the closest point to p being selected if there are
//     multiple points, p must be strictly outside the polygon
//   pointTangents(p): returns a pair of iterators, with the first iterator
//     being the left tangent point and the second iterator being the right
//     tangent if p is considered to be below the polygon; all points strictly
//     between the tangent iterator are strictly within the tangent lines,
//     while all other points are on or outside the tangent lines,
//     p must be strictly outside the polygon
//   circleTangents(c, inner): returns a pair of iterators, with the first
//     iterator being the left inner or outer tangent point and the second
//     iterator being the right tangent if c is considered to be below the
//     polygon; all points strictly between the tangent iterator are strictly
//     within the tangent lines, while all other points are on or outside the
//     tangent lines (same iterator means all points are inside or on the
//     tangent lines), the circle must be strictly and completely outside the
//     polygon
//   hullTangents(hull, inner): returns a vector of pair of iterators, with the
//     first element in each pair being the iterator in the first polygon,
//     and the second element being the iteraetor in the second polygon; first
//     point in each pair is the left inner/outer tangent point of the first
//     polygon if the second polygon is considered to be below the
//     first polygon
//   closestPt(p): returns the closest point on the edge of the polygon to
//     the point p where p is strictly outside the polygon
//   addPoint(p): adds the point p and removes any points that are no longer
//     in the convex hull; returns true if p is in the hull, false otherwise;
//     Angle::pivot is set to (0, 0)
//   extremeVertex(dir): returns an iterator pointing to an extreme vertex
//     in the direction dir (a vertex that is the furthest point in that
//     direction, selecting the rightmost vertex if there are multiple);
//     Angle::pivot is set to (0, 0)
//   lineIntersection(l): if l does not intersect with this convex hull, an
//     empty vector is returned; if there is one point of intersection, a
//     vector containing the point of intersection is returned; if there is a
//     segment of intersection, a vector containing two endpoints of the
//     line segment intersection is returned
//   halfPlaneIntersection(l): intersects this convex hull with the half-plane
//     specified by the left side of l (including l itself); resulting hull
//     is strictly convex; Angle::pivot is set to (0, 0)
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(1)
//   clear: O(N)
//   isIn, singlePointTangent, pointTangents, circleTangents,
//     closestPt, extremeVertex, lineIntersection: O(log N)
//   hullTangents: O(log N log M)
//   addPoint, halfPlaneIntersection: O(log N + K) for K removed points
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
//   https://codeforces.com/problemsets/acmsguru/problem/99999/277
//   https://open.kattis.com/problems/bigbrother
//   https://open.kattis.com/problems/marshlandrescues
//   https://dmoj.ca/problem/ccoprep3p3
//   https://www.acmicpc.net/problem/4225
struct IncrementalConvexHull : public set<Ray> {
  using iter = set<Ray>::iterator; T a2; IncrementalConvexHull() : a2(0) {}
  iter mod(iter it) const { return !empty() && it == end() ? begin() : it; }
  iter prv(iter it) const { return prev(it == begin() ? end() : it); }
  iter nxt(iter it) const { return mod(next(it)); }
  void erase(iter it) {
    iter a = prv(it), b = nxt(it); set<Ray>::erase(it); if (empty()) return;
    a->nxt = &*b;
  }
  void emplace(pt p, Line l) {
    iter it = set<Ray>::emplace(p, l).first;
    prv(it)->nxt = &*it; it->nxt = &*nxt(it);
  }
  iter rem(iter it) {
    iter a = prv(it), b = nxt(it); a2 += cross(a->p, b->p);
    a2 -= cross(a->p, it->p) + cross(it->p, b->p); erase(it); return b;
  }
  void clear() { set<Ray>::clear(); a2 = 0; }
  int isIn(pt p) const {
    if (empty()) return -1;
    pt a = begin()->p, b = prev(end())->p; if (onSeg(p, a, b)) return 0;
    auto it = lower_bound(IsInCmp(a, b, p)); pt q = nxt(it)->p;
    if (onSeg(p, it->p, q)) return 0;
    vector<pt> P{a, it->p, q}; sort(P.begin(), P.end());
    P.erase(unique(P.begin(), P.end()), P.end());
    if (int(P.size()) < 3) return onSeg(p, P[0], P.back()) ? 0 : -1;
    return ccw(it->p, q, p) >= 0 ? 1 : -1;
  }
  iter singlePointTangent(pt p, bool left) const {
    pt a = begin()->p, b = prev(end())->p; int o = ccw(p, a, b);
    bool farSide = o ? o < 0 : lt(distSq(p, b), distSq(p, a));
    return mod(lower_bound(PointTangentCmp(a, b, p, left, farSide)));
  }
  pair<iter, iter> pointTangents(pt p) const {
    return make_pair(singlePointTangent(p, 1), singlePointTangent(p, 0));
  }
  pair<iter, iter> circleTangents(Circle c, bool inner) const {
    pair<iter, iter> ret; pt a = begin()->p, b = prev(end())->p;
    for (int h = 0; h < 2; h++) {
      vector<pair<pt, pt>> t;
      assert(circleCircleTangent(Circle(a, 0), c, inner, t) == 1);
      pt q = t[h].second; int o = ccw(q, a, b);
      bool farSide = o ? o < 0 : lt(distSq(q, b), distSq(q, a));
      (h ? ret.second : ret.first) = mod(lower_bound(CircleTangentCmp(
          a, b, c, inner, h, farSide)));
    }
    return ret;
  }
  pt closestPt(pt p) const {
    auto f = [&] (iter a, iter b) {
      if (a == b) return a->p;
      if (a != (b = prv(b)))
        a = mod(lower_bound(ClosestPointCmp(
            p, Angle(a->l.v), Angle(prv(b)->l.v))));
      return closestPtOnSeg(p, a->p, a->nxt->p);
    };
    iter a, b; tie(a, b) = pointTangents(p); if (*b < *a) {
      pt q = f(a, prev(end())), r = f(begin(), b);
      pt s = distSq(p, q) < distSq(p, r) ? q : r;
      pt t = closestPtOnSeg(p, prev(end())->p, begin()->p);
      return distSq(p, s) < distSq(p, t) ? s : t;
    } else return f(a, b);
  }
  vector<pair<iter, iter>> hullTangents(const IncrementalConvexHull &hull,
                                        bool inner) const;
  bool addPoint(pt p) {
    if (empty()) { emplace(p, Line(p, p)); return true; }
    if (isIn(p) >= 0) return false;
    iter l, r; tie(l, r) = pointTangents(p); l = prv(l);
    if (ccw(l->p, nxt(l)->p, p) > 0) l = nxt(l);
    if (ccw(p, r->p, nxt(r)->p) == 0) r = nxt(r);
    for (l = nxt(l); l != r; l = rem(l));
    l = prv(l); pt a = l->p, b = r->p;
    a2 += cross(a, p) + cross(p, b) - cross(a, b);
    erase(l); emplace(a, Line(a, p)); emplace(p, Line(p, b)); return true;
  }
  iter extremeVertex(pt dir) const {
    return mod(lower_bound(Ray(pt(0, 0), Line(perp(dir), 0))));
  }
  vector<pt> lineIntersection(Line l) const {
    vector<pt> ret;
    auto check = [&] (iter a) {
      pt p = a->p, q = a->nxt->p; if (p != q) {
        pt r; lineLineIntersection(l, a->l, r);
        if (onSeg(r, p, q)) ret.push_back(r);
      }
      if (l.onLeft(p) == 0) ret.push_back(p);
      if (l.onLeft(q) == 0) ret.push_back(q);
    };
    auto f = [&] (iter a, iter b) {
      if (a == b) check(a);
      else check(mod(lower_bound(LineIntersectionCmp(
          l, Angle(a->l.v), Angle(prv(b)->l.v)))));
    };
    Line l2(-l.v, -l.c); iter a = mod(lower_bound(Ray(pt(0, 0), l)));
    iter b = mod(lower_bound(Ray(pt(0, 0), l2)));
    if (*b < *a) { std::swap(a, b); std::swap(l, l2); }
    f(a, b); std::swap(l, l2); f(b, prev(end())); f(begin(), a);
    sort(ret.begin(), ret.end());
    ret.erase(unique(ret.begin(), ret.end()), ret.end()); return ret;
  }
  void halfPlaneIntersection(Line l) {
    if (empty()) return;
    iter b = mod(lower_bound(Ray(pt(0, 0), l)));
    if (l.onLeft(b->p) >= 0) return;
    iter a = prv(b); while (a != b && l.onLeft(a->p) < 0) a = prv(a);
    if (a == b) { clear(); return; }
    iter c = nxt(a); while (c != b) c = rem(c);
    c = nxt(b); while (l.onLeft(c->p) < 0) c = nxt(b = rem(b));
    pt p, q, r = a->p, s = c->p; Line la = a->l, lb = b->l;
    if (lineLineIntersection(la, l, p) != 1
        || lineLineIntersection(l, lb, q) != 1) {
      assert(lineLineIntersection(la, lb, p) == 1); q = p;
    }
    a2 -= cross(r, b->p) + cross(b->p, s); erase(a); erase(b);
    a2 += cross(r, p) + cross(p, q) + cross(q, s);
    if (r != p) emplace(r, la);
    if (q != s) emplace(q, lb);
    if ((r == p && q == s) || p != q) emplace(p, l);
  }
};

// Helper struct for hullTangents
struct HullTangentCmp : public Ray {
  pt q; const IncrementalConvexHull &hull; bool inner, h, farSide;
  HullTangentCmp(pt p, pt q, const IncrementalConvexHull &hull,
                 bool inner, bool h, bool farSide)
      : Ray(p), q(q), hull(hull), inner(inner), h(h), farSide(farSide) {}
  bool cmp(const Ray &o) const override {
    if (farSide == h && o.p == p) return true;
    if (farSide != h && o.p == q) return false;
    pt q = hull.singlePointTangent(o.p, inner ^ h)->p;
    if (ccw(q, p, o.p) == (h ? 1 : -1)) return farSide == h;
    else return (ccw(o.p, o.nxt->p, q) < 0) == h;
  }
};

vector<pair<IncrementalConvexHull::iter, IncrementalConvexHull::iter>>
    IncrementalConvexHull::hullTangents(const IncrementalConvexHull &hull,
                                        bool inner) const {
  vector<pair<iter, iter>> ret(2); pt a = begin()->p, b = prev(end())->p;
  for (int h = 0; h < 2; h++) {
    pt q = hull.singlePointTangent(a, inner ^ h)->p; int o = ccw(q, a, b);
    bool farSide = o ? o < 0 : lt(distSq(q, b), distSq(q, a));
    ret[h].first = mod(lower_bound(HullTangentCmp(
          a, b, hull, inner, h, farSide)));
    ret[h].second = hull.singlePointTangent(ret[h].first->p, inner ^ h);
  }
  return ret;
}


text:
#pragma once
#include <bits/stdc++.h>
#include "../../search/BinarySearch.h"
using namespace std;

// Supports adding lines in the form f(x) = mx + b and finding
//   the maximum value of f(x) at any given x; this version only supports
//   adding lines in sorted order
// Template Arguments:
//   T: the type of the slope (m) and intercept (b) of the line, as well as
//     the type of the function argument (x), must be able to store m * b
//   Cmp: the comparator to compare two f(x) values,
//       convention is same as std::priority_queue in STL
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
//   REVERSE: boolean to indicate whether the lines are added in reverse order
//     of the comparator
// Constructor Arguments:
//   cmp: an instance of the Cmp struct
// Functions:
//   addLine(m, b): adds a line in the form f(x) = mx + b to the set of lines,
//     lines must be added in the order of slope sorted by Cmp, or in reverse
//     order if REVERSE is true
//   getMax(x): finds the maximum x value (based on the comparator) for all
//     inserted lines
//   getMaxMonoInc(x): finds the maximum x value (based on the comparator)
//     for all inserted lines, where all queries have non decreasing x values
//     (based on the standard < operator)
//   getMaxMonoDec(x): finds the maximum x value (based on the comparator)
//     for all inserted lines, where all queries have non increasing x values
//     (based on the standard < operator)
//   size(): returns the number of lines in the convex hull
//   reserve(N): reserves space for N lines in the convex hull
// In practice, has a moderate constant
// Time Complexity:
//   addLine, getMaxMonoInc, getMaxMonoDec: O(1) amortized
//   getMax: O(log N) for N lines in the convex hull
//   size: O(1)
//   reserve: O(N)
// Memory Complexity: O(N) for N lines in the convex hull
// Tested:
//   https://www.spoj.com/problems/CHTPRAC/
//   https://atcoder.jp/contests/dp/tasks/dp_z
template <class T, class Cmp = less<T>, const bool REVERSE = false>
struct ConvexHullTrick {
  struct Line {
    T m, b; Line(T m, T b) : m(m), b(b) {}
    T eval(T x) const { return m * x + b; }
  };
  vector<Line> L; Cmp cmp; int front, back;
  ConvexHullTrick(Cmp cmp = Cmp()) : cmp(cmp), front(0), back(0) {}
  int size() const { return L.size(); }
  void addLine(T m, T b) {
    auto ccw = [&] {
      T c1 = (L.back().m - m) * (L[size() - 2].b - b);
      T c2 = (L.back().b - b) * (L[size() - 2].m - m);
      return REVERSE ? c1 >= c2 : c1 <= c2;
    };
    while (!L.empty() && !cmp(m, L.back().m) && !cmp(L.back().m, m)
           && !cmp(b, L.back().b))
      L.pop_back();
    while (size() >= 2 && ccw()) L.pop_back();
    if (size() == back) back++;
    L.emplace_back(m, b);
    front = min(front, size() - 1); back = min(back, size());
  }
  T moveFront(T x) {
    while (front + 1 < size() && !cmp(L[front + 1].eval(x), L[front].eval(x)))
      front++;
    return L[front].eval(x);
  }
  T moveBack(T x) {
    while (back - 2 >= 0 && !cmp(L[back - 2].eval(x), L[back - 1].eval(x)))
      back--;
    return L[back - 1].eval(x);
  }
  T getMaxMonoInc(T x) { return REVERSE ? moveBack(x) : moveFront(x); }
  T getMaxMonoDec(T x) { return REVERSE ? moveFront(x) : moveBack(x); }
  T getMax(T x) const {
    return L[bsearch<FIRST>(0, size() - 1, [&] (int i) {
      return cmp(L[i + 1].eval(x), L[i].eval(x));
    })].eval(x);
  }
  void reserve(int N) { L.reserve(N); }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
#include "Angle.h"
#include "Line.h"
#include "../../datastructures/IntervalUnion.h"
using namespace std;

// Functions for a 2D circle
struct Circle {
  pt o; T r; Circle(T r = 0) : o(0, 0), r(r) {}
  Circle(pt o, T r) : o(o), r(r) {}
  // 1 if p is inside this circle, 0 if on this circle,
  //   -1 if outside this circle
  int contains(pt p) const { return sgn(r * r - distSq(o, p)); }
  // 1 if c is strictly inside this circle, 0 if inside and touching this
  //   circle, -1 otherwise
  int contains(Circle c) const {
    T dr = r - c.r; return lt(dr, 0) ? -1 : sgn(dr * dr - distSq(o, c.o));
  }
  // 1 if c is strictly outside this circle, 0 if outside and touching this
  //   circle, -1 otherwise
  int disjoint(Circle c) const {
    T sr = r + c.r; return sgn(distSq(o, c.o) - sr * sr);
  }
  pt proj(pt p) const { return o + (p - o) * r / dist(o, p); }
  pt inv(pt p) const { return o + (p - o) * r * r / distSq(o, p); }
};

// Determine the intersection of a circle and a line
// Function Arguments:
//   c: the circle
//   l: the line
// Return Value: the points of intersection, if any, of the circle and
//   the line, guaranteed to be sorted based on projection on the line
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/noi05p6
vector<pt> circleLineIntersection(Circle c, Line l) {
  vector<pt> ret; T h2 = c.r * c.r - l.distSq(c.o); pt p = l.proj(c.o);
  if (eq(h2, 0)) ret.push_back(p);
  else if (lt(0, h2)) {
    pt h = l.v * sqrt(h2) / abs(l.v);
    ret.push_back(p - h); ret.push_back(p + h);
  }
  return ret;
}

// Determine the intersection of a circle and a line segment
// Function Arguments:
//   c: the circle
//   a: the first point on the line segment
//   b: the second point on the line segment
// Return Value: the points of intersection, if any, of the circle and
//   the line segment, guaranteed to be sorted based on projection on the line
//   from a to b
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/noi05p6
vector<pt> circleSegIntersection(Circle c, pt a, pt b) {
  vector<pt> ret; if (a == b) { if (c.contains(a) == 0) ret.push_back(a); }
  else {
    Line l(a, b); for (auto &&p : circleLineIntersection(c, l))
      if (l.cmpProj(a, p) <= 0 && l.cmpProj(p, b) <= 0) ret.push_back(p);
  }
  return ret;
}

// Determine the area of the intersection of a circle and a half-plane defined
//   by the left side of a line
// Function Arguments:
//   c: the circle
//   l: the line with the half-plane defined by the left side
// Return Value: the area of the intersection of the circle and the half-plane
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/noi05p6
T circleHalfPlaneIntersectionArea(Circle c, Line l) {
  T h2 = c.r * c.r - l.distSq(c.o), ret = 0; if (!lt(h2, 0)) {
    pt p = l.proj(c.o), h = l.v * sqrt(max(h2, T(0))) / abs(l.v);
    pt a = p - h, b = p + h; T theta = abs(ang(a, c.o, b));
    ret = c.r * c.r * (theta - sin(theta)) / 2;
  }
  if (l.onLeft(c.o) > 0) ret = acos(T(-1)) * c.r * c.r - ret;
  return ret;
}

// Determine the intersection of two circles
// Function Arguments:
//   c1: the first circle
//   c2: the second circle
//   res: the points of intersection, if any, of the two circles;
//     the first point is guaranteed to not be on the left side of the
//     line from c1.o to c2.o
// Return Value: 0 if no intersection, 2 if identical circles, 1 otherwise
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://codeforces.com/contest/420/problem/E
//   https://open.kattis.com/problems/drawingcircles
//   https://dmoj.ca/problem/noi05p6
int circleCircleIntersection(Circle c1, Circle c2, vector<pt> &res) {
  pt d = c2.o - c1.o; T d2 = norm(d);
  if (eq(d2, 0)) return eq(c1.r, c2.r) ? 2 : 0;
  T pd = (d2 + c1.r * c1.r - c2.r * c2.r) / 2;
  T h2 = c1.r * c1.r - pd * pd / d2; pt p = c1.o + d * pd / d2;
  if (eq(h2, 0)) res.push_back(p);
  else if (lt(0, h2)) {
    pt h = perp(d) * sqrt(h2 / d2); res.push_back(p - h); res.push_back(p + h);
  }
  return !res.empty();
}

// Determine the area of the intersection of two circles
// Function Arguments:
//   c1: the first circle
//   c2: the second circle
// Return Value: the area of the intersection of the two circles
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://codeforces.com/contest/600/problem/D
T circleCircleIntersectionArea(Circle c1, Circle c2) {
  T d = dist(c1.o, c2.o); if (!lt(d, c1.r + c2.r)) return 0;
  if (!lt(c2.r, d + c1.r)) return acos(T(-1)) * c1.r * c1.r;
  if (!lt(c1.r, d + c2.r)) return acos(T(-1)) * c2.r * c2.r;
  auto A = [&] (T r1, T r2) {
    T a = (d * d + r1 * r1 - r2 * r2) / (2 * d * r1);
    T theta = 2 * acos(max(T(-1), min(T(1), a)));
    return r1 * r1 * (theta - sin(theta)) / 2;
  };
  return A(c1.r, c2.r) + A(c2.r, c1.r);
}

// Determine the tangents of two circles
// Function Arguments:
//   c1: the first circle
//   c2: the second circle
//   inner: whether to find the inner or outer tangents
//   res: a vector of pairs of size 2 of the tangents, with each pair
//     representing a point on the first circle and the second circle;
//     the first point in each pair is guaranteed to not be on the left side of
//     the line from c1.o to c2.o
// Return Value: 0 if no tangents, 2 if identical circles, 1 otherwise
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/nccc7s4
//   https://dmoj.ca/problem/noi05p6
int circleCircleTangent(Circle c1, Circle c2, bool inner,
                        vector<pair<pt, pt>> &res) {
  pt d = c2.o - c1.o; T r2 = inner ? -c2.r : c2.r, dr = c1.r - r2;
  T d2 = norm(d), h2 = d2 - dr * dr;
  if (eq(d2, 0) || lt(h2, 0)) return eq(h2, 0) ? 2 : 0;
  for (T sign : {T(-1), T(1)}) {
    pt v = (d * dr + perp(d) * sqrt(max(h2, T(0))) * sign) / d2;
    res.emplace_back(c1.o + v * c1.r, c2.o + v * r2);
  }
  return 1;
}

// Determines the circumcircle from 3 non-collinear points
// Function Arguments:
//   a: the first point
//   b: the second point
//   c: the third point
// Return Value: the circumcircle of the 3 points
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://www.spoj.com/problems/QCJ4/
Circle circumcircle(pt a, pt b, pt c) {
  b -= a; c -= a;
  pt ret = b * c * (conj(c) - conj(b)) / (b * conj(c) - conj(b) * c);
  return Circle(a + ret, abs(ret));
}

// Computes the area of union of multiple circles
// Function Arguments:
//   circles: a vector of the circles
// Return Value: the area of the union of all the circles
// Time Complexity: O(N^2 log N)
// Memory Complexity: O(N)
// Tested:
//   https://www.spoj.com/problems/CIRU/
T circleUnionArea(const vector<Circle> &circles) {
  int n = circles.size(); T ret = 0; for (int i = 0; i < n; i++) {
    vector<pair<Angle, Angle>> intervals; Angle::setPivot(circles[i].o);
    bool inside = false; for (int j = 0; j < n; j++) if (i != j) {
      int o = circles[j].contains(circles[i]);
      if (o > 0 || (o == 0 && (lt(circles[i].r, circles[j].r) || j < i))) {
        inside = true; break;
      }
      vector<pt> p; circleCircleIntersection(circles[i], circles[j], p);
      if (int(p.size()) == 2) {
        Angle a(p[0]), b(p[1]); if (a < b) intervals.emplace_back(a, b);
        else {
          intervals.emplace_back(a, Angle(circles[i].o));
          Angle c(circles[i].o - pt(circles[i].r, 0));
          intervals.emplace_back(c, b);
        }
      }
    }
    if (inside) continue;
    if (intervals.empty()) ret += acos(T(-1)) * circles[i].r * circles[i].r;
    else {
      intervalUnion(intervals); if (intervals.back().second == circles[i].o) {
        intervals.front().first = intervals.back().first; intervals.pop_back();
      }
      for (int j = 0, k = int(intervals.size()); j < k; j++) {
        pt a = intervals[j].second.p;
        pt b = intervals[j + 1 == k ? 0 : j + 1].first.p;
        ret += cross(a, b) / 2;
        ret += circleHalfPlaneIntersectionArea(circles[i], Line(b, a));
      }
    }
  }
  return ret;
}


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Helper struct for IncrementalConvexHullTrick
template <class T> struct CHTLine {
  bool isQuery; T m, b; mutable T x;
  CHTLine(T m, T b) : isQuery(false), m(m), b(b), x(T()) {}
  CHTLine(T x) : isQuery(true), m(T()), b(T()), x(x) {}
};

// Supports adding lines in the form f(x) = mx + b and finding
//   the maximum value of f(x) at any given x; this version allows for
//   updates and queries in arbitrary order
// Template Arguments:
//   T: the type of the slope (m) and intercept (b) of the line, as well as
//     the type of the function argument (x)
//   Cmp: the comparator to compare two f(x) values,
//       convention is same as std::priority_queue in STL
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
// Constructor Arguments:
//   cmp: an instance of the Cmp struct
//   INF: a value for positive infinity, must be negatable
// Functions (in addition to std::multiset):
//   addLine(m, b): adds a line in the form f(x) = mx + b to the set of lines
//   getMax(x): finds the maximum value of f(x) (based on the comparator)
//     for all inserted lines
// In practice, has a moderate constant, performance compared to
//   IncrementalConvexHullTrickSqrtBuffer and SparseLiChao can vary, slower
//   than LiChao
// Time Complexity:
//   constructor: O(1)
//   addLine, getMax: O(log(N)) amortized for N lines in the convex hull
// Memory Complexity: O(N) for N lines in the convex hull
// Tested:
//   https://judge.yosupo.jp/problem/line_add_get_min
//   https://open.kattis.com/problems/longestlife
//   https://www.spoj.com/problems/CHTPRAC/
#define FUN function<bool(const CHTLine<T> &, const CHTLine<T> &)>
template <class T, class Cmp = less<T>>
struct IncrementalConvexHullTrick : public multiset<CHTLine<T>, FUN> {
  using MS = multiset<CHTLine<T>, FUN>; using iter = typename MS::iterator;
  using MS::begin; using MS::end; using MS::emplace; using MS::erase;
  using MS::lower_bound; Cmp cmp; T INF;
  FUN makeFun(Cmp cmp) {
    return [=] (const CHTLine<T> &a, const CHTLine<T> &b) {
      return a.isQuery || b.isQuery ? a.x < b.x : cmp(a.m, b.m);
    };
  }
  IncrementalConvexHullTrick(Cmp cmp = Cmp(), T INF = numeric_limits<T>::max())
      : MS(makeFun(cmp)), cmp(cmp), INF(INF) {}
  template <const bool _ = is_integral<T>::value
                || is_same<__int128_t, T>::value>
  typename enable_if<!_, T>::type div(T a, T b) { return a / b; }
  template <const bool _ = is_integral<T>::value
                || is_same<__int128_t, T>::value>
  typename enable_if<_, T>::type div(T a, T b) {
    return a / b - T((a < T()) != (b < T()) && T() != a % b);
  }
  bool intersect(iter x, iter y) {
    if (y == end()) { x->x = INF; return false; }
    if (!cmp(x->m, y->m) && !cmp(y->m, x->m))
      x->x = cmp(y->b, x->b) ? INF : -INF;
    else x->x = div(y->b - x->b, x->m - y->m);
    return x->x >= y->x;
  }
  void addLine(T m, T b) {
    auto z = emplace(m, b), y = z++, x = y;
    while (intersect(y, z)) z = erase(z);
    if (x != begin() && intersect(--x, y)) intersect(x, y = erase(y));
    while ((y = x) != begin() && (--x)->x >= y->x) intersect(x, erase(y));
  }
  T getMax(T x) const {
    auto l = *lower_bound(CHTLine<T>(x)); return l.m * x + l.b;
  }
};
#undef FUN


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
#include "Angle.h"
#include "Line.h"
using namespace std;

// Helper struct for polygonTriangulation
struct Chain {
  static T X; mutable Line l;
  Chain(T Y) : l(0, 1, Y) {}
  Chain(pt p, pt q) : l(p, q) {}
  bool operator < (const Chain &c) const {
    return lt((l.c + l.v.y * X) * c.l.v.x, (c.l.c + c.l.v.y * X) * l.v.x);
  }
};

// Helper struct for polygonTriangulation
struct ChainInfo {
  int chainId, last, link;
  ChainInfo(int chainId, int last)
      : chainId(chainId), last(last), link(last) {}
};

T Chain::X = 0;

// Triangulates a polygon such that all N - 2 triangles are disjoint, with the
//   union of the triangles being equal to the original polygon
// Function Arguments:
//   P: the points of the simple polygon in ccw order, no two points can
//     be identical
// Return Value: a vector of arrays of 3 points, representing a triangle in
//   the triangulation, points given in ccw order
// In practice, has a moderate constant
// Time Complexity: O(N log N)
// Memory Complexity: O(N)
// Tested:
//   Fuzz Tested
vector<array<pt, 3>> polygonTriangulation(const vector<pt> &P) {
  vector<array<pt, 3>> ret;
  auto monotone = [&] (const vector<int> &ind) {
    auto add = [&] (int i, int j, int k) {
      ret.push_back(array<pt, 3>{P[ind[i]], P[ind[j]], P[ind[k]]});
    };
    auto pccw = [&] (int i, int j, int k) {
      return ccw(P[ind[i]], P[ind[j]], P[ind[k]]);
    };
    int st = min_element(ind.begin(), ind.end(), [&] (int i, int j) {
      return P[i].x < P[j].x;
    }) - ind.begin();
    auto prv = [&] (int i) { return i == 0 ? int(ind.size()) - 1 : i - 1; };
    auto nxt = [&] (int i) { return i + 1 == int(ind.size()) ? 0 : i + 1; };
    deque<int> lo{st}, hi{st}; for (int k = 1; k < int(ind.size()); k++) {
      int a = nxt(lo.back()), b = prv(hi.back());
      if (P[ind[a]].x < P[ind[b]].x) {
        while (int(lo.size()) >= 2
            && pccw(lo[lo.size() - 2], lo.back(), a) == 1) {
          add(lo[lo.size() - 2], lo.back(), a); lo.pop_back();
        }
        while (int(hi.size()) >= 2 && pccw(a, hi[1], hi[0]) == 1) {
          add(a, hi[1], hi[0]); hi.pop_front(); lo.front() = hi.front();
        }
        lo.push_back(a);
      } else {
        while (int(hi.size()) >= 2
            && pccw(b, hi.back(), hi[hi.size() - 2]) == 1) {
          add(b, hi.back(), hi[hi.size() - 2]); hi.pop_back();
        }
        while (int(lo.size()) >= 2 && pccw(lo[0], lo[1], b) == 1) {
          add(lo[0], lo[1], b); lo.pop_front(); hi.front() = lo.front();
        }
        hi.push_back(b);
      }
    }
  };
  int n = P.size(); vector<int> ord(n); iota(ord.begin(), ord.end(), 0);
  sort(ord.begin(), ord.end(), [&] (int a, int b) { return P[a].x < P[b].x; });
  auto isAdj = [&] (int i, int j) {
    int d = abs(i - j); return d == 1 || d == n - 1;
  };
  auto prv = [&] (int i) { return i == 0 ? n - 1 : i - 1; };
  auto nxt = [&] (int i) { return i + 1 == n ? 0 : i + 1; };
  vector<vector<pair<int, int>>> G(n); int curEdge = 0, curChain = 0;
  auto addBiEdge = [&] (int i, int j) {
    if (i == -1 || j == -1 || isAdj(i, j)) return;
    G[i].emplace_back(j, curEdge++); G[j].emplace_back(i, curEdge++);
  };
  for (int i = 0; i < n; i++) G[i].emplace_back(nxt(i), curEdge++);
  multimap<Chain, ChainInfo> lo, hi; for (int l = 0, r = 0; l < n; l = r) {
    for (r = l + 1; r < n && !lt(P[ord[l]].x, P[ord[r]].x); r++);
    sort(ord.begin() + l, ord.begin() + r, [&] (int a, int b) {
      return P[a].y < P[b].y;
    });
    Chain::X = P[ord[l]].x; for (int d = l, u = l; d < r; d = u) {
      for (u = d + 1; u < r && isAdj(ord[u - 1], ord[u]); u++);
      auto it1 = lo.upper_bound(Chain(P[ord[u - 1]].y));
      auto it2 = hi.lower_bound(Chain(P[ord[d]].y));
      if (it1 == lo.begin() || it2 == hi.end()
          || ((--it1)->second.chainId != it2->second.chainId
                && (nxt(it1->second.last) != ord[u - 1]
                      || prv(it2->second.last) != ord[d]))) {
        lo.emplace(Chain(P[ord[d]], P[nxt(ord[d])]),
                   ChainInfo(curChain, ord[d]));
        hi.emplace(Chain(P[ord[u - 1]], P[prv(ord[u - 1])]),
                   ChainInfo(curChain++, ord[u - 1]));
      } else if (it1->second.chainId != it2->second.chainId) {
        auto it3 = next(it2), it4 = prev(it1);
        addBiEdge(ord[u - 1], it1->second.link);
        addBiEdge(ord[d], it2->second.link);
        addBiEdge(ord[u - 1], it3->second.link);
        addBiEdge(ord[d], it4->second.link);
        lo.erase(it1); hi.erase(it2);
        it3->second.chainId = it4->second.chainId;
        it3->second.link = -1; it4->second.link = ord[u - 1];
      } else {
        bool conLo = nxt(it1->second.last) == ord[d];
        bool conHi = prv(it2->second.last) == ord[u - 1];
        if (conLo && conHi) {
          addBiEdge(ord[d], it1->second.link);
          addBiEdge(ord[d], it2->second.link);
          lo.erase(it1); hi.erase(it2);
        } else if (conLo) {
          addBiEdge(ord[d], it1->second.link);
          addBiEdge(ord[d], it2->second.link);
          it1->first.l = Line(P[ord[u - 1]], P[nxt(ord[u - 1])]);
          it1->second.last = it1->second.link = ord[u - 1];
          it2->second.link = -1;
        } else if (conHi) {
          addBiEdge(ord[d], it1->second.link);
          addBiEdge(ord[d], it2->second.link);
          it2->first.l = Line(P[ord[d]], P[prv(ord[d])]);
          it2->second.last = it2->second.link = ord[d];
          it1->second.link = -1;
        } else if (nxt(it1->second.last) == ord[u - 1]) {
          addBiEdge(ord[u - 1], it1->second.link);
          addBiEdge(ord[u - 1], it2->second.link);
          it1->first.l = Line(P[ord[d]], P[nxt(ord[d])]);
          it1->second.last = ord[d]; it1->second.link = ord[u - 1];
          it2->second.link = -1;
        } else if (prv(it2->second.last) == ord[d]) {
          addBiEdge(ord[d], it1->second.link);
          addBiEdge(ord[d], it2->second.link);
          it2->first.l = Line(P[ord[u - 1]], P[prv(ord[u - 1])]);
          it2->second.last = it2->second.link = ord[u - 1];
          it1->second.link = -1;
        } else {
          addBiEdge(ord[d], it1->second.link);
          addBiEdge(ord[d], it2->second.link);
          hi.emplace(Chain(P[ord[d]], P[prv(ord[d])]),
                     ChainInfo(it1->second.chainId, ord[d]));
          it1->second.link = it2->second.link = -1;
          it2->second.chainId = curChain;
          lo.emplace(Chain(P[ord[u - 1]], P[nxt(ord[u - 1])]),
                     ChainInfo(curChain++, ord[u - 1]));
        }
      }
    }
  }
  auto cmpAng = [&] (const pair<int, int> &a, const pair<int, int> &b) {
    return Angle(P[a.first]) < Angle(P[b.first]);
  };
  auto eqAng = [&] (const pair<int, int> &a, const pair<int, int> &b) {
    return Angle(P[a.first]) == Angle(P[b.first]);
  };
  for (int i = 0; i < n; i++) {
    Angle::setPivot(P[i]); sort(G[i].begin(), G[i].end(), cmpAng);
    G[i].erase(unique(G[i].begin(), G[i].end(), eqAng), G[i].end());
  }
  vector<bool> vis(curEdge, false);
  for (int i = 0; i < n; i++) for (auto &&e : G[i]) if (!vis[e.second]) {
    vector<int> ind; int cur = i, nxt = e.first, eid = e.second;
    while (!vis[eid]) {
      ind.push_back(cur); vis[eid] = true; Angle::setPivot(P[nxt]);
      auto it = lower_bound(G[nxt].begin(), G[nxt].end(),
                            make_pair(cur, -1), cmpAng);
      it = prev(it == G[nxt].begin() ? G[nxt].end() : it);
      cur = nxt; nxt = it->first; eid = it->second;
    }
    monotone(ind);
  }
  return ret;
}


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
using namespace std;

// Computes the convex hull of a set of N points
// Function Arguments:
//   P: the vector of points
// Return Value: a vector of points in the convex hull in ccw order
// In practice, has a small constant
// Time Complexity: O(N log N)
// Memory Complexity: O(N)
// Tested:
//   https://open.kattis.com/problems/convexhull
vector<pt> convexHull(vector<pt> P) {
  vector<pt> hull; sort(P.begin(), P.end()); for (int h = 0; h < 2; h++) {
    int st = hull.size(); for (auto &&p : P) {
      while (int(hull.size()) >= st + 2
          && ccw(hull[hull.size() - 2], hull.back(), p) <= 0)
        hull.pop_back();
      hull.push_back(p);
    }
    hull.pop_back(); reverse(P.begin(), P.end());
  }
  if (int(hull.size()) == 2 && hull[0] == hull[1]) hull.pop_back();
  if (hull.empty() && !P.empty()) hull.push_back(P[0]);
  return hull;
}


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
using namespace std;

// Computes the closest pair of points out of a set of N points
// Constructor Arguments:
//   P: the vector of points
// Fields:
//   best1: the first point in the pair of closest points
//   best2: the second point in the pair of closest points
//   bestDistSq: the distance squared between best1 and best2
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(N log N)
// Memory Complexity: O(N)
// Tested:
//   https://open.kattis.com/problems/closestpair2
struct ClosestPair {
  static bool yOrdLt(pt p, pt q) { return lt(p.y, q.y); }
  pt best1, best2; T bestDistSq;
  void closest(vector<pt> &P, int lo, int hi) {
    if (hi <= lo) return;
    int mid = lo + (hi - lo) / 2; pt median = P[mid];
    closest(P, lo, mid); closest(P, mid + 1, hi); vector<pt> aux;
    merge(P.begin() + lo, P.begin() + mid + 1, P.begin() + mid + 1,
          P.begin() + hi + 1, back_inserter(aux), yOrdLt);
    copy(aux.begin(), aux.end(), P.begin() + lo);
    for (int i = lo, k = 0; i <= hi; i++) {
      T dx = P[i].x - median.x, dx2 = dx * dx; if (lt(dx2, bestDistSq)) {
        for (int j = k - 1; j >= 0; j--) {
          T dy = P[i].y - aux[j].y, dy2 = dy * dy;
          if (!lt(dy2, bestDistSq)) break;
          T dSq = distSq(P[i], aux[j]); if (lt(dSq, bestDistSq)) {
            bestDistSq = dSq; best1 = P[i]; best2 = aux[j];
          }
        }
        aux[k++] = P[i];
      }
    }
  }
  ClosestPair(vector<pt> P) : bestDistSq(numeric_limits<T>::max()) {
    sort(P.begin(), P.end()); closest(P, 0, int(P.size()) - 1);
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
#include "ConvexHull.h"
using namespace std;

// Computes the farthest pair of points out of a set of N points
// Constructor Arguments:
//   P: the vector of points
// Fields:
//   best1: the first point in the pair of farthest points
//   best2: the second point in the pair of farthest points
//   bestDistSq: the distance squared between best1 and best2
//   hull: the points in the convex hull of P
// In practice, has a small constant
// Time Complexity:
//   constructor: O(N log N)
// Memory Complexity: O(N)
// Tested:
//   Fuzz Tested
struct FarthestPair {
  pt best1, best2; T bestDistSq; vector<pt> hull;
  FarthestPair(const vector<pt> &P) : bestDistSq(0), hull(convexHull(P)) {
    int H = hull.size(); pt o(0, 0);
    for (int i = 0, j = H < 2 ? 0 : 1; i < j; i++) for (;; j = (j + 1) % H) {
      T dSq = distSq(hull[i], hull[j]); if (lt(bestDistSq, dSq)) {
        bestDistSq = dSq; best1 = hull[i]; best2 = hull[j];
      }
      pt a = hull[i + 1] - hull[i], b = hull[(j + 1) % H] - hull[j];
      if (ccw(o, a, b) <= 0) break;
    }
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
using namespace std;

// Functions for a 2D line
struct Line {
  pt v; T c;
  // ax + by = c, left side is ax + by > c
  Line(T a = 0, T b = 0, T c = 0) : v(b, -a), c(c) {}
  // direction vector v with offset c
  Line(pt v, T c) : v(v), c(c) {}
  // points p and q
  Line(pt p, pt q) : v(q - p), c(cross(v, p)) {}
  T eval(pt p) const { return cross(v, p) - c; }
  // sign of onLeft, dist: 1 if left of line, 0 if on line, -1 if right of line
  int onLeft(pt p) const { return sgn(eval(p)); }
  T dist(pt p) const { return eval(p) / abs(v); }
  T distSq(pt p) const { T e = eval(p); return e * e / norm(v); }
  // rotated 90 degrees ccw
  Line perpThrough(pt p) const { return Line(p, p + perp(v)); }
  Line translate(pt p) const { return Line(v, c + cross(v, p)); }
  Line shiftLeft(T d) const { return Line(v, c + d * abs(v)); }
  pt proj(pt p) const { return p - perp(v) * eval(p) / norm(v); }
  pt refl(pt p) const { return p - perp(v) * T(2) * eval(p) / norm(v); }
  // compares points by orthogonal projection (3 way comparison)
  int cmpProj(pt p, pt q) const { return sgn(dot(v, p) - dot(v, q)); }
};

// Bisector of 2 lines
// Function Arguments:
//   l1: the first line
//   l2: the second line
//   interior: whether the interior or exterior bisector is used
// Return Value: a line representing the bisector
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/secret
Line bisector(Line l1, Line l2, bool interior) {
  T s = interior ? 1 : -1;
  return Line(l2.v / abs(l2.v) + l1.v / abs(l1.v) * s,
              l2.c / abs(l2.v) + l1.c / abs(l1.v) * s);
}

// Intersection of 2 lines
// Function Arguments:
//   l1: the first line
//   l2: the second line
//   res: a reference to a point to store the intersection if it exists
// Return Value: 0 if no intersection, 1 if proper intersection, 2 otherwise
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/nccc7s5
int lineLineIntersection(Line l1, Line l2, pt &res) {
  T d = cross(l1.v, l2.v);
  if (eq(d, 0)) return l2.v * l1.c == l1.v * l2.c ? 2 : 0;
  res = (l2.v * l1.c - l1.v * l2.c) / d; return 1;
}

// Determines if a point is on a line segment
// Function Arguments:
//   p: the point to check if on the line segment
//   a: one endpoint of the line segment
//   b: the other endpoint of the line segment
// Return Value: true if p is on the line segment a-b, false otherwise
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/segmentintersection
bool onSeg(pt p, pt a, pt b) {
  return !ccw(p, a, b) && !lt(0, dot(a - p, b - p));
}

// Determine if two line segments intersect
// Function Arguments:
//   a: one endpoint of the first line segment
//   b: the other endpoint of the first line segment
//   p: one endpoint of the second line segment
//   q: the other endpoint of the second line segment
// Return Value: 0 if no intersection, 1 if proper intersection (a single
//   point and not an endpoint), 2 otherwise
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/segmentintersection
int segSegIntersects(pt a, pt b, pt p, pt q) {
  if (ccw(a, b, p) * ccw(a, b, q) < 0 && ccw(p, q, a) * ccw(p, q, b) < 0)
    return 1;
  if (onSeg(p, a, b) || onSeg(q, a, b) || onSeg(a, p, q) || onSeg(b, p, q))
    return 2;
  return 0;
}

// Determine the intersection of two line segments
// Function Arguments:
//   a: one endpoint of the first line segment
//   b: the other endpoint of the first line segment
//   p: one endpoint of the second line segment
//   q: the other endpoint of the second line segment
// Return Value: if the line segments do not intersect, an empty vector
//   of points; if the line segments intersect at a point, a vector containing
//   the point of intersection; if the line segments have a line segment of
//   intersection, a vector containing the two endpoints of the
//   line segment intersection (it can return more if there is precision error)
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/segmentintersection
vector<pt> segSegIntersection(pt a, pt b, pt p, pt q) {
  int intersects = segSegIntersects(a, b, p, q);
  if (!intersects) return vector<pt>();
  if (intersects == 1) {
    T c1 = cross(p - a, b - a), c2 = cross(q - a, b - a);
    return vector<pt>{(c1 * q - c2 * p) / (c1 - c2)};
  }
  vector<pt> ret; if (onSeg(p, a, b)) ret.push_back(p);
  if (onSeg(q, a, b)) ret.push_back(q);
  if (onSeg(a, p, q)) ret.push_back(a);
  if (onSeg(b, p, q)) ret.push_back(b);
  sort(ret.begin(), ret.end());
  ret.erase(unique(ret.begin(), ret.end()), ret.end()); return ret;
}

// Finds the closest point on a line segment to another point
// Function Arguments
//   p: the reference point
//   a: one endpoint of the line segment
//   b: the other endpoint of the line segment
// Return Value: the closest point to p on the line segment a-b
// Time Complexity: O(1)
// Memory Complexity: O(1)
pt closestPtOnSeg(pt p, pt a, pt b) {
  if (a == b) return a;
  T d = distSq(a, b), t = min(d, max(T(0), dot(p - a, b - a)));
  return a + (b - a) * t / d;
}

// Finds the distance to the closest point on a line segment to another point
// Function Arguments
//   p: the reference point
//   a: one endpoint of the line segment
//   b: the other endpoint of the line segment
// Return Value: the distance to the closest point to p on the line segment a-b
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/segmentdistance
T ptSegDist(pt p, pt a, pt b) {
  if (a == b) return dist(a, p);
  T d = distSq(a, b), t = min(d, max(T(0), dot(p - a, b - a)));
  return abs((p - a) * d - (b - a) * t) / d;
}

// Finds the closest distance between two line segments
// Function Arguments
//   a: one endpoint of the first line segment
//   b: the other endpoint of the first line segment
//   p: one endpoint of the second line segment
//   q: the other endpoint of the second line segment
// Return Value: the closest distance between the two line
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/segmentdistance
T segSegDist(pt a, pt b, pt p, pt q) {
  return segSegIntersects(a, b, p, q) > 0
      ? 0
      : min({ptSegDist(p, a, b), ptSegDist(q, a, b),
             ptSegDist(a, p, q), ptSegDist(b, p, q)});
}


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
using namespace std;

// Kd Tree for randomized 2d points
// Constructor Arguments:
//   xmin: the minimum x value
//   ymin: the minimum y value
//   xmax: the maximum x value
//   ymax: the maximum y value
//   P: the points
// Functions:
//   empty(): returns true if the tree is empty, false otherwise
//   size(): the number of points in the tree
//   insert(p): inserts p into the tree
//   contains(p): returns true if the tree contains p, false otherwise
//   range(rect): returns a vector of the points in the rectangle rect
//   nearest(p, nearest): sets nearest to the nearest point in the tree to p,
//     returns true if a point exists, false otherwise (tree is empty)
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(N) for N randomized points
//   empty, size: O(1)
//   insert, contains: O(log N) for N randomized points
//   range: O(sqrt N + K) for N randomized points and K points in the rectangle
//   nearest: O(log N) for N randomized points
// Memory Complexity: O(N)
// Tested:
//   https://open.kattis.com/problems/closestpair1
struct KdTree {
  struct Rectangle {
    T xmin, ymin, xmax, ymax;
    Rectangle(T xmin = 0, T ymin = 0, T xmax = 0, T ymax = 0)
        : xmin(xmin), ymin(ymin), xmax(xmax), ymax(ymax) {
      assert(xmin <= xmax); assert(ymin <= ymax);
    }
    bool intersects(const Rectangle &that) const {
      return !lt(xmax, that.xmin) && !lt(ymax, that.ymin)
          && !lt(that.xmax, xmin) && !lt(that.ymax, ymin);
    }
    bool contains(pt p) const {
      return !lt(p.x, xmin) && !lt(p.y, ymin)
          && !lt(xmax, p.x) && !lt(ymax, p.y);
    }
    T distSq(pt p) const {
      T dx = 0, dy = 0;
      if (p.x < xmin) dx = p.x - xmin;
      else if (p.x > xmax) dx = p.x - xmax;
      if (p.y < ymin) dy = p.y - ymin;
      else if (p.y > ymax) dy = p.y - ymax;
      return dx * dx + dy * dy;
    }
  };
  struct Node {
    pt p; Rectangle r; Node *lu, *rd;
    Node(pt p, const Rectangle &r) : p(p), r(r), lu(nullptr), rd(nullptr) {}
  };
  deque<Node> TR;
  static bool xOrdLt(pt p, pt q) { return lt(p.x, q.x); }
  static bool yOrdLt(pt p, pt q) { return lt(p.y, q.y); }
  Node *makeNode(pt p, const Rectangle &r) {
    TR.emplace_back(p, r); return &TR.back();
  }
  T XMIN, YMIN, XMAX, YMAX; int cnt; Node *root;
  template <class It>
  Node *build(Node *n, It st, int lo, int hi, bool partition,
              T xmin, T ymin, T xmax, T ymax) {
    if (lo > hi) return nullptr;
    int mid = lo + (hi - lo) / 2; if (partition)
      nth_element(st + lo, st + mid, st + hi + 1, xOrdLt);
    else nth_element(st + lo, st + mid, st + hi + 1, yOrdLt);
    pt p = *(st + mid); n = makeNode(p, Rectangle(xmin, ymin, xmax, ymax));
    if (partition) {
      n->lu = build(n->lu, st, lo, mid - 1, !partition,
                    xmin, ymin, n->p.x, ymax);
      n->rd = build(n->rd, st, mid + 1, hi, !partition,
                    n->p.x, ymin, xmax, ymax);
    } else {
      n->lu = build(n->lu, st, lo, mid - 1, !partition,
                    xmin, ymin, xmax, n->p.y);
      n->rd = build(n->rd, st, mid + 1, hi, !partition,
                    xmin, n->p.y, xmax, ymax);
    }
    return n;
  }
  Node *insert(Node *n, pt p, bool partition,
               T xmin, T ymin, T xmax, T ymax) {
    if (!n) { cnt++; return makeNode(p, Rectangle(xmin, ymin, xmax, ymax)); }
    if (n->p == p) return n;
    if (partition) {
      if (xOrdLt(p, n->p)) n->lu = insert(n->lu, p, !partition,
                                          xmin, ymin, n->p.x, ymax);
      else n->rd = insert(n->rd, p, !partition,
                          n->p.x, ymin, xmax, ymax);
    } else {
      if (yOrdLt(p, n->p)) n->lu = insert(n->lu, p, !partition,
                                          xmin, ymin, xmax, n->p.y);
      else n->rd = insert(n->rd, p, !partition,
                          xmin, n->p.y, xmax, ymax);
    }
    return n;
  }
  bool contains(Node *n, pt p, bool partition) {
    if (!n) return false;
    if (n->p == p) return true;
    if (partition) {
      if (xOrdLt(p, n->p)) return contains(n->lu, p, !partition);
      else return contains(n->rd, p, !partition);
    } else {
      if (yOrdLt(p, n->p)) return contains(n->lu, p, !partition);
      else return contains(n->rd, p, !partition);
    }
  }
  void range(Node *n, vector<pt> &ret, const Rectangle &rect) {
    if (!n || !rect.intersects(n->r)) return;
    if (rect.contains(n->p)) ret.push_back(n->p);
    range(n->lu, ret, rect); range(n->rd, ret, rect);
  }
  bool findNearest(Node *n, pt p, bool hasNearest, pt &nearest) {
    if (!n || (hasNearest && lt(distSq(nearest, p), n->r.distSq(p))))
      return hasNearest;
    if (!hasNearest || lt(distSq(n->p, p), distSq(nearest, p))) {
      hasNearest = true; nearest = n->p;
    }
    if (n->lu && n->lu->r.contains(p)) {
      hasNearest = findNearest(n->lu, p, hasNearest, nearest);
      hasNearest = findNearest(n->rd, p, hasNearest, nearest);
    } else {
      hasNearest = findNearest(n->rd, p, hasNearest, nearest);
      hasNearest = findNearest(n->lu, p, hasNearest, nearest);
    }
    return hasNearest;
  }
  KdTree(T xmin, T ymin, T xmax, T ymax)
      : XMIN(xmin), YMIN(ymin), XMAX(xmax), YMAX(ymax),
        cnt(0), root(nullptr) {}
    KdTree(T xmin, T ymin, T xmax, T ymax, vector<pt> P)
      : XMIN(xmin), YMIN(ymin), XMAX(xmax), YMAX(ymax), cnt(P.size()) {
    root = build(root, P.begin(), 0, cnt - 1, true, XMIN, YMIN, XMAX, YMAX);
  }
  bool empty() { return cnt == 0; }
  int size() { return cnt; }
  void insert(pt p) { root = insert(root, p, true, XMIN, YMIN, XMAX, YMAX); }
  bool contains(pt p) { return contains(root, p, true); }
  vector<pt> range(const Rectangle &rect) {
    vector<pt> ret; range(root, ret, rect); return ret;
  }
  bool findNearest(pt p, pt &nearest) {
    return findNearest(root, p, false, nearest);
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
using namespace std;

// Struct to compare angles in the range [-PI, PI) around a pivot based
//   on arg(p - pivot)
// Points directly to the left have an angle of -PI, points equal to the pivot
//   are placed after all points
// Constructor Arguments:
//   p: the point
// Fields:
//   static pivot: the pivot point
// Functions:
//   static setPivot(p): sets p as the pivot point
//   <, <=, >, >=, ==, !=: compares 2 angles
//   +, +=, -, -=: arithmetic operators
// Time Complexity:
//   constructor, setPivot, <, <=, >, >=, ==, !=, +, +=, -, -=: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://judge.yosupo.jp/problem/sort_points_by_argument
#define OP(op, body) Angle operator op (Angle a) const { return body; } \
  Angle &operator op##= (Angle a) { return *this = *this op a; }
#define CMP(op, body) bool operator op (Angle a) const { return body; }
struct Angle {
  static pt pivot; static void setPivot(pt p) { pivot = p; }
  pt p; Angle(pt p = pt(0, 0)) : p(p) {}
  int half() const {
    if (eq(p.x, pivot.x) && eq(p.y, pivot.y)) return 2;
    return int(!lt(p.y, pivot.y) && (!eq(p.y, pivot.y) || !lt(p.x, pivot.x)));
  }
  bool operator < (Angle a) const {
    int h = half() - a.half(); return h == 0 ? ccw(pivot, p, a.p) > 0 : h < 0;
  }
  CMP(<=, !(a < *this)) CMP(>, a < *this) CMP(>=, !(*this < a))
  CMP(==, !(*this < a) && !(a < *this)) CMP(!=, *this < a || a < *this)
  Angle operator + () const { return *this; }
  Angle operator - () const { return Angle(pivot + conj(p - pivot)); }
  OP(+, Angle(pivot + (p - pivot) * (a.p - pivot))) OP(-, *this + (-a))
};
#undef OP
#undef CMP

pt Angle::pivot = pt(0, 0);


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
#include "Circle.h"
#include "DelaunayTriangulation.h"
using namespace std;

// Computes the Voronoi Diagram and Delaunay Triangulation of a set of
//   distinct points
// If all points are collinear, there is no triangulation or diagram
// If there are 4 or more points on the same circle, the triangulation is
//   ambiguous, otherwise it is unique
// Each circumcircle does not completely contain any of the input points
// The Voronoi Diagram is represented as a graph with the first tri.size()
//   vertices being the circumcenters of each triangle in the triangulation
//   and additional vertices being added to represent infinty points in a
//   certain direction
// An edge is added between (v, w) if either triangles
//   v and w share an edge, or if v is a triangle that has an edge that is
//   not shared with any other triangle and rayDir[w - tri.size()] is the
//   direction of the ray from the circumcenter of v passing perpendicular
//   to the edge of the triangle of v that is not shared
// Constructor Arguments:
//   P: the distinct points
// Fields:
//   tri: a vector of arrays of 3 points, representing a triangle in
//     the triangulation, points given in ccw order
//   circumcircles: the circumcircles of each triangle
//   rayDir: for edges in the Voronoi Diagram that are infinity, this vector
//     contains the direction vector of those edges
//   G: the adjacency list of each vertex in the Voronoi Diagram
// In practice, has a moderate constant
// Time Complexity:
//   construction: O(N log N)
// Memory Complexity: O(N log N)
// Tested:
//   https://open.kattis.com/problems/pandapreserve
struct VoronoiDiagram {
  vector<array<pt, 3>> tri; vector<Circle> circumcircles; vector<pt> rayDir;
  vector<vector<int>> G;
  VoronoiDiagram(const vector<pt> &P)
      : tri(DelaunayTriangulation(P).tri), G(tri.size()) {
    map<pair<pt, pt>, int> seen; circumcircles.reserve(tri.size());
    for (auto &&t : tri)
      circumcircles.push_back(circumcircle(t[0], t[1], t[2]));
    for (int h = 0; h < 2; h++) for (int i = 0; i < int(tri.size()); i++)
      for (int k = 0; k < 3; k++) {
        pt p1 = tri[i][k], p2 = tri[i][(k + 1) % 3];
        auto it = seen.find(make_pair(p1, p2)); if (h == 0) {
          if (it == seen.end()) it = seen.find(make_pair(p2, p1));
          if (it == seen.end()) seen[make_pair(p1, p2)] = i;
          else {
            int j = it->second; G[i].push_back(j); G[j].push_back(i);
            seen.erase(it);
          }
        } else if (it != seen.end()) {
          int j = G.size(); G.emplace_back(); rayDir.push_back(perp(p1 - p2));
          G[i].push_back(j); G[j].push_back(i);
        }
      }
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
using namespace std;

// Affine Transformations in 2D
// Functions:
//   prependMatrix(m2, b2): sets m = m2 * m and b = m2 * b + b2
//   transform(t): applies the AffineTransformation t to this
//   scale(p): scales the x and y coordinates by p.x and p.y respectively
//   translate(p): translates the point by p
//   rotate(theta): rotates the point theta radians around the origin
//   reflect(dir): reflects the point across the line passing through the
//     origin with direction dir
//   project(dir): projects the point onto the line passing through the
//     origin with direction dir
//   applyTransform(p): applies the transformation to the point p
//   inverse(): returns the inverse of this transformation,
//     determinant of m must nonzero
// Time Complexity:
//   constructor, prependMatrix, transform, scale, translate, rotate, reflect,
//     project, applyTransform, inverse: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/tree3
struct AffineTransform {
  array<array<T, 2>, 2> m; array<T, 2> b;
  AffineTransform() {
    for (int i = 0; i < 2; i++) {
      b[i] = T(0); for (int j = 0; j < 2; j++) m[i][j] = i == j ? T(1) : T(0);
    }
  }
  void prependMatrix(const array<array<T, 2>, 2> &m2, const array<T, 2> &b2) {
    array<array<T, 2>, 2> resm; array<T, 2> resb;
    for (int i = 0; i < 2; i++) for (int j = 0; j < 2; j++) {
      resm[i][j] = T(0);
      for (int k = 0; k < 2; k++) resm[i][j] += m2[i][k] * m[k][j];
    }
    for (int i = 0; i < 2; i++) {
      resb[i] = b2[i]; for (int j = 0; j < 2; j++) resb[i] += m2[i][j] * b[j];
    }
    m = resm; b = resb;
  }
  void transform(const AffineTransform &t) { prependMatrix(t.m, t.b); }
  void scale(pt p) {
    prependMatrix({array<T, 2>{p.x, T(0)}, {T(0), p.y}}, {T(0), T(0)});
  }
  void translate(pt p) {
    prependMatrix({array<T, 2>{T(1), T(0)}, {T(0), T(1)}}, {p.x, p.y});
  }
  void rotate(T theta) {
    T c = cos(theta), s = sin(theta);
    prependMatrix({array<T, 2>{c, -s}, {s, c}}, {T(0), T(0)});
  }
  void reflect(pt dir) {
    T a = dir.x * dir.x - dir.y * dir.y, b = T(2) * dir.x * dir.y;
    T n = norm(dir); a /= n; b /= n;
    prependMatrix({array<T, 2>{a, b}, {b, -a}}, {T(0), T(0)});
  }
  void project(pt dir) {
    T a = dir.x * dir.x, b = dir.x * dir.y, c = dir.y * dir.y;
    T n = norm(dir); a /= n; b /= n; c /= n;
    prependMatrix({array<T, 2>{a, b}, {b, c}}, {T(0), T(0)});
  }
  pt applyTransform(pt p) {
    return pt(m[0][0] * p.x + m[0][1] * p.y + b[0],
              m[1][0] * p.x + m[1][1] * p.y + b[1]);
  }
  AffineTransform inverse() const {
    AffineTransform ret;
    ret.prependMatrix({array<T, 2>{T(1), T(0)}, {T(0), T(1)}}, {-b[0], -b[1]});
    T det = m[0][0] * m[1][1] - m[0][1] * m[1][0]; assert(!eq(det, 0));
    ret.prependMatrix({array<T, 2>{m[1][1] / det, -m[0][1] / det},
                                  {-m[1][0] / det, m[0][0] / det}},
                      {T(0), T(0)});
    return ret;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "Point.h"
using namespace std;

// Helper struct to maintain upper hull
struct DecrementalUpperHull {
  struct Link; using ptr = Link *;
  struct Link {
    int id; pt p; ptr prv, nxt;
    Link(int id, pt p) : id(id), p(p), prv(nullptr), nxt(nullptr) {}
  };
  struct Node {
    ptr chain, back, tangent;
    Node() : chain(nullptr), back(nullptr), tangent(nullptr) {}
  };
  int N; vector<Node> TR; vector<Link> links;
  template <class F, class G>
  pair<ptr, ptr> findBridge(ptr l, ptr r, F f, G g) {
    while (f(l) || f(r)) {
      if (!f(r) || (f(l) && g(pt(0, 0), f(l)->p - l->p, f(r)->p - r->p))) {
        if (g(l->p, f(l)->p, r->p)) l = f(l);
        else break;
      } else {
        if (!g(l->p, r->p, f(r)->p)) r = f(r);
        else break;
      }
    }
    return make_pair(l, r);
  }
  void fixChain(int x, ptr l, ptr r, bool rev) {
    if (rev) {
      tie(r, l) = findBridge(r, l, [&] (ptr q) { return q->prv; },
                             [&] (pt a, pt b, pt c) {
                               return ccw(a, b, c) >= 0;
                             });
    } else {
      tie(l, r) = findBridge(l, r, [&] (ptr q) { return q->nxt; },
                             [&] (pt a, pt b, pt c) {
                               return ccw(a, b, c) <= 0;
                             });
    }
    TR[x].tangent = l;
    TR[x].chain = TR[x * 2].chain; TR[x].back = TR[x * 2 + 1].back;
    TR[x * 2].chain = l->nxt; TR[x * 2 + 1].back = r->prv;
    if (l->nxt) l->nxt->prv = nullptr;
    else TR[x * 2].chain = nullptr;
    if (r->prv) r->prv->nxt = nullptr;
    else TR[x * 2 + 1].chain = nullptr;
    l->nxt = r; r->prv = l;
  }
  void build(int x, int tl, int tr) {
    if (tl == tr) { TR[x].chain = TR[x].back = &links[tl]; return; }
    int m = tl + (tr - tl) / 2;
    build(x * 2, tl, m); build(x * 2 + 1, m + 1, tr);
    fixChain(x, TR[x * 2].chain, TR[x * 2 + 1].chain, false);
  }
  void rob(int x, int y) {
    TR[x].chain = TR[y].chain; TR[x].back = TR[y].back;
    TR[y].chain = TR[y].back = nullptr;
  }
  void rem(int x, int tl, int tr, int i) {
    if (i < tl || tr < i) return;
    int m = tl + (tr - tl) / 2, y = x * 2 + int(i > m); if (!TR[x].tangent) {
      TR[y].chain = TR[x].chain; TR[y].back = TR[x].back;
      if (i <= m) rem(x * 2, tl, m, i);
      else rem(x * 2 + 1, m + 1, tr, i);
      rob(x, y); return;
    }
    ptr l = TR[x].tangent, r = l->nxt; l->nxt = TR[x * 2].chain;
    if (TR[x * 2].chain) TR[x * 2].chain->prv = l;
    else TR[x * 2].back = l;
    TR[x * 2].chain = TR[x].chain; r->prv = TR[x * 2 + 1].back;
    if (TR[x * 2 + 1].back) TR[x * 2 + 1].back->nxt = r;
    else TR[x * 2 + 1].chain = r;
    TR[x * 2 + 1].back = TR[x].back;
    if (TR[y].chain == TR[y].back && TR[y].chain->id == i) {
      TR[y].chain = TR[y].back = nullptr;
      rob(x, y ^ 1); TR[x].tangent = nullptr; return;
    }
    if (i <= m) {
      if (l->id == i) l = l->nxt;
      rem(x * 2, tl, m, i); if (!l) l = TR[x * 2].back;
    } else {
      if (r->id == i) r = r->prv;
      rem(x * 2 + 1, m + 1, tr, i); if (!r) r = TR[x * 2 + 1].chain;
    }
    fixChain(x, l, r, i <= m);
  }
  void rem(int i) {
    if (TR[1].chain == TR[1].back) {
      TR[1].chain = TR[1].back = nullptr; return;
    }
    rem(1, 0, N - 1, i);
  }
  DecrementalUpperHull(const vector<pt> &P)
      : N(P.size()), TR(N == 0 ? 0 : 1 << __lg(N * 4 - 1)) {
    links.reserve(N); for (int i = 0; i < N; i++) links.emplace_back(i, P[i]);
    build(1, 0, N - 1);
  }
  vector<int> getHull() {
    vector<int> ret;
    for (ptr x = TR[1].chain; x; x = x->nxt) ret.push_back(x->id);
    return ret;
  }
};

// Computes the convex layer each point is on
// Function Arguments:
//   P: a vector of points
// Return Value: a vector of integers representing the convex layer each
//   point is on, with 0 being the outermost layer, and increasing when
//   moving inwards
// In practice, has a small constant
// Time Complexity: O(N (log N)^2)
// Memory Complexity: O(N)
// Tested:
//   https://judge.yosupo.jp/problem/convex_layers
vector<int> convexLayers(const vector<pt> &P) {
  if (P.empty()) return vector<int>();
  int N = P.size(); vector<int> ind(N); iota(ind.begin(), ind.end(), 0);
  sort(ind.begin(), ind.end(), [&] (int i, int j) { return P[i] < P[j]; });
  vector<pt> tmp(N); for (int i = 0; i < N; i++) tmp[i] = P[ind[i]];
  DecrementalUpperHull up(tmp);
  for (int i = 0; i < N; i++) tmp[i] = P[ind[N - i - 1]] * T(-1);
  DecrementalUpperHull down(move(tmp)); vector<int> ret(N, -1);
  for (int layer = 0, done = 0; done < N; layer++) {
    vector<int> hull; for (int i : up.getHull()) hull.push_back(i);
    for (int i : down.getHull()) hull.push_back(N - i - 1);
    for (int i : hull) if (ret[ind[i]] == -1) {
      ret[ind[i]] = layer; done++; up.rem(i); down.rem(N - i - 1);
    }
  }
  return ret;
}