二维计算几何模板

#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 )旋转