思路:考试的时候我非常地**，写了圆并，然后还TM写了半平面交和三角剖分，虽然只有30分。。但是看在我写了500行的份上还是挂着吧。。
1 #include<cstdio>
2 #include<cmath>
3 #include<cstring>
4 #include<iostream>
5 #include<algorithm>
6 const double Pi=acos(-1);
7 const double eps=1e-10;
8 int n,m,cnt,tot,num;
9 bool pd[200005],check;
10 int read(){
11 int t=0,f=1;char ch=getchar();
12 while (ch<‘0‘||ch>‘9‘) {if (ch==‘-‘) f=-1;ch=getchar();}
13 while (‘0‘<=ch&&ch<=‘9‘) {t=t*10+ch-‘0‘;ch=getchar();}
14 return t*f;
15 }
16 struct node{
17 double l,r;
18 }L[200005];
19 struct Point{
20 double x,y;
21 Point(){}
22 Point(double x0,double y0):x(x0),y(y0){}
23 }p[200005],c[200005],d[200005];
24 struct Line{
25 Point s,e;
26 double slop;
27 Line(){}
28 Line(Point s0,Point e0):s(s0),e(e0){}
29 }l[200005],q[200005];
30 struct Circle{
31 Point p;
32 double r;
33 Circle(){}
34 Circle(double x0,double y0,double r0):p(Point(x0,y0)),r(r0){}
35 }a[200005];
36 bool cmp2(node a,node b){
37 if (fabs(a.l-b.l)<eps) return a.r<b.r;
38 else return (a.l<b.l);
39 }
40 int sgn(double x){if (x>eps) return 1;if (x<-eps) return -1;return 0;}
41 double operator *(Point p1,Point p2){return p1.x*p2.y-p1.y*p2.x;}
42 double operator /(Point p1,Point p2){return p1.x*p2.x+p1.y*p2.y;}
43 Point operator -(Point p1,Point p2){return Point(p1.x-p2.x,p1.y-p2.y);}
44 Point operator +(Point p1,Point p2){return Point(p1.x+p2.x,p1.y+p2.y);}
45 Point operator *(Point p1,double x){return Point(p1.x*x,p1.y*x);}
46 Point operator /(Point p1,double x){return Point(p1.x/x,p1.y/x);}
47 double sqr(double x){return x*x;}
48 double llen(Point p1){return sqrt(sqr(p1.x)+sqr(p1.y)); }
49 double dis(Point p1,Point p2){return llen(p1-p2);}
50 double dis(Point p1){return llen(p1);}
51 Point evector(Point p1){double t=llen(p1);if (t<eps) return Point(0,0);else return p1/t;}
52 bool cmp3(Line p1,Line p2){
53 if (p1.slop!=p2.slop) return p1.slop<p2.slop;
54 else return (p1.e-p1.s)*(p2.e-p1.s)>0;
55 }
56 bool cmp(Point p1,Point p2){
57 double t=(p1-p[1])*(p2-p[1]);
58 if (fabs(t)<eps) return dis(p[1],p1)<dis(p[1],p2);
59 return t>0;
60 }
61 bool conclude(Circle p1,Circle p2){
62 double Len=dis(p1.p,p2.p);
63 if (fabs(p1.r-p2.r)>=Len) return 1;
64 return 0;
65 }
66 bool intersect(Circle p1,Circle p2){
67 double Len=dis(p1.p,p2.p);
68 if (fabs(p1.r-p2.r)<=Len&&Len<=p1.r+p2.r) return 1;
69 return 0;
70 }
71 double dist_line(Line p){
72 double A,B,C,dist;
73 A=p.s.y-p.e.y;
74 B=p.s.x-p.e.x;
75 C=p.s.x*p.e.y-p.s.y*p.e.x;
76 dist=fabs(C)/sqrt(sqr(A)+sqr(B));
77 return dist;
78 }
79 double dist_line(Point p1,Line p){
80 p.s.x-=p1.x;
81 p.e.x-=p1.x;
82 p.s.y-=p1.y;
83 p.e.y-=p1.y;
84 return dist_line(p);
85 }
86 double get_cos(double a,double b,double c){
87 return (b*b+c*c-a*a)/(2*b*c);
88 }
89 double get_angle(Point p1,Point p2){
90 if (!sgn(dis(p1))||!sgn(dis(p2))) return 0.0;
91 double A,B,C;
92 A=dis(p1);
93 B=dis(p2);
94 C=dis(p1,p2);
95 if (C<=eps) return 0.0;
96 return fabs(acos(get_cos(C,A,B)));
97 }
98 Point get_point(Point p){
99 double T=sqr(p.x)+sqr(p.y);
100 return Point(sgn(p.x)*sqrt(sqr(p.x)/T),sgn(p.y)*sqrt(sqr(p.y)/T));
101 }
102 bool bigconclude(Circle p1,Circle p2,Point p){
103 Point e1=p2.p-p1.p;
104 Point e2=p-p1.p;
105 if (sgn(e1/e2)<0) return 1;
106 else return 0;
107 }
108 double S(Point p1,Point p2,Point p3){
109 return fabs((p2-p1)*(p3-p1))/2;
110 }
111 double work(Point p1,Point p2,double R){
112 Point O=Point(0,0);
113 double f=sgn(p1*p2),res=0;
114 if (!sgn(f)||!sgn(dis(p1))||!sgn(dis(p2))) return 0.0;
115 double l=dist_line(Line(p1,p2));
116 double a=dis(p1);
117 double b=dis(p2);
118 double c=dis(p1,p2);
119 if (a<=R&&b<=R){
120 return fabs(p1*p2)/2.0*f;
121 }
122 if (a>=R&&b>=R&&l>=R){
123 double ang=get_angle(p1,p2);
124 return fabs((ang/(2.0))*(R*R))*f;
125 }
126 if (a>=R&&b>=R&&l<=R&&(get_cos(a,b,c)<=0||get_cos(b,a,c)<=0)){
127 double ang=get_angle(p1,p2);
128 return fabs((ang/(2.0))*(R*R))*f;
129 }
130 if (a>=R&&b>=R&&l<=R&&(get_cos(a,b,c)>0&&get_cos(b,a,c)>0)){
131 double dist=dist_line(Line(p1,p2));
132 double len=sqrt(sqr(R)-sqr(dist))*2.0;
133 double ang1=get_angle(p1,p2);
134 double cos2=get_cos(len,R,R);
135 res+=fabs(len*dist/2.0);
136 double ang2=ang1-acos(cos2);
137 res+=fabs((ang2/(2))*(R*R));
138 return res*f;
139 }
140 if ((a>=R&&b<R)||(a<R&&b>=R)){
141 if (b>a) std::swap(a,b),std::swap(p1,p2);
142 double T=sqr(p1.x-p2.x)+sqr(p1.y-p2.y);
143 Point u=Point(sgn(p1.x-p2.x)*sqrt(sqr(p1.x-p2.x)/T),sgn(p1.y-p2.y)*sqrt(sqr(p1.y-p2.y)/T));
144 double dist=dist_line(Line(p1,p2));
145 double len=sqrt(R*R-dist*dist);
146 double len2=sqrt(sqr(dis(p2))-sqr(dist));
147 if (fabs(dis(p2+u*len2)-dist)<=eps) len+=len2;
148 else len-=len2;
149 Point p=p2+u*len;
150 res+=S(O,p2,p);
151 double ang=get_angle(p1,p);
152 res+=fabs((ang/2.0)*R*R);
153 return res*f;
154 }
155 return 0;
156 }
157 Point inter(Line p1,Line p2){
158 double t1=(p2.e-p1.s)*(p1.e-p1.s);
159 double t2=(p1.e-p1.s)*(p2.s-p1.s);
160 if (fabs(t1+t2)<eps) check=1;
161 double k=t1/(t1+t2);
162 Point p;
163 p.x=p2.e.x+(p2.s.x-p2.e.x)*k;
164 p.y=p2.e.y+(p2.s.y-p2.e.y)*k;
165 return p;
166 }
167 bool jud(Line p1,Line p2,Line p3){
168 Point p=inter(p1,p2);
169 return (p3.e-p3.s)*(p-p3.s)<0;
170 }
171 double inter(Circle P){
172 double res=0;p[n+1]=p[1];
173 for (int i=1;i<=n+1;i++)
174 c[i].x=p[i].x-P.p.x,c[i].y=p[i].y-P.p.y;
175 for (int i=1;i<=n;i++)
176 res+=work(c[i],c[i+1],P.r);
177 return res;
178 }
179 void hpi(){
180 int L=1,R=0;tot=0;
181 for (int i=1;i<=cnt;i++){
182 if (l[i].slop!=l[i-1].slop) tot++;
183 l[tot]=l[i];
184 }
185 q[++R]=l[1];q[++R]=l[2];
186 int all=tot;
187 for (int i=3;i<=all;i++){
188 while (L<R&&jud(q[R],q[R-1],l[i])) R--;
189 while (L<R&&jud(q[L],q[L+1],l[i])) L++;
190 q[++R]=l[i];
191 }
192 while (L<R&&jud(q[R],q[R-1],q[L])) R--;
193 while (L<R&&jud(q[L],q[L+1],q[R])) L++;
194 tot=0;
195 for (int i=L;i<R;i++)
196 d[++tot]=inter(q[i],q[i+1]);
197 }
198 bool judge(Point p1,Point p2){
199 return (fabs(p1.x-p2.x)<eps&&fabs(p1.y-p2.y)<eps);
200 }
201 double gworking1(Point c1,Point c2,Point p1,Point p2){
202 double rev=sgn(c1*c2);
203 if (fabs(rev)<eps) return 0;
204 num=0;
205 Point O=Point(0,0);
206 if ((c1*c2)<0) std::swap(c1,c2);
207 l[++num].s=O,l[num].e=c1;if (judge(l[num].s,l[num].e)) num--;
208 l[++num].s=c1,l[num].e=c2;if (judge(l[num].s,l[num].e)) num--;
209 l[++num].s=c2,l[num].e=O;if (judge(l[num].s,l[num].e)) num--;
210 if ((p1*p2)<0) std::swap(p1,p2);
211 l[++num].s=O,l[num].e=p1;if (judge(l[num].s,l[num].e)) num--;
212 l[++num].s=p1,l[num].e=p2;if (judge(l[num].s,l[num].e)) num--;
213 l[++num].s=p2,l[num].e=O;if (judge(l[num].s,l[num].e)) num--;
214 for (int i=1;i<=num;i++)
215 l[i].slop=atan2(l[i].e.y-l[i].s.y,l[i].e.x-l[i].s.x);
216 std::sort(l+1,l+1+num,cmp3);
217 check=0;
218 hpi();
219 if (check) return 0;
220 d[tot+1]=d[1];
221 double res=0;
222 for (int i=1;i<=tot;i++)
223 res+=(d[i]*d[i+1])/2.0;
224 return fabs(res)*rev;
225 }
226 double Gworking1(double a1,double a2,Circle w){
227 double r=w.r;
228 Point p1=Point(r*cos(a1),r*sin(a1));
229 Point p2=Point(r*cos(a2),r*sin(a2));
230 p[n+1]=p[1];
231 for (int i=1;i<=n+1;i++)
232 c[i]=Point(p[i].x-w.p.x,p[i].y-w.p.y);
233 double res=0;
234 for (int i=1;i<=n;i++)
235 res+=gworking1(c[i],c[i+1],p1,p2);
236 return res;
237 }
238 Point inter(Circle p,Line l,Point p1){
239 double Dist=dist_line(l);
240 double dis1=llen(p1);
241 double len=sqrt(sqr(dis1)-sqr(Dist));
242 Point e=evector(l.e-l.s);
243 Point sx=p1;
244 double R=p.r;
245 double len2=sqrt(sqr(R)-sqr(Dist));
246 sx=sx-e*len;
247 sx=sx+e*len2;
248 return sx;
249 }
250 double gworking2(Point c1,Point c2,Point p1,Point p2,double R){
251 double rev=sgn(c1*c2),res=0;Point O=Point(0,0);
252 if ((p1*p2)<0) std::swap(p1,p2);
253 if ((c1*c2)<0) std::swap(c1,c2);
254 if (fabs(rev)<eps) return 0;
255 if (!sgn(dis(p1))||!sgn(dis(p2))) return 0.0;
256 double l=dist_line(Line(p1,p2));
257 double a=dis(p1);
258 double b=dis(p2);
259 double c=dis(p1,p2);
260 if (a<=R&&b<=R){
261 res=fabs(p1*p2)/2.0;
262 }
263 else
264 if (a>=R&&b>=R&&l>=R){
265 double ang=get_angle(p1,p2);
266 res=abs((ang/(2.0))*(R*R));
267 }
268 else
269 if (a>=R&&b>=R&&l<=R&&(get_cos(a,b,c)<=0||get_cos(b,a,c)<=0)){
270 double ang=get_angle(p1,p2);
271 res=fabs((ang/(2.0))*(R*R));
272 }
273 else
274 if (a>=R&&b>=R&&l<=R&&(get_cos(a,b,c)>0&&get_cos(b,a,c)>0)){
275 double dist=dist_line(Line(p1,p2));
276 double len=sqrt(sqr(R)-sqr(dist))*2.0;
277 double ang1=get_angle(p1,p2);
278 double cos2=get_cos(len,R,R);
279 res+=fabs(len*dist/2.0);
280 double ang2=ang1-acos(cos2);
281 res+=fabs((ang2/(2))*(R*R));
282 }
283 else
284 if ((a>=R&&b<R)||(a<R&&b>=R)){
285 if (b>a) std::swap(a,b),std::swap(p1,p2);
286 double T=sqr(p1.x-p2.x)+sqr(p1.y-p2.y);
287 Point u=Point(sgn(p1.x-p2.x)*sqrt(sqr(p1.x-p2.x)/T),sgn(p1.y-p2.y)*sqrt(sqr(p1.y-p2.y)/T));
288 double dist=dist_line(Line(p1,p2));
289 double len=sqrt(R*R-dist*dist);
290 double len2=sqrt(sqr(dis(p2))-sqr(dist));
291 if (fabs(dis(p2+u*len2)-dist)<=eps) len+=len2;
292 else len-=len2;
293 Point p=p2+u*len;
294 res+=S(O,p2,p);
295 double ang=get_angle(p1,p);
296 res+=fabs((ang/2.0)*R*R);
297 }
298 int in1=((sgn(c1*p1)*sgn(p1*c2))>=0);
299 int in2=((sgn(c1*p2)*sgn(p2*c2))>=0);
300 if (in1&&in2) return res*rev;
301 if (!in1){
302 if (dis(p1)<=R){
303 Point sb=inter(Line(O,c1),Line(p1,p2));
304 double Di=dis(sb);
305 if (Di<=R) res-=S(O,sb,p1);
306 else{
307 Point sx=inter(Circle(0,0,R),Line(p1,p2),p1);
308 double ang=get_angle(c1,sx);
309 res-=fabs((ang/2.0)*R*R);
310 res-=S(O,sx,p1);
311 }
312 }else{
313 Point sb=inter(Line(O,c1),Line(p1,p2));
314 double Di=dis(sb);
315 if (Di>=R){
316 double ang=get_angle(c1,p1);
317 res-=fabs((ang/2.0)*R*R);
318 }else{
319 Point sx=inter(Circle(0,0,R),Line(p1,p2),p1);
320 double ang2=get_angle(p1,sx);
321 res-=fabs((ang2/2.0)*R*R);
322 res-=S(sb,sx,O);
323 }
324 }
325 }
326 if (!in2){
327 if (dis(p2)<=R){
328 Point sb=inter(Line(O,c2),Line(p1,p2));
329 double Di=dis(sb);
330 if (Di<=R) res-=S(O,sb,p2);
331 else{
332 Point sx=inter(Circle(0,0,R),Line(p2,p1),p2);
333 double ang=get_angle(c2,sx);
334 res-=fabs((ang/2.0)*R*R);
335 res-=S(O,sx,p2);
336 }
337 }else{
338 Point sb=inter(Line(O,c2),Line(p1,p2));
339 double Di=dis(sb);
340 if (Di>=R){
341 double ang=get_angle(c2,p2);
342 res-=fabs((ang/2.0)*R*R);
343 }else{
344 Point sx=inter(Circle(0,0,R),Line(p2,p1),p2);
345 double ang2=get_angle(p2,sx);
346 res-=fabs((ang2/2.0)*R*R);
347 res-=S(sb,sx,O);
348 }
349 }
350 }
351 return res*rev;
352 }
353 double Gworking2(double a1,double a2,Circle w){
354 double res=0,r=w.r;
355 Point p1=Point(r*cos(a1),r*sin(a1));
356 Point p2=Point(r*cos(a2),r*sin(a2));
357 p[n+1]=p[1];
358 for (int i=1;i<=n+1;i++)
359 c[i]=Point(p[i].x-w.p.x,p[i].y-w.p.y);
360 for (int i=1;i<=n;i++)
361 res+=gworking2(c[i],c[i+1],p1,p2,r);
362 return res;
363 }
364 void graham(){
365 int k=1;
366 for (int i=2;i<=n;i++)
367 if (p[i].y<p[k].y||(p[i].x<p[k].x&&p[i].y==p[k].y)) k=i;
368 std::sort(p+2,p+1+n,cmp);
369 int top=1;c[top]=p[1];
370 for (int i=2;i<=n;i++){
371 while (top>1&&(c[top]-c[top-1])*(p[i]-c[top])<=0) top--;
372 c[++top]=p[i];
373 }
374 n=top;
375 for (int i=1;i<=top;i++)
376 p[i]=c[i];
377 }
378 void Work(Circle p1,Point a1,Point a2){
379 Point O=p1.p;double r=p1.r;
380 a1.x-=O.x;a2.x-=O.x;
381 a1.y-=O.y;a2.y-=O.y;
382 double x1=atan2(a1.y,a1.x),x2=atan2(a2.y,a2.x);
383 if (x1<0) x1+=2*Pi;
384 if (x2<0) x2+=2*Pi;
385 if (x1<=x2){
386 cnt++;
387 L[cnt].l=x1;
388 L[cnt].r=x2;
389 }else{
390 cnt++;
391 L[cnt].l=x1;
392 L[cnt].r=2*Pi;
393 cnt++;
394 L[cnt].l=0.0;
395 L[cnt].r=x2;
396 }
397 }
398 void add(Circle p1,Circle p2){
399 double x1=p1.p.x,x2=p2.p.x,y1=p1.p.y,y2=p2.p.y,r1=p1.r,r2=p2.r;
400 double A=-2*(x1+x2),B=-2*(y1+y2),C=sqr(x1)+sqr(x2)+sqr(y1)+sqr(y2)-sqr(r1)-sqr(r2);
401 Point S,E;
402 if (fabs(A)<eps){
403 S=Point(1,(C)/B);
404 E=Point(2,(C)/B);
405 }else
406 if (fabs(B)<eps){
407 S=Point((C)/A,1);
408 E=Point((C)/A,2);
409 }else{
410 S=Point(1,(C-A)/B);
411 E=Point(2,(C-2*A)/B);
412 }
413 Line l=Line(S,E);
414 double Dist1=dist_line(p1.p,l);
415 double Dist2=llen(p1.p-S);
416 Point e=evector(E-S);
417 double Dist3=sqrt(sqr(Dist2)-sqr(Dist1));
418 e=e*Dist3;
419 Point a1,a2;
420 if (fabs(dis(S+e,p1.p)-Dist1)<eps){
421 S=S+e;
422 }
423 else S=S-e;
424 e=e/Dist3;
425 double Dist4=sqrt(sqr(p1.r)-sqr(llen(S-p1.p)));
426 e=e*Dist4;
427 a1=S+e;
428 a2=S-e;
429 if ((a1-p1.p)*(a2-p1.p)<0) std::swap(a1,a2);
430 //a1放在a2的顺时针方向、
431 if (bigconclude(p1,p2,a1)) std::swap(a1,a2);
432 //a1到a2必须是逆时针。
433 Work(p1,a1,a2);
434 }
435 int main(){
436 freopen("aerolite.in","r",stdin);
437 freopen("aerolite.out","w",stdout);
438 scanf("%d",&m);
439 for (int i=1;i<=m;i++){
440 scanf("%lf%lf%lf",&a[i].p.x,&a[i].p.y,&a[i].r);
441 }
442 scanf("%d",&n);
443 for (int i=1;i<=n;i++){
444 scanf("%lf%lf",&p[i].x,&p[i].y);
445 }
446 graham();
447 for (int i=1;i<=m;i++) pd[i]=1;
448 for (int i=1;i<=m;i++){
449 for (int j=1;j<=m;j++)
450 if (i!=j){
451 if (conclude(a[i],a[j])) {pd[i]=0;break;}
452 }
453 }
454 double ans=0;
455 for (int i=1;i<=m;i++)
456 if (pd[i]){
457 bool mark=0;
458 for (int j=1;j<=m;j++)
459 if (i!=j&&pd[j]){
460 if (intersect(a[i],a[j])) {mark=1;break;}
461 }
462 if (!mark) ans+=inter(a[i]),pd[i]=0;
463 }
464 int j=0;
465 for (int i=1;i<=m;i++)
466 if (pd[i]) a[++j]=a[i];
467 for (int i=1;i<=m;i++) pd[i]=1;
468 for (int i=1;i<=m;i++){
469 cnt=0;
470 for (int j=1;j<=m;j++)
471 if (i!=j){
472 add(a[i],a[j]);
473 }
474 std::sort(L+1,L+1+cnt,cmp2);
475 L[0].r=0.0;
476 L[cnt+1].l=2*Pi;
477 double Reach=0.0,Last=0.0;
478 int ww=0;
479 for (int j=1;j<=cnt;j++)
480 if (L[j].l>Reach){
481 if (j==126) ww++;
482 ans+=Gworking1(Last,Reach,a[i]);
483 Last=L[j].l;
484 ans+=Gworking2(Reach,L[j].l,a[i]);
485 Reach=L[j].r;
486 }else{
487 Reach=L[j].r;
488 }
489 if (fabs(Reach-2*Pi)<eps){
490 ans+=Gworking1(Last,Reach,a[i]);
491 }else
492 if (fabs(Reach-2*Pi)>=eps){
493 ans+=Gworking1(Last,Reach,a[i]);
494 ans+=Gworking2(Reach,2*Pi,a[i]);
495 }
496 }
497 printf("%.8lf\n",ans);
498 return 0;
499 }
时间： 2024-11-05 13:32:10