[半平面交]小凸想跑步 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=0Line(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)));elsehp.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;
}

[半平面交]小凸想跑步 LibreOJ2008相关推荐

「SCOI2015」小凸想跑步解题报告
「SCOI2015」小凸想跑步最开始以为和多边形的重心有关,后来发现多边形的重心没啥好玩的性质实际上你把面积小于的不等式列出来,发现是一次的,那么就可以半平面交了 Code: #include & ...
[Scoi2015]小凸玩矩阵
bzoj 4443: [Scoi2015]小凸玩矩阵 http://www.lydsy.com/JudgeOnline/problem.php?id=4443 Time Limit: 10 Sec ...
LibreOJ #2006. 「SCOI2015」小凸玩矩阵二分答案+二分匹配
#2006. 「SCOI2015」小凸玩矩阵内存限制:256 MiB时间限制:1000 ms标准输入输出题目类型:传统评测方式:文本比较上传者: 匿名提交提交记录统计讨论测试数据题目描述小 ...
bzoj 4443 [Scoi2015]小凸玩矩阵网络流,二分
[Scoi2015]小凸玩矩阵 Time Limit: 10 Sec Memory Limit: 128 MB Submit: 1564 Solved: 734 [Submit][Status][ ...
[SCOI2015]小凸玩矩阵（匈牙利+二分）
description 题目描述小凸和小方是好朋友,小方给小凸一个 N×M(N≤M)的矩阵 A,要求小凸从其中选出 N 个数,其中任意两个数字不能在同一行或同一列,现小凸想知道选出来的 N 个数中第 ...
【bzoj4443】[Scoi2015]小凸玩矩阵二分+二分图最大匹配
题目描述小凸和小方是好朋友,小方给小凸一个N*M(N<=M)的矩阵A,要求小秃从其中选出N个数,其中任意两个数字不能在同一行或同一列,现小凸想知道选出来的N个数中第K大的数字的最小值是多少. ...
bzoj 4443: [Scoi2015]小凸玩矩阵（二分+二分匹配）
4443: [Scoi2015]小凸玩矩阵 Time Limit: 10 Sec Memory Limit: 128 MB Submit: 1002 Solved: 505 [Submit][St ...
bzoj-4433 小凸玩矩阵（二分图，二分+匈牙利）
4443: [Scoi2015]小凸玩矩阵 Description 小凸和小方是好朋友,小方给小凸一个N*M(N<=M)的矩阵A,要求小秃从其中选出N个数,其中任意两个数字不能在同一行或同一列, ...
UVA1396 Most Distant Point from the Sea（AM - ICPC - Tokyo - 2007）（计算几何，半平面交 + 二分答案）
整理的算法模板合集: ACM模板题目传送门见<训练指南>P279 很明显就是一个二分答案,它问的是最远的点,直接枚举因为这里都是double类型的数所以有无限个点,我们可以直接二分. ...

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

[半平面交]小凸想跑步 LibreOJ2008相关推荐

最新文章

热门文章