2D Geometry

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using ld = long double;
const ld EPS = 1e-9L;

using T = ld;
using pt = complex<T>;
#define x real()
#define y imag()

int sgn(T val) { return (T(0) < val) - (val < T(0)); }
T sq(pt p) { return p.x * p.x + p.y * p.y; }
pt perp(pt p) { return {-p.y, p.x}; }
T dot(pt v, pt w) { return v.x * w.x + v.y * w.y; }
T cross(pt v, pt w) { return v.x * w.y - v.y * w.x; }
T orient(pt a, pt b, pt c) { return cross(b - a, c - a); }

struct line {
    pt v; T c;
    line(pt v, T c) : v(v), c(c) {}
    line(T a, T b, T _c) : v(b, -a), c(_c) {}
    line(pt p, pt q) : v(q - p), c(cross(v, p)) {}
    T side(pt p) { return cross(v, p) - c; }
    ld dist(pt p) { return fabsl(side(p)) / abs(v); }
    pt proj(pt p) { return p - perp(v) * side(p) / sq(v); }
    pt refl(pt p) { return p - perp(v) * (T)2.0L * side(p) / sq(v); }
};

bool lineIntersect(line l1, line l2, pt &out) {
    T d = cross(l1.v, l2.v);
    if (fabsl(d) <= EPS) return false;
    out = (l2.v * l1.c - l1.v * l2.c) / d;
    return true;
}

bool onSegment(pt a, pt b, pt p) {
    return fabsl(orient(a, b, p)) <= EPS && dot(a - p, b - p) <= EPS;
}

bool properInter(pt a, pt b, pt c, pt d, pt &out) {
    T oa = orient(c, d, a), ob = orient(c, d, b);
    T oc = orient(a, b, c), od = orient(a, b, d);
    if (sgn(oa) * sgn(ob) < 0 && sgn(oc) * sgn(od) < 0) {
        out = (a * ob - b * oa) / (ob - oa);
        return true;
    }
    return false;
}

ld segPoint(pt a, pt b, pt p) {
    if (a != b && dot(p - a, b - a) >= 0 && dot(p - b, a - b) >= 0)
        return line(a, b).dist(p);
    return min(abs(p - a), abs(p - b));
}

bool isConvex(vector<pt> &p) {
    bool hasPos = false, hasNeg = false;
    int n = p.size();
    for (int i = 0; i < n; i++) {
        int o = sgn(orient(p[i], p[(i+1)%n], p[(i+2)%n]));
        if (o > 0) hasPos = true;
        if (o < 0) hasNeg = true;
    }
    return !(hasPos && hasNeg);
}

ld areaTriangle(pt a, pt b, pt c) {
    return fabsl(cross(b - a, c - a)) / 2.0L;
}

ld areaPolygon(vector<pt> &p) {
    ld area = 0;
    int n = p.size();
    for (int i = 0; i < n; i++)
        area += cross(p[i], p[(i+1)%n]);
    return fabsl(area) / 2.0L;
}

bool inPolygon(vector<pt> &p, pt a, bool strict = true) {
    int cnt = 0, n = p.size();
    for (int i = 0; i < n; i++) {
        if (onSegment(p[i], p[(i+1)%n], a)) return !strict;
        bool above = (p[(i+1)%n].y >= a.y) - (p[i].y >= a.y);
        cnt += above * orient(a, p[i], p[(i+1)%n]) > 0;
    }
    return cnt & 1;
}