计算几何

#include <bits/stdc++.h>

using namespace std;

const int N = 2e5 + 5;

#define eps 1e-6

template<typename T> inline T abs_(T a) { return a < 0 ? -a : a; }
template<typename T> inline T min_(T a, T b) { return a < b ? a : b; }
template<typename T> inline T max_(T a, T b) { return a > b ? a : b; }

struct node {
	double x, y;
	inline node(double X, double Y) : x(X), y(Y) { }
	inline node() { x = y = 0; }
	inline bool operator > (const node &b) const { return abs_(y - b.y) < eps ? x > b.x : y > b.y; }
	inline bool operator < (const node &b) const { return abs_(y - b.y) < eps ? x < b.x : y < b.y; }
	inline bool operator == (const node &b) const { return abs_(x - b.x) < eps && abs_(y - b.y) < eps; }
	inline double mo() { return x * x + y * y; }
	inline double len() { return sqrt(mo()); }
};

inline bool equ(double a, double b) { return abs_(a - b) < eps; }
inline int dcmp(double x) { return abs_(x) < eps ? 0 : (x > 0 ? 1 : -1); }
inline double dist(node a, node b) { return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y)); }
inline node vec(node a, node b) { return node(b.x - a.x, b.y - a.y); }
inline double operator * (node a, node b) { return a.x * b.x + a.y * b.y; }
inline node operator * (double k, node b) { return node(k * b.x, k * b.y); }
inline double operator ^ (node a, node b) { return a.x * b.y - a.y * b.x; }
inline node operator + (node a, node b) { return node(a.x + b.x, a.y + b.y); }
inline node operator - (node a, node b) { return node(a.x - b.x, a.y - b.y); }
inline int perp(node a, node b) { return dcmp(a * b); }
inline bool isinter(node a, node b, node p, node q) { return dcmp(((a - b) ^ (p - b)) * ((a - b) ^ (q - b))) < 0 && dcmp(((p - q) ^ (a - q)) * ((p - q) ^ (b - q))) < 0; }
inline node inter(node a, node b, node c, node d) { return a + (((a - c) ^ (d - c)) / ((d - c) ^ (b - a))) * (b - a); }
inline double dis(node a, node p, node q) { return p == q ? 0 : abs_((a - p) ^ (a - q)) / (p - q).len(); }
inline node foot(node c, node a, node b) { return a + (abs_(((c - a) * (c - b)) / (b - a).mo()) * (b - a)); }
inline node rotate(node a, double t) { return node(a.x * cos(t) - a.y * sin(t), a.y * cos(t) + a.x * sin(t)); }
inline double triarea(node a, node b, node c) { return abs_((a - b) ^ (a - c)) / 2.0; }
inline bool cmp(node a, node b) { return dcmp(a ^ b) > 0 || (equ(a ^ b, 0) && (a.mo() < b.mo())); }

#define poly vector<node>

inline double area(vector<node> A) {
	int n = A.size();
	if (n <= 2) return 0;
	double res = A[n - 1] ^ A[0];
	for (int i = 0; i < n - 1; ++i)
		res += A[i] ^ A[i + 1];
	return res / 2.0;
}

inline double length(vector<node> A) {
	int n = A.size(); double res = dist(A[0], A[n - 1]);
	for (int i = 0; i < n - 1; ++i)
		res += dist(A[i], A[i + 1]);
	return res;
}

namespace Conv {

int stk[N], top;

inline poly convex(poly a) {
	sort(a.begin(), a.end()); top = 0;
	poly res; int n = a.size(); node p = a[0];
	for (int i = 0; i < n; ++i) a[i] = a[i] - p;
	poly :: iterator tmp = a.begin();
	sort(++tmp, a.end(), cmp);
	stk[++top] = 0;
	for (int i = 1; i < n; ++i) {
		while (top > 1 && ((a[i] - a[stk[top - 1]]) ^ (a[stk[top]] - a[stk[top - 1]])) >= 0) --top;
		stk[++top] = i;
	}
	for (int i = 1; i <= top; ++i) res.push_back(a[stk[i]] + p);
	return res;
}

inline poly minkowski(poly a, poly b) {
	int n = a.size(), m = b.size();
	poly A, B;
	for (int i = 0; i < n; ++i) A.push_back(vec(a[i], a[(i + 1) % n]));
	for (int i = 0; i < m; ++i) B.push_back(vec(b[i], b[(i + 1) % m]));
	poly res; int i = 0, j = 0, cnt = 0; res.push_back(a[0] + b[0]);
	while (i < n && j < m) {
		res.push_back(res[cnt] + (dcmp(A[i] ^ B[j]) >= 0 ? A[i++] : B[j++]));
		++cnt;
	}
	while (i < n) {
		res.push_back(res[cnt] + A[i++]);
		++cnt;
	}
	while (j < m) {
		res.push_back(res[cnt] + B[j++]);
		++cnt;
	} return convex(res);
}

}

namespace Cartr {

poly cv;

inline double cartr(poly a) {
	cv = Conv :: convex(a); int n = cv.size();
	double ans = (cv[1] - cv[0]).mo(); int t = 0;
	for (int i = 0; i < n; ++i) {
		while (dis(cv[t], cv[i], cv[(i + 1) % n]) < dis(cv[(t + 1) % n], cv[i], cv[(i + 1) % n]))
			t = (t + 1) % n;
		ans = max_(ans, max_((cv[i] - cv[t]).mo(), (cv[(i + 1) % n] - cv[t]).mo()));
	} return ans;
}

}

node t;

inline bool checkin(poly &a, node x) {
	t = a[0];
	if (dcmp((x - t) ^ (a[0] - t)) > 0 || dcmp((*a.rbegin() - t) ^ (x - t)) > 0) return 0;
	poly :: iterator p = lower_bound(a.begin(), a.end(), x, [](node a, node b){return cmp(a - t, b - t);});
	--p; node i = *p, j; ++p;
	if (p == a.end()) j = a[0];
	else j = *p;
	return dcmp((x - i) ^ (j - i)) <= 0;
}

struct line {
	node a, b;
	inline line() { }
	inline line(double x1, double y1, double x2, double y2) : a(node(x1, y1)), b(node(x2, y2)) { }
	inline line(node A, node B) : a(A), b(B) { }
	inline bool operator < (const line &bb) {
		double t1 = atan2(b.y - a.y, b.x - a.x), t2 = atan2(bb.b.y - bb.a.y, bb.b.x - bb.a.x);
		if (abs_(t1 - t2) > eps) return t1 < t2;
		return ((bb.a - a) ^ (bb.b - a)) > eps;
	}
	inline bool operator > (const line &bb) {
		double t1 = atan2(b.y - a.y, b.x - a.x), t2 = atan2(bb.b.y - bb.a.y, bb.b.x - bb.a.x);
		if (abs_(t1 - t2) > eps) return t1 > t2;
		return ((bb.a - a) ^ (bb.b - a)) < -eps;
	}
	inline bool operator == (const line &bb) {
		return a == bb.a && b == bb.b;
	}
	inline double angle() {
		return atan2(b.y - a.y, b.x - a.x);
	}
};

inline node inter(line a, line b) {
	return inter(a.a, a.b, b.a, b.b);
}

namespace HP {
	
line q[N];
int l, r;
node p[N];
	
inline poly polyAnd(vector<line> &s) {
	int n = s.size(); sort(s.begin(), s.end());
	l = 1, r = 0;
	for (int i = 0; i < n; ++i) {
		if (i && equ(s[i].angle(), s[i - 1].angle())) continue;
		while (r > l && dcmp((s[i].b - p[r]) ^ (s[i].a - p[r])) > 0) --r;
		while (r > l && dcmp((s[i].b - p[l + 1]) ^ (s[i].a - p[l + 1])) > 0) ++l;
		q[++r] = s[i];
		if (r > l) p[r] = inter(q[r], q[r - 1]);
	}
	while (r > l && ((q[l].b - p[r]) ^ (q[l].a - p[r])) > eps) --r;
	p[r + 1] = inter(q[l], q[r]);
	++r; poly res;
	for (int i = l + 1; i <= r; ++i) res.push_back(p[i]);
	return res;
}
}

const double PI = acos(-1); 

struct circ {
	node o;
	double r;
	inline circ() { }
	inline circ(node O, double R) : o(O), r(R) { }
	inline circ(node a, node b, node c) {
		double A1 = 2 * (b.x - a.x), B1 = 2 * (b.y - a.y), C1 = b.mo() - a.mo();
		double A2 = 2 * (c.x - b.x), B2 = 2 * (c.y - b.y), C2 = c.mo() - b.mo();
		o = node((C1 * B2 - C2 * B1) / (A1 * B2 - A2 * B1), (A1 * C2 - A2 * C1) / (A1 * B2 - A2 * B1));
		r = (a - o).len();
	}
	inline bool cover(const circ &p) { return (o - p.o).len() <= r - p.r; }
	inline bool cover(node p) { return dcmp(r - (p - o).len()) >= 0; }
	inline bool area() { return PI * r * r; }
};

inline circ minCover(poly a) {
	mt19937 rnd(time(0));
	shuffle(a.begin(), a.end(), rnd);
	circ res(node(0, 0), 0); int n = a.size();
	for (int i = 0; i < n; ++i) {
		if (!res.cover(a[i])) {
			res = circ(a[i], 0);
			for (int j = 0; j < i; ++j) {
				if (!res.cover(a[j])) {
					res = circ(0.5 * (a[i] + a[j]), (a[i] - a[j]).len() / 2);
					for (int k = 0; k < j; ++k) {
						if (!res.cover(a[k]))
							res = circ(a[i], a[j], a[k]);
					}
				}
			}
		}
	} return res;
}

李超树：

template<typename T> inline T abs_(T a) { return a < 0 ? -a : a; }
template<typename T> inline T min_(T a, T b) { return a < b ? a : b; }
template<typename T> inline T max_(T a, T b) { return a > b ? a : b; }

struct line {
	int x1, y1, x2, y2;
	
	inline double f(int x) {
		if (x1 == x2) return max_(y1, y2);
		return 1.0 * (y2 - y1) / (x2 - x1) * (x - x1) + y1;
	}
} p[100005];

int lc[N << 1 | 1], rc[N << 1 | 1], id[N << 1 | 1], ndnum = 0, rt = 0;

inline void insert(int &now, int l, int r, int x) {
	if (!now) now = ++ndnum;
	if (p[x].x1 <= l && r <= p[x].x2) {
		if (!id[now]) {
			id[now] = x;
			return ;
		}
		if (p[id[now]].f(l) >= p[x].f(l) && p[id[now]].f(r) >= p[x].f(r)) return ;
		if (p[id[now]].f(l) < p[x].f(l) && p[id[now]].f(r) < p[x].f(r)) { id[now] = x; return ; }
		int mid = l + r >> 1; if (p[id[now]].f(mid) < p[x].f(mid)) swap(id[now], x);
		if (p[id[now]].f(l) <= p[x].f(l)) insert(lc[now], l, mid, x);
		else insert(rc[now], mid + 1, r, x);
		return ;
	} int mid = l + r >> 1;
	if (p[x].x1 <= mid) insert(lc[now], l, mid, x);
	if (p[x].x2 > mid) insert(rc[now], mid + 1, r, x);
}

inline bool equ(double a, double b) {
	return abs_(a - b) < 1e-6;
}

inline int query(int now, int l, int r, int q) {
	if (!now) return 0;
	if (l >= r) return id[now];
	int res = 0, mid = l + r >> 1;
	if (q <= mid) res = query(lc[now], l, mid, q);
	if (q > mid) res = query(rc[now], mid + 1, r, q);
	if (!id[now]) return res;
	if (!res) return id[now];
	double a = p[id[now]].f(q), b = p[res].f(q);
	if (equ(a, b)) return min_(res, id[now]);
	return a < b ? res : id[now];
}