NOI十连测 第六测 T3
思路:考试的时候我非常地**，写了圆并，然后还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 }
