二维计算几何模板

#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <set>
#include <map>
#include <vector>
#include <queue>
#include <ctime>

using namespace std;

#define LL long long
#define inf 0x3f3f3f3f
#define MP make_pair
#define Pi acos(-1.0)
#define eps 1e-8
#define N 410
#define M 200020

#define sqr(x) ((x) * (x))

double torad(double deg) {
    return deg / 180 * Pi;
}
//经纬度（角度）转化为空间坐标( lat 纬度，lng 经度 )
void get_coord(double r, double lat, double lng, double &x, double &y, double &z) {
    lat = torad(lat);
    lng = torad(lng);
    x = r * cos(lat) * cos(lng);
    y = r * cos(lat) * sin(lng);
    z = r * sin(lat);
}
int dcmp(double x) {
    if (fabs(x) < eps) return 0;
    return x < 0 ? -1 : 1;
}
struct point {
    double x, y;
    point(double x = 0, double y = 0) : x(x), y(y) {}
    point operator - (const point &b) const {
        return point(x - b.x, y - b.y);
    }
    point operator + (const point &b) const {
        return point(x + b.x, y + b.y);
    }
    point operator * (const double &k) const {
        return point(x * k, y * k);
    }
    point operator / (const double &k) const {
        return point(x / k, y / k);
    }
    bool operator == (const point &b) const {
        return dcmp(x - b.x) == 0 && dcmp(y - b.y) == 0;
    }
    double ang() {
        return atan2(y, x);
    }
    double len() {
        return sqrt(x * x + y * y);
    }
    void input() {
        scanf("%lf%lf", &x, &y);
    }
};
double dot(point a, point b) {
    return a.x * b.x + a.y * b.y;
}
double cross(point a, point b) {
    return a.x * b.y - a.y * b.x;
}
// 向量旋转
point rotate(point a, double rad) {
    return point(a.x * cos(rad) - a.y * sin(rad), a.x * sin(rad) + a.y * cos(rad));
}
// 向量的单位法线
point normal(point a) {
    double L = a.len();
    return point(-a.y / L, a.x / L);
}
// 点到直线的距离
double point_to_line(point p, point a, point b) {
    point v1 = b - a, v2 = p - a;
    return fabs(cross(v1, v2)) / v1.len();
}
// 点到线段的距离
double point_to_segment(point p, point a, point b) {
    if (a == b) return (p - a).len();
    point v1 = b - a, v2 = p - a, v3 = p - b;
    if (dcmp(dot(v1, v2)) < 0) return v2.len();
    if (dcmp(dot(v1, v3)) > 0) return v3.len();
    return fabs(cross(v1, v2)) / v1.len();
}
// 点到直线上的投影
point get_line_projection(point p, point a, point b) {
    point v = b - a;
    return a + v * (dot(v, p - a) / dot(v, v));
}
// 判断点是是否在一条线段上
bool on_segment(point p, point a1, point a2) {
    if (p == a1 || p == a2) return 1;    //端点
    return dcmp(cross(a1 - p, a2 - p)) == 0 && dcmp(dot(a1 - p, a2 - p)) < 0;
}
// 两直线平行
bool parallel(point a1, point a2, point b1, point b2) {
    return dcmp(cross(a2 - a1, b2 - b1)) == 0;
}
// 线段相交
bool segment_intersection(point a1, point a2, point b1, point b2) {
    if (on_segment(a1, b1, b2) || on_segment(a2, b1, b2)) return 1;  //端点
    if (on_segment(b1, a1, a2) || on_segment(b2, a1, a2)) return 1;  //端点
    double c1 = cross(a2 - a1, b1 - a1), c2 = cross(a2 - a1, b2 - a1);
    double c3 = cross(b2 - b1, a1 - b1), c4 = cross(b2 - b1, a2 - b1);
    return dcmp(c1) * dcmp(c2) < 0 && dcmp(c3) * dcmp(c4) < 0;
}
// 多边形的有向面积
double polygon_area(point *p, int n) {
    double area = 0;
    for (int i = 1; i < n - 1; ++i)
        area += cross(p[i] - p[0], p[i + 1] - p[0]);
    return area / 2;
}
// 点在多边形内判定
int point_in_polygon(point p, point *poly, int n) {
    int wn = 0;
    poly[n] = poly[0];
    for (int i = 0; i < n; ++i) {
        if (on_segment(p, poly[i], poly[i + 1])) return -1; //在边界上
        int k = dcmp(cross(poly[i + 1] - poly[i], p - poly[i]));
        int d1 = dcmp(poly[i].y - p.y);
        int d2 = dcmp(poly[i + 1].y - p.y);
        if (k > 0 && d1 <= 0 && d2 > 0) wn++;
        if (k < 0 && d2 <= 0 && d1 > 0) wn--;
    }
    if (wn != 0) return 1;   //内部
    return 0;    //外部
}
// 多边形重心
point masscenter(point *p, int n) {
    point ret = point(0, 0);
    double sum = polygon_area(p, n);
    if (dcmp(sum) == 0) return ret;
    p[n] = p[0];
    for (int i = 0; i < n; ++i)
        ret = ret + (p[i] + p[i + 1]) * cross(p[i + 1], p[i]);
    return ret / sum / 6.0;
}
struct Line {
    point p, v;
    double ang;
    Line() {}
    Line(point p, point v) : p(p), v(v) {
        ang = atan2(v.y, v.x);
    }
    bool operator < (const Line &b) const {
        return ang < b.ang;
    }
};
// 点p在有向直线L的左边（线上不算）
bool onleft(Line L, point p) {
    return cross(L.v, p - L.p) > 0;
}
// 二直线交点
point get_line_intersection(Line a, Line b) {
    point u = a.p - b.p;
    double t = cross(b.v, u) / cross(a.v, b.v);
    return a.p + a.v * t;
}
int half_plane_intersection(Line *L, int n, point *poly) {
    sort(L, L + n);
    int first, last;
    point *p = new point[n];
    Line *q = new Line[n];
    q[first = last = 0] = L[0];
    for (int i = 1; i < n; ++i) {
        while (first < last && !onleft(L[i], p[last - 1])) last--;
        while (first < last && !onleft(L[i], p[first])) first++;
        q[++last] = L[i];
        if (dcmp(cross(q[last].v, q[last - 1].v)) == 0) {
            last--;
            if (onleft(q[last], L[i].p)) q[last] = L[i];
        }
        if (first < last)
            p[last - 1] = get_line_intersection(q[last - 1], q[last]);
    }
    while (first < last && !onleft(q[first], p[last - 1])) last--;
    if (last - first <= 1) return 0;
    p[last] = get_line_intersection(q[last], q[first]);
    int m = 0;
    for (int i = first; i <= last; ++i) poly[m++] = p[i];
    return m;
}
// 旋转卡壳求最远点对
double rotate_calipers(point *ch, int n) {
    int j = 1;
    double ret = 0;
    ch[n] = ch[0];
    for (int i = 0; i < n; ++i) {
        while (fabs(cross(ch[i + 1] - ch[i], ch[j + 1] - ch[i])) > fabs(cross(ch[i + 1] - ch[i], ch[j] - ch[i])))
            j = (j + 1) % n;
        ret = max(ret, max((ch[i] - ch[j]).len(), (ch[i + 1] - ch[j]).len()));
    }
    return ret;
}
// 求平面点集最小的三角形(先按照x轴排序)
point p[N], q[N];
bool cmpY(const point &a, const point &b) {
    return a.y < b.y || (a.y == b.y && a.x < b.x);
}
double solve(int L, int R) {
    if (R - L + 1 <= 6) {
        double ret = 1e20;
        for (int i = L; i <= R; ++i) {
            for (int j = i + 1; j <= R; ++j) {
                double len1 = (p[i] - p[j]).len(), len2;
                if (len1 > ret / 2) continue;
                for (int k = j + 1; k <= R; ++k) {
                    len2 = (p[i] - p[k]).len() + (p[j] - p[k]).len();
                    ret = min(ret, len1 + len2);
                }
            }
        }
        return ret;
    }
    int m = (L + R) >> 1, cnt = 0;
    double d = min(solve(L, m), solve(m + 1, R));
    for (int i = L; i <= R; ++i)
        if (fabs(p[m].x - p[i].x) <= d)
            q[cnt++] = p[i];
    sort(q, q + cnt, cmpY);
    for (int i = 0; i < cnt; ++i) {
        for (int j = i + 1; j < cnt && j < i + 7; ++j) {
            double len1 = (q[i] - q[j]).len(), len2;
            if (len1 > d / 2) continue;
            for (int k = j + 1; k < cnt && k < j + 7; ++k) {
                len2 = (q[i] - q[k]).len() + (q[j] - q[k]).len();
                d = min(d, len1 + len2);
            }
        }
    }
    return d;
}
// 平面最近点对(先按照x轴排序)
double solve(int L, int R) {
    if (R - L <= 3) {
        double ret = 1e20;
        for (int i = L; i <= R; ++i)
            for (int j = i + 1; j <= R; ++j)
                ret = min(ret, (p[i] - p[j]).len());
        return ret;
    }
    int m = (L + R) >> 1, cnt = 0;
    double ret = min(solve(L, m), solve(m + 1, R));
    for (int i = L; i <= R; ++i)
        if (fabs(p[m].x - p[i].x) <= ret)
            q[cnt++] = p[i];
    sort(q, q + cnt, cmpY);
    for (int i = 0; i < cnt; ++i)
        for (int j = i + 1; j < cnt && (q[j].y - q[i].y) <= ret; ++j)
            ret = min(ret, (q[j] - q[i]).len());
    return ret;
}
struct circle {
    point c;
    double r;
    circle() {}
    circle(point c, double r) : c(c), r(r) {}
    point getPoint(double a) {
        return point(c.x + cos(a) * r, c.y + sin(a) * r);
    }
};

//点与圆的切点(点在圆外)
int get_tangents(point p, circle c, point &a, point &b) {
    point u = p - c.c;
    double len = u.len();
    double ang = acos(c.r / len);
    double bas = atan2(u.y, u.x);
    a = c.getPoint(bas + ang);
    b = c.getPoint(bas - ang);
    return 2;
}
//圆的公切线 返回切线条数
int get_tangents(circle A, circle B, point *a, point *b) {
    int cnt = 0;
    if (A.r < B.r) swap(A, B), swap(a, b);
    double d2 = (A.c.x - B.c.x) * (A.c.x - B.c.x) + (A.c.y - B.c.y) * (A.c.y - B.c.y);
    double rcha = A.r - B.r;
    double rsum = A.r + B.r;
    if (dcmp(d2 - rcha * rcha) < 0) return 0; //内含
    double bas = atan2(B.c.y - A.c.y, B.c.x - A.c.x);
    if (dcmp(d2) == 0 && dcmp(A.r - B.r) == 0) return -1; //无数条切线
    if (dcmp(d2 - rcha * rcha) == 0) {                //内含，一条切线
        a[cnt] = A.getPoint(bas);
        b[cnt++] = B.getPoint(bas);
        return 1;
    }
    double ang = acos((A.r - B.r) / sqrt(d2));
    a[cnt] = A.getPoint(bas + ang); b[cnt++] = B.getPoint(bas + ang);
    a[cnt] = A.getPoint(bas - ang); b[cnt++] = B.getPoint(bas - ang);
    if (dcmp(d2 - rsum * rsum) == 0) {
        a[cnt] = A.getPoint(bas); b[cnt++] = B.getPoint(bas + Pi);
    }
    else if (dcmp(d2 - rsum * rsum) > 0) {
        double ang = acos((A.r + B.r) / sqrt(d2));
        a[cnt] = A.getPoint(bas + ang); b[cnt++] = B.getPoint(Pi + bas + ang);
        a[cnt] = A.getPoint(bas - ang); b[cnt++] = B.getPoint(Pi + bas - ang);
    }
    return cnt;
}
//圆与线段的交点
void circle_cross_line(point a, point b, point o, double r, point *ret, int &num) {
    double x0 = o.x, y0 = o.y, x1 = a.x, y1 = a.y, x2 = b.x, y2 = b.y;
    double dx = x2 - x1, dy = y2 - y1;
    double A = dx * dx + dy * dy, B = 2 * dx * (x1 - x0) + 2 * dy * (y1 - y0);
    double C = (x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0) - r * r;
    double delta = B * B - 4 * A * C;
    num = 0;
    if (dcmp(delta) >= 0) {
        double t1 = (-B - sqrt(delta)) / (2 * A);
        double t2 = (-B + sqrt(delta)) / (2 * A);
        if (dcmp(t1 - 1) <= 0 && dcmp(t1) >= 0)
            ret[num++] = point(x1 + t1 * dx, y1 + t1 * dy);
        if (dcmp(t2 - 1) <= 0 && dcmp(t2) >= 0)
            ret[num++] = point(x1 + t2 * dx, y1 + t2 * dy);
    }
}
// 两圆相交的面积
double cal_area_sec(double ra, double d, double rb) {
    double cosh = (ra * ra + d - rb * rb) / (ra * sqrt(d) * 2);
    double sinh = sqrt(1 - cosh * cosh);
    double area_san = ra * (cosh * ra) * sinh;
    double ang = acos(cosh);
    double area_sec = ra * ra * ang;
    return area_sec - area_san;
}
double circle_cross(point a, double ra, point b, double rb) {
    point p = b - a;
    double d = p.x * p.x + p.y * p.y;
    if (dcmp(sqrt(d) - ra - rb) >= 0) return 0;      //相离
    if (dcmp(fabs(ra - rb) - sqrt(d)) >= 0 || dcmp(d) == 0)  //内含或相切
        return Pi * min(ra, rb) * min(ra, rb);
    return cal_area_sec(ra, d, rb) + cal_area_sec(rb, d, ra);
}
// 两圆相交求交点
bool circle_cross(point a, double ra, point b, double rb, point &v1, point &v2) {
    double d = (a - b).len();
    if (dcmp(d - ra - rb) >= 0 || dcmp(fabs(ra - rb) - d) >= 0) return 1;
    double da = (ra * ra + d * d - rb * rb) / (2 * ra * d);
    double aa = atan2((b - a).y, (b - a).x);
    double rad = acos(da);
    v1 = point(a.x + cos(aa - rad) * ra, a.y + sin(aa - rad) * ra);
    v2 = point(a.x + cos(aa + rad) * ra, a.y + sin(aa + rad) * ra);
    return 0;
}
// 圆心在原点，圆上两点a，b之间的扇形面积
double sector_area(point a, point b, double r) {
    double theta = atan2(a.y, a.x) - atan2(b.y, b.x);
    while (theta <= 0) theta += 2 * Pi;
    while (theta > 2 * Pi) theta -= 2 * Pi;
    theta = min(theta, 2 * Pi - theta);
    return r * r * theta / 2;
}
// 圆和多边形交的面积
double cal(point a, point b, double r) {
    point p[3];
    int num = 0;
    bool ina = (a.len() - r) < 0, inb = (b.len() - r) < 0;
    if (ina) {
        if (inb)
            return fabs(cross(a, b)) / 2.0;
        else {
            circle_cross_line(a, b, point(0, 0), r, p, num);
            return sector_area(b, p[0], r) + fabs(cross(a, p[0])) / 2.0;
        }
    }
    else {
        if (inb) {
            circle_cross_line(a, b, point(0, 0), r, p, num);
            return sector_area(p[0], a, r) + fabs(cross(p[0], b)) / 2.0;
        }
        else {
            circle_cross_line(a, b, point(0, 0), r, p, num);
            if (num == 2)
                return sector_area(a, p[0], r) + sector_area(p[1], b, r) + fabs(cross(p[0], p[1])) / 2.0;
            else
                return sector_area(a, b, r);
        }
    }
}
double circle_poly_area(point *res, int n, double r) { //圆心在原点,半径为r
    double ret = 0;
    res[n] = res[0];
    for (int i = 0; i < n; ++i) {
        int sgn = dcmp(cross(res[i], res[i + 1]));
        if (sgn != 0)
            ret += sgn * cal(res[i], res[i + 1], r);
    }
    return ret;
}
// 三点求圆心
void circle_center(point p0, point p1, point p2, point &cp) {
    double a1 = p1.x - p0.x, b1 = p1.y - p0.y, c1 = (a1 * a1 + b1 * b1) / 2;
    double a2 = p2.x - p0.x, b2 = p2.y - p0.y, c2 = (a2 * a2 + b2 * b2) / 2;
    double d = a1 * b2 - a2 * b1;
    cp.x = p0.x + (c1 * b2 - c2 * b1) / d;
    cp.y = p0.y + (a1 * c2 - a2 * c1) / d;
}
// 两点求圆心
void circle_center(point p0, point p1, point &cp) {
    cp.x = (p0.x + p1.x) / 2;
    cp.y = (p0.y + p1.y) / 2;
}
// 点是否在圆内（包括边界）
bool point_in_circle(point p, point cp, double r) {
    return dcmp((p - cp).len() - r) <= 0;
}
// 最小圆覆盖
void min_circle_cover(point *a, int n, point &cp, double &r) {
    r = 0, cp = a[0];
    for (int i = 1; i < n; ++i)
        if (!point_in_circle(a[i], cp, r)) {
            cp = a[i], r = 0;
            for (int j = 0; j < i; ++j)
                if (!point_in_circle(a[j], cp, r)) {
                    circle_center(a[i], a[j], cp);
                    r = (a[j] - cp).len();
                    for (int k = 0; k < j; ++k)
                        if (!point_in_circle(a[k], cp, r)) {
                            circle_center(a[i], a[j], a[k], cp);
                            r = (a[k] - cp).len();
                        }
                }
        }
}


// 圆面积并
struct arc {
    point u, v;
    double st, ed;
    arc(point u = point(0, 0), point v = point(0, 0), double st = 0, double ed = 0) : u(u), v(v), st(st), ed(ed) {}
    bool operator < (const arc &b) const {
        return dcmp(st - b.st) < 0 || (dcmp(st - b.st) == 0 && dcmp(ed - b.ed) < 0);
    }
}A[N << 1];
struct circle {
    point p;
    double r;
    bool operator < (const circle &b) const {
        return r > b.r;
    }
    void input() {
        p.input(); scanf("%lf", &r);
    }
}cir[N], cirt[N];
int n, m;
void prepare() {
    for (int i = 1; i <= n; ++i) cir[i].input();
    sort(cir + 1, cir + n + 1);
    m = 0;
    for (int i = 1; i <= n; ++i) {
        if (dcmp(cir[i].r) == 0) break;
        bool ok = 0;
        for (int j = 0; j < m; ++j) {
            if (dcmp((cir[i].p - cirt[j].p).len() - (cirt[j].r - cir[i].r)) <= 0) {
                ok = 1; break;
            }
        }
        if (ok) continue;
        cirt[m++] = cir[i];
    }
    n = m;
    memcpy(cir, cirt, sizeof cirt);
}
void solve() {
    double ans = 0;
    for (int i = 0; i < n; ++i) {
        m = 0;
        for (int j = 0; j < n; ++j) {
            if (i == j) continue;
            if ((cir[i].p - cir[j].p).len() >= (cir[i].r + cir[j].r)) continue;
            point v1, v2;
            circle_cross(cir[i].p, cir[i].r, cir[j].p, cir[j].r, v1, v2);
            double st = atan2(v1.y - cir[i].p.y, v1.x - cir[i].p.x), ed = atan2(v2.y - cir[i].p.y, v2.x - cir[i].p.x);
            if (st <= ed) {
                A[m++] = arc(v1, v2, st, ed);
            }
            else {
                A[m++] = arc(v1, point(cir[i].p.x - cir[i].r, cir[i].p.y), st, Pi);
                A[m++] = arc(point(cir[i].p.x - cir[i].r, cir[i].p.y), v2, -Pi, ed);
            }
        }
        sort(A, A + m);
        if (m > 0) {
            double L = A[0].st, R = A[0].ed;
            point p = A[0].v, q;
            for (int j = 1; j < m; ++j) {
                if (A[j].st > R) {
                    double ang = A[j].st - R;
                    ans += 0.5 * cir[i].r * cir[i].r * (ang - sin(ang));
                    q = A[j].u;
                    ans += cross(p, q) / 2;
                    L = A[j].st, R = A[j].ed, p = A[j].v;

                }
                else if (A[j].ed > R) {
                    R = A[j].ed, p = A[j].v;
                }
            }
            double ang = 2 * Pi - R + A[0].st;
            ans += 0.5 * cir[i].r * cir[i].r * (ang - sin(ang));
            q = A[0].u;
            ans += cross(p, q) / 2;
        }
        else ans += Pi * cir[i].r * cir[i].r;
    }
    printf("%.3lf\n", ans);
}
// 圆的k次面积并
struct Circle {
    point c;
    double r, ang;
    int d;
    Circle() {}
    Circle(point c, double ang = 0, int d = 0) :c(c), ang(ang), d(d) {}
    void input() {
        scanf("%lf%lf", &c.x, &c.y); d = 1;
        scanf("%lf", &r);
    }
}cir[N], tp[N << 1];
bool circmp(const Circle &a, const Circle &b) {
    return dcmp(a.r - b.r) < 0;
}
bool cmp(const Circle &a, const Circle &b) {
    return dcmp(a.ang - b.ang) < 0 || (dcmp(a.ang - b.ang) == 0 && a.d > b.d);
}
double calc(Circle o, Circle a, Circle b) {
    double ans = (b.ang - a.ang) * sqr(o.r)
        - cross(a.c - o.c, b.c - o.c) + cross(a.c - point(0, 0), b.c - point(0, 0));
    return ans / 2;
}
double area[N];
void CirUnion(Circle cir[], int n) {
    Circle res1, res2;
    sort(cir, cir + n, circmp);
    for (int i = 0; i < n; ++i)
        for (int j = i + 1; j < n; ++j)
            if (dcmp((cir[i].c - cir[j].c).len() + cir[i].r - cir[j].r) <= 0)
                cir[i].d++;
    for (int i = 0; i < n; ++i) {
        int tn = 0, cnt = 0;
        for (int j = 0; j < n; ++j) {
            if (i == j) continue;
            if (circle_cross(cir[i].c, cir[i].r, cir[j].c, cir[j].r, res1.c, res2.c))
                continue;
            res1.ang = atan2(res1.c.y - cir[i].c.y, res1.c.x - cir[i].c.x);
            res2.ang = atan2(res2.c.y - cir[i].c.y, res2.c.x - cir[i].c.x);
            res1.d = 1, tp[tn++] = res1;
            res2.d = -1, tp[tn++] = res2;
            if (dcmp(res1.ang - res2.ang) > 0) cnt++;
        }
        tp[tn++] = Circle(point(cir[i].c.x - cir[i].r, cir[i].c.y), Pi, -cnt);
        tp[tn++] = Circle(point(cir[i].c.x - cir[i].r, cir[i].c.y), -Pi, cnt);
        sort(tp, tp + tn, cmp);
        int p, s = cir[i].d + tp[0].d;
        for (int j = 1; j < tn; ++j) {
            p = s; s += tp[j].d;
            area[p] += calc(cir[i], tp[j - 1], tp[j]);
        }
    }
}
/*
void solve() {
    int n;
    scanf("%d", &n);
    for (int i = 0; i < n; ++i)
        cir[i].input();
    for (int i = 0; i <= n; ++i)
        area[i] = 0;
    CirUnion(cir, n);
    for (int i = 1; i <= n; ++i)
        area[i] -= area[i + 1];
    for (int i = 1; i <= n; ++i)
        printf("[%d] = %.3lf\n", i, area[i]);
}
*/
// 多边形面积并
struct polygon {
    point p[N];    //N比较小(N<=50)
    int sz;
    void input() {
        for (int i = 0; i < sz; ++i)
            p[i].input();
        if (dcmp(polygon_area(p, sz)) < 0)
            reverse(p, p + sz);
        p[sz] = p[0];
    }
}g[5];
pair<double, int> C[100020];
double segP(point a, point b, point c) {
    if (dcmp(b.x - c.x))
        return (a.x - b.x) / (c.x - b.x);
    return (a.y - b.y) / (c.y - b.y);
}

double polyUnion(int n) {   //n是多边形的数目
    double sum = 0;
    for (int i = 0; i < n; ++i)
        for (int ii = 0; ii < g[i].sz; ++ii) {
            int tot = 0;
            C[tot++] = MP(0, 0);
            C[tot++] = MP(1, 0);
            for (int j = 0; j < n; ++j) if (i != j)
                for (int jj = 0; jj < g[j].sz; ++jj) {
                    int d1 = dcmp(cross(g[i].p[ii + 1] - g[i].p[ii], g[j].p[jj] - g[i].p[ii]));
                    int d2 = dcmp(cross(g[i].p[ii + 1] - g[i].p[ii], g[j].p[jj + 1] - g[i].p[ii]));
                    if (!d1 && !d2) {
                        point t1 = g[j].p[jj + 1] - g[j].p[jj];
                        point t2 = g[i].p[ii + 1] - g[i].p[ii];
                        if (dcmp(dot(t1, t2)) > 0 && j < i) {
                            C[tot++] = MP(segP(g[j].p[jj], g[i].p[ii], g[i].p[ii + 1]), 1);
                            C[tot++] = MP(segP(g[j].p[jj + 1], g[i].p[ii], g[i].p[ii + 1]), -1);
                        }
                    }
                    else if (d1 >= 0 && d2 < 0 || d1 < 0 && d2 >= 0) {
                        double d3 = cross(g[j].p[jj + 1] - g[j].p[jj], g[i].p[ii] - g[j].p[jj]);
                        double d4 = cross(g[j].p[jj + 1] - g[j].p[jj], g[i].p[ii + 1] - g[j].p[jj]);
                        if (d2 < 0)
                            C[tot++] = MP(d3 / (d3 - d4), 1);
                        else C[tot++] = MP(d3 / (d3 - d4), -1);
                    }
                }
            sort(C, C + tot);
            double cur = min(max(C[0].first, 0.0), 1.0);
            int sgn = C[0].second;
            double s = 0;
            for (int j = 1; j < tot; ++j) {
                double nxt = min(max(C[j].first, 0.0), 1.0);
                if (!sgn) s += nxt - cur;
                sgn += C[j].second;
                cur = nxt;
            }
            sum += cross(g[i].p[ii], g[i].p[ii + 1]) * s;
        }
    return fabs(sum) / 2;
}

// 判断射线与线段是否相交（s1，e1射线；s2，e2线段）
bool intersect(point s1, point e1, point s2, point e2) {
    double a, t1, t2;
    a = cross(s2 - s1, s2 - e1) - cross(e2 - s1, e2 - e1);
    if (fabs(a) < eps) return false;
    t1 = cross(s2 - s1, s2 - e1) / a;
    t2 = cross(s2 - s1, s2 - e2) / a;
    return (t1 > -eps && t1 < 1 + eps && t2 > -eps);
}

/*
//三角剖分
struct point {
    double x, y;
    int id;
    struct Edge *e;
    bool operator < (const point & p) const {
        return dcmp(x - p.x) != 0 ? x < p.x : dcmp(y - p.y) < 0;
    }
    bool operator == (const point & p) const {
        return dcmp(x - p.x) == 0 && dcmp(y - p.y) == 0;
    }
    void input(int i) {
        id = i;
        scanf("%lf%lf", &x, &y);
    }
}pnt[N];
double cross(point & o, point & a, point & b) {
    return (a.x - o.x) * (b.y - o.y) - (b.x - o.x) * (a.y - o.y);
}
double dot(point & o, point & a, point & b) {
    return (a.x - o.x) * (b.x - o.x) + (a.y - o.y) * (b.y - o.y);
}
double dis(point a, point b) {
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
struct Edge {
    point *o, *d;
    Edge *on, *op, *dn, *dp;
};
#define Op(e,p) ((e)->o==p?(e)->d:(e)->o)
#define Next(e,p) ((e)->o==p?(e)->on:(e)->dn)
#define Prev(e,p) ((e)->o==p?(e)->op:(e)->dp)
struct Delaunay {
    void solve(point * ps, int n) { //点集需要 sort 和 unique
        sort(ps, ps + n);
        edge_num = 0;
        rubb = NULL;
        for (int i = 0; i < n; i++) ps[i].e = NULL;
        Edge* l_cw, *r_ccw;
        divide(ps, 0, n, l_cw, r_ccw);
    }
    Edge es[M], *rubb;
    int edge_num;
    Edge *make_edge(point &u, point &v) {
        Edge * e;
        if (rubb == NULL) {
            e = es + edge_num++;
        }
        else {
            e = rubb;
            rubb = rubb->dn;
        }
        e->on = e->op = e->dn = e->dp = e;
        e->o = &u; e->d = &v;
        if (u.e == NULL) u.e = e;
        if (v.e == NULL) v.e = e;
        return e;
    }
    void delete_edge(Edge *e) {
        point *u = e->o, *v = e->d;
        if (u->e == e) u->e = e->on;
        if (v->e == e) v->e = e->dn;
        Prev(e->on, u) = e->op;
        Next(e->op, u) = e->on;
        Prev(e->dn, v) = e->dp;
        Next(e->dp, v) = e->dn;
        e->dn = rubb;
        rubb = e;
    }
    void splice(Edge *a, Edge *b, point *v) {
        Edge *n;
        n = Next(a, v); Next(a, v) = b;
        Prev(n, v) = b;
        Next(b, v) = n; Prev(b, v) = a;
    }
    Edge *join(Edge *a, point *u, Edge *b, point *v, int s) {
        Edge *e = make_edge(*u, *v);
        if (s == 0) {
            splice(Prev(a, u), e, u);
            splice(b, e, v);
        }
        else {
            splice(a, e, u);
            splice(Prev(b, v), e, v);
        }
        return e;
    }
    void lower_tangent(Edge * & l, Edge * & r, point * & s, point * & u) {
        point *dl = Op(l, s), *dr = Op(r, u);
        while (1) {
            if (dcmp(cross((*s), (*dl), (*u))) > 0) {
                l = Prev(l, dl); s = dl; dl = Op(l, s);
            }
            else if (dcmp(cross((*u), (*dr), (*s))) < 0) {
                r = Next(r, dr); u = dr; dr = Op(r, u);
            }
            else break;
        }
    }
    void merge(Edge *r_cw_l, point *s, Edge *l_ccw_r, point *u, Edge
        **l_tangent) {
        Edge *b, *lc, *rc;
        point *dlc, *drc;
        double crc, clc;
        lower_tangent(r_cw_l, l_ccw_r, s, u);
        b = join(r_cw_l, s, l_ccw_r, u, 1);
        *l_tangent = b;
        do {
            lc = Next(b, s); rc = Prev(b, u); dlc = Op(lc, s); drc = Op(rc, u);
            double cplc = cross(*dlc, *s, *u);
            double cprc = cross(*drc, *s, *u);
            bool alc = dcmp(cplc)>0, arc = dcmp(cprc)>0;
            if (!alc && !arc) break;
            if (alc) {
                clc = dot(*dlc, *s, *u) / cplc;
                do {
                    Edge * next = Next(lc, s);
                    point & dest = *Op(next, s);
                    double cpn = cross(dest, *s, *u);
                    if (dcmp(cpn) <= 0) break;
                    double cn = dot(dest, *s, *u) / cpn;
                    if (dcmp(cn - clc)>0) break;
                    delete_edge(lc);
                    lc = next;
                    clc = cn;
                } while (1);
            }
            if (arc) {
                crc = (double)dot(*drc, *s, *u) / cprc;
                do {
                    Edge * prev = Prev(rc, u);
                    point & dest = *Op(prev, u);
                    double cpp = cross(dest, *s, *u);
                    if (dcmp(cpp) <= 0) break;
                    double cp = dot(dest, *s, *u) / cpp;
                    if (dcmp(cp - crc) > 0) break;
                    delete_edge(rc);
                    rc = prev;
                    crc = cp;
                } while (1);
            }
            dlc = Op(lc, s); drc = Op(rc, u);
            if (!alc || (alc && arc && dcmp(crc - clc) < 0)) {
                b = join(b, s, rc, drc, 1);
                u = drc;
            }
            else {
                b = join(lc, dlc, b, u, 1);
                s = dlc;
            }
        } while (1);
    }
    void divide(point *p, int l, int r, Edge * & l_ccw, Edge * & r_cw) {
        int n = r - l;
        Edge *l_ccw_l, *r_cw_l, *l_ccw_r, *r_cw_r, *l_tangent, *c;
        if (n == 2) {
            l_ccw = r_cw = make_edge(p[l], p[l + 1]);
        }
        else if (n == 3) {
            Edge * a = make_edge(p[l], p[l + 1]), *b = make_edge(p[l + 1], p[l + 2]);
            splice(a, b, &p[l + 1]);
            double c_p = cross(p[l], p[l + 1], p[l + 2]);
            if (dcmp(c_p)>0) {
                c = join(a, &p[l], b, &p[l + 2], 1); l_ccw = a; r_cw = b;
            }
            else if (dcmp(c_p) < 0) {
                c = join(a, &p[l], b, &p[l + 2], 0); l_ccw = c; r_cw = c;
            }
            else {
                l_ccw = a; r_cw = b;
            }
        }
        else if (n > 3) {
            int split = (l + r) / 2;
            divide(p, l, split, l_ccw_l, r_cw_l);
            divide(p, split, r, l_ccw_r, r_cw_r);
            merge(r_cw_l, &p[split - 1], l_ccw_r, &p[split], &l_tangent);
            if (l_tangent->o == &p[l]) l_ccw_l = l_tangent;
            if (l_tangent->d == &p[r - 1]) r_cw_r = l_tangent;
            l_ccw = l_ccw_l; r_cw = r_cw_r;
        }
    }
} de;
void getEdge(int &k, int n) {
    k = 0;
    Edge *st, *cur;
    point *u, *v;
    for (int i = 0; i < n; ++i) {
        u = &pnt[i];
        st = cur = u->e;
        do {
            v = Op(cur, u);
            if (u < v)
                addEdge(k, u->id, v->id, dis(*u, *v));
        } while ((cur = Next(cur, u)) != st);
    }

}
void enum_triangle(point *ps, int n) {
    Edge *e_start, *e, *nxt;
    point *u, *v, *w;
    for (int i = 0; i < n; i++) {
        u = &ps[i];
        e_start = e = u->e;
        do {
            v = Op(e, u);
            if (u < v) {
                nxt = Next(e, u);
                w = Op(nxt, u);
                if (u < w && Next(nxt, w) == Prev(e, v)) {
                    // now, (u v w) is a triangle!!!!!!
                    // 这时，uvw 的外接圆是空的（不包含 ps 中的其他点），如果要求最大空圆，则计算 uvw 的外接圆就可以！
                }
            }
        } while ((e = Next(e, u)) != e_start);
    }
}

*/
// 任意四面体体积公式
// 已知任意四面体P-ABC，记PA=a，PB=b，PC=c，cos角APB=x，cos角BPC=y，cos角CPA=z，
// 则：V=1/6*abc*sqrt(1+2xyz-x^2-y^2-z^2)
//Pick定理：2S = 2a + b - 2，其中a表示多边形内部的点数，b表示多边形边界上的点数，s表示多边形的面积

//三维坐标绕 ( x, y, z )旋转