【LuoguP4081】[SCOI2015]小凸想跑步
题意
给你一个凸多边形,求出在其内部选择一个点,这个点与最开始输入的两个点形成的三角形是以该点对凸多边形三角剖分的三角形中面积最小的一个三角形的概率。
Sol
答案就是 可行域面积与该凸多边形面积之比。
通过数学方法列出第一个三角形和其他三角形面积关系的式子,解出来发现都是一个半平面,所以我们要做的就是快速求解半平面交。
把所有要加入的直线用向量表示 , 按照极角排序 ( 用 atan2() ) , 然后依次加入直线。
维护一个双端队列 , 每次加入一条直线时判断最左最右的交点和当前直线的位置关系,如果点在不在加入向量的左侧那么要把最左或最右的直线弹掉。
加入完后还要判一下最右边交点和最左边的向量是否满足左右关系。
求完半平面交之后算个面积除一下就行了。
#include<bits/stdc++.h>
using namespace std;
#define Set(a,b) memset(a,b,sizeof(a))
template<class T>inline void init(T&x){
x=0;char ch=getchar();bool t=0;
for(;ch>'9'||ch<'0';ch=getchar()) if(ch=='-') t=1;
for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+(ch-48);
if(t) x=-x;return;
}typedef long long ll;
typedef double R;
#define Sqr(a) ((a)*(a))
const R PI=acos(-1),eps=1e-9,INF=1e9;
const int N=2e5+10;
struct point{
R x,y;
point(R _x=0.0,R _y=0.0) {x=_x,y=_y;}
inline R operator *(const point b){return x*b.y-y*b.x;}
inline point operator *(const R d){return point(x*d,y*d);}
inline point operator /(const R d){return point(x/d,y/d);}
inline point operator +(const point b){return point(x+b.x,y+b.y);}
inline point operator -(const point b){return point(x-b.x,y-b.y);}
inline void output(){printf("( %lf , %lf )",x,y);return;}
inline R len(){return sqrt(Sqr(x)+Sqr(y));}
}P[N];
struct line{
point A,B;R ang;
line(point _A=point(),point _B=point()){A=_A,B=_B;ang=atan2(B.y,B.x);}
inline bool operator <(const line b)const{return ang<b.ang;}
}L[N];int cnt=0;
int n;
inline int fcmp(R r){if(r>eps) return 1;else if(r<-eps) return -1;return 0;}
inline bool Left(point A,point B){return fcmp(A*B)>0;}
inline bool isleft(point A,line B){return fcmp(B.B*(A-B.A))>0;}
point Cr[N];R S=0;
inline point Inter(line L1,line L2){return L1.A+L1.B*((L2.B*(L2.A-L1.A))/(L2.B*L1.B));}
inline void HalfPlane(){
int l,r=1;sort(L+1,L+1+cnt);// 按照斜率(极角)排序
for(int i=2;i<=cnt;++i) {
if(fcmp(L[i].ang-L[r].ang)) L[++r]=L[i];
else if(isleft(L[i].A,L[r])) L[r]=L[i];
}
cnt=r,l=r=1;
for(int i=2;i<=cnt;++i) {// 类似维护凸包
while(l<r&&!isleft(Cr[r],L[i])) --r;
while(l<r&&!isleft(Cr[l+1],L[i])) ++l;// 交点必须在内部, 否则上一个向量没有用
L[++r]=L[i];
if(l<r) Cr[r]=Inter(L[r],L[r-1]);
}
while(l<r&&!isleft(Cr[r],L[l])) --r;
Cr[l]=Cr[r+1]=Inter(L[l],L[r]);
R area=0;
for(int i=l+2;i<=r;++i) area+=(Cr[i-1]-Cr[l])*(Cr[i]-Cr[l]);
printf("%.4lf\n",area/S);
}
int main()
{
init(n);S=0;
for(int i=1;i<=n;++i) {scanf("%lf %lf",&P[i].x,&P[i].y);}
for(int i=3;i<=n;++i) S+=(P[i-1]-P[1])*(P[i]-P[1]);
P[n+1]=P[1];
for(int i=1;i<=n;++i) {
point Q=P[i+1]-P[i];
L[++cnt]=line(P[i],Q);
}
R kx=P[1].y-P[2].y,ky=P[2].x-P[1].x,re=P[1].x*P[2].y-P[2].x*P[1].y;
for(int i=2;i<=n;++i){
R kxx=P[i].y-P[i+1].y,kyy=P[i+1].x-P[i].x,ree=P[i].x*P[i+1].y-P[i+1].x*P[i].y;
R kY=ky-kyy,kX=kxx-kx,Kre=ree-re;
if(kY) L[++cnt]=line(point(0,Kre/kY),point(-kY,-kX));
else if(kX) L[++cnt]=line(point(-Kre/kX,0),point(0,-kX));
}HalfPlane();
return 0;
}