POJ - 3525 Most Distant Point from the Sea

题目链接：https://vjudge.net/problem/POJ-3525

题意：求半平面交内的最大圆半径。

思路：二分半径，每次让所有的边推进mid，然后再跑半平面交，看是否有解，然后再不断更新即可。然后就只剩下考虑如何将每条边推进R了。

可以把直线表示成Ax + By + C = 0后，把直线向内推进d时，只需要写：C=C+R*sqrt(A*A+B*B)。此时的Ax + By + C = 0就是推进R后的函数了。

在vj上用G++交一直答案错误，用C++交就过了。又是这样。。。。。。。。

用C++交重载减号编译过不了，要重新写个函数改下代码才能过。。。。。。。。

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
const int maxn = 1e4;
const double EPS = 1e-7;
struct point//点
{
    double x,y;
    point friend operator -(point A,point B)
    {
        return {A.x-B.x,A.y-B.y};
    }
};
struct line//线
{
    point x,y;
};
point p[maxn],q[maxn];
line L[maxn],que[maxn],R[maxn];
double chaj(point A,point B)//差积
{
    return A.x*B.y-A.y*B.x;
}
bool fanx(int N)//判断是否为逆序
{
    double sum=0;
    for(int i=2; i<N; i++)
        sum+=chaj(p[i]-p[1],p[i+1]-p[1]);
    return sum<0;//正序
}
double getAngle(point A)//获得极角角度
{
    return atan2(A.y,A.x);
}
bool cmp(line A,line B)//根据极角排序，极角相同，左边的排到前面
{
    double sum1=getAngle(A.y-A.x);
    double sum2=getAngle(B.y-B.x);
    if (fabs(sum1 - sum2) < EPS)
        return chaj(A.y-A.x,B.y-A.x)<=0;
    return sum1<sum2;
}
point getIntersectPoint(line A, line B)//得到两条线的交点
{
    double a1 = A.x.y - A.y.y, b1 = A.y.x - A.x.x, c1 = A.x.x * A.y.y - A.y.x * A.x.y;
    double a2 = B.x.y - B.y.y, b2 = B.y.x - B.x.x, c2 = B.x.x * B.y.y - B.y.x * B.x.y;
    point PO;
    PO.x= (c1*b2-c2*b1)/(a2*b1-a1*b2);
    PO.y= (a2*c1-a1*c2)/(a1*b2-a2*b1);
    return PO;
}
line getLine(line Q,double R)//直线向左平移R,Ax+By+C=0
{
    double A=Q.y.y-Q.x.y;
    double B=Q.x.x-Q.y.x;
    double C=Q.x.y*Q.y.x-Q.y.y*Q.x.x;
    C=C+R*sqrt(A*A+B*B);
    line P;
    if(fabs(Q.x.x-Q.y.x)>0)//线不是竖直的
    {
        P.x.x=Q.x.x;
        P.x.y=-(A*Q.x.x+C)/B;
        P.y.x=Q.y.x;
        P.y.y=-(A*Q.y.x+C)/B;
    }
    else
    {
        P.x.y=Q.x.y;
        P.x.x=-(B*Q.x.y+C)/A;
        P.y.y=Q.y.y;
        P.y.x=-(B*Q.y.y+C)/A;
    }
    return P;
}
bool onRight(line A, line B, line C)//判断 b,c 的交点是否在 a 的右边
{
    point o = getIntersectPoint(B, C);
    if( chaj( A.y - A.x,o - A.x )<0 )
        return true;
    return false;
}
bool HalfPlaneIntersection(int cnt)
{
    int head=0,tail=0;
    for(int i=1; i<=cnt; i++)
    {
        //判断新加入直线产生的影响
        while(tail-head>1 && onRight(L[i], que[tail - 1], que[tail - 2]))
            tail--;
        while(tail-head>1 && onRight(L[i], que[head], que[head + 1]))
            head++;
        que[tail++]=L[i];
    }
    //最后判断最先加入的直线和最后的直线的影响
    while(tail-head>1 && onRight(que[head], que[tail-1], que[tail-2]))
        tail--;
    while(tail-head>1 && onRight(que[tail-1], que[head], que[head+1]))
        head++;
    //半平面交存在que中， head到tail-1
    if (tail-head >= 3)
        return true;
    return false;
}
int main()
{
    int N;
    while(~scanf("%d",&N))
    {
        if(N==0)
            break;
        for(int i=1; i<=N; i++)
            cin>>p[i].x>>p[i].y;
        for(int i=1; i<N; i++)
        {
            R[i].x=p[i];
            R[i].y=p[i+1];
        }
        R[N].x=p[N];
        R[N].y=p[1];
        sort(R+1,R+N+1,cmp);
        int cnt=1;//去重后的个数
        for(int i=2; i<=N; i++)
        {
            if (fabs(getAngle(R[i].y-R[i].x) - getAngle(R[i-1].y-R[i-1].x)) < EPS)
                continue;
            R[++cnt]=R[i];
        }
        double l=0,r=10000000;//二分左右
        while(r-l>EPS)
        {
            double mid=(l+r)/2.0;
            for(int i=1;i<=cnt;i++)
                L[i]=getLine(R[i],mid);
            if(HalfPlaneIntersection(cnt))
                l=mid;
            else
                r=mid;
        }
        printf("%.6lf\n",r);
    }
}