[半平面交]小凸想跑步 LibreOJ2008

题目描述

小凸晚上喜欢到操场跑步，今天他跑完两圈之后，他玩起了这样一个游戏。

操场是个凸 n 边形，n 个顶点按照逆时针从 0∼n−1 编号。现在小凸随机站在操场中的某个位置，标记为 P 点。将 P 点与 n 个顶点各连一条边，形成 n 个三角形。如果这时 P 点，0 号点，1 号点形成的三角形的面积是 n 个三角形中最小的一个，小凸则认为这是一次正确站位。

现在小凸想知道他一次站位正确的概率是多少。

输入格式

第一行包含 1 个整数 n，表示操场的顶点数和游戏的次数。
接下来有 n 行，每行包含两个整数 Xi​、Yi​ 表示顶点的坐标。
输入保证按逆时针顺序输入点，所有点保证构成一个 n 多边形。所有点保证不存在三点共线。

输出格式

输出一个数，正确站位的概率，保留 4 位小数。

样例

Inputcopy	Outputcopy
5 1 8 0 7 0 0 8 0 8 8	0.6316

数据范围与提示

3≤N≤105,−109≤X,Y≤109

题意： 给出一个凸多边形，然后向该多边形内随机撒点p，将点p与多边形n个顶点连线构成n个三角形，问第0个、第1个顶点和点p构成的三角形是n个三角形中最小三角形的概率。

分析： 类似几何概型，最终概率就是一个面积比值。只要求出能使点p(x, y)满足条件的那个区域面积就行了。现在来分析下如何确定这块区域，要使题目中说的那个三角形是最小的那就需要它的面积小于或等于其余n-1个三角形面积，这样就得到了n-1个不等式，另外还需要保证取点是在n边形内部取，又需要n个不等式，这样加起来是2*n-1个不等式，化简一下得到2*n-1个半平面，求一下半平面交得到面积。

具体代码如下：

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <vector>
#define double long double
using namespace std;
const double eps = 1e-8;
const double inf = 1e20;
const double pi = acos(-1.0);
const int maxp = 200010;
//`Compares a double to zero`
int sgn(double x)
{
	if(fabs(x) < eps)return 0;
	if(x < 0)return -1;
	else return 1;
}
//square of a double
inline double sqr(double x){return x*x;}
struct Point
{
	double x, y;
	Point(){}
	Point(double _x,double _y){x = _x, y = _y;}
	void input(){scanf("%Lf%Lf",&x,&y);}
	void output(){printf("%.2f %.2f\n",x,y);}
	bool operator == (Point b)const{return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;}
	//第一关键字为x，第二关键字为y 
	bool operator < (Point b)const{return sgn(x-b.x)== 0?sgn(y-b.y)<0:x<b.x;}
	Point operator -(const Point &b)const{return Point(x-b.x,y-b.y);}
	//叉积
	double operator ^(const Point &b)const{return x*b.y - y*b.x;}
	//点积
	double operator *(const Point &b)const{return x*b.x + y*b.y;}
	//返回长度
	double len(){return hypot(x,y);/*库函数*/}
	//返回长度的平方
	double len2(){return x*x + y*y;}
	//返回两点的距离
	double distance(Point p){return hypot(x-p.x,y-p.y);}
	Point operator +(const Point &b)const{return Point(x+b.x,y+b.y);}
	Point operator *(const double &k)const{return Point(x*k,y*k);}
	Point operator /(const double &k)const{return Point(x/k,y/k);}
	//`计算pa  和  pb 的夹角`
	//`就是求这个点看a,b 所成的夹角`
	//`测试 LightOJ1203`
	double rad(Point a,Point b)
	{
		Point p = *this;
		return fabs(atan2( fabs((a-p)^(b-p)),(a-p)*(b-p) ));
	}
	//`化为长度为r的向量`
	Point trunc(double r)
	{
		double l = len();
		if(!sgn(l))return *this;
		r /= l;
		return Point(x*r,y*r);
	}
	//`逆时针旋转90度`
	Point rotleft(){return Point(-y,x);}
	//`顺时针旋转90度`
	Point rotright(){return Point(y,-x);}
	//`绕着p点逆时针旋转angle`
	Point rotate(Point p,double angle)
	{
		Point v = (*this) - p;
		double c = cos(angle), s = sin(angle);
		return Point(p.x + v.x*c - v.y*s,p.y + v.x*s + v.y*c);
	}
}; 
struct Line
{
	Point s,e;
	Line(){}
	Line(Point _s,Point _e){s = _s, e = _e;}
	bool operator ==(Line v){return (s == v.s)&&(e == v.e);}
	//`根据一个点和倾斜角angle确定直线,0<=angle<pi`
	Line(Point p,double angle)
	{
		s = p;
		if(sgn(angle-pi/2) == 0){e = (s + Point(0,1));}
		else{e = (s + Point(1,tan(angle)));}
	}
	//ax+by+c=0
	Line(double a,double b,double c)
	{
		if(sgn(a) == 0)	s = Point(0,-c/b), e = Point(1,-c/b);
		else if(sgn(b) == 0) s = Point(-c/a,0), e = Point(-c/a,1);
		else s = Point(0,-c/b), e = Point(1,(-c-a)/b);
	}
	void input()
	{
		s.input();
		e.input();
	}
	void adjust(){if(e < s)swap(s,e);}
	//求线段长度
	double length(){return s.distance(e);}
	//`返回直线倾斜角 0<=angle<pi`
	double angle()
	{
		double k = atan2(e.y-s.y,e.x-s.x);
		if(sgn(k) < 0)k += pi;
		if(sgn(k-pi) == 0)k -= pi;
		return k;
	}
	//`点和直线关系`
	//`1  在左侧`
	//`2  在右侧`
	//`3  在直线上`
	int relation(Point p)
	{
		int c = sgn((p-s)^(e-s));
		if(c < 0)return 1;
		else if(c > 0)return 2;
		else return 3;
	}
	// 点在线段上的判断
	bool pointonseg(Point p){return sgn((p-s)^(e-s)) == 0 && sgn((p-s)*(p-e)) <= 0;}
	//`两向量平行(对应直线平行或重合)`
	bool parallel(Line v){return sgn((e-s)^(v.e-v.s)) == 0;/*两向量叉积为0*/ }
	//`两线段相交判断`
	//`2 规范相交`
	//`1 非规范相交`
	//`0 不相交`
	int segcrossseg(Line v)
	{
		int d1 = sgn((e-s)^(v.s-s));
		int d2 = sgn((e-s)^(v.e-s));
		int d3 = sgn((v.e-v.s)^(s-v.s));
		int d4 = sgn((v.e-v.s)^(e-v.s));
		if( (d1^d2)==-2 && (d3^d4)==-2 )return 2;//如果线段两端点在另一条线段两侧 
		return (d1==0 && sgn((v.s-s)*(v.s-e))<=0) ||//枚举交点 
			(d2==0 && sgn((v.e-s)*(v.e-e))<=0) ||
			(d3==0 && sgn((s-v.s)*(s-v.e))<=0) ||
			(d4==0 && sgn((e-v.s)*(e-v.e))<=0);
	}
	//`直线和线段相交判断`
	//`-*this line   -v seg`
	//`2 规范相交`
	//`1 非规范相交`
	//`0 不相交`
	int linecrossseg(Line v)
	{
		int d1 = sgn((e-s)^(v.s-s));
		int d2 = sgn((e-s)^(v.e-s));
		if((d1^d2)==-2) return 2;
		return (d1==0||d2==0);
	}
	//`两直线关系`
	//`0 平行`
	//`1 重合`
	//`2 相交`
	int linecrossline(Line v)
	{
		if((*this).parallel(v))//此时平行或者重合 
			return v.relation(s)==3;//如果当前直线起点在另一条直线上 
		return 2;
	}
	//`求两直线的交点`
	//`要保证两直线不平行或重合`
	Point crosspoint(Line v)//同底三角形面积比，以v为底  
	{
		double a1 = (v.e-v.s)^(s-v.s);
		double a2 = (v.e-v.s)^(e-v.s);
		return Point((s.x*a2-e.x*a1)/(a2-a1),(s.y*a2-e.y*a1)/(a2-a1));
	}
	//点到直线的距离
	double dispointtoline(Point p){return fabs((p-s)^(e-s))/length();}
	//点到线段的距离
	double dispointtoseg(Point p)
	{
		if(sgn((p-s)*(e-s))<0 || sgn((p-e)*(s-e))<0)
			return min(p.distance(s),p.distance(e));
		return dispointtoline(p);
	}
	//`返回线段到线段的距离`
	//`前提是两线段不相交，相交距离就是0了`
	double dissegtoseg(Line v){return min(min(dispointtoseg(v.s),dispointtoseg(v.e)),min(v.dispointtoseg(s),v.dispointtoseg(e)));}
	//`返回点p在直线上的投影`利用点乘的定义 
	Point lineprog(Point p){return s + ( ((e-s)*((e-s)*(p-s)))/((e-s).len2()) );}
	//`返回点p关于直线的对称点`
	Point symmetrypoint(Point p)
	{
		Point q = lineprog(p);
		return Point(2*q.x-p.x,2*q.y-p.y);
	}
};
struct polygon
{
	int n;
	Point p[maxp];
	Line l[maxp];
	void input(int _n){
		n = _n;
		for(int i = 0;i < n;i++)
			p[i].input();
	}
	void add(Point q){
		p[n++] = q;
	}
	//`得到面积`，点要相邻
	double getarea()
	{
		double sum = 0;
		for(int i = 0;i < n;i++)
			sum += (p[i]^p[(i+1)%n]);
		return fabs(sum)/2;
	}
	//`得到方向`
	//` 1 表示逆时针，0表示顺时针`
	bool getdir(){
		double sum = 0;
		for(int i = 0;i < n;i++)
			sum += (p[i]^p[(i+1)%n]);
		if(sgn(sum) > 0)return 1;
		return 0;
	}
};
//`AB X AC`
double cross(Point A,Point B,Point C){
	return (B-A)^(C-A);
}
//`AB*AC`
double dot(Point A,Point B,Point C){
	return (B-A)*(C-A);
}
//`半平面交`
//`测试 POJ3335 POJ1474 POJ1279`
//***************************
struct halfplane:public Line
{
	double angle;
	halfplane(){}
	//`表示向量s->e逆时针(左侧)的半平面`
	halfplane(Point _s,Point _e)
	{
		s = _s;
		e = _e;
	}
	halfplane(Line v)
	{
		s = v.s;
		e = v.e;
	}
	void calcangle(){angle = atan2(e.y-s.y,e.x-s.x);}
	bool operator <(const halfplane &b)const{return angle < b.angle;}
};

struct halfplanes
{
	int n;//需要输入 
	halfplane hp[maxp];//需要输入，且封闭区域都在向量逆时针方向 
	Point p[maxp];
	int que[maxp];
	int st,ed;//队列的头尾指针，且下标从0开始，指向元素就是头和尾 
	void push(halfplane tmp){hp[n++] = tmp;}
	//去重
	void unique()
	{
		int m = 1;
		for(int i = 1;i < n;i++)
		{
			if(sgn(hp[i].angle-hp[i-1].angle) != 0)
				hp[m++] = hp[i];
			//去除极角相同的情况下，位置在右边(沿向量方向)的边 
			else if(sgn( (hp[m-1].e-hp[m-1].s)^(hp[i].s-hp[m-1].s) ) > 0)
				hp[m-1] = hp[i];
		}
		n = m;
	}
	bool halfplaneinsert()//如果半平面交不存在或者不封闭，返回false 
	{
		for(int i = 0;i < n;i++)hp[i].calcangle();
		sort(hp,hp+n);//先对倾斜角排序 
		unique();
		que[st=0] = 0;
		que[ed=1] = 1;
		p[1] = hp[0].crosspoint(hp[1]);
		for(int i = 2;i < n;i++){
			while(st<ed && sgn((hp[i].e-hp[i].s)^(p[ed]-hp[i].s))<0)ed--;
			while(st<ed && sgn((hp[i].e-hp[i].s)^(p[st+1]-hp[i].s))<0)st++;
			que[++ed] = i;
			if(hp[i].parallel(hp[que[ed-1]]))return false;
			p[ed]=hp[i].crosspoint(hp[que[ed-1]]);
		}
		while(st<ed && sgn((hp[que[st]].e-hp[que[st]].s)^(p[ed]-hp[que[st]].s))<0)ed--;
		while(st<ed && sgn((hp[que[ed]].e-hp[que[ed]].s)^(p[st+1]-hp[que[ed]].s))<0)st++;
		if(st+1>=ed)return false;//最后剩下小于三条直线，表明半平面交不存在 
		return true;
	}
	//`得到最后半平面交得到的凸多边形`
	//`需要先调用halfplaneinsert() 且返回true`
	void getconvex(polygon &con)
	{
		p[st] = hp[que[st]].crosspoint(hp[que[ed]]);
		con.n = ed-st+1;
		for(int j = st,i = 0;j <= ed;i++,j++)
			con.p[i] = p[j];
	}
};

polygon con, final;
halfplanes hp;

signed main()
{
	int n;
	cin >> n;
	con.input(n);
	con.add(con.p[0]);
	//ax + by + c <= 0 
	for(int i = 1; i < n; i++)
	{
		double a = con.p[0].y-con.p[1].y+con.p[i+1].y-con.p[i].y;
		double b = -(con.p[0].x-con.p[1].x+con.p[i+1].x-con.p[i].x);
		double c = con.p[0].x*con.p[1].y-con.p[1].x*con.p[0].y+con.p[i+1].x*con.p[i].y-con.p[i].x*con.p[i+1].y;
		if(b == 0)
		{
			if(a > 0)
				hp.push(halfplane(Point(-c/a, 0), Point(-c/a, 1)));
			else if(a < 0)
				hp.push(halfplane(Point(-c/a, 1), Point(-c/a, 0)));
		}
		else if(b < 0)
			hp.push(halfplane(Point(0, -c/b), Point(1, (-a-c)/b)));
		else
			hp.push(halfplane(Point(1, (-a-c)/b), Point(0, -c/b)));
	}
	for(int i = 0; i < n; i++)
		hp.push(halfplane(con.p[i], con.p[i+1]));
	hp.halfplaneinsert();
	hp.getconvex(final);
	printf("%.4Lf\n", final.getarea()/con.getarea()+eps);
	return 0;
}