BZOJ2411 : 黄牌

将裁判作为原点，求出原点到每个圆的切点。

将这些切点以及矩形的顶点极角排序，用堆维护最靠近裁判的圆。

对于一段相邻的极角区间，如果没有圆，那么对答案的贡献是三角形的面积。

否则求出两条射线与圆的交点$A,B$，则对答案的贡献是三角形$OAB$的面积减去三角形$OAB$与圆的交的面积。

时间复杂度$O(nlog n)$。

注意因为精度问题，求射线与圆的交点时可能会被判断为不存在交点，此时应该直接用切点。

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<queue>
#include<vector>
using namespace std;
typedef long double ld;
typedef pair<int,int>PI;
const int N=1005;
const ld pi=acos(-1.0),eps=1e-9;
int n,m,h,w,w0,i,radius[N],ce,in[N],d[N];ld ans;
priority_queue<PI,vector<PI>,greater<PI> >q;
struct E{ld x;int p;E(){}E(ld _x,int _p,int _t){x=_x,p=_p<<3|_t;}}e[N<<1];
inline bool cmp(const E&a,const E&b){return a.x<b.x;}
inline void add(int x){in[x]=1;q.push(PI(d[x],x));}
inline int sgn(ld x){
  if(x<-eps)return -1;
  if(x>eps)return 1;
  return 0;
}
inline bool Quadratic(ld A,ld B,ld C,ld *t0,ld *t1){
  ld discrim=B*B-4.f*A*C;
  if(discrim<0.)return 0;
  ld rootDiscrim=sqrt(discrim);
  ld q;
  if(B<0)q=-.5f*(B-rootDiscrim);
  else   q=-.5f*(B+rootDiscrim);
  *t0=q/A;*t1=C/q;
  if(*t0>*t1)swap(*t0,*t1);
  return 1;
}
struct P{
  ld x,y;
  P(){x=y=0;}
  P(ld _x,ld _y){x=_x,y=_y;}
  P operator+(const P&b)const{return P(x+b.x,y+b.y);}
  P operator-(const P&b)const{return P(x-b.x,y-b.y);}
  P operator*(ld b)const{return P(x*b,y*b);}
  P operator/(ld b)const{return P(x/b,y/b);}
  ld operator*(const P&b)const{return x*b.x+y*b.y;}
  ld len()const{return hypot(x,y);}
  ld len_sqr()const{return x*x+y*y;}
  P rotate(ld c)const{return P(x*cos(c)-y*sin(c),x*sin(c)+y*cos(c));}
  P trunc(ld l)const{return(*this)*l/len();}
}O,left[N],right[N];
inline ld cross(const P&a,const P&b){return a.x*b.y-a.y*b.x;}
inline ld get_angle(const P&a,const P&b){return fabs(atan2(fabs(cross(a,b)),a*b));}
inline P lerp(const P&a,const P&b,ld t){return a*(1-t)+b*t;}
struct circle{
  P c;ld r;
  circle(){c=P(0,0),r=0;}
  circle(P _c,ld _r){c=_c,r=_r;}
}cir[N];
inline void getTangents(const P&p,const circle&C,P&A,P&B){
  P u=C.c-p;
  ld dist=u.len(),ang=asin(C.r/dist);
  u=u.trunc(sqrt(u.len_sqr()-C.r*C.r));
  A=u.rotate(-ang)+p;
  B=u.rotate(ang)+p;
}
inline bool circle_line_intersection(const circle&c,const P&a,const P&b,ld *t0,ld *t1){
  P d=b-a;ld A=d*d;
  ld B=d*(a-c.c)*2.0;
  ld C=(a-c.c).len_sqr()-c.r*c.r;
  return Quadratic(A,B,C,t0,t1);
}
inline ld circle_triangle_intersection_area(const circle&c,const P&a,const P&b){
  if(sgn(cross(a-c.c,b-c.c))==0)return 0;
  static P q[5];int len=0;ld t0,t1;q[len++]=a;
  if(circle_line_intersection(c,a,b,&t0,&t1)){
    if(0<=t0&&t0<=1)q[len++]=lerp(a,b,t0);
    if(0<=t1&&t1<=1)q[len++]=lerp(a,b,t1);
  }
  q[len++]=b;ld s=0;
  for(int i=1;i<len;i++){
    P z=(q[i-1]+q[i])/2;
    if((z-c.c).len_sqr()<=c.r*c.r)
      s+=fabs(cross(q[i-1]-c.c,q[i]-c.c))/2;
    else
      s+=c.r*c.r*get_angle(q[i-1]-c.c,q[i]-c.c)/2;
  }
  return s;
}
inline ld circle_polygon_intersection_area(const circle&c,const P&A,const P&B){
  static P v[9];
  int n=3;
  v[0]=O;
  v[1]=A;
  v[2]=B;
  ld s=0;
  for(int i=0;i<n;i++){
    int j=(i+1)%n;
    s+=circle_triangle_intersection_area(c,v[i],v[j])*sgn(cross(v[i]-c.c,v[j]-c.c));
  }
  return fabs(s);
}
inline P line_intersection(const P&a,const P&b,const P&p,const P&q){
  ld U=cross(p-a,q-p),D=cross(b-a,q-p);
  return a+(b-a)*(U/D);
}
inline ld triangle_area(const P&A,const P&B){
  return fabs(cross(O,A)+cross(A,B)+cross(B,O))/2;
}
inline void solve(ld L,ld R){
  if(L+eps>R)return;
  P A=P(cos(L),sin(L))+O;
  P B=P(cos(R),sin(R))+O;
  while(!q.empty()&&!in[q.top().second])q.pop();
  if(q.empty()){
    P C,D;
    if(m==0)C=P(w,0),D=P(w,h);
    else if(m==1)C=P(0,h),D=P(w,h);
    else C=P(0,0),D=P(0,h);
    ans+=triangle_area(line_intersection(O,A,C,D),line_intersection(O,B,C,D));
    return;
  }
  int x=q.top().second;
  ld t0,t1;P _A,_B;
  if(circle_line_intersection(cir[x],O,A,&t0,&t1))_A=lerp(O,A,t0);else _A=left[x];
  if(circle_line_intersection(cir[x],O,B,&t0,&t1))_B=lerp(O,B,t0);else _B=right[x];
  ans+=triangle_area(_A,_B)-circle_polygon_intersection_area(cir[x],_A,_B);
}
int main(){
  scanf("%d%d%d%d",&n,&w,&h,&w0);
  O=P(w0,0);
  e[++ce]=E(atan2(h,w-w0),0,1);
  e[++ce]=E(atan2(h,-w0),0,2);
  for(i=1;i<=n;i++){
    int x,y;
    scanf("%d%d%d",&x,&y,&radius[i]);
    cir[i]=circle(P(x,y),radius[i]);
    int w=(x-w0)*(x-w0)+y*y;
    d[i]=w-radius[i]*radius[i];
    ld t=atan2(y,x-w0),o=asin(radius[i]/sqrt(w)),l=t-o,r=t+o;
    e[++ce]=E(l,i,4);
    e[++ce]=E(r,i,5);
    P A,B;
    getTangents(O,cir[i],A,B);
    left[i]=A;
    right[i]=B;
  }
  sort(e+1,e+ce+1,cmp);
  for(i=1;i<=ce;i++){
    solve(e[i-1].x,e[i].x);
    if((e[i].p&7)==4)add(e[i].p>>3);
    else if((e[i].p&7)==5)in[e[i].p>>3]=0;
    else m=e[i].p&7;
  }
  solve(e[ce].x,pi);
  ans/=h*w;
  return printf("%.8f",(double)ans),0;
}