bzoj 2178
这题调精度真痛苦啊（向管理员要了数据才调出来）。
用的是hwd在WC2015上讲的方法，考虑将原图分割，根据每个圆的左右边界和圆与圆交点的横坐标来分割，这样原图就被分成很多竖着的长条，并且每一条中间都没有交点，这样就有一个性质：每一条都是"弓形-梯形-弓形 弓形-梯形-弓形..."的形式，然后从一个方向开始，记录当前进入的圆的数量，每当出来就就计算面积。
1 /************************************************************** 2 Problem: 2178 3 User: idy002 4 Language: C++ 5 Result: Accepted 6 Time:1248 ms 7 Memory:48560 kb 8 ****************************************************************/ 9 10 #include <cstdio> 11 #include <cstring> 12 #include <cmath> 13 #include <algorithm> 14 #define eps 1e-11 15 #define N 2010 16 using namespace std; 17 18 19 inline int sg( long double x ) { return (x>-eps)-(x<eps); } 20 struct Vector { 21 long double x, y; 22 Vector(){} 23 Vector( long double x, long double y ):x(x),y(y){} 24 Vector operator+( const Vector & b ) const { return Vector(x+b.x,y+b.y); } 25 Vector operator-( const Vector & b ) const { return Vector(x-b.x,y-b.y); } 26 Vector operator*( long double b ) const { return Vector(x*b,y*b); } 27 Vector operator/( long double b ) const { return Vector(x/b,y/b); } 28 long double operator^( const Vector & b ) const { return x*b.y-y*b.x; } 29 long double operator&( const Vector & b ) const { return x*b.x+y*b.y; } 30 long double ang() { return atan2(y,x); } 31 long double len() { return sqrt(x*x+y*y); } 32 long double len2() { return x*x+y*y; } 33 }; 34 typedef Vector Point; 35 struct Circle { 36 Point c; 37 long double r; 38 Circle(){} 39 Circle( Point c, long double r ):c(c),r(r){} 40 inline bool in( Point &p ) { 41 Vector vv=p-c; 42 return (vv.x*vv.x+vv.y*vv.y) < r*r*0.999; 43 } 44 Point pt( long double ang ) const { 45 return c+Vector(cos(ang),sin(ang))*r; 46 } 47 }; 48 struct Arc { 49 int cid, type; 50 long double al, ar; 51 Point pa, pb; 52 Arc(){} 53 Arc( int cid, int type, const Point &a, const Point &b ); 54 long double area(); 55 }; 56 57 int n; 58 bool del[N], has[N]; 59 long double spt[N*N]; int stot; 60 Circle cir[N]; 61 Arc arc[N+N]; 62 63 bool cmp_arc( const Arc &a, const Arc &b ) { 64 if( sg(a.pa.y-b.pa.y)!=0 ) return sg(a.pa.y-b.pa.y)<0; 65 if( sg(a.pb.y-b.pb.y)!=0 ) return sg(a.pb.y-b.pb.y)<0; 66 if( a.type!=b.type ) return a.type>b.type; 67 if( a.type==1 ) { 68 return cir[a.cid].r < cir[b.cid].r; 69 } else { 70 return cir[a.cid].r > cir[b.cid].r; 71 } 72 } 73 bool cmp_cir_r( const Circle &a, const Circle &b ) { return a.r<b.r; } 74 Arc::Arc( int cid, int type, const Point &a, const Point &b ) { 75 this->cid = cid; 76 this->type = type; 77 pa = a; 78 pb = b; 79 this->al = (a-cir[cid].c).ang(); 80 this->ar = (b-cir[cid].c).ang(); 81 if( pa.x>pb.x ) swap(pa,pb); 82 } 83 long double Arc::area() { 84 long double da = ar-al; 85 while( sg(da)<=0 ) da+=M_PI+M_PI; 86 while( sg(da-M_PI-M_PI)>0 ) da-=M_PI+M_PI; 87 return (da-sin(da))*cir[cid].r*cir[cid].r/2.0; 88 } 89 int ccinter( int ca, int cb, Point *p ) { 90 if( cir[ca].r<cir[cb].r ) swap(ca,cb); 91 long double cd = (cir[ca].c - cir[cb].c).len2(); 92 long double s1 = cir[ca].r-cir[cb].r; 93 long double s2 = cir[ca].r+cir[cb].r; 94 s1=s1*s1, s2=s2*s2; 95 96 if( sg(cd-s1)<0 || sg(cd-s2)>0 ) return 0; 97 if( sg(cd-s1)==0 || sg(cd-s2)==0 ) { 98 long double ang = (cir[cb].c-cir[ca].c).ang(); 99 p[0] = cir[ca].pt(ang); 100 return 1; 101 } 102 long double r1 = cir[ca].r, r2 = cir[cb].r; 103 long double base = (cir[cb].c-cir[ca].c).ang(); 104 long double delta = acos( (r1*r1+cd-r2*r2)/(2.0*r1*sqrt(cd)) ); 105 p[0] = cir[ca].pt( base-delta ); 106 p[1] = cir[ca].pt( base+delta ); 107 return 2; 108 } 109 bool clinter( int c, long double xl, long double xr, Arc *arc ) { 110 long double xlb = cir[c].c.x-cir[c].r; 111 long double xrb = cir[c].c.x+cir[c].r; 112 if( sg(xlb-xr)>=0 || sg(xrb-xl)<=0 ) return false; 113 Point p = cir[c].c; 114 long double r = cir[c].r; 115 long double d1, d2; 116 long double s1 = r*r-(p.x-xl)*(p.x-xl); 117 long double s2 = r*r-(p.x-xr)*(p.x-xr); 118 if( s1<0.0 ) s1=0.0; 119 if( s2<0.0 ) s2=0.0; 120 d1 = sqrt( s1 ); 121 d2 = sqrt( s2 ); 122 arc[0] = Arc( c, -1, Point(xr,p.y+d2), Point(xl,p.y+d1) ); 123 arc[1] = Arc( c, +1, Point(xl,p.y-d1), Point(xr,p.y-d2) ); 124 return true; 125 } 126 long double area( long double lf, long double rg ) { 127 int m=0; 128 for( int i=0; i<n; i++ ) 129 if( clinter(i,lf,rg,arc+m) ) 130 m += 2; 131 sort( arc, arc+m, cmp_arc ); 132 133 int cc = 0, top; 134 long double rt = 0.0; 135 for( int i=0; i<m; i++ ) { 136 if( cc==0 ) 137 top = i; 138 cc += arc[i].type; 139 if( cc==0 ) { 140 rt += (rg-lf)*((arc[i].pa+arc[i].pb).y-(arc[top].pa+arc[top].pb).y)/2.0; 141 rt += arc[top].area()+arc[i].area(); 142 } 143 } 144 return rt; 145 } 146 bool cmp_eq( long double x, long double y ) { 147 return sg(x-y)==0; 148 } 149 void clean() { 150 int j=0; 151 for( int i=0; i<n; i++ ) 152 if( !del[i] ) cir[j++]=cir[i]; 153 n = j; 154 } 155 long double area() { 156 Point ip[2]; 157 int ic; 158 long double rt = 0.0; 159 for( int i=0; i<n; i++ ) { 160 for( int j=i+1; j<n; j++ ) { 161 ic = ccinter( i, j, ip ); 162 if( ic<=1 ) continue; 163 for( int k=0; k<ic; k++ ) { 164 int q=n; 165 for( q=0; q<n; q++ ) { 166 if( q==i || q==j ) continue; 167 if( cir[q].in(ip[k]) ) 168 break; 169 } 170 if( q==n ) spt[stot++]=ip[k].x; 171 } 172 } 173 } 174 for( int i=0; i<n; i++ ) { 175 spt[stot++] = cir[i].c.x-cir[i].r; 176 spt[stot++] = cir[i].c.x+cir[i].r; 177 } 178 sort( spt, spt+stot ); 179 stot = unique( spt, spt+stot, cmp_eq ) - spt; 180 181 for( int i=1; i<stot; i++ ) 182 rt += area( spt[i-1], spt[i] ); 183 return rt; 184 } 185 void init() { 186 sort( cir, cir+n, cmp_cir_r ); 187 for( int i=0; i<n; i++ ) 188 for( int j=i+1; j<n; j++ ) { 189 long double dx = cir[i].c.x-cir[j].c.x; 190 long double dy = cir[i].c.y-cir[j].c.y; 191 long double dij = dx*dx+dy*dy; 192 long double cc = cir[j].r-cir[i].r; 193 cc = cc*cc; 194 if( dij<cc ) del[i]=true; 195 } 196 clean(); 197 } 198 int main() { 199 scanf( "%d", &n ); 200 for( int i=0; i<n; i++ ) 201 scanf( "%Lf%Lf%Lf", &cir[i].c.x, &cir[i].c.y, &cir[i].r ); 202 init(); 203 printf( "%.3Lf\n", area() ); 204 }
