传送门
Solution
考虑求每个点的贡献
等价于一个以OA长为半径的圆心为原点的圆在多边形内的弧对应的角度/(2pi)
求弧度可以利用三角剖分
在原点的点要特判,采用射线法就可以了
Code
#include <bits/stdc++.h>
#define reg register
#define ll long long
#define db double
using namespace std;
int read()
{
int x=0,f=1;char ch=getchar();
while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
while(ch<='9'&&ch>='0'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return x*f;
}
const db eps=1e-6,Pi=acos(-1.);
const int N=1005;
int sign(db x){return ((x>eps)-(x<-eps));}
struct P{
db x,y;
P(db x=0,db y=0):x(x),y(y){}
P operator+(P o){return P(x+o.x,y+o.y);}
P operator-(P o){return P(x-o.x,y-o.y);}
P operator *(db o){return P(x*o,y*o);}
P rev(){return P(-x,-y);}
db operator *(P o){return x*o.x+y*o.y;}
db operator ^(P o){return x*o.y-y*o.x;}
db mo(){return sqrt(x*x+y*y);}
db thi(P o){return acos(((*this)*o)/(mo()*o.mo()));}
}s[N],p[N];
db calc(P a,P b,db R)
{
db ras=0;
db f=sign(a^b);
if(sign(max(a.mo(),b.mo())-R)==-1) return 0;
db thi=a.thi(b);
db h=fabs((a^b)/(a-b).mo());
if(sign(h-R)>-1) return thi*f;
db a0=asin(h/R),a1=a.rev().thi(b-a),a2=b.rev().thi(a-b);
if(sign(a0-a1)==1) ras+=a0-a1;
if(sign(a0-a2)==1) ras+=a0-a2;
return min(ras,thi)*f;
}
int chk(P x,P y)
{
if(!sign(x.y)) return sign(x.x)==1;
if(!sign(y.y)) return 0;
if(x.y>y.y) swap(x,y);
if(sign(x.y*y.y)==1) return 0;
P xxx=x+(y-x)*(x.y/(y.y-x.y));
return sign(xxx.x)==1;
}
int main()
{
reg int i,j,n,m;
db ans=0.;
n=read(),m=read();
for(i=1;i<=n;++i) p[i].x=read(),p[i].y=read();
for(i=1;i<=m;++i) s[i].x=read(),s[i].y=read();
s[m+1]=s[1];
for(i=1;i<=m;++i)
{
if(!sign(s[i]^s[i+1])) continue;
for(j=1;j<=n;++j)
{
if(!sign(p[j].mo())) continue;
ans+=calc(s[i],s[i+1],p[j].mo());
}
}
for(i=1;i<=n;++i)
{
if(sign(p[i].mo())) continue;
int flag=1,cnt=0;
for(j=1;j<=m;++j)
{
if(sign(s[j]^s[j+1])) cnt+=chk(s[j],s[j+1]);
else if(sign(s[j]*s[j+1])<1){flag=0;break;}
}
if(flag&&cnt%2==1) ans+=Pi*2.;
}
ans/=(Pi*2.);
return 0*printf("%.5lf
",ans);
}
Blog来自PaperCloud,未经允许,请勿转载,TKS!