Description
Updog is watching a plane object with a telescope. The field of vision in the telescope can be described as a circle. The center is at the origin and the radius is R. The object can be seen as a simple polygon of N vertexes. Updog wants to know the area of the part of the object that can be seen in the telescope.
Input
The input will contain several test cases. For each case:
The first line contains only one real number R.
The second line contains an integer N. The following N lines contain two real numbers xi and yi each, which describe the coordinates of a vertex. Two vertexes in adjacent lines are also adjacent on the polygon.
Constraints: 3 ≤ N ≤50, 0.1 ≤ R ≤1000, -1000 ≤ xi, yi ≤ 1000
Output
For each test case, output one real number on a separate line, which is the area of the part that can be seen. The result should be rounded to two decimal places.
题目大意：给一个圆心在原点的圆，半径为r。再给一个简单多边形，问这个简单多边形和圆的公共面积是多少。
思路：跟POJ2986一个思路，把简单多边形的每一个点和原点连线，就把这个多边形和圆的交变成了多个三角形与圆的交，根据有向面积的思路，加加减减就可以得到公共面积。
PS：所谓简单多边形，是指所有边都不相交的多边形。
PS2：之前求直线与圆的交点的时候，没留意到-EPS<delta<0，然后要求sqrt(delta)的情况……WA了很多很多次……计算几何真是可怕……
代码（0MS）：
1 #include <cstdio> 2 #include <cstring> 3 #include <iostream> 4 #include <algorithm> 5 #include <cmath> 6 using namespace std; 7 #define sqr(x) ((x) * (x)) 8 9 const int MAXN = 55; 10 const double EPS = 1e-10; 11 const double PI = acos(-1.0);//3.14159265358979323846 12 const double INF = 1; 13 14 inline int sgn(double x) { 15 return (x > EPS) - (x < -EPS); 16 } 17 18 struct Point { 19 double x, y, ag; 20 Point() {} 21 Point(double x, double y): x(x), y(y) {} 22 void read() { 23 scanf("%lf%lf", &x, &y); 24 } 25 bool operator == (const Point &rhs) const { 26 return sgn(x - rhs.x) == 0 && sgn(y - rhs.y) == 0; 27 } 28 bool operator < (const Point &rhs) const { 29 if(y != rhs.y) return y < rhs.y; 30 return x < rhs.x; 31 } 32 Point operator + (const Point &rhs) const { 33 return Point(x + rhs.x, y + rhs.y); 34 } 35 Point operator - (const Point &rhs) const { 36 return Point(x - rhs.x, y - rhs.y); 37 } 38 Point operator * (const double &b) const { 39 return Point(x * b, y * b); 40 } 41 Point operator / (const double &b) const { 42 return Point(x / b, y / b); 43 } 44 double operator * (const Point &rhs) const { 45 return x * rhs.x + y * rhs.y; 46 } 47 double length() { 48 return sqrt(x * x + y * y); 49 } 50 double angle() { 51 return atan2(y, x); 52 } 53 Point unit() { 54 return *this / length(); 55 } 56 void makeAg() { 57 ag = atan2(y, x); 58 } 59 void print() { 60 printf("%.10f %.10f\n", x, y); 61 } 62 }; 63 typedef Point Vector; 64 65 double dist(const Point &a, const Point &b) { 66 return (a - b).length(); 67 } 68 69 double cross(const Point &a, const Point &b) { 70 return a.x * b.y - a.y * b.x; 71 } 72 //ret >= 0 means turn right 73 double cross(const Point &sp, const Point &ed, const Point &op) { 74 return cross(sp - op, ed - op); 75 } 76 77 double area(const Point& a, const Point &b, const Point &c) { 78 return fabs(cross(a - c, b - c)) / 2; 79 } 80 //counter-clockwise 81 Point rotate(const Point &p, double angle, const Point &o = Point(0, 0)) { 82 Point t = p - o; 83 double x = t.x * cos(angle) - t.y * sin(angle); 84 double y = t.y * cos(angle) + t.x * sin(angle); 85 return Point(x, y) + o; 86 } 87 88 double cosIncludeAngle(const Point &a, const Point &b, const Point &o) { 89 Point p1 = a - o, p2 = b - o; 90 return (p1 * p2) / (p1.length() * p2.length()); 91 } 92 93 double includedAngle(const Point &a, const Point &b, const Point &o) { 94 return acos(cosIncludeAngle(a, b, o)); 95 /* 96 double ret = abs((a - o).angle() - (b - o).angle()); 97 if(sgn(ret - PI) > 0) ret = 2 * PI - ret; 98 return ret; 99 */ 100 } 101 102 struct Seg { 103 Point st, ed; 104 double ag; 105 Seg() {} 106 Seg(Point st, Point ed): st(st), ed(ed) {} 107 void read() { 108 st.read(); ed.read(); 109 } 110 void makeAg() { 111 ag = atan2(ed.y - st.y, ed.x - st.x); 112 } 113 }; 114 typedef Seg Line; 115 116 //ax + by + c > 0 117 Line buildLine(double a, double b, double c) { 118 if(sgn(a) == 0 && sgn(b) == 0) return Line(Point(sgn(c) > 0 ? -1 : 1, INF), Point(0, INF)); 119 if(sgn(a) == 0) return Line(Point(sgn(b), -c/b), Point(0, -c/b)); 120 if(sgn(b) == 0) return Line(Point(-c/a, 0), Point(-c/a, sgn(a))); 121 if(b < 0) return Line(Point(0, -c/b), Point(1, -(a + c) / b)); 122 else return Line(Point(1, -(a + c) / b), Point(0, -c/b)); 123 } 124 125 void moveRight(Line &v, double r) { 126 double dx = v.ed.x - v.st.x, dy = v.ed.y - v.st.y; 127 dx = dx / dist(v.st, v.ed) * r; 128 dy = dy / dist(v.st, v.ed) * r; 129 v.st.x += dy; v.ed.x += dy; 130 v.st.y -= dx; v.ed.y -= dx; 131 } 132 133 bool isOnSeg(const Seg &s, const Point &p) { 134 return (p == s.st || p == s.ed) || 135 (((p.x - s.st.x) * (p.x - s.ed.x) < 0 || 136 (p.y - s.st.y) * (p.y - s.ed.y) < 0) && 137 sgn(cross(s.ed, p, s.st)) == 0); 138 } 139 140 bool isInSegRec(const Seg &s, const Point &p) { 141 return sgn(min(s.st.x, s.ed.x) - p.x) <= 0 && sgn(p.x - max(s.st.x, s.ed.x)) <= 0 142 && sgn(min(s.st.y, s.ed.y) - p.y) <= 0 && sgn(p.y - max(s.st.y, s.ed.y)) <= 0; 143 } 144 145 bool isIntersected(const Point &s1, const Point &e1, const Point &s2, const Point &e2) { 146 return (max(s1.x, e1.x) >= min(s2.x, e2.x)) && 147 (max(s2.x, e2.x) >= min(s1.x, e1.x)) && 148 (max(s1.y, e1.y) >= min(s2.y, e2.y)) && 149 (max(s2.y, e2.y) >= min(s1.y, e1.y)) && 150 (cross(s2, e1, s1) * cross(e1, e2, s1) >= 0) && 151 (cross(s1, e2, s2) * cross(e2, e1, s2) >= 0); 152 } 153 154 bool isIntersected(const Seg &a, const Seg &b) { 155 return isIntersected(a.st, a.ed, b.st, b.ed); 156 } 157 158 bool isParallel(const Seg &a, const Seg &b) { 159 return sgn(cross(a.ed - a.st, b.ed - b.st)) == 0; 160 } 161 162 //return Ax + By + C =0 's A, B, C 163 void Coefficient(const Line &L, double &A, double &B, double &C) { 164 A = L.ed.y - L.st.y; 165 B = L.st.x - L.ed.x; 166 C = L.ed.x * L.st.y - L.st.x * L.ed.y; 167 } 168 //point of intersection 169 Point operator * (const Line &a, const Line &b) { 170 double A1, B1, C1; 171 double A2, B2, C2; 172 Coefficient(a, A1, B1, C1); 173 Coefficient(b, A2, B2, C2); 174 Point I; 175 I.x = - (B2 * C1 - B1 * C2) / (A1 * B2 - A2 * B1); 176 I.y = (A2 * C1 - A1 * C2) / (A1 * B2 - A2 * B1); 177 return I; 178 } 179 180 bool isEqual(const Line &a, const Line &b) { 181 double A1, B1, C1; 182 double A2, B2, C2; 183 Coefficient(a, A1, B1, C1); 184 Coefficient(b, A2, B2, C2); 185 return sgn(A1 * B2 - A2 * B1) == 0 && sgn(A1 * C2 - A2 * C1) == 0 && sgn(B1 * C2 - B2 * C1) == 0; 186 } 187 188 double Point_to_Line(const Point &p, const Line &L) { 189 return fabs(cross(p, L.st, L.ed)/dist(L.st, L.ed)); 190 } 191 192 double Point_to_Seg(const Point &p, const Seg &L) { 193 if(sgn((L.ed - L.st) * (p - L.st)) < 0) return dist(p, L.st); 194 if(sgn((L.st - L.ed) * (p - L.ed)) < 0) return dist(p, L.ed); 195 return Point_to_Line(p, L); 196 } 197 198 double Seg_to_Seg(const Seg &a, const Seg &b) { 199 double ans1 = min(Point_to_Seg(a.st, b), Point_to_Seg(a.ed, b)); 200 double ans2 = min(Point_to_Seg(b.st, a), Point_to_Seg(b.ed, a)); 201 return min(ans1, ans2); 202 } 203 204 struct Circle { 205 Point c; 206 double r; 207 Circle() {} 208 Circle(Point c, double r): c(c), r(r) {} 209 void read() { 210 c.read(); 211 scanf("%lf", &r); 212 } 213 double area() const { 214 return PI * r * r; 215 } 216 bool contain(const Circle &rhs) const { 217 return sgn(dist(c, rhs.c) + rhs.r - r) <= 0; 218 } 219 bool contain(const Point &p) const { 220 return sgn(dist(c, p) - r) <= 0; 221 } 222 bool intersect(const Circle &rhs) const { 223 return sgn(dist(c, rhs.c) - r - rhs.r) < 0; 224 } 225 bool tangency(const Circle &rhs) const { 226 return sgn(dist(c, rhs.c) - r - rhs.r) == 0; 227 } 228 Point pos(double angle) const { 229 Point p = Point(c.x + r, c.y); 230 return rotate(p, angle, c); 231 } 232 }; 233 234 double CommonArea(const Circle &A, const Circle &B) { 235 double area = 0.0; 236 const Circle & M = (A.r > B.r) ? A : B; 237 const Circle & N = (A.r > B.r) ? B : A; 238 double D = dist(M.c, N.c); 239 if((D < M.r + N.r) && (D > M.r - N.r)) { 240 double cosM = (M.r * M.r + D * D - N.r * N.r) / (2.0 * M.r * D); 241 double cosN = (N.r * N.r + D * D - M.r * M.r) / (2.0 * N.r * D); 242 double alpha = 2 * acos(cosM); 243 double beta = 2 * acos(cosN); 244 double TM = 0.5 * M.r * M.r * (alpha - sin(alpha)); 245 double TN = 0.5 * N.r * N.r * (beta - sin(beta)); 246 area = TM + TN; 247 } 248 else if(D <= M.r - N.r) { 249 area = N.area(); 250 } 251 return area; 252 } 253 254 int intersection(const Seg &s, const Circle &cir, Point &p1, Point &p2) { 255 double angle = cosIncludeAngle(s.ed, cir.c, s.st); 256 //double angle1 = cos(includedAngle(s.ed, cir.c, s.st)); 257 double B = dist(cir.c, s.st); 258 double a = 1, b = -2 * B * angle, c = sqr(B) - sqr(cir.r); 259 double delta = sqr(b) - 4 * a * c; 260 if(sgn(delta) < 0) return 0; 261 if(sgn(delta) == 0) delta = 0; 262 double x1 = (-b - sqrt(delta)) / (2 * a), x2 = (-b + sqrt(delta)) / (2 * a); 263 Vector v = (s.ed - s.st).unit(); 264 p1 = s.st + v * x1; 265 p2 = s.st + v * x2; 266 return 1 + sgn(delta); 267 } 268 269 double CommonArea(const Circle &cir, Point p1, Point p2) { 270 if(p1 == cir.c || p2 == cir.c) return 0; 271 if(cir.contain(p1) && cir.contain(p2)) { 272 return area(cir.c, p1, p2); 273 } else if(!cir.contain(p1) && !cir.contain(p2)) { 274 Point q1, q2; 275 int t = intersection(Line(p1, p2), cir, q1, q2); 276 if(t == 0) { 277 double angle = includedAngle(p1, p2, cir.c); 278 return 0.5 * sqr(cir.r) * angle; 279 } else { 280 double angle1 = includedAngle(p1, p2, cir.c); 281 double angle2 = includedAngle(q1, q2, cir.c); 282 if(isInSegRec(Seg(p1, p2), q1))return 0.5 * sqr(cir.r) * (angle1 - angle2 + sin(angle2)); 283 else return 0.5 * sqr(cir.r) * angle1; 284 } 285 } else { 286 if(cir.contain(p2)) swap(p1, p2); 287 Point q1, q2; 288 intersection(Line(p1, p2), cir, q1, q2); 289 double angle = includedAngle(q2, p2, cir.c); 290 double a = area(cir.c, p1, q2); 291 double b = 0.5 * sqr(cir.r) * angle; 292 return a + b; 293 } 294 } 295 296 struct Triangle { 297 Point p[3]; 298 Triangle() {} 299 Triangle(Point *t) { 300 for(int i = 0; i < 3; ++i) p[i] = t[i]; 301 } 302 void read() { 303 for(int i = 0; i < 3; ++i) p[i].read(); 304 } 305 double area() const { 306 return ::area(p[0], p[1], p[2]); 307 } 308 Point& operator[] (int i) { 309 return p[i]; 310 } 311 }; 312 313 double CommonArea(Triangle tir, const Circle &cir) { 314 double ret = 0; 315 ret += sgn(cross(tir[0], cir.c, tir[1])) * CommonArea(cir, tir[0], tir[1]); 316 ret += sgn(cross(tir[1], cir.c, tir[2])) * CommonArea(cir, tir[1], tir[2]); 317 ret += sgn(cross(tir[2], cir.c, tir[0])) * CommonArea(cir, tir[2], tir[0]); 318 return abs(ret); 319 } 320 321 struct Poly { 322 int n; 323 Point p[MAXN];//p[n] = p[0] 324 void init(Point *pp, int nn) { 325 n = nn; 326 for(int i = 0; i < n; ++i) p[i] = pp[i]; 327 p[n] = p[0]; 328 } 329 double area() { 330 if(n < 3) return 0; 331 double s = p[0].y * (p[n - 1].x - p[1].x); 332 for(int i = 1; i < n; ++i) 333 s += p[i].y * (p[i - 1].x - p[i + 1].x); 334 return s / 2; 335 } 336 }; 337 //the convex hull is clockwise 338 void Graham_scan(Point *p, int n, int *stk, int &top) {//stk[0] = stk[top] 339 sort(p, p + n); 340 top = 1; 341 stk[0] = 0; stk[1] = 1; 342 for(int i = 2; i < n; ++i) { 343 while(top && cross(p[i], p[stk[top]], p[stk[top - 1]]) <= 0) --top; 344 stk[++top] = i; 345 } 346 int len = top; 347 stk[++top] = n - 2; 348 for(int i = n - 3; i >= 0; --i) { 349 while(top != len && cross(p[i], p[stk[top]], p[stk[top - 1]]) <= 0) --top; 350 stk[++top] = i; 351 } 352 } 353 //use for half_planes_cross 354 bool cmpAg(const Line &a, const Line &b) { 355 if(sgn(a.ag - b.ag) == 0) 356 return sgn(cross(b.ed, a.st, b.st)) < 0; 357 return a.ag < b.ag; 358 } 359 //clockwise, plane is on the right 360 bool half_planes_cross(Line *v, int vn, Poly &res, Line *deq) { 361 int i, n; 362 sort(v, v + vn, cmpAg); 363 for(i = n = 1; i < vn; ++i) { 364 if(sgn(v[i].ag - v[i-1].ag) == 0) continue; 365 v[n++] = v[i]; 366 } 367 int head = 0, tail = 1; 368 deq[0] = v[0], deq[1] = v[1]; 369 for(i = 2; i < n; ++i) { 370 if(isParallel(deq[tail - 1], deq[tail]) || isParallel(deq[head], deq[head + 1])) 371 return false; 372 while(head < tail && sgn(cross(v[i].ed, deq[tail - 1] * deq[tail], v[i].st)) > 0) 373 --tail; 374 while(head < tail && sgn(cross(v[i].ed, deq[head] * deq[head + 1], v[i].st)) > 0) 375 ++head; 376 deq[++tail] = v[i]; 377 } 378 while(head < tail && sgn(cross(deq[head].ed, deq[tail - 1] * deq[tail], deq[head].st)) > 0) 379 --tail; 380 while(head < tail && sgn(cross(deq[tail].ed, deq[head] * deq[head + 1], deq[tail].st)) > 0) 381 ++head; 382 if(tail <= head + 1) return false; 383 res.n = 0; 384 for(i = head; i < tail; ++i) 385 res.p[res.n++] = deq[i] * deq[i + 1]; 386 res.p[res.n++] = deq[head] * deq[tail]; 387 res.n = unique(res.p, res.p + res.n) - res.p; 388 res.p[res.n] = res.p[0]; 389 return true; 390 } 391 392 //ix and jx is the points whose distance is return, res.p[n - 1] = res.p[0], res must be clockwise 393 double dia_rotating_calipers(Poly &res, int &ix, int &jx) { 394 double dia = 0; 395 int q = 1; 396 for(int i = 0; i < res.n - 1; ++i) { 397 while(sgn(cross(res.p[i], res.p[q + 1], res.p[i + 1]) - cross(res.p[i], res.p[q], res.p[i + 1])) > 0) 398 q = (q + 1) % (res.n - 1); 399 if(sgn(dist(res.p[i], res.p[q]) - dia) > 0) { 400 dia = dist(res.p[i], res.p[q]); 401 ix = i; jx = q; 402 } 403 if(sgn(dist(res.p[i + 1], res.p[q]) - dia) > 0) { 404 dia = dist(res.p[i + 1], res.p[q]); 405 ix = i + 1; jx = q; 406 } 407 } 408 return dia; 409 } 410 //a and b must be clockwise, find the minimum distance between two convex hull 411 double half_rotating_calipers(Poly &a, Poly &b) { 412 int sa = 0, sb = 0; 413 for(int i = 0; i < a.n; ++i) if(sgn(a.p[i].y - a.p[sa].y) < 0) sa = i; 414 for(int i = 0; i < b.n; ++i) if(sgn(b.p[i].y - b.p[sb].y) < 0) sb = i; 415 double tmp, ans = dist(a.p[0], b.p[0]); 416 for(int i = 0; i < a.n; ++i) { 417 while(sgn(tmp = cross(a.p[sa], a.p[sa + 1], b.p[sb + 1]) - cross(a.p[sa], a.p[sa + 1], b.p[sb])) > 0) 418 sb = (sb + 1) % (b.n - 1); 419 if(sgn(tmp) < 0) ans = min(ans, Point_to_Seg(b.p[sb], Seg(a.p[sa], a.p[sa + 1]))); 420 else ans = min(ans, Seg_to_Seg(Seg(a.p[sa], a.p[sa + 1]), Seg(b.p[sb], b.p[sb + 1]))); 421 sa = (sa + 1) % (a.n - 1); 422 } 423 return ans; 424 } 425 426 double rotating_calipers(Poly &a, Poly &b) { 427 return min(half_rotating_calipers(a, b), half_rotating_calipers(b, a)); 428 } 429 430 /*******************************************************************************************/ 431 432 Point p[MAXN]; 433 Circle cir; 434 double r; 435 int n; 436 437 int main() { 438 while(scanf("%lf", &r) != EOF) { 439 scanf("%d", &n); 440 for(int i = 0; i < n; ++i) p[i].read(); 441 p[n] = p[0]; 442 cir = Circle(Point(0, 0), r); 443 double ans = 0; 444 for(int i = 0; i < n; ++i) 445 ans += sgn(cross(p[i], Point(0, 0), p[i + 1])) * CommonArea(cir, p[i], p[i + 1]); 446 printf("%.2f\n", fabs(ans)); 447 } 448 }