【loj6437】 【PKUSC2018】 PKUSC 计算几何
题目大意：给你一个m个点的简单多边形。对于每个点i∈[1,n]，作一个以O点为原点且过点i的圆，求该圆在多边形内的圆弧长度/圆长。
其中n≤200，m≤500。
我们将n个点分开处理。
首先，我们要判断需处理的圆，是否被包含在多边形内，或者圆把多边形包含了。
我们显然可以从原点出发，向上作一条x=0的射线，判断该射线与多边形有多少个交点。
显然，若交点数量为奇数个，那么该点就在多边形内，否则在多边形外。
若圆与多边形存在交点，我们对多边形的每条边，求出其与圆，有多少个交点（0个，1个，2个，其实1个点可以当2个去处理）
这里就是简单地推个式子就可以求出，详见代码。
求出这些交点后，进行极角排序。
对于两个排序后相邻的点，我们要判断，由这两个点所构成的弧，是否被多边形包含。
有一种很暴力的思路，就是求出这两个点在弧上的中点，然后作一条从该点出发，向上作一条平行于x轴的射线，求这条射线与多边形的相交次数。
显然，若交点数量为奇数个，那么该点就在多边形内，否则在多边形外（和上面判断圆是否被包含的部分相同）。
对于每个圆，最多会出现2m个交点，最后单次判断两点所构乘的圆弧是否被多边形包含，需要O(m)的时间，则总时间复杂度为O(n*m^2)。
下面是一些需要注意的细节：
1，上述判断奇数偶数个点时，所做射线可能会经过多边形两条边的交点。考虑到输入的点全部都是整点，故对构成多边形的点集偏移2eps即可。
2，由于点要被多边形严格包含，那么在判断时，范围要向两边缩减eps。
思路简单，代码复杂，qwq....
（听说还卡精度）。。。
1 #include<bits/stdc++.h> 2 #define M 505 3 #define DB long double 4 #define eps 1e-6 5 #define PI 3.14159265358979323846 6 using namespace std; 7 8 bool zero(DB x){return fabs(x)<eps;} 9 struct pt{ 10 DB x,y; 11 pt(){x=y=0;} 12 pt(DB xx,DB yy){x=xx; y=yy;} 13 friend pt operator -(pt a,pt b){return pt(a.x-b.x,a.y-b.y);} 14 friend bool operator <(pt a,pt b){return atan2(a.y,a.x)<atan2(b.y,b.x);} 15 DB mo(){return x*x+y*y;} 16 }; 17 18 struct line{ 19 DB a,b,c; 20 line(){a=b=c=0;} 21 line(DB aa,DB bb,DB cc){a=aa; b=bb; c=cc;} 22 line(pt A,pt B){ 23 DB x1=A.x,y1=A.y; 24 DB x2=B.x,y2=B.y; 25 a=y1-y2; b=x2-x1; 26 c=-a*x2-b*y2; 27 } 28 }; 29 30 pt a[M],b[M],c[M]; int n,m; 31 int type[M]={0}; 32 DB sita[M]={0}; 33 DB solve(pt hh){ 34 memset(c,0,sizeof(c)); 35 memset(type,0,sizeof(type)); 36 int cnt=0,ok=1; 37 DB R=hh.x*hh.x+hh.y*hh.y,r=sqrt(R); 38 for(int i=1;i<=m;i++){ 39 line p=line(b[i],b[i+1]); 40 if(fabs(p.c/sqrt(p.a*p.a+p.b*p.b))-eps<r) ok=0; 41 DB xl=min(b[i].x,b[i+1].x),xr=max(b[i].x,b[i+1].x); 42 DB yl=min(b[i].y,b[i+1].y),yr=max(b[i].y,b[i+1].y); 43 xl-=eps; yl-=eps; xr+=eps; yr+=eps; 44 DB x1,x2,y1,y2; bool is0=0; 45 if(zero(p.a)){ 46 y1=y2=-p.c/p.b; 47 DB delta=R-y1*y1; 48 if(delta<-eps) continue; 49 if(zero(delta)) is0=1; 50 x1=-sqrt(delta); x2=sqrt(delta); 51 }else{ 52 DB A=p.b*p.b+p.a*p.a; 53 DB B=p.b*p.c*2; 54 DB C=p.c*p.c-R*p.a*p.a; 55 DB delta=B*B-4*A*C; 56 if(delta<-eps) continue; 57 if(zero(delta)) is0=1; 58 y1=(-B+sqrt(delta))/(2*A); 59 y2=(-B-sqrt(delta))/(2*A); 60 x1=(-p.b*y1-p.c)/p.a; 61 x2=(-p.b*y2-p.c)/p.a; 62 } 63 int ok1,ok2; 64 if(xl<=x1&&x1<=xr&&yl<=y1&&y1<=yr){ c[++cnt]=pt(x1,y1); if(is0) continue;} 65 if(xl<=x2&&x2<=xr&&yl<=y2&&y2<=yr) c[++cnt]=pt(x2,y2); 66 } 67 if(ok){ 68 //return 1; 69 c[cnt=1]=pt(0,0); 70 }else{ 71 sort(c+1,c+cnt+1); 72 for(int i=1;i<=cnt;i++) sita[i]=atan2(c[i].y,c[i].x); 73 } 74 c[cnt+1]=c[1]; sita[cnt+1]=sita[1]+2*PI; 75 DB ans=0; 76 for(int i=1;i<=cnt;i++){ 77 double SITA=(sita[i]+sita[i+1])/2; 78 pt hh=pt(r*cos(SITA)+eps*2,r*sin(SITA)); 79 int cnt=0; 80 for(int j=1;j<=m;j++){ 81 double l=b[j].x,r=b[j+1].x; 82 if(zero(b[j].x-b[j+1].x)) continue; 83 line p=line(b[j],b[j+1]); 84 double y=-(p.a*hh.x+p.c)/p.b; 85 if(l>r) swap(l,r); 86 l+=eps; r-=eps; 87 if(l<=hh.x&&hh.x<=r&&y+eps>hh.y) cnt++; 88 } 89 if(cnt&1) ans+=sita[i+1]-sita[i]; 90 } 91 ans/=(2.*PI); 92 return ans; 93 } 94 95 96 int main(){ 97 // freopen("in.txt","r",stdin); 98 // freopen("out.txt","w",stdout); 99 scanf("%d%d",&n,&m); 100 for(int i=1;i<=n;i++) cin>>a[i].x>>a[i].y; 101 for(int i=1;i<=m;i++) cin>>b[i].x>>b[i].y; 102 b[m+1]=b[1]; 103 DB ans=0; 104 for(int i=1;i<=n;i++){ 105 ans+=solve(a[i]); 106 //printf("%.5Lf\n",ans); 107 } 108 double res=ans; 109 printf("%.5lf\n",res); 110 }
