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); }