【洛谷5467】[PKUSC2018] PKUSC（计算几何）

点此看题面

• 给定\(n\)个点和一个\(m\)个点构成的多边形。
• 随机将多边形旋转一个角度，求落在多边形内部的点数的期望。
• \(n\le200,m\le500\)

---

能自己把这道题的做法口胡出来，但具体实现中一些地方还是参考了一下题解。

不过计算几何应该都是口胡容易代码难吧。。。

解题思路

考虑求多边形内点数的期望，由于点与点之间是独立的，所以这就相当于是求每个点落在多边形内的概率之和。

而一个点落在多边形内的概率，就是\([0,2\pi)\)内所有能使它落在多边形内的弧度总数除以\(2\pi\)。

转多边形实在是太麻烦了，因此我们改为考虑点的旋转，而这实际上就是一个圆。

那么我们的问题就变成了以原点为圆心、以该点到原点的距离为半径画一个圆，圆上在多边形内部的弧对应的角度总和。

这只要先求出圆与多边形上每一条边的交点，然后分别考虑每相邻两个交点之间的弧是否在多边形内部即可，而这也等价于这条弧的中点是否在多边形内部。

口胡到这里其实已经结束了，接下来具体介绍几种函数的实现。

求线段与圆的交点

注意这道题中的圆必然以原点为圆心。

我们先求出原点到线段的垂足，如果原点到该垂足的距离\(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\n",ans/(2*Pi)),0;//求所有点概率和
}