uva 11978 Fukushima Nuclear Blast (二分+多边形与圆交)
二分+多边形与圆交就可以了
1 //#include<bits/stdc++.h> 2 #include<iostream> 3 #include<cstdio> 4 #include<cmath> 5 #include<cstring> 6 #include<vector> 7 #include<algorithm> 8 using namespace std; 9 typedef long long LL; 10 11 const double PI = acos(-1.0); 12 const double EPS = 1e-8; 13 14 inline int sgn(double x) { 15 return (x > EPS) - (x < -EPS); 16 } 17 18 struct Point { 19 double x, y; 20 Point() {} 21 Point(double x, double y): x(x), y(y) {} 22 void read() { 23 scanf("%lf%lf", &x, &y); 24 } 25 double angle() { 26 return atan2(y, x); 27 } 28 Point operator + (const Point &rhs) const { 29 return Point(x + rhs.x, y + rhs.y); 30 } 31 Point operator - (const Point &rhs) const { 32 return Point(x - rhs.x, y - rhs.y); 33 } 34 Point operator * (double t) const { 35 return Point(x * t, y * t); 36 } 37 Point operator / (double t) const { 38 return Point(x / t, y / t); 39 } 40 double operator *(const Point &b)const 41 { 42 return x*b.x + y*b.y; 43 } 44 double length() const { 45 return sqrt(x * x + y * y); 46 } 47 Point unit() const { //单位向量 48 double l = length(); 49 return Point(x / l, y / l); 50 } 51 }; 52 double cross(const Point &a, const Point &b) { 53 return a.x * b.y - a.y * b.x; 54 } 55 double sqr(double x) { 56 return x * x; 57 } 58 double dist(const Point &p1, const Point &p2) { 59 return (p1 - p2).length(); 60 } 61 double sdist(Point a,Point b){ 62 return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y); 63 } 64 //向量 op 逆时针旋转 angle 65 Point rotate(const Point &p, double angle, const Point &o = Point(0, 0)) { 66 Point t = p - o; 67 double x = t.x * cos(angle) - t.y * sin(angle); 68 double y = t.y * cos(angle) + t.x * sin(angle); 69 return Point(x, y) + o; 70 } 71 Point line_inter(Point A,Point B,Point C,Point D){ //直线相交交点 72 Point ans; 73 double a1=A.y-B.y; 74 double b1=B.x-A.x; 75 double c1=A.x*B.y-B.x*A.y; 76 77 double a2=C.y-D.y; 78 double b2=D.x-C.x; 79 double c2=C.x*D.y-D.x*C.y; 80 81 ans.x=(b1*c2-b2*c1)/(a1*b2-a2*b1); 82 ans.y=(a2*c1-a1*c2)/(a1*b2-a2*b1); 83 return ans; 84 } 85 Point p_to_seg(Point p,Point a,Point b){ //点到线段的最近点 86 Point tmp=p; 87 tmp.x+=a.y-b.y; 88 tmp.y+=b.x-a.x; 89 if(cross(a-p,tmp-p)*cross(b-p,tmp-p)>0) return dist(p,a)<dist(p,b)?a:b; 90 return line_inter(p,tmp,a,b); 91 } 92 void line_circle(Point c,double r,Point a,Point b,Point &p1,Point &p2){ 93 Point tmp=c; 94 double t; 95 tmp.x+=(a.y-b.y);//求垂直于ab的直线 96 tmp.y+=(b.x-a.x); 97 tmp=line_inter(tmp,c,a,b); 98 t=sqrt(sqr(r)-sqr( dist(c,tmp)))/dist(a,b); //比例 99 p1.x=tmp.x+(b.x-a.x)*t; 100 p1.y=tmp.y+(b.y-a.y)*t; 101 p2.x=tmp.x-(b.x-a.x)*t; 102 p2.y=tmp.y-(b.y-a.y)*t; 103 } 104 struct Region { 105 double st, ed; 106 Region() {} 107 Region(double st, double ed): st(st), ed(ed) {} 108 bool operator < (const Region &rhs) const { 109 if(sgn(st - rhs.st)) return st < rhs.st; 110 return ed < rhs.ed; 111 } 112 }; 113 struct Circle { 114 Point c; 115 double r; 116 vector<Region> reg; 117 Circle() {} 118 Circle(Point c, double r): c(c), r(r) {} 119 void read() { 120 c.read(); 121 scanf("%lf", &r); 122 } 123 void add(const Region &r) { 124 reg.push_back(r); 125 } 126 bool contain(const Circle &cir) const { 127 return sgn(dist(cir.c, c) + cir.r - r) <= 0; 128 } 129 bool intersect(const Circle &cir) const { 130 return sgn(dist(cir.c, c) - cir.r - r) < 0; 131 } 132 }; 133 void intersection(const Circle &cir1, const Circle &cir2, Point &p1, Point &p2) { //两圆相交 交点 134 double l = dist(cir1.c, cir2.c); //两圆心的距离 135 double d = (sqr(l) - sqr(cir2.r) + sqr(cir1.r)) / (2 * l); //cir1圆心到交点直线的距离 136 double d2 = sqrt(sqr(cir1.r) - sqr(d)); //交点到 两圆心所在直线的距离 137 Point mid = cir1.c + (cir2.c - cir1.c).unit() * d; 138 Point v = rotate(cir2.c - cir1.c, PI / 2).unit() * d2; 139 p1 = mid + v, p2 = mid - v; 140 } 141 Point calc(const Circle &cir, double angle) { 142 Point p = Point(cir.c.x + cir.r, cir.c.y); 143 return rotate(p, angle, cir.c); 144 } 145 const int MAXN = 1010; 146 Circle cir[MAXN],cir2[MAXN]; 147 bool del[MAXN]; 148 int n; 149 double get_area(Circle* cir,int n) { //多个圆的相交面积 150 double ans = 0; 151 memset(del,0,sizeof(del)); 152 for(int i = 0; i < n; ++i) { 153 for(int j = 0; j < n; ++j) if(!del[j]) { //删除被包含的圆 154 if(i == j) continue; 155 if(cir[j].contain(cir[i])) { 156 del[i] = true; 157 break; 158 } 159 } 160 } 161 for(int i = 0; i < n; ++i) if(!del[i]) { 162 Circle &mc = cir[i]; 163 Point p1, p2; 164 bool flag = false; 165 for(int j = 0; j < n; ++j) if(!del[j]) { 166 if(i == j) continue; 167 if(!mc.intersect(cir[j])) continue; 168 flag = true; 169 intersection(mc, cir[j], p1, p2); //求出两圆的交点 170 double rs = (p2 - mc.c).angle(), rt = (p1 - mc.c).angle(); 171 if(sgn(rs) < 0) rs += 2 * PI; 172 if(sgn(rt) < 0) rt += 2 * PI; 173 if(sgn(rs - rt) > 0) mc.add(Region(rs, PI * 2)), mc.add(Region(0, rt)); //添加相交区域 174 else mc.add(Region(rs, rt)); 175 } 176 if(!flag) { 177 ans += PI * sqr(mc.r); 178 continue; 179 } 180 sort(mc.reg.begin(), mc.reg.end()); //对相交区域进行排序 181 int cnt = 1; 182 for(int j = 1; j < int(mc.reg.size()); ++j) { 183 if(sgn(mc.reg[cnt - 1].ed - mc.reg[j].st) >= 0) { //如果有区域可以合并，则合并 184 mc.reg[cnt - 1].ed = max(mc.reg[cnt - 1].ed, mc.reg[j].ed); 185 } else mc.reg[cnt++] = mc.reg[j]; 186 } 187 mc.add(Region()); 188 mc.reg[cnt] = mc.reg[0]; 189 for(int j = 0; j < cnt; ++j) { 190 p1 = calc(mc, mc.reg[j].ed); 191 p2 = calc(mc, mc.reg[j + 1].st); 192 ans += cross(p1, p2) / 2; // 193 double angle = mc.reg[j + 1].st - mc.reg[j].ed; 194 if(sgn(angle) < 0) angle += 2 * PI; 195 ans += 0.5 * sqr(mc.r) * (angle - sin(angle)); //弧所对应的的面积 196 } 197 } 198 return ans; 199 } 200 double two_cir(Circle t1,Circle t2){ //两个圆的相交面积 201 if(t1.contain(t2)||t2.contain(t1)) return PI * sqr(min(t2.r,t1.r)); 202 if(!t1.intersect(t2)) return 0; 203 double ans=0,len=dist(t1.c,t2.c); 204 double x=(sqr(t1.r)+sqr(len)-sqr(t2.r))/(2*len); 205 double angle1=acos(x/t1.r),angle2=acos((len-x)/t2.r); 206 ans=sqr(t1.r)*angle1+sqr(t2.r)*angle2-len*t1.r*sin(angle1); // 两个扇形 减去一个四边形面积 207 return ans; 208 } 209 double triangle_circle(Point a,Point b,Point o,double r){ 210 double sign=1.0; 211 double ans=0; 212 Point p1,p2; 213 a=a-o;b=b-o; 214 o=Point(0,0); 215 if(sgn(cross(a,b))==0) return 0.0; 216 if(sdist(a,o)>sdist(b,o)){ 217 swap(a,b); 218 sign=-1.0; 219 } 220 if(sdist(a,o)<r*r+EPS){ //a 在内， b 在外 221 if(sdist(b,o)<r*r+EPS) return cross(a,b)/2.0*sign; 222 line_circle(o,r,a,b,p1,p2); 223 if(dist(p1,b)>dist(p2,b)) swap(p1,p2); 224 double ans1=fabs(cross(a,p1)); 225 double ans2=acos(p1*b/p1.length()/b.length())*r*r; 226 ans=(ans1+ans2)/2.0; 227 if(cross(a,b)<EPS&&sign>0.0||cross(a,b)>EPS&&sign<0.0) return -ans; 228 return ans; 229 } 230 Point tmp=p_to_seg(o,a,b); 231 if(sdist(o,tmp)>r*r-EPS){ //a,b所在直线与圆没有交点 232 double angle=a*b/a.length()/b.length(); 233 if(angle>1.0) angle=1; 234 if(angle<-1.0) angle=-1.0; 235 ans=acos(a*b/a.length()/b.length())*r*r/2.0; 236 if(cross(a,b)<EPS&&sign>0.0||cross(a,b)>EPS&&sign<0.0) return -ans; 237 return ans; 238 } 239 line_circle(o,r,a,b,p1,p2); 240 double cm=r/(dist(a,o)-r); 241 Point m=Point((o.x+cm*a.x)/(1+cm),(o.y+cm*a.y)/(1+cm) ); 242 double cn=r/(dist(o,b)-r); 243 Point n=Point( (o.x+cn*b.x)/(1+cn) , (o.y+cn*b.y)/(1+cn) ); 244 double angle= m*n/m.length()/n.length(); 245 if(angle>1.0) angle=1; 246 if(angle<-1.0) angle=-1.0; 247 double ans1 = acos(m*n/m.length()/n.length())*r*r; 248 angle=p1*p2/p1.length()/p2.length() ; 249 if(angle>1.0) angle=1; 250 if(angle<-1.0) angle=-1.0; 251 double ans2 = acos(p1*p2/p1.length()/p2.length() )*r*r-fabs(cross(p1,p2)); 252 ans=(ans1-ans2)/2.0; 253 if(cross(a,b)<EPS&&sign>0.0||cross(a,b)>EPS&&sign<0.0) return -ans; 254 return ans; 255 } 256 Point pt[5100]; 257 int main() { 258 #ifndef ONLINE_JUDGE 259 freopen("input.txt","r",stdin); 260 #endif // ONLINE_JUDGE 261 Point a,b,c; 262 double x1,y1,h,x2,y2,r,p; 263 int t,cas=0;cin>>t; 264 while(t--){ 265 int n;cin>>n; 266 for(int i=0;i<n;i++) pt[i].read(); 267 c.read(); 268 cin>>p; 269 pt[n]=pt[0]; 270 double aLL=0,ans; 271 for(int i=0;i<n;i++) aLL+=cross(pt[i],pt[i+1])/2.0; 272 double l=0,r=1000000,mid; 273 while(r-l>1e-3){ 274 mid=(l+r)/2; 275 double sum=0; 276 for(int i=0;i<n;i++) sum+=triangle_circle(pt[i],pt[i+1],c,mid); 277 sum=fabs(sum); 278 if(sum>=aLL*p/100) { 279 ans=mid; 280 r=mid; 281 } 282 else l=mid; 283 } 284 printf("Case %d: %.0f\n",++cas,ans); 285 } 286 }
