BZOJ4445 SCOI2015小凸想跑步（半平面交）

最新推荐文章于 2019-09-21 03:12:45 发布

weixin_30698297

最新推荐文章于 2019-09-21 03:12:45 发布

阅读量127

点赞数

原文链接：http://www.cnblogs.com/Gloid/p/10283928.html

版权

本文深入探讨了几何算法中一个特定问题的解决方案：如何计算一个随机点落在给定多边形内的概率。通过使用半平面交算法和凸包的概念，文章详细解释了数学推导过程，并提供了完整的C++代码实现，展示了如何处理点与多边形的关系，为理解和解决类似几何问题提供了实用的方法。

摘要由CSDN通过智能技术生成

　　考虑怎样的点满足条件。设其为(x_p,y_p)，则要满足(x₀-x_p,y₀-y_p)×(x₁-x_p,y₁-y_p)<=(x_i-x_p,y_i-y_p)×(x_i+1-x_p,y_i+1-y_p)对任意i成立。拆开式子，有(x₀-x_p)*(y₁-y_p)-(y₀-y_p)*(x₁-x_p)<=(x_i-x_p)*(y_i+1-y_p)-(y_i-y_p)*(x_i+1-x_p)，也即x₀y₁-x₀y_p-x_py₁-y₀x₁+y₀x_p+y_px₁<=x_iy_i+1-x_iy_p-x_py_i+1-y_ix_i+1+y_ix_p+y_px_i+1。移项，得(y₀-y₁+y_i+1-y_i)x_p+(x₁-x₀+x_i-x_i+1)y_p<=x_iy_i+1-y_ix_i+1-x₀y₁+y₀x₁。一长串乱七八糟的限制都是对(x_p,y_p)这个二元组的，半平面交即可。当然再套上一个凸包自身的限制。

#include<iostream> 
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;
#define ll long long
#define N 100010
#define vector point
#define double long double
char getc(){char c=getchar();while ((c<'A'||c>'Z')&&(c<'a'||c>'z')&&(c<'0'||c>'9')) c=getchar();return c;}
int gcd(int n,int m){return m==0?n:gcd(m,n%m);}
int read()
{
    int x=0,f=1;char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
    while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return x*f;
}
const double eps=1E-8;
struct point
{
    double x,y;
    vector operator +(const vector&a) const
    {
        return (vector){x+a.x,y+a.y};
    }
    vector operator -(const vector&a) const
    {
        return (vector){x-a.x,y-a.y};
    }
    double operator *(const vector&a) const
    {
        return x*a.y-y*a.x;
    }
    vector operator *(const double&a) const
    {
        return (vector){x*a,y*a};
    }
}a[N];
struct line
{
    point a;vector p;
    bool operator <(const line&a) const
    {
        return atan2(p.x,p.y)>atan2(a.p.x,a.p.y);
    }
}b[N<<1];
int n,m,head,tail;
point p[N<<1];
line q[N<<1];
double area; 
bool onright(point a,line p)
{
    return (a-p.a)*p.p>=0;
}
point cross(line a,line b)
{
    return a.a+a.p*((b.p*(b.a-a.a))/(b.p*a.p));
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("bzoj4445.in","r",stdin);
    freopen("bzoj4445.out","w",stdout);
    const char LL[]="%I64d\n";
#else
    const char LL[]="%lld\n";
#endif
    n=read();
    for (int i=0;i<n;i++) a[i].x=read(),a[i].y=read();
    for (int i=0;i<n;i++) area+=a[i]*a[(i+1)%n];
    for (int i=0;i<n;i++) m++,b[m].a=a[i],b[m].p=a[(i+1)%n]-a[i];
    for (int i=1;i<n;i++)
    {
        double A=(a[0].y-a[1].y+a[(i+1)%n].y-a[i].y),B=(a[1].x-a[0].x+a[i].x-a[(i+1)%n].x),C=(a[i]*a[(i+1)%n]-a[0]*a[1]);
        m++;b[m].p=(vector){-B,A};
        if (fabs(B)>0.5) b[m].a=(point){0,C/B};
        else b[m].a=(point){C/A,0};
    }
    sort(b+1,b+m+1);
    head=tail=1;q[1]=b[1];
    for (int i=2;i<=m;i++)
    {
        while (head<tail&&onright(p[tail],b[i])) tail--;
        while (head<tail&&onright(p[head+1],b[i])) head++;
        q[++tail]=b[i];
        if (fabs(q[tail].p*q[tail-1].p)<eps)
        {
            tail--;
            if (onright(q[tail].a,b[i])) q[tail]=b[i];
        }
        if (head<tail) p[tail]=cross(q[tail],q[tail-1]);
    }
    while (head<tail&&onright(p[tail],q[head])) tail--;
    p[head]=cross(q[head],q[tail]);
    double area2=0;
    for (int i=head;i<tail;i++)
    area2+=p[i]*p[i+1];area2+=p[tail]*p[head];
#undef double
    printf("%.4f",(double)(area2/area));
    return 0;
}