LA 3890 半平面交

二分查询答案，判断每一个新形成的向量合在一块能否形成半平面交

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
#define N 110
#define eps 1e-7

int dcmp(double x)
{
    if(fabs(x)<eps) return 0;
    else return x<0?-1:1;
}

struct Point{
    double x , y;
    Point(double x=0 , double y=0):x(x),y(y){}
}po[N] , poly[N];

typedef Point Vector;
Vector vec[N]; //记录每条边上对应的法向量

Vector operator+(Vector a , Vector b) { return Vector(a.x+b.x , a.y+b.y); }
Vector operator-(Point a , Point b) { return Vector(a.x-b.x , a.y-b.y); }
Vector operator*(Vector a , double b) { return Vector(a.x*b , a.y*b); }
Vector operator/(Vector a , double b) { return Vector(a.x/b , a.y/b); }
bool operator==(const Point &a , const Point &b) { return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0; }

double Dot(Vector a  , Vector b) { return a.x*b.x+a.y*b.y; }
double Cross(Vector a , Vector b) { return a.x*b.y - b.x*a.y; }
double Length(Vector a) { return sqrt(Dot(a,a)); }
double Angle(Vector A , Vector B) { return acos(Dot(A,B)) / Length(A) / Length(B); }
double Area2(Point A , Point B , Point C) { return Cross(B-A , C-A); }

Vector Rotate(Vector A , double rad) { return Vector(A.x*cos(rad)-A.y*sin(rad) , A.x*sin(rad)+A.y*cos(rad)); }

Vector Normal(Vector a)
{
    double l = Length(a);
    return Vector(-a.y/l , a.x/l);
}

struct Line{
    Point P;
    Vector v;
    double ang;
    Line(){}
    Line(Point P , Vector v):P(P),v(v){ang = atan2(v.y,v.x);}
    bool operator<(const Line &m)const {
        return ang<m.ang;
    }
}line[N] , L[N];

bool OnLeft(Line L , Point P)
{
    return Cross(L.v , P-L.P) > 0;
}

Point GetIntersection(Line a , Line b)
{
    Vector u = a.P-b.P;
    double t = Cross(b.v , u) / Cross(a.v , b.v);
    return a.P+a.v*t;
}

/***半平面交的主过程，返回形成半平面交点的个数，无法形成就返回0***/
int HalfplaneIntersection(Line *L , int n , Point *poly)
{
    sort(L , L+n);
    int first , last;
    Point *p = new Point[n];
    Line *q = new Line[n];
    q[first=last=0] = L[0];
    for(int i=1 ; i<n ; i++)
    {
        while(first<last && !OnLeft(L[i] , p[last-1])) last--;
        while(first<last && !OnLeft(L[i] , p[first])) first++;
        q[++last] = L[i];
        if(fabs(Cross(q[last].v , q[last-1].v)) < eps)
        {
            last--;
            if(OnLeft(q[last] , L[i].P)) q[last]=L[i];
        }
        if(first<last) p[last-1] = GetIntersection(q[last-1] , q[last]);
    }
    while(first < last && !OnLeft(q[first] , p[last-1])) last--;
    //删除无用平面
    if(last-first<=1) return 0;
    p[last] = GetIntersection(q[last] , q[first]);

    //从deque复制到输出中
    int m=0;
    for(int i=first ; i<=last ; i++) poly[m++] = p[i];
    return m;
}

int main()
{
   // freopen("a.in" , "r" , stdin);
    int n;
    while(scanf("%d" , &n) , n)
    {
        for(int i=0 ; i<n ; i++)
            scanf("%lf%lf" , &po[i].x , &po[i].y);

        for(int i=0 ; i<n ; i++) vec[i] = Normal(po[(i+1)%n]-po[i]);
        double l=0 , r=20000;
        while(r-l>eps){
            double m=(l+r)/2;
            for(int i=0 ; i<n ; i++) L[i] = Line(po[i]+vec[i]*m , po[(i+1)%n]-po[i]);
            if(HalfplaneIntersection(L , n , poly)>0) l=m;
            else r=m;
        }
        printf("%.6f
" , l);
    }
    return 0;
}