- 给定(n)个点和一个(m)个点构成的多边形。
- 随机将多边形旋转一个角度,求落在多边形内部的点数的期望。
- (nle200,mle500)
能自己把这道题的做法口胡出来,但具体实现中一些地方还是参考了一下题解。
不过计算几何应该都是口胡容易代码难吧。。。
解题思路
考虑求多边形内点数的期望,由于点与点之间是独立的,所以这就相当于是求每个点落在多边形内的概率之和。
而一个点落在多边形内的概率,就是([0,2pi))内所有能使它落在多边形内的弧度总数除以(2pi)。
转多边形实在是太麻烦了,因此我们改为考虑点的旋转,而这实际上就是一个圆。
那么我们的问题就变成了以原点为圆心、以该点到原点的距离为半径画一个圆,圆上在多边形内部的弧对应的角度总和。
这只要先求出圆与多边形上每一条边的交点,然后分别考虑每相邻两个交点之间的弧是否在多边形内部即可,而这也等价于这条弧的中点是否在多边形内部。
口胡到这里其实已经结束了,接下来具体介绍几种函数的实现。
求线段与圆的交点
注意这道题中的圆必然以原点为圆心。
我们先求出原点到线段的垂足,如果原点到该垂足的距离(l)已经超过半径(r)说明二者相离,无交点。
否则,我们分别按该线段的两种方向,给这个垂足加上一个长度为(sqrt{r^2-l^2})的向量,然后看看得到的这个新点是不是在线段上就好了。
判断点是否在多边形内
这倒是一个我之前就写过的东西,但现在发现有更简单的写法。
考虑从一个点随便向一个方向引出一条线,那么与多边形的边相交次数的奇偶性就意味着它是否在多边形中。
但是,一个交点同时存在于两条边中,所以当这条线经过交点的时候就可能会出错。
以前我傻乎乎地直接向右引一条线,每次要讨论是否与端点相交,有些麻烦。
现在发现,其实只要随便向一个刁钻的角度引一条线,就不用考虑端点的问题了。
代码:(O(nm^2))
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 200
#define M 500
#define DB double
#define eps 1e-8
using namespace std;
int n,m,cnt,x[N+5],y[N+5];DB Pi=acos(-1);
I int dcmp(Con DB& x) {return fabs(x)<eps?0:(x<0?-1:1);}
struct P
{
#define CP Con P&
DB x,y,ang;I P() {x=y=ang=0;}I P(Con DB& a,Con DB& b):x(a),y(b){ang=0;}
I P operator + (CP o) Con {return P(x+o.x,y+o.y);}
I P operator - (CP o) Con {return P(x-o.x,y-o.y);}
I P operator * (Con DB& o) Con {return P(x*o,y*o);}
I P operator / (Con DB& o) Con {return P(x/o,y/o);}
I bool operator == (CP o) Con {return !dcmp(ang-o.ang);}
I bool operator < (CP o) Con {return dcmp(ang-o.ang)<0;}
I DB operator * (CP o) Con {return x*o.x+y*o.y;}
I DB operator ^ (CP o) Con {return x*o.y-y*o.x;}
I DB Len() {return sqrt(x*x+y*y);}
}p[M+5],f[2*M+5];
struct S
{
#define CS Con S&
P A,B;I S() {A=B=P();}I S(CP a,CP b):A(a),B(b){}
};
I P Foot(CS s) {return P(-(s.B-s.A).y,(s.B-s.A).x);}//垂线上的一点
I bool P_on_S(CP p,CS s)//判断这条直线上的一点是否在线段上
{
if(dcmp(p.x-min(s.A.x,s.B.x))<0||dcmp(p.x-max(s.A.x,s.B.x))>0) return 0;
if(dcmp(p.y-min(s.A.y,s.B.y))<0||dcmp(p.y-max(s.A.y,s.B.y))>0) return 0;return 1;
}
I P L_with_L(CS X,CS Y)//求两条直线交点
{
DB s1=(X.B-X.A)^(Y.A-X.A),s2=(X.B-X.A)^(Y.B-X.A);
return Y.A+(Y.B-Y.A)*s1/(s1-s2);//利用面积比得出长度比
}
I void S_with_C(CS s,Con DB& r)//求线段与圆的交点
{
P t=L_with_L(S(P(),Foot(s)),s);if(t.Len()>r) return;//求出垂足,如果距离超过r说明相离
P g=(s.B-s.A)/(s.B-s.A).Len()*sqrt(r*r-t.Len()*t.Len());//求出用于变化的向量
P_on_S(t+g,s)&&(f[++cnt]=t+g,0),P_on_S(t-g,s)&&(f[++cnt]=t-g,0);//两种方向
}
I bool P_in_G(CP x)//判断点是否在多边形内
{
RI i,res=0;P t;S s(x,P(x.x+0.413,x.y+0.521));//向一个刁钻的角度引一条线
for(i=1;i<=m;++i) t=L_with_L(s,S(p[i],p[i+1])),//求出与边所在直线的交点
res^=P_on_S(t,S(p[i],p[i+1]))&&dcmp(t.x-x.x)>=0;return res;//判断是否与边有交
}
I DB Calc(CI x,CI y)//求解点(x,y)的答案
{
RI i;DB r=P(x,y).Len();for(cnt=0,i=1;i<=m;++i) S_with_C(S(p[i],p[i+1]),r);//求出圆与多边形每条边的交点
if(!cnt) return P_in_G(P(r,0))?2*Pi:0;for(i=1;i<=cnt;++i) f[i].ang=atan2(f[i].x,f[i].y);//特判没有交点,判断是否完全在多边形内部
sort(f+1,f+cnt+1),cnt=unique(f+1,f+cnt+1)-f-1,(f[cnt+1]=f[1]).ang+=2*Pi;
DB res=0;P t;for(i=1;i<=cnt;++i) t=Foot(S(f[i],f[i+1])),//检验相邻两交点之间的弧
P_in_G(t/t.Len()*r)&&(res+=f[i+1].ang-f[i].ang);return res;//取中点判断是否在多边形内部
}
int main()
{
RI i;for(scanf("%d%d",&n,&m),i=1;i<=n;++i) scanf("%d%d",x+i,y+i);
for(i=1;i<=m;++i) scanf("%lf%lf",&p[i].x,&p[i].y);p[m+1]=p[1];
DB ans=0;for(i=1;i<=n;++i) ans+=Calc(x[i],y[i]);return printf("%.5lf
",ans/(2*Pi)),0;//求所有点概率和
}