牛客网暑期ACM多校训练营（第三场）J 多边形与圆相交的面积

  1 #include <bits/stdc++.h>
  2 #define LL long long
  3 #define PI 3.1415926535897932384626
  4 #define maxn 1000
  5 #define EXIT exit(0);
  6 #define DEBUG puts("Here is a BUG");
  7 #define CLEAR(name, init) memset(name, init, sizeof(name))
  8 const double eps = 1e-9;
  9 const int MAXN = (int)1e9 + 5;
 10 using namespace std;
 11 #define Vector Point
 12 int dcmp(double x) { return fabs(x) < eps ? 0 : (x < 0 ? -1 : 1); }
 13 struct Point {
 14     double x, y;
 15  
 16     Point(const Point& rhs): x(rhs.x), y(rhs.y) { } //拷贝构造函数
 17     Point(double x = 0.0, double y = 0.0): x(x), y(y) { }   //构造函数
 18  
 19     friend istream& operator >> (istream& in, Point& P) { return in >> P.x >> P.y; }
 20     friend ostream& operator << (ostream& out, const Point& P) { return out << P.x << ‘ ‘ << P.y; }
 21  
 22     friend Vector operator + (const Vector& A, const Vector& B) { return Vector(A.x+B.x, A.y+B.y); }
 23     friend Vector operator - (const Point& A, const Point& B) { return Vector(A.x-B.x, A.y-B.y); }
 24     friend Vector operator * (const Vector& A, const double& p) { return Vector(A.x*p, A.y*p); }
 25     friend Vector operator / (const Vector& A, const double& p) { return Vector(A.x/p, A.y/p); }
 26     friend bool operator == (const Point& A, const Point& B) { return dcmp(A.x-B.x) == 0 && dcmp(A.y-B.y) == 0; }
 27     friend bool operator < (const Point& A, const Point& B) { return A.x < B.x || (A.x == B.x && A.y < B.y); }
 28  
 29     void in(void) { scanf("%lf%lf", &x, &y); }
 30     void out(void) { printf("%lf %lf", x, y); }
 31 };
 32  
 33 template <class T> T sqr(T x) { return x * x;}
 34 double Dot(const Vector& A, const Vector& B) { return A.x*B.x + A.y*B.y; }  //点积
 35 double Length(const Vector& A){ return sqrt(Dot(A, A)); }
 36 double Angle(const Vector& A, const Vector& B) { return acos(Dot(A, B)/Length(A)/Length(B)); }  //向量夹角
 37 double Cross(const Vector& A, const Vector& B) { return A.x*B.y - A.y*B.x; }    //叉积
 38 double Area(const Point& A, const Point& B, const Point& C) { return fabs(Cross(B-A, C-A)); }
 39 Vector normal(Vector x) { return Point(-x.y, x.x) / Length(x);}
 40 double angle(Vector x) { return atan2(x.y, x.x);}
 41  
 42 Vector vecunit(Vector x){ return x / Length(x);} //单位向量
 43 struct Circle {
 44     Point c;    //圆心
 45     double r;   //半径
 46  
 47     Circle() { }
 48     Circle(const Circle& rhs): c(rhs.c), r(rhs.r) { }
 49     Circle(const Point& c, const double& r): c(c), r(r) { }
 50  
 51     Point point(double ang) const { return Point(c.x + cos(ang)*r, c.y + sin(ang)*r); } //圆心角所对应的点
 52     double area(void) const { return PI * r * r; }
 53 };
 54 struct Line {
 55     Point P;    //直线上一点
 56     Vector dir; //方向向量(半平面交中该向量左侧表示相应的半平面)
 57     double ang; //极角，即从x正半轴旋转到向量dir所需要的角（弧度）
 58  
 59     Line() { }  //构造函数
 60     Line(const Line& L): P(L.P), dir(L.dir), ang(L.ang) { }
 61     Line(const Point& P, const Vector& dir): P(P), dir(dir) { ang = atan2(dir.y, dir.x); }
 62  
 63     bool operator < (const Line& L) const { //极角排序
 64         return ang < L.ang;
 65     }
 66  
 67     Point point(double t) { return P + dir*t; }
 68 };
 69  
 70 bool InCircle(Point x, Circle c) { return dcmp(c.r*c.r - Length(c.c - x)*Length(c.c - x)) >= 0;}
 71 Point GetIntersection(Line a, Line b) //线段交点
 72 {
 73     Vector u = a.P-b.P;
 74     double t = Cross(b.dir, u) / Cross(a.dir, b.dir);
 75     return a.P + a.dir*t;
 76 }
 77  
 78 bool OnSegment(Point p, Point a1, Point a2)
 79 {
 80     return dcmp(Cross(a1-p, a2-p)) == 0 && dcmp(Dot(a1-p, a2-p)) < 0;
 81 }
 82 int getSegCircleIntersection(Line L, Circle C, Point* sol)
 83 {
 84     Vector nor = normal(L.dir);
 85     Line pl = Line(C.c, nor);
 86     Point ip = GetIntersection(pl, L);
 87     double dis = Length(ip - C.c);
 88     if (dcmp(dis - C.r) > 0) return 0;
 89     Point dxy = vecunit(L.dir) * sqrt(C.r*C.r - dis*dis);
 90     int ret = 0;
 91     sol[ret] = ip + dxy;
 92     if (OnSegment(sol[ret], L.P, L.point(1))) ret++;
 93     sol[ret] = ip - dxy;
 94     if (OnSegment(sol[ret], L.P, L.point(1))) ret++;
 95     return ret;
 96 }
 97  
 98 double SegCircleArea(Circle C, Point a, Point b) //线段切割圆
 99 {
100     double a1 = angle(a - C.c);
101     double a2 = angle(b - C.c);
102     double da = fabs(a1 - a2);
103     if (da > PI) da = PI * 2.0 - da;
104     return dcmp(Cross(b - C.c, a - C.c)) * da * sqr(C.r) / 2.0;
105 }
106  
107 double PolyCiclrArea(Circle C, Point *p, int n)//多边形与圆相交面积
108 {
109     double ret = 0.0;
110     Point sol[2];
111     p[n] = p[0];
112     for(int i=0;i<n;i++)
113     {
114         double t1, t2;
115         int cnt = getSegCircleIntersection(Line(p[i], p[i+1]-p[i]), C, sol);
116         if (cnt == 0)
117         {
118             if (!InCircle(p[i], C) || !InCircle(p[i+1], C)) ret += SegCircleArea(C, p[i], p[i+1]);
119             else ret += Cross(p[i+1] - C.c, p[i] - C.c) / 2.0;
120         }
121         if (cnt == 1)
122         {
123             if (InCircle(p[i], C) && !InCircle(p[i+1], C)) ret += Cross(sol[0] - C.c, p[i] - C.c) / 2.0, ret += SegCircleArea(C, sol[0], p[i+1]);
124             else ret += SegCircleArea(C, p[i], sol[0]), ret += Cross(p[i+1] - C.c, sol[0] - C.c) / 2.0;
125         }
126         if (cnt == 2)
127         {
128             if ((p[i] < p[i + 1]) ^ (sol[0] < sol[1])) swap(sol[0], sol[1]);
129             ret += SegCircleArea(C, p[i], sol[0]);
130             ret += Cross(sol[1] - C.c, sol[0] - C.c) / 2.0;
131             ret += SegCircleArea(C, sol[1], p[i+1]);
132         }
133     }
134     return fabs(ret);
135 }
136 double PolygonArea(Point *po, int n) {
137     double area = 0.0;
138     for(int i = 1; i < n-1; i++) {
139         area += Cross(po[i]-po[0], po[i+1]-po[0]);
140     }
141     return area * 0.5;
142 }
143 Point a[210], b;
144 double p, q;
145 int main(){
146     int n, m;
147     scanf("%d", &n);
148     for(int i=0;i<n;i++) a[i].in();
149     double ar = fabs(PolygonArea(a, n));
150     scanf("%d", &m);
151     for(int i=0;i<m;i++){
152         b.in();
153         scanf("%lf%lf", &p, &q);
154         double l = 0, r = 10000, num = 50, mid, aa = 1-p/q;
155         while(num--){
156             mid = (l+r)/2;
157             Circle yuan(b, mid);
158             if(PolyCiclrArea(yuan, a, n)/ar < aa) l = mid;
159             else r = mid;
160         }
161         printf("%.10f\n", mid);
162     }
163     return 0;
164 }