Harry Potter and J.K.Rowling
题意：半平面交 ＋ 圆与多边形的面积并．
转换成圆与三角形的面积并，分四种情况讨论一下．
1 #include <bits/stdc++.h> 2 using namespace std; 3 const double pi = acos(-1.0); 4 const int maxn = 2010; 5 const int inf = 0x3f3f3f3f; 6 const double eps = 1e-8; 7 struct Point { 8 double x,y; 9 Point (double x = 0, double y = 0) : x(x), y(y) {} 10 }; 11 typedef Point Vector; 12 struct Circle{ 13 Point o; 14 double r; 15 Circle(){} 16 Circle(Point o, double r) : o(o), r(r){} 17 //Point point(double a){ 18 // return Point(o.x + cos(a)*r, o.y + sin(a)*r); 19 //} 20 }; 21 Vector operator + (Vector a, Vector b) { 22 return Vector (a.x + b.x, a.y + b.y); 23 } 24 Vector operator * (Vector a, double s) { 25 return Vector (a.x * s, a.y * s); 26 } 27 Vector operator / (Vector a, double p) { 28 return Vector (a.x / p, a.y / p); 29 } 30 Vector operator - (Point a, Point b) { 31 return Vector (a.x - b.x, a.y - b.y); 32 } 33 bool operator < (Point a, Point b) { 34 return a.x < b.x || (a.x == b.x && a.y < b.y); 35 } 36 int dcmp (double x) { 37 if(fabs(x) < eps) return 0; 38 return x < 0 ? -1 : 1; 39 } 40 bool operator == (const Point &a, const Point &b) { 41 return dcmp(a.x - b.x) == 0 && dcmp(a.y - b.y) == 0; 42 } 43 double Dot(Vector a, Vector b) { 44 return a.x * b.x + a.y * b.y; 45 } 46 double Angle(Vector a){ 47 return atan2(a.y, a.x); 48 } 49 double Length(Vector a) { 50 return sqrt(Dot(a, a)); 51 } 52 double Cross (Vector a, Vector b) { 53 return a.x * b.y - a.y * b.x; 54 } 55 //半平面交 56 struct Line { 57 Point p; 58 Vector v; 59 double rad; 60 Line () {} 61 Line (Point p, Vector v) : p(p), v(v) { 62 rad = atan2(v.y,v.x); 63 } 64 bool operator < (const Line &L) const { 65 return rad < L.rad; 66 } 67 }; 68 bool OnLeft(Line L, Point p) { 69 return Cross(L.v, p - L.p) > 0; 70 } 71 Point GetLineIntersection (Line a, Line b) { 72 Vector u = a.p - b.p; 73 double t = Cross(b.v, u) / Cross(a.v, b.v); 74 return a.p + a.v*t; 75 } 76 77 int HalfplaneIntersection(Line *L, int n,Point *poly) { 78 sort(L, L+n); 79 int first,last; 80 Point *p = new Point[n]; 81 Line *q = new Line[n]; //双端队列 82 q[first = last = 0] = L[0]; 83 for(int i = 1; i < n; i++) { 84 while(first < last && !OnLeft(L[i], p[last-1])) last--; //去尾 85 while(first < last && !OnLeft(L[i], p[first])) first++; 86 q[++last] = L[i]; 87 if(dcmp(Cross(q[last].v, q[last-1].v)) == 0) { 88 last--; 89 if(OnLeft(q[last], L[i].p)) q[last] = L[i]; 90 } 91 if(first < last) p[last-1] = GetLineIntersection (q[last-1],q[last]); 92 } 93 while(first < last && !OnLeft(q[first], p[last-1])) last--; //删除无用平面 94 if(last - first <= 1) { 95 delete []p; 96 delete []q; 97 return 0; //空集 98 } 99 p[last] = GetLineIntersection (q[last], q[first]); 100 int m = 0; 101 for(int i = first; i <= last; i++) poly[m++] = p[i]; 102 delete []p; 103 delete []q; 104 return m; 105 } 106 Point p[maxn][2], poly[maxn]; 107 Line line[maxn]; 108 Circle cr; 109 Point goal; 110 //线段AB与圆的交点 111 void LineIntersectionCircle(Point A, Point B, Circle cr, Point *p, int &num){ 112 double x0 = cr.o.x, y0 = cr.o.y; 113 double x1 = A.x, y1 = A.y, x2 = B.x, y2 = B.y; 114 double dx = x2 - x1, dy = y2 - y1; 115 double a = dx * dx + dy * dy; 116 double b = 2 * dx *(x1 - x0) + 2 * dy * (y1 -y0); 117 double c = (x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0) - cr.r *cr.r; 118 double delta = b*b - 4*a*c; 119 num = 0; 120 if(dcmp(delta) >= 0){ 121 double t1 = (-b - sqrt(delta)) / (2*a); 122 double t2 = (-b + sqrt(delta)) / (2*a); 123 if(dcmp(t1-1) <= 0 && dcmp(t1) >= 0) { 124 p[num++] = Point(x1 + t1*dx, y1 + t1*dy); 125 } 126 if(dcmp(t2-1) <= 0 && dcmp(t2) >= 0) { 127 p[num++] = Point(x1 + t2*dx, y1 + t2*dy); 128 } 129 } 130 } 131 //扇形面积 132 double SectorArea(Point a, Point b, Circle cr){ 133 double theta = Angle(a - cr.o) - Angle(b - cr.o); 134 while(theta <= 0) theta += 2*pi; 135 while(theta > 2*pi) theta -= 2*pi; 136 theta = min(theta, 2*pi - theta); 137 return cr.r * cr.r * theta / 2; 138 } 139 140 double cal(Point a, Point b, Circle cr){ 141 Point tp[3]; 142 int num = 0; 143 Vector ta = a - cr.o; 144 Vector tb = b - cr.o; 145 bool ina = (Length(ta) - cr.r) < 0; 146 bool inb = (Length(tb) - cr.r) < 0; 147 if(ina){ 148 if(inb){ 149 return fabs(Cross(ta, tb))/2; 150 } else{ 151 LineIntersectionCircle(a, b, cr, tp, num); 152 return SectorArea(b, tp[0], cr) + fabs(Cross(ta, tp[0] - cr.o))/2; 153 } 154 } else{ 155 if(inb){ 156 LineIntersectionCircle(a, b, cr, tp, num); 157 return SectorArea(a, tp[0], cr) + fabs(Cross(tb, tp[0] - cr.o))/2; 158 } else { 159 LineIntersectionCircle(a, b, cr, tp, num); 160 if(num == 2) { 161 return SectorArea(a, tp[0], cr) + SectorArea(b, tp[1], cr) + fabs(Cross(tp[0]-cr.o, tp[1]-cr.o))/2; 162 } else { 163 return SectorArea(a, b, cr); 164 } 165 } 166 } 167 } 168 //圆与多边形的面积并 169 double CirclePolyArea(Point *p, int n, Circle cr){ 170 double res = 0; 171 p[n] = p[0]; 172 for(int i = 0; i < n; i++){ 173 int sgn = dcmp(Cross(p[i] - cr.o, p[i+1] - cr.o)); 174 if(sgn){ 175 res += sgn * cal(p[i], p[i+1], cr); 176 } 177 } 178 return res; 179 } 180 181 182 int main(){ 183 //freopen("in.txt", "r", stdin); 184 int t, kase = 0; 185 scanf("%d", &t); 186 while(t--){ 187 int n; 188 double r; 189 scanf("%lf %d", &r, &n); 190 cr.r = r; 191 cr.o = Point(0, 0); 192 for(int i = 0; i < n; i++){ 193 for(int j = 0; j < 2; j++){ 194 scanf("%lf %lf", &p[i][j].x, &p[i][j].y); 195 } 196 } 197 scanf("%lf %lf", &goal.x, &goal.y); 198 for(int i = 0; i < n; i++){ 199 Line lp = Line(p[i][0], p[i][0] - p[i][1]); 200 if(OnLeft(lp, goal)){ 201 line[i] = lp; 202 }else { 203 line[i] = Line(p[i][0], p[i][1] - p[i][0]); 204 } 205 } 206 Point p1 = Point(-11111, -11111), p2 = Point(-11111, 11111); 207 Point p3 = Point(11111, 11111), p4 = Point(11111, -11111); 208 line[n] = Line(p1, p1-p2); 209 line[n+1] = Line(p2, p2-p3); 210 line[n+2] = Line(p3, p3-p4); 211 line[n+3] = Line(p4, p4-p1); 212 int sz = HalfplaneIntersection(line, n+4, poly); 213 double area = CirclePolyArea(poly, sz, cr); 214 printf("Case %d: %.5lf", ++kase, area/pi/r/r*100); 215 puts("%"); 216 } 217 return 0; 218 }