【BZOJ4445】[SCOI2015]小凸想跑步(半平面交)
【BZOJ4445】[SCOI2015]小凸想跑步(半平面交)
题面
题解
首先把点给设出来,\(A(x_a,y_a),B(x_b,y_b),C(x_c,y_c),D(x_d,y_d),P(x,y)\)
然后我们考虑\(S_\Delta ABP<S_\Delta CDP\)什么情况下满足。
根据点积来求面积,得到:
\[(x_a-x,y_a-y)\times(x_b-x,y_b-y)<(x_c-x,y_c-y)\times(x_d-x,y_d-y)
\]
这个东西左边拆开之后得到:
\[\begin{aligned}
&\ \ \ \ (x_a-x)(y_b-y)-(x_b-x)(y_a-y)\\
&=x(y_a-y_b)-y(x_a-x_b)+x_ay_b-x_by_a
\end{aligned}\]
右侧类似。然后就可以移项,得到:
\[x(y_a-y_b-y_c+y_d)-y(x_a-x_b-x_c+x_d)+x_ay_b-x_by_a-x_cy_d+x_dy_c<0
\]
那么对于相邻的两个点和\(0,1\)两个点进行一次比较,每次都可以得到一个半平面,最后求解这个半平面交就是结果了。
#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
#define MAX 100100
const double eps=1e-7;
inline int read()
{
int x=0;bool t=false;char ch=getchar();
while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
if(ch=='-')t=true,ch=getchar();
while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
return t?-x:x;
}
struct Node{double x,y;}p[MAX],Qp[MAX<<1];
Node operator+(Node a,Node b){return (Node){a.x+b.x,a.y+b.y};}
Node operator-(Node a,Node b){return (Node){a.x-b.x,a.y-b.y};}
Node operator*(Node a,double b){return (Node){a.x*b,a.y*b};}
double operator*(Node a,Node b){return a.x*b.x+a.y*b.y;}
double cross(Node a,Node b){return a.x*b.y-a.y*b.x;}
struct Line{Node p,v;double alpha;}S[MAX<<1],Q[MAX<<1];int tot;
bool cmp(Line a,Line b){return a.alpha<b.alpha;}
Node Intersection(Line a,Line b)
{
Node c=b.p-a.p;
double t=cross(b.v,c)/cross(b.v,a.v);
return a.p+a.v*t;
}
int zero(double x)
{
if(fabs(x)<eps)return 0;
return x>0?1:-1;
}
bool Left(Node a,Line b){return zero(cross(b.v,a-b.p))>0;}
int n;double area;
int l,r;
bool HalfPlaneIntersection()
{
for(int i=1;i<=tot;++i)S[i].alpha=atan2(S[i].v.y,S[i].v.x);
sort(&S[1],&S[tot+1],cmp);
Q[l=r=1]=S[1];
for(int i=2;i<=tot;++i)
{
while(l<r&&!Left(Qp[r-1],S[i]))--r;
while(l<r&&!Left(Qp[l],S[i]))++l;
if(zero(cross(Q[r].v,S[i].v))==0)
Q[r]=Left(Q[r].p,S[i])?Q[r]:S[i];
else Q[++r]=S[i];
if(l<r)Qp[r-1]=Intersection(Q[r],Q[r-1]);
}
while(l<r&&!Left(Qp[r-1],Q[l]))--r;
return (r-l)>1;
}
int main()
{
n=read();
for(int i=0;i<n;++i)p[i].x=read(),p[i].y=read();p[n]=p[0];
for(int i=0;i<n;++i)S[++tot]=(Line){p[i],p[i+1]-p[i]};
for(int i=1;i<n;++i)area+=fabs(cross(p[i]-p[0],p[i+1]-p[0]));
for(int i=1;i<n;++i)
{
double a=p[0].y-p[1].y-p[i].y+p[i+1].y;
double b=-p[0].x+p[1].x+p[i].x-p[i+1].x;
double c=cross(p[0],p[1])-cross(p[i],p[i+1]);
Line d;
if(fabs(a)>eps)d=(Line){(Node){-c/a,0},(Node){-b,a}};
else d=(Line){(Node){0,-c/b},(Node){-b,a}};
S[++tot]=d;
}
HalfPlaneIntersection();
double ans=0;Qp[r]=Intersection(Q[l],Q[r]);
for(int i=l+1;i<r;++i)ans+=fabs(cross(Qp[i]-Qp[l],Qp[i+1]-Qp[l]));
printf("%.4lf\n",ans/area);
return 0;
}