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

Source: Internet
Author: User
Tags acos

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)

Contact Us

The content source of this page is from Internet, which doesn't represent Alibaba Cloud's opinion; products and services mentioned on that page don't have any relationship with Alibaba Cloud. If the content of the page makes you feel confusing, please write us an email, we will handle the problem within 5 days after receiving your email.

If you find any instances of plagiarism from the community, please send an email to: info-contact@alibabacloud.com and provide relevant evidence. A staff member will contact you within 5 working days.

A Free Trial That Lets You Build Big!

Start building with 50+ products and up to 12 months usage for Elastic Compute Service

  • Sales Support

    1 on 1 presale consultation

  • After-Sales Support

    24/7 Technical Support 6 Free Tickets per Quarter Faster Response

  • Alibaba Cloud offers highly flexible support services tailored to meet your exact needs.