SCOI2015 小凸想跑步 半平面交 向量法

SCOI2015 小凸想跑步

内存限制：256 MiB 时间限制：1000 ms
问题描述

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

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

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

输入格式

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

输出格式

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

样例输入 1

5
1 8
0 7
0 0
8 0
8 8

样例输出 1

0.6316

样例输入 2

4
0 0
7 0
5 5
-1 4

样例输出 2

0.2475

数据范围

3≤N≤10^​5​​,−10^​9​​≤X,Y≤10^​9​​

---

按照题意说的去做，列出一些不等式即可。发现P点可以取的范围是若干半平面的交，答案就是两块面积之比。由于数据范围限制，采用O(nlogn)的半平面交。

注意两点：
1.题目给出了P点在凸多边形内，这个限制要加入到用来求答案的半平面里。
2.如果用直线的一般式，就会出现一些较繁琐的特判。为避免这个采用纯向量推导。

下面是推导过程,方法来自Sparrow：

$由题意得，(A_1-A_0)\times(P-A_1)\leq(A_{i+1}-A_i)\times(P-A_i)$
$打开括号，整理得P\times(A_{i+1}-A_i+A_0-A_1)-(A_1\times A_0-A_{i+1}\times A_i)\leq 0$
$注意到我们想要的“点P在半平面l内的限制”的形式是(P-l.p)\times l.v\leq0$
$展开得P\times l.v-l.p\times l.v\leq0，这与上面的式子形式相似$
$因此考虑将上面的约束转化为点P在半平面l内。那么l.v=A_{i+1}-A_i+A_0-A_1$
$考虑如何求l.p，设k=A_1\times A_0-A_{i+1}\times A_i，那么有l.p\times l.v=k$
$之前已经构造出了l.v，k也是已知的，现在也采取构造法弄出l.p$
$设ver为与l.v垂直、模长与l.v相等、且在l.v的顺时针方向的向量$
$那么有 \frac{k}{|ver|^2}ver\times l.v=k恒成立$
$所以，l.p=\frac{k}{|ver|^2}ver$

这样做就避免了特判。

---

下面的代码给整体下标都加了一。

#include<stdio.h>
#include<algorithm>
#include<cstring>
#include<cmath>
#define MAXN 200005
#define db double
using namespace std;

int N;
db S1,S2;

struct Point{
    db x,y;
}A[MAXN];

typedef Point Vector;

struct Line{
    Vector v;
    Point p;
    db ang;
    void GetAng(){
        ang=atan2(v.y,v.x);
    }
}L[MAXN];int tot;

Vector operator+(Vector a,Vector b){
    a.x+=b.x;a.y+=b.y;
    return a;
}

Vector operator-(Vector a,Vector b){
    a.x-=b.x;a.y-=b.y;
    return a;
}

Vector operator*(db k,Vector a){
    a.x*=k;a.y*=k;
    return a;
}

db operator*(Vector a,Vector b){
    return a.x*b.y-a.y*b.x;
}

bool operator<(Line a,Line b){
    if(a.ang==b.ang)
        return (a.p-b.p)*b.v>0;
    return a.ang<b.ang;
}

db sqlen(Vector a){
    return a.x*a.x+a.y*a.y;
}

bool onleft(Line b,Point a){
    return (a-b.p)*b.v<=0;
}

Point Insc(Line A,Line B){
    db k=(A.p-B.p)*A.v/(B.v*A.v);
    return k*B.v+B.p;
}

Line Q[MAXN];
Point P[MAXN];int cnt;

void HPI(){
    int i,t=0,j,head,tail;
    sort(L+1,L+tot+1);
    for(i=1;i<=tot;i++){
        if(i==1||L[i-1].ang!=L[i].ang)t++;
        L[t]=L[i];
    }
    head=1;tail=1;Q[tail]=L[1];
    for(i=2;i<=t;i++){
        while(head<tail&&!onleft(L[i],Insc(Q[tail],Q[tail-1])))tail--;
        while(head<tail&&!onleft(L[i],Insc(Q[head],Q[head+1])))head++;
        Q[++tail]=L[i];
    }
    while(head<tail&&!onleft(Q[head],Insc(Q[tail],Q[tail-1])))tail--;
    if(head==tail)return;
    Q[tail+1]=Q[head];
    for(i=head;i<=tail;i++)P[++cnt]=Insc(Q[i],Q[i+1]);
}

int main(){
    int i;
    scanf("%d",&N);
    for(i=1;i<=N;i++)scanf("%lf%lf",&A[i].x,&A[i].y);A[N+1]=A[1];
    for(i=1;i<=N;i++)S1+=A[i]*A[i+1];

    for(i=1;i<=N;i++)L[++tot].v=A[i+1]-A[i],L[tot].p=A[i];
    Vector ver;db k;
    for(i=2;i<=N;i++){
        tot++;
        L[tot].v=A[i+1]-A[i]-A[2]+A[1];
        ver.x=L[tot].v.y;ver.y=-L[tot].v.x;
        k=A[i+1]*A[i]-A[2]*A[1];
        L[tot].p=-k/sqlen(ver)*ver;
    }
    for(i=1;i<=tot;i++)L[i].GetAng();
    HPI();
    if(cnt)for(i=1,P[cnt+1]=P[1];i<=cnt;i++)S2+=P[i]*P[i+1];
    printf("%.4lf\n",S2/S1);
}