bzoj 4445 小凸想跑步 - 半平面交
假设点的位置是$(x, y)$,那么可以用叉积来计算三角形的面积。
这样可以列出$n - 1$个不等式。
将每个化成形如$ax + by + c \leqslant 0$的形式。
然后分类讨论($b = 0的时候需要特殊处理$)将它转换成二维平面上的半平面。
接着做半平面交,算面积就好了。
Code
1 /** 2 * bzoj 3 * Problem#4445 4 * Accepted 5 * Time: 288ms 6 * Memory: 20068k 7 */ 8 #include <algorithm> 9 #include <iostream> 10 #include <cassert> 11 #include <cstring> 12 #include <cstdio> 13 #include <cmath> 14 using namespace std; 15 typedef bool boolean; 16 17 const double eps = 1e-8; 18 #define check(_s) cerr << _s << endl; 19 20 inline int dcmp(double x) { 21 if (fabs(x) <= eps) return 0; 22 return (x < 0) ? (-1) : (1); 23 } 24 25 typedef class Vector { 26 public: 27 double x, y; 28 29 Vector(double x = 0.0, double y = 0.0):x(x), y(y) { } 30 }Point, Vector; 31 32 ostream& operator << (ostream& os, Point a) { 33 os << "(" << a.x << ", " << a.y << ")"; 34 return os; 35 } 36 37 Vector operator + (Vector a, Vector b) { 38 return Vector(a.x + b.x, a.y + b.y); 39 } 40 41 Vector operator - (Vector a, Vector b) { 42 return Vector(a.x - b.x, a.y - b.y); 43 } 44 45 Vector operator * (double x, Vector a) { 46 return Vector(a.x * x, a.y * x); 47 } 48 49 double cross(Vector a, Vector b) { 50 return a.x * b.y - a.y * b.x; 51 } 52 53 double dot(Vector a, Vector b) { 54 return a.x * b.x + a.y * b.y; 55 } 56 57 Point getLineIntersection(Point A, Vector v, Point B, Vector u) { 58 double t = (cross(B - A, u) / cross(v, u)); 59 return A + t * v; 60 } 61 62 typedef class Line { 63 public: 64 Point p; 65 Vector v; 66 double ang; 67 68 Line() { } 69 Line(Point p, Vector v):p(p), v(v) { 70 ang = atan2(v.y, v.x); 71 } 72 73 boolean operator < (Line b) const { 74 // return dcmp(cross(v, b.v)) > 0; 75 return ang < b.ang; 76 } 77 }Line; 78 79 const int N = 1e5 + 5; 80 81 int n; 82 int st = 1, ed = 0, top; 83 Line *ls, *qs; 84 Point *ps; 85 86 inline void init() { 87 scanf("%d", &n); 88 ps = new Point[(n << 1) + 5]; 89 for (int i = 0; i < n; i++) 90 scanf("%lf%lf", &ps[i].x, &ps[i].y); 91 } 92 93 inline void build() { 94 int d; 95 Point cur, nxt, vec; 96 ls = new Line[(n << 1) + 5]; 97 vec = ps[1] - ps[0]; 98 ls[++top] = Line(ps[0], vec); 99 double bx = vec.y, by = vec.x, bc = cross(vec, ps[1]), nx, ny, nc; 100 for (int i = 1; i < n; i++) { 101 cur = ps[i], nxt = ps[(i + 1) % n], vec = nxt - cur; 102 nx = bx - vec.y, ny = by - vec.x, nc = bc - cross(vec, cur); 103 d = dcmp(ny); 104 if (!d) { 105 d = dcmp(nx); 106 if (!d) continue; 107 ls[++top] = Line(Point(-nc / nx, 0), Vector(0, -d)); 108 109 } else { 110 nx /= ny, nc /= ny; 111 ls[++top] = Line(Point(0, nc), Vector(-d, -nx * d)); 112 } 113 ls[++top] = Line(cur, vec); 114 } 115 } 116 117 boolean isCrossingOut(Line b, Line cur, Line a) { 118 assert (dcmp(cross(cur.v, a.v))); 119 // if (!dcmp(cross(cur.v, a.v))) return true; 120 Point p = getLineIntersection(cur.p, cur.v, a.p, a.v); 121 return dcmp(cross(p - b.p, b.v)) > 0; 122 } 123 124 inline void HalfPlaneIntersection(int n) { 125 sort(ls + 1, ls + n + 1); 126 qs = new Line[n + 5]; 127 for (int i = 1; i <= n; i++) { 128 while (st < ed && isCrossingOut(ls[i], qs[ed], qs[ed - 1])) ed--; 129 while (st < ed && isCrossingOut(ls[i], qs[st], qs[st + 1])) st++; 130 qs[++ed] = ls[i]; 131 if (st < ed && !dcmp(cross(qs[ed].v, qs[ed - 1].v))) { 132 ed--; 133 if (dcmp(cross(qs[ed].p - ls[i].v, ls[i].v)) < 0) 134 qs[ed] = ls[i]; 135 } 136 } 137 while (st < ed && isCrossingOut(qs[st], qs[ed], qs[ed - 1])) ed--; 138 } 139 140 inline double PolygonArea(Point *ps, int n) { 141 double rt = cross(ps[n], ps[1]); 142 for (int i = 1; i < n; i++) 143 rt += cross(ps[i], ps[i + 1]); 144 return fabs(rt) / 2; 145 } 146 147 inline void solve() { 148 double all = PolygonArea(ps - 1, n); 149 build(); 150 HalfPlaneIntersection(top); 151 // assert(st < ed - 1); 152 if (st >= ed - 1) { 153 puts("0.0000"); 154 return; 155 } 156 for (int i = st; i < ed; i++) 157 ps[i] = getLineIntersection(qs[i].p, qs[i].v, qs[i + 1].p, qs[i + 1].v); 158 ps[ed] = getLineIntersection(qs[ed].p, qs[ed].v, qs[st].p, qs[st].v); 159 printf("%.4f", PolygonArea(ps + (st - 1), ed - st + 1) / all); 160 } 161 162 int main() { 163 init(); 164 solve(); 165 return 0; 166 }