Don't Burst the Balloon (AOJ1342)
円とかで削って残る点があるか→交点すべてをゆるい条件で調べる
containCCは細かいことを考えると結構面倒
#include <bits/stdc++.h> #define rep(i,n) for(int i=0;i<(int)(n);i++) #define rep1(i,n) for(int i=1;i<=(int)(n);i++) #define all(c) c.begin(),c.end() #define pb push_back #define fs first #define sc second #define show(x) cout << #x << " = " << x << endl #define chmin(x,y) x=min(x,y) #define chmax(x,y) x=max(x,y) #define X real() #define Y imag() using namespace std; typedef double D; typedef complex<D> P; typedef pair<P,P> L; typedef vector<P> Pol; struct C{P p;D r;}; D eps=1e-9; bool eq(D a,D b){return abs(a-b)<eps;} int large(D a,D b){return eq(a,b)?0:(a>b?1:-1);} bool containCC(C a,C b){ return abs(a.p-b.p)+b.r<a.r+eps; } D dot(P a,P b){return real(conj(a)*b);} P perp(L l,P p){ D t=dot(p-l.fs,l.fs-l.sc)/norm(l.fs-l.sc); return l.fs+t*(l.fs-l.sc); } D dLP(L l,P p){return abs(perp(l,p)-p);} Pol intCC(C a,C b){ int x=large(abs(a.p-b.p),a.r+b.r); Pol ret; if(x==1) return ret; if(x==0){ ret.pb((a.r*b.p+b.r*a.p)/(a.r+b.r)); return ret; } if(containCC(a,b)||containCC(b,a)) return ret; D d=abs(a.p-b.p); D theta=acos((a.r*a.r+d*d-b.r*b.r)/(2.0*a.r*d)); ret.pb(a.p+(b.p-a.p)/d*polar(a.r,theta)); ret.pb(a.p+(b.p-a.p)/d*polar(a.r,-theta)); return ret; } Pol intCL(C c,L l){ int x=large(dLP(l,c.p),c.r); Pol ret; if(x==1) return ret; if(x==0){ ret.pb(perp(l,c.p)); return ret; } D d=dLP(l,c.p); if(d<eps){ ret.pb(c.p+(l.fs-l.sc)/abs(l.fs-l.sc)*c.r); ret.pb(c.p-(l.fs-l.sc)/abs(l.fs-l.sc)*c.r); return ret; } D theta=acos(d/c.r); P p=perp(l,c.p); ret.pb(c.p+(p-c.p)/d*polar(c.r,theta)); ret.pb(c.p+(p-c.p)/d*polar(c.r,-theta)); return ret; } P p[10]; int N; D h[10],w; D dis(D r,D h){ if(r<=h) return r; return sqrt(2*r*h-h*h); } void append(Pol& pol,Pol a){ for(P p:a) pol.pb(p); } bool can(D r){ D e=dis(r,w); if(e>=50) return 0; vector<L> ls; vector<P> ps; ps.pb(P(e,e)); ps.pb(P(e,100-e)); ps.pb(P(100-e,100-e)); ps.pb(P(100-e,e)); rep(i,4) ls.pb(L(ps[i],ps[(i+1)%4])); vector<C> cs; rep(i,N) cs.pb(C{p[i],dis(r,h[i])}); rep(i,N) rep(j,i) append(ps,intCC(cs[i],cs[j])); rep(i,N) rep(j,4) append(ps,intCL(cs[i],ls[j])); for(P p:ps){ if(e-eps<p.X&&p.X<100-e+eps&&e-eps<p.Y&&p.Y<100-e+eps){ bool ok=1; rep(i,N){ if(abs(p-cs[i].p)+eps<cs[i].r) ok=0; } if(ok) return 1; } } return 0; } int main(){ while(true){ cin>>N>>w; if(N==0) break; rep(i,N){ int x,y; cin>>x>>y>>h[i]; p[i]=P(x,y); } D ub=1e5,lb=0; rep(tt,50){ D r=(ub+lb)/2; if(can(r)) lb=r; else ub=r; } printf("%.12f\n",ub); } }