zoj 1670 Jewels from Heaven

题意：三个人，在给定正方形内，求第一个人拿到珠宝的概率。珠宝随机出现在正方形内。

思路：中垂线+半平面相交。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<string>
#include<cmath>
#include<vector>
using namespace std;
const int maxn=1e5+7;
const double eps=1e-8;
const double pi=acos(-1);

double 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) { }
};

typedef Point Vector;

Vector operator + (const Point& A, const Point& B)
{
    return Vector(A.x+B.x, A.y+B.y);
}

Vector operator - (const Point& A, const Point& B)
{
    return Vector(A.x-B.x, A.y-B.y);
}

Vector operator * (const Point& A, double v)
{
    return Vector(A.x*v, A.y*v);
}

Vector operator / (const Point& A, double v)
{
    return Vector(A.x/v, A.y/v);
}

double Cross(const Vector& A, const Vector& B)
{
    return A.x*B.y - A.y*B.x;
}

double Dot(const Vector& A, const Vector& B)
{
    return A.x*B.x + A.y*B.y;
}

double Length(const Vector& A)
{
    return sqrt(Dot(A,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));
}

bool operator < (const Point& p1, const Point& p2)
{
    return p1.x < p2.x || (p1.x == p2.x && p1.y < p2.y);
}

bool operator == (const Point& p1, const Point& p2)
{
    return p1.x == p2.x && p1.y == p2.y;
}

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& L) const
    {
        return ang < L.ang;
    }
};

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


Point GetLineIntersection(const Line& a, const 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;
}

const double INF = 1e8;

Point ansPoly[maxn];
int HalfplaneIntersection(vector<Line> L)     //L为切割平面的直线集合，求半平面交，返回点的个数，点存在anspoly数组中
{
    int n = L.size();
    sort(L.begin(), L.end()); // 按极角排序
    int first, last;         // 双端队列的第一个元素和最后一个元素的下标
    vector<Point> p(n);      // p[i]为q[i]和q[i+1]的交点
    vector<Line> q(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] = GetLineIntersection(q[last-1], q[last]);
    }
    while(first < last && !OnLeft(q[first], p[last-1])) last--; //
    if(last - first <= 1) return 0; //
    p[last] = GetLineIntersection(q[last], q[first]); //
    // 从deque复制到输出中
    int index=0;
    for(int i = first; i <= last; i++) ansPoly[index++]=p[i];
    return index;
}

double PolygonArea(int n,Point *p)
{
    double area=0;
    for(int i=1; i<n-1; i++)
        area+=Cross(p[i]-p[0],p[i+1]-p[0]);
    return area/2;
}

Point p[5];
int main()
{
//  freopen("in.txt","r",stdin);
    while(cin>>p[0].x>>p[0].y>>p[1].x>>p[1].y>>p[2].x>>p[2].y)
    {
        if(p[0].x==p[1].x&&p[0].y==p[1].y&&p[0].x==0)break;
//        Point zd1,zd2;
        vector<Line>  vec;
        vec.push_back(Line(Point(0,0),Point(1,0)));
        vec.push_back(Line(Point(10000,0),Point(0,1)));
        vec.push_back(Line(Point(10000,10000),Point(-1,0)));
        vec.push_back(Line(Point(0,10000),Point(0,-1)));
        Vector v=(p[1]-p[0]);
        vec.push_back(Line((p[1]+p[0])*0.5,Normal(v)));
        v=(p[2]-p[0]);
        vec.push_back(Line((p[2]+p[0])*0.5,Normal(v)));
        int m=HalfplaneIntersection(vec);
        double ans=PolygonArea(m,ansPoly);
        printf("%.3f
",ans/(1.0*10000*10000));
    }
    return 0;
}