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
",ans);
498     return 0;
499 }
NOI十连测 第六测 T3

NOI十连测第六测 T3
