标签：黄牌 return BZOJ2411 ld t0 t1 const circle

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

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

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

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

时间复杂度$O(n\log 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;
}

标签：黄牌,return,BZOJ2411,ld,t0,t1,const,circle
来源： https://www.cnblogs.com/clrs97/p/12163880.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。