UVA 11978-fukushima Nuclear Blast (two points + geometry)

Topic Link: UVA 11978-fukushima nuclear Blast

Two points, circle and polygon area intersection. Splits the polygon into an n-part triangle (each two points and the center of the circle) to calculate the direction area. The area of each triangle and circle is discussed in four categories.

#include <cstdio> #include <cstring> #include <cmath> #include <vector> #include <algorithm >using namespace Std;const Double pi = 4 * atan (1); const double EPS = 1e-10;inline int dcmp (double x) {if (Fabs (x) &L T EPS) return 0; else return x < 0? -1:1; }inline double getdistance (double x, double y) {return sqrt (x * x + y * y);} struct point {double x, y; Point (Double x = 0, double y = 0): x (x), Y (y) {}void read () {scanf ("%lf%lf", &x, &y);} void Write () {printf ("%lf%lf", x, y);} BOOL operator = = (CONST point& u) Const {return dcmp (x-u.x) = = 0 && dcmp (y-u.y) = = 0;} BOOL operator! = (const point& u) Const {return! ( *this = = u); }bool operator < (const point& u) Const {return x < u.x | | (x = = u.x && y < u.y); }bool operator > (const point& u) Const {return u < *this;} BOOL operator <= (const point& u) Const {return *this < u | | *this = u;} BOOL operator >= (const point&u) Const {return *this > U | | *this = = u;} Point operator + (const point& u) {return point (x + u.x, y + u.y);} Point operator-(const point& u) {return point (x-u.x, y-u.y);} Point operator * (const double u) {return point (x * u, y * u);} Point operator/(const double u) {return point (x/u, y/u);} Double operator * (const point& u) {return x*u.y-y*u.x;}}; typedef point VECTOR;STRUCT Line {Double A, B, C; Line (Double A = 0, double b = 0, double c = 0): A (a), B (b), C (c) {}};struct Circle {point o;double R; Circle () {}circle (Point O, double r = 0): O (O), R (r) {}void read () {o.read (), scanf ("%lf", &r);} Point Point (Double rad) {return point (o.x + cos (rad) *r, O.y + sin (rad) *r);} Double Getarea (double rad) {return rad * R * R/2;}}; Namespace punctual {Double getdistance (point A, point B) {double x=a.x-b.x, y=a.y-b.y; return sqrt (x*x + y*y);}}; namespace Vectorial {/* dot product: The product of two vector lengths multiplied by the cosine of their angle, the angle is greater than 90 degrees the time product is negative */double Getdot (vector A, vector b) {returN a.x * b.x + a.y * B.Y;  }/* Cross product: The cross product is equal to two vectors consisting of a triangle with a direction area of twice times, crosses (V, w) =-cross (W, v) */double Getcross (vector A, vector b) {return a.x * B.Y-A.Y * b.x; }double GetLength (Vector a) {return sqrt (Getdot (A, a));} Double Getplength (Vector a) {return Getdot (A, a);} Double Getangle (Vector u) {return atan2 (u.y, u.x);} Double Getangle (vector A, vector b) {return ACOs (Getdot (A, B)/GetLength (a)/GetLength (b));} Vector Rotate (vector A, double rad) {return vector (A.x*cos (RAD)-a.y*sin (RAD), A.x*sin (RAD) +a.y*cos (RAD));} /* unit normal */vector Getnormal (Vector a) {double L = getlength (a); return Vector (-a.y/l, a.x/l);}}; Namespace Linear {using namespace vectorial; Line GetLine (double x1, double y1, double x2, double y2) {return line (y2-y1, x1-x2, y1* (x2-x1)-x1* (y2-y1));} Line GetLine (double A, double b, point u) {return line (A, B, u.y * b-u.x * a);} BOOL Getintersection (line p, line Q, point& o) {if (Fabs (P.A * q.b-q.a * p.b) < EPS) return false;o.x = (Q.C * p. B-P.C * q.b)/(P.A * q.b-q.a * p.b); o.y = (Q.C * p.a-p.c * q.a)/(P.B * q.a-q.b * p.a); return true;} /* Intersection of linear PV and line QW */bool getintersection (point P, Vector v, point Q, Vector W, point& o) {if (dcmp (Getcross (V, w)) = = 0 ) return false; Vector u = p-q;double k = Getcross (w, u)/Getcross (V, w); o = p + v * K;return true;} /* points P to line ab distance */double getdistancetoline (point P, dot a, bit b) {return fabs (Getcross (b-a, p-a)/GetLength (B-A));} Double getdistancetosegment (point P, point A, point B) {if (a = = b) return getlength (P-A); Vector v1 = b-a, v2 = p-a, V3 = p-b;if (dcmp (Getdot, V1) < 0) return v2 (getlength); else if (V2 (dcmp (Getdot, V 3)) > 0) return GetLength (V3), Else return fabs (Getcross (v1, v2)/GetLength (v1));} /* point P on the line ab projection */point getpointtoline (points P, a, point B) {Vector v = b-a; return a+v* (Getdot (V, p-a)/Getdot (V, v)); }/* to determine if a segment has an intersection */bool haveintersection (Point A1, Point A2, point B1, point B2) {double C1=getcross (a2-a1, b1-a1), c2=getc Ross (A2-A1,B2-A1), C3=getcross (B2-B1, A1-B1), C4=getcross (B2-B1,A2-B1), return dcmp (C1) *dcmp (C2) < 0 && dcmp (C3) *dcmp ( C4) < 0;} /* Determine if the point is */bool onsegment on a line segment {return dcmp (Getcross (A-p, b-p)) = = 0 && dcmp (Getdot (A- P, b-p)) < 0; }}namespace Triangular {using namespace vectorial;double getangle (double A, double b, double c) {return ACOs (a*a+b*b-c* c)/(2*A*B)); }double Getarea (Double A, double b, double c) {double s = (a+b+c)/2; return sqrt (s* (s-a) * (s-b) * (S-C));} Double Getarea (Double A, double h) {return a * H/2;} Double Getarea (Point A, point B, point C) {return fabs (Getcross (b-a, c-a))/2;} Double Getdirarea (Point A, point B, point C) {return Getcross (b-a, c-a)/2;}};  Namespace polygonal {using namespace vectorial;using namespace linear;double getarea (point* p, int n) {double ret = 0;for (int i = 1; i < n-1; i++) ret + = Getcross (P[i]-p[0], p[i+1]-p[0]); return fabs (ret)/2;} /* Convex package */int getconvexhull (point* p, int n, point* ch) {sort (p, p + N); int m = 0;for (int i = 0; i < n; i++) {/* can be collinear *///while (M > 1 && dcmp (Getcross (ch[m-1 ]-ch[m-2], p[i]-ch[m-1])) < 0) M--;while (M > 1 && dcmp (Getcross (Ch[m-1]-ch[m-2], p[i]-ch[m-1])) <= 0) m -;ch[m++] = P[i];} int k = m;for (int i = n-2; I >= 0; i--) {/* can be collinear *///while (M > K && dcmp (Getcross (ch[m-1]-ch[m-2), p[i]-ch[ M-2])) < 0) M--;while (M > K && dcmp (Getcross (Ch[m-1]-ch[m-2], p[i]-ch[m-2]) <= 0) m--;ch[m++] = P[i];} if (n > 1) m--;return m;} int Ispointinpolygon (Point O, point* p, int n) {int wn = 0;for (int i = 0; i < n; i++) {int J = (i + 1)% n;if (onsegm ENT (o, p[i], p[j])) return 0; the boundary int k = dcmp (Getcross (P[j]-p[i], o-p[i])); int D1 = DCMP (p[i].y-o.y); int d2 = DCMP (P[J].Y-O.Y); if (k > 0 &am p;& D1 <= 0 && d2 > 0) wn++;if (k < 0 && D2 <= 0 && d1 > 0) wn--;} Return WN? -1:1;}}; Namespace Circular {using namespace linear;using namespace VecTorial;using namespace triangular;/* line and original intersection */int getlinecircleintersection (point P, point Q, Circle O, double& t1, double& T2, vector<point>& sol) {vector v = q-p;/* to be emptied before use Sol *///sol.clear ();d ouble a = v.x, B = p.x-o.o . x, C = v.y, d = p.y-o.o.y;double e = a*a+c*c, F = (a*b+c*d), g = b*b+d*d-o.r*o.r;double Delta = f*f-4*e*g;if (DCMP ( Delta) < 0) return 0;if (dcmp (delta) = = 0) {T1 = t2 =-F/(2 * e); Sol.push_back (P + v * t1); return 1;} T1 = (-f-sqrt (delta))/(2 * e); Sol.push_back (p + v * t1); t2 = (-f + sqrt (delta))/(2 * e); Sol.push_back (p + v * T2); return 2;} /* Circle intersection */int getcirclecircleintersection (Circle O1, Circle O2, vector<point>& sol) {Double d = getlength (O1. O-O2.O); if (dcmp (d) = = 0) {if (dcmp (O1.R-O2.R) = = 0) Return-1;return 0;} if (dcmp (O1.R + o2.r-d) < 0) return 0;if (DCMP (fabs)-D) > 0) return O1.R-O2.R a = 0;double (Getangle) ;d ouble da = ACOs ((O1.R*O1.R + D*D-O2.R*O2.R)/(2*o1.r*d)); Point P1 = O1. Point (A-da), p2 = O1.point (A+da), Sol.push_back (p1), if (P1 = = p2) return 1;sol.push_back (P2); return 2;} /* tangents */int gettangents (point P, Circle O, vector* v) {Vector U = o.o-p;double d = getlength (u); if (d < O.R) R Eturn 0;else if (dcmp (D-O.R) = = 0) {V[0] = rotate (u, PI/2); return 1;} else {double ang = asin (O.R/D); v[0] = rotate (u ,-ang); v[1] = rotate (U, ANG); return 2;}} /* A[i] and b[i] are the tangency points of the section I tangent on O1 and O2 respectively */int gettangents (Circle O1, Circle O2, point* A, point* b) {int cnt = 0;IF (O1.R < O2.R) {swap (O1, O2); swap (A, b);} Double D2 = GetLength (O1.O-O2.O); D2 = D2 * d2;double rdif = O1.R-O2.R, rsum = O1.R + o2.r;if (D2 < RDIF * Rdif) return 0;if (dcmp (d2) = = 0 && DCMP (O1.R-O2.R) = = 0) return-1;double base = Getangle (O2.O-O1.O); if (dcmp (D2-RDIF * rdif) = = 0) {a[cnt] = O1.point ( Base); B[CNT] = O2.point (base); Cnt++;return CNT;} Double ang = ACOs ((O1.R-O2.R)/sqrt (D2)); a[cnt] = O1.point (Base+ang); B[CNT] = O2.point (Base+ang); CNT++;A[CNT] = O1.poiNT (Base-ang); B[CNT] = O2.point (Base-ang); Cnt++;if (dcmp (d2-rsum * rsum) = = 0) {a[cnt] = O1.point (base); b[cnt] = O2.point (base); cnt++;} else if (D2 > Rsum * rsum) {Double ang = ACOs ((O1.R + O2.R)/sqrt (D2)); a[cnt] = O1.point (Base+ang); b[cnt] = O2.point (Base+ang); cnt++;a[cnt ] = O1.point (Base-ang); B[CNT] = O2.point (Base-ang); cnt++;} return CNT;}  /* Three points to determine the tangent circle */circle circumscribedcircle (Point P1, Spot, point P3) {Double Bx = p2.x-p1.x, by = P2.Y-P1.Y;  Double Cx = p3.x-p1.x, Cy = p3.y-p1.y;  Double D = 2 * (Bx * cy-by * Cx);  Double CX = (Cy * (BX * bx + by *) – by * (CX * CX + cy * cy))/D + p1.x;  Double cy = (BX * (CX * CX + cy * CY)-CX * (BX * bx + by *))/D + p1.y;  Point P = Point (CX, CY);  Return Circle (P, GetLength (p1-p));  }/* three to determine inscribed circle */circle inscribedcircle (Point P1, points P2, dot p3) {Double A = GetLength (P2-P3);  Double b = getlength (P3-P1);  Double c = getlength (P1-P2);  Point P = (P1 * A + P2 * b + p3 * c)/(A + B + c); ReTurn Circle (P, Getdistancetoline (P, p1, p2));  }/* Triangle vertex is center */double getpublicareatotriangle (Circle O, point A, point B) {if (dcmp ((A-O.O) * (b-o.o)) = = 0) return 0;int sig = 1;double da = getplength (o.o-a), db = Getplength (O.o-b), if (dcmp (DA-DB) > 0) {swap (DA, db), swap (A, b); sig =-1;} Double T1, t2;vector<point> sol;int n = getlinecircleintersection (A, B, O, T1, T2, Sol); if (dcmp (DA-O.R*O.R) <= 0 {if (dcmp (DB-O.R*O.R) <= 0) return Getdirarea (O.O, A, b) * Sig;int k = 0;if (Getplength (sol[0]-b) > Getplength (sol[ 1]-b)) K = 1;double ret = Getarea (o.o, A, sol[k]) + O.getarea (Getangle (SOL[K]-O.O, b-o.o));d ouble tmp = (a-o.o) * (B-O.O); re TURN ret * SIG * DCMP (TMP);}  Double d = getdistancetosegment (o.o, A, B), if (dcmp (D-O.R) >= 0) {DOUBLE ret = O.getarea (Getangle (A-O.O, b-o.o));d ouble TMP = (A-O.O) * (B-O.O); return ret * sig * DCMP (TMP);} Double K1 = o.r/getlength (a-o.o), K2 = O.r/getlength (B-O.O); Point P = o.o + (a-o.o) * K1, q = o.o + (b-o.o) * k2;double Ret1 = O. Getarea (Getangle (P-O.O, q-o.o));d ouble Ret2 = O.getarea (Getangle (SOL[0]-O.O, SOL[1]-O.O))-Getarea (O.O, sol[0], sol[1 ]);d ouble ret = (ret1-ret2), TMP = (A-O.O) * (B-O.O); return ret * sig * DCMP (TMP);} Double Getpublicareatopolygon (Circle O, point* p, int n) {if (dcmp (O.R) = = 0) return 0;double area = 0;for (int i = 0; i < n; i++) {int u = (i + 1)% n;area + = Getpublicareatotriangle (O, P[i], p[u]);} return Fabs (area);}}; Using namespace Vectorial;using namespace circular;using namespace Polygonal;const int maxn = 5005;const double inf = 1e4; int N; Point P[MAXN], C;int main () {int cas;scanf ('%d ', &cas); for (int kcas = 1; kcas <= cas; kcas++) {Double x, y, K;sca NF ("%d", &n), for (int i = 0; i < N; i++) P[i].read () scanf ("%lf%lf%lf", &x, &y, &k);d ouble S = Getarea ( P, N) * k/100;double L = 0, R = Inf;while (dcmp (r-l) > 0) {Double mid = (L + r)/2;double area = Getpublicareatopo Lygon (Circle (Point (x, Y), mid), P, N), if (dcmp (area-s) >= 0) r = Mid;elsel =Mid;} Double ans = (L + R)/2;printf ("Case%d:%.0lf\n", kcas, ans);} return 0;}

Copyright NOTICE: This article for Bo Master original article, without Bo Master permission not reproduced.

UVA 11978-fukushima Nuclear Blast (two points + geometry)