计算几何板子

struct node{
    double x,y;
};
node a,b,c;
//求两个点之间的长度
double len(node a,node b) {
    double tmp = sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
    return tmp;
}
//给出三个点，求三角形的面积  海伦公式: p=(a+b+c)/2，S=sqrt(p(p-a)(p-b)(p-c))
double Area(node a,node b,node c){
    double lena = len(a,b);
    double lenb = len(b,c);
    double lenc = len(a,c);

    double p = (lena + lenb + lenc) / 2.0;
    double S = sqrt(p * (p - lena) * (p - lenb) * (p - lenc));
    return S;
}
//三角形求每条边对应的圆心角
void Ran() {
    double lena = len(a,b);
    double lenb = len(b,c);
    double lenc = len(a,c);
    double A = acos((lenb * lenb + lenc * lenc - lena * lena) / (2 * lenb * lenc));
    double B = acos((lena * lena + lenc * lenc - lenb * lenb) / (2 * lena * lenc));
    double C = acos((lena * lena + lenb * lenb - lenc * lenc) / (2 * lena * lenb));
}
//求外接圆半径r = a * b * c / 4S 
double R(node a,node b,node c) {
    double lena = len(a,b);
    double lenb = len(b,c);
    double lenc = len(a,c);
    double S = Area(a,b,c);
    double R = lena * lenb * lenc / (4.0 * S);
}

 

const double eps = 1e-8, pi = acos(-1);
const int N = 1e5 + 10;
struct point {
    double x, y;

    point() {}

    point(double _x, double _y) {
        x = _x;
        y = _y;
    }

    point operator-(const point &p) const {
        return point(x - p.x, y - p.y);
    }

    point operator+(const point &p) const {
        return point(x + p.x, y + p.y);
    }

    double operator*(const point &p) const {
        return x * p.y - y * p.x;
    }

    double operator/(const point &p) const {
        return x * p.x + y * p.y;
    }
}a[N], b[N << 1];

struct line {
    point p1, p2;
    double ang;

    line() {}

    line(point p01, point p02) {
        p1 = p01;
        p2 = p02;
    }

    void getang() {
        ang = atan2(p2.y - p1.y, p2.x - p1.x);
    }
}c[N], q[N];

struct circle
{
    double x, y, r;
}cc;

struct point3 {
    double x, y, z;

    point3() {}

    point3(double _x, double _y, double _z) {
        x = _x;
        y = _y;
        z = _z;
    }

    point3 operator-(const point3 &p) const {
        return point3(x - p.x, y - p.y, z - p.z);
    }

    point3 operator+(const point3 &p) const {
        return point3(x + p.x, y + p.y, z - p.z);
    }

    point3 operator*(const point3 &p) const {
        return point3(y * p.z - z * p.y, z * p.x - x * p.z, x * p.y - y * p.x);
    }

    double operator/(const point3 &p) const {
        return x * p.x + y * p.y + z * p.z;
    }
};

int n;

bool cmp(const point &aa, const point &bb) {
    return aa.x < bb.x || (aa.x == bb.x && aa.y < bb.y);
}

bool cmpl(const line &aa, const line &bb) {
    if (abs(aa.ang - bb.ang) > eps)
        return aa.ang < bb.ang;
    return (bb.p2 - aa.p1) * (bb.p1 - bb.p2) > 0;
}

int sgn(double x) {
    if (x > eps)
        return 1;
    if (x < -eps)
        return -1;
    return 0;
}

double getdist2(const point &aa) {
    return (aa.x * aa.x + aa.y * aa.y);
}

bool online(const point &aa, const point &bb, const point &cc) {     //三点共线判断
    if (!(min(bb.x, cc.x) <= aa.x && aa.x <= max(bb.x, cc.x)))
        return 0;
    if (!(min(bb.y, cc.y) <= aa.y && aa.y <= max(bb.y, cc.y)))
        return 0;
    return !sgn((bb - aa) * (cc - aa));
}

bool checkline(const point &aa, const point &bb, const point &cc, const point &dd) {      //线段是否相交
    int c1 = sgn((bb - aa) * (cc - aa)), c2 = sgn((bb - aa) * (dd - aa)),
            c3 = sgn((dd - cc) * (aa - cc)), c4 = sgn((dd - cc) * (bb - cc));
    if (c1 * c2 < 0 && c3 * c4 < 0)
        return 1;//规范相交
    if (c1 == 0 && c2 == 0) {
        if (min(aa.x, bb.x) > max(cc.x, dd.x) || min(cc.x, dd.x) > max(aa.x, bb.x)
            || min(aa.y, bb.y) > max(cc.y, dd.y) || min(cc.y, dd.y) > max(aa.y, bb.y))
            return 0;
        return 1;//共线有重叠
    }
    if (online(cc, aa, bb) || online(dd, aa, bb)
        || online(aa, cc, dd) || online(bb, cc, dd))
        return 1;//端点相交
    return 0;
}

point getpoint(const point &aa, const point &bb,
               const point &cc, const point &dd) {      //两直线交点
    double a1, b1, c1, a2, b2, c2;
    point re;
    a1 = aa.y - bb.y;
    b1 = bb.x - aa.x;
    c1 = aa * bb;
    a2 = cc.y - dd.y;
    b2 = dd.x - cc.x;
    c2 = cc * dd;
    //以下为交点横纵坐标
    re.x = (c1 * b2 - c2 * b1) / (a2 * b1 - a1 * b2);
    re.y = (a2 * c1 - a1 * c2) / (a1 * b2 - a2 * b1);
    return re;
}

bool protrusion(point aa[]){    //判断多边形凹(凸)
    int flag = sgn((aa[1] - aa[0]) * (aa[2] - aa[0]));
    for (int i = 1; i < n; ++i)
        if (sgn((aa[(i + 1) % n] - aa[i]) * (aa[(i + 2) % n] - aa[i])) != flag)
            return 1;
    return 0;
}
bool inpolygon(const point &aa) {   //点在凸多边形内(外)
    int flag = sgn((a[0] - aa) * (a[1] - aa));
    for (int i = 1; i < n; ++i)
        if (sgn((a[i] - aa) * (a[(i + 1) % n] - aa)) != flag)
            return 0;
    return 1;
}

void transline(line &aa, double dist) {     //逆时针方向平移线段
    double d = sqrt((aa.p1.x - aa.p2.x) * (aa.p1.x - aa.p2.x) +
                    (aa.p1.y - aa.p2.y) * (aa.p1.y - aa.p2.y));
    point ta;
    ta.x = aa.p1.x + dist / d * (aa.p1.y - aa.p2.y);
    ta.y = aa.p1.y - dist / d * (aa.p1.x - aa.p2.x);
    aa.p2 = ta + aa.p2 - aa.p1;
    aa.p1 = ta;
}

void convexhull(point aa[], point bb[]) {   //凸包
    int len = 0;
    sort(aa, aa + n, cmp);
    bb[len++] = aa[0];//bb数组要开2倍(防止出现直线)
    bb[len++] = aa[1];
    for (int i = 2; i < n; ++i) {
        while (len > 1 && (aa[i] - bb[len - 2]) * (bb[len - 1] - bb[len - 2]) > 0)
            //若严格则加上等于
            --len;
        bb[len++] = aa[i];
    }
    int t = len;
    for (int i = n - 2; i >= 0; --i) {
        while (len > t && (aa[i] - bb[len - 2]) * (bb[len - 1] - bb[len - 2]) > 0)
            //同上
            --len;
        bb[len++] = aa[i];
    }
    --len;
    n = len;
}
bool checkout(const line &aa, const line &bb, const line &cc) {     //检查交点是否在向量顺时针侧
    point p = getpoint(aa.p1, aa.p2, bb.p1, bb.p2);
    return (cc.p1 - p) * (cc.p2 - p) < 0;//如果不允许共线或算面积 则此处不取等
}

double halfplane(point aa[], line bb[]) {   //半平面交
    sort(bb, bb + n, cmpl);
    int n2 = 1;
    for (int i = 1; i < n; ++i) {
        if (bb[i].ang - bb[i - 1].ang > eps)
            ++n2;
        bb[n2 - 1] = bb[i];
    }
    n = n2;
    int front = 0, tail = 0;
    q[tail++] = bb[0], q[tail++] = bb[1];
    for (int i = 2; i < n; ++i) {
        while (front + 1 < tail && checkout(q[tail - 2], q[tail - 1], bb[i]))
            --tail;
        while (front + 1 < tail && checkout(q[front], q[front + 1], bb[i]))
            ++front;
        q[tail++] = bb[i];
    }
    while (front + 2 < tail && checkout(q[tail - 2], q[tail - 1], q[front]))
        --tail;
    while (front + 2 < tail && checkout(q[front], q[front + 1], q[tail - 1]))
        ++front;
    if (front + 2 >= tail)
        return 0;
    int j = 0;
    for (int i = front; i < tail; ++i, ++j) {
        aa[j] = getpoint(q[i].p1, q[i].p2, q[(i != tail - 1 ? i + 1 : front)].p1,
                         q[(i != tail - 1 ? i + 1 : front)].p2);
    }
    double re = 0;
    for (int i = 1; i < j - 1; ++i)
        re += (aa[i] - aa[0]) * (aa[i + 1] - aa[0]);
    return abs(re * 0.5);
}
double calipers(point aa[]) {   //旋转卡壳
    double re = 0;
    aa[n] = aa[0];
    int now = 1;
    for (int i = 0; i < n; ++i) {
        while ((aa[i + 1] - aa[i]) * (aa[now + 1] - aa[i]) >
               (aa[i + 1] - aa[i]) * (aa[now] - aa[i]))
            now = (now == n - 1 ? 0 : now + 1);
        re = max(re, getdist2(aa[now] - aa[i]));
    }
    return re;
}

line bisector(const point &aa, const point &bb) {       //中垂线(正方形)
    double mx, my;
    mx = (aa.x + bb.x) / 2;
    my = (aa.y + bb.y) / 2;
    line cc;
    cc.p1.x = mx - (aa.y - bb.y) / 2;
    cc.p1.y = my + (aa.x - bb.x) / 2;
    cc.p2 = aa + bb - cc.p1;
    cc.getang();
    return cc;
}

point rotate(const point &p, double cost, double sint) {    //逆时针向量旋转
    double x = p.x, y = p.y;
    return point(x * cost - y * sint, x * sint + y * cost);
}
void getpoint(circle c1, circle c2) {   //已确保两圆有交点时求出两圆交点
    long double dab = sqrt(getdist2(point(c1.x, c1.y) -
                                    point(c2.x, c2.y)));
    if (c1.r > c2.r)
        swap(c1, c2);
    long double cost = (c1.r * c1.r + dab * dab - c2.r * c2.r) /
                       (c1.r * dab * 2);
    long double sint = sqrt(1 - cost * cost);
    point re = rotate(point(c2.x, c2.y) - point(c1.x, c1.y), cost, sint);
    re.x = c1.x + re.x * (c1.r / dab);
    re.y = c1.y + re.y * (c1.r / dab);
    point re2 = rotate(point(c2.x, c2.y) - point(c1.x, c1.y), cost, -sint);
    re2.x = c1.x + re2.x * (c1.r / dab);
    re2.y = c1.y + re2.y * (c1.r / dab);
}

double mycos(double B, double C, double A) {   //余弦定理 给定边长
    return (B * B + C * C - A * A) / (B * C * 2);
}
double mycos2(const point &aa, const point &bb, const point &cc) {//余弦定理 给定点坐标
    double C2 = getdist2(aa - bb), A2 = getdist2(bb - cc), B2 = getdist2(cc - aa);
    return (B2 + C2 - A2) / (sqrt(B2 * C2) * 2);
}
point mirror_point(const point &aa, const point &bb, const point &cc) {     //轴对称点 也可用来求垂足
    double cost, sint;
    cost = mycos2(bb, cc, aa);
    sint = sqrt(1 - cost * cost);
    if ((cc - bb) * (aa - bb) > 0)
        sint = -sint;
    point re;
    re = rotate(aa - bb, cost, sint) + bb;
    re = rotate(re - bb, cost, sint) + bb;
    return re;
}
double area(const point &aa, const point &bb, const point &cc) {//求三角形面积
    return abs((bb - aa) * (cc - aa) / 2);
}
void max_triangle(const point aa[]) {   //求凸包上最大三角形
    double ans = 0;
    int p1, p2, p3;
    for (int i = 0; i < n; ++i) {
        int j = (i + 1) % n;
        int k = (j + 1) % n;
        while (k != i && area(aa[i], aa[j], aa[k]) < area(aa[i], aa[j], aa[(k + 1) % n]))
            k = (k + 1) % n;
        if (k == i)
            continue;
        int kk = (k + 1) % n;
        while (j != kk && k != i) {
            if (ans < area(aa[i], aa[j], aa[k])) {
                ans = area(aa[i], aa[j], aa[k]);//三角形面积
                p1 = i;
                p2 = j;
                p3 = k;
            }
            while (k != i && area(aa[i], aa[j], aa[k]) < area(aa[i], aa[j], aa[(k + 1) % n]))
                k = (k + 1) % n;
            j = (j + 1) % n;
        }
    }
    point q1, q2, q3;
    q1 = b[p2] + b[p3] - b[p1];
    q2 = b[p1] + b[p3] - b[p2];
    q3 = b[p1] + b[p2] - b[p3];
    //三角形上三个点
}

void get_panel(const point3 &p1, const point3 &p2, const point3 &p3,
               double &A, double &B, double &C, double &D) {   //由三点确定一个平面方程
    A = ((p2.y - p1.y) * (p3.z - p1.z) - (p2.z - p1.z) * (p3.y - p1.y));
    B = ((p2.z - p1.z) * (p3.x - p1.x) - (p2.x - p1.x) * (p3.z - p1.z));
    C = ((p2.x - p1.x) * (p3.y - p1.y) - (p2.y - p1.y) * (p3.x - p1.x));
    D = (-(A * p1.x + B * p1.y + C * p1.z));
}

double dist_panel(const point3 &pt, double A, double B, double C, double D) {   //点到平面距离
    return abs(A * pt.x + B * pt.y + C * pt.z + D) /
           sqrt(A * A + B * B + C * C);
}