Newcoder 141 J.Distance to Work（二分+计算几何）

Description

给出$n$个顶点简单多边形每点坐标，$m$次查询，每次给出圆心坐标，求圆的半径使得该多边形不在圆内部分面积占全部面积的$\frac{P}{Q}$

Input

第一行输入一整数$n$表示简单多边形点数，之后顺序输入$n$个顶点的横纵坐标$x_i,y_i$，之后输入一整数$m$表示查询，每组查询输入四个整数$x_0,y_0,P,Q$

$(3\le n\le 200,1\le m\le 200,1\le P<Q\le 200,|x_0|,|y_0|,|x_i|,|y_i|\le 10^3$

Output

对于每组查询，输出圆的半径

Sample Input

4
0 0
2 0
2 2
0 2
1
1 1 1 2

Sample Output

0.797884560809

Solution

二分半径求出圆与多边形相交面积判断是否合法即可

Code

#include<cstdio>
#include<cmath>
#include<iostream>
#include<iomanip>
using namespace std;
#define maxn 250
typedef long double ld;
const ld eps=1e-10;
struct node
{
    ld x,y;
}p[maxn],q,pp[maxn];
int sgn(ld x)
{
    if(fabs(x)<eps)return 0;
    if(x>eps)return 1;
    return -1;
}
ld multi(node a,node b)//向量叉乘
{
    return a.x*b.y-a.y*b.x;
}
int get_sign(node a,node b)//判断向量叉乘符号
{
    if(multi(a,b)>0) return 1;
    return -1;
}
ld dis(node a,node b)//两点间距离
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
ld dis_line(node a,node b)//求原点到直线ab的距离
{
    return fabs(multi(a,b))/dis(a,b);
}
ld get_angle(ld a,ld b,ld c)//求边c对角的值
{
    return (a*a+b*b-c*c)/(2.0*a*b);
}
node get_point(node a,ld r)//求a点与原点连线与圆的交点坐标
{
    node ans;
    ld k;
    if(a.x!=0)//斜率存在
    {
        k=a.y/a.x;
        ans.x=fabs(r)/sqrt(1.0+k*k);
        if(a.x<0) ans.x=-ans.x;//判断ans.x的符号
        ans.y=k*ans.x;
    }
    else//斜率不存在
    {
        ans.x=0;
        if(a.y>0) ans.y=r;//判断ans.y的符号
        else ans.y=-r;
    }
    return ans;
}
ld get_area(node a,node b,node o,ld r)//求a,b,o三点所形成的三角形与圆相交的有向面积
{
    int sign=get_sign(a,b);//判断有向面积的符号
    ld ans=0;
    ld oa=dis(o,a),ob=dis(o,b),ab=dis(a,b);
    ld l=dis_line(a,b);
    if(oa==0||ob==0)//有一边边长为0则面积为0
        return 0;
    //第一种情况,两点都在圆内
    if(oa<=r&&ob<=r)
    {
        ans=fabs(multi(a,b))/2.0;//相交面积即为三角形oab的面积
        return sign*ans;
    }
    //第二种情况,两点都在圆外且圆心到直线ab距离大于半径
    else if(oa>=r&&ob>=r&&l>=r)
    {
        node t1=get_point(a,r);//oa与圆的交点
        node t2=get_point(b,r);//ob与圆的交点
        ld d=dis(t1,t2);
        ld ang=acos(get_angle(r,r,d));//角t1ot2
        ans=fabs(ang*r*r/2.0);//相交面积即为扇形ot1t2的面积
        return sign*ans;
    }
    //第三种情况,两点都在圆外且圆心到直线距离小于半径但三角形oab有一底角是钝角即直线ab与圆没有交点
    else if(oa>=r&&ob>=r&&l<=r&&(get_angle(ab,oa,ob)<=0||get_angle(ab,ob,oa)<=0))
    {
        node t1=get_point(a,r);//oa与圆的交点
        node t2=get_point(b,r);//ob与圆的交点
        ld dist=dis(t1,t2);//线段t1t2长度
        ld ang=acos(get_angle(r,r,dist));//角t1ot2
        ans=fabs(ang*r*r/2.0);//相交面积即为扇形ot1t2的面积
        return sign*ans;
    }
    //第四种情况,两点都在圆外且原点到直线ab距离小于半径且三角形oab底角均为锐角即直线ab与圆有两个交点
    else if(oa>=r&&ob>=r&&l<=r&&get_angle(ab,oa,ob)>0&&get_angle(ab,ob,oa)>0)
    {
        node c,d;//c,d为直线ab与圆的两个交点
        if(a.x!=b.x)//直线ab斜率存在
        {
            ld k=(a.y-b.y)/(a.x-b.x);//直线ab斜率
            ld h=a.y-k*a.x;//直线ab截距
            //解一元二次方程求圆与直线ab的交点
            ld a0=1.0+k*k;
            ld b0=2.0*k*h;
            ld c0=h*h-r*r;
            c.x=(-b0+sqrt(b0*b0-4.0*a0*c0))/(2.0*a0);
            c.y=k*c.x+h;
            d.x=(-b0-sqrt(b0*b0-4.0*a0*c0))/(2.0*a0);
            d.y=k*d.x+h;
        }
        else//直线ab斜率不存在
        {
            c.x=a.x;
            d.x=a.x;
            c.y=sqrt(r*r-a.x*a.x);
            d.y=-sqrt(r*r-a.x*a.x);
        }
        node t1=get_point(a,r);//oa与圆的交点
        node t2=get_point(b,r);//ob与圆的交点
        ld d1=dis(c,d);//线段cd长度
        ld d2=dis(t1,t2);//线段t1t2长度
        ld ang1=acos(get_angle(r,r,d1));//角cod
        ld ang2=acos(get_angle(r,r,d2));//角t1ot2
        ld s1=fabs(ang1*r*r/2.0);//小扇形ocd面积
        ld s2=fabs(ang2*r*r/2.0);//小大扇形ot1t2面积
        ld s3=fabs(multi(c,d))/2.0;//三角形ocd面积
        ans=s2+s3-s1;//相交面积即为三角形ocd面积加上两个扇形面积之差
        return sign*ans;
    }
    //第五种情况,a点在圆外,b点在圆内
    else if(oa>=r&&ob<=r)
    {
        node c,d,e;//c,d两点为直线ab与圆的两个交点,e点为介于ab两点之间的那个交点
        if(a.x!=b.x)//直线ab斜率存在
        {
            ld k=(a.y-b.y)/(a.x-b.x);//直线ab斜率
            ld h=a.y-k*a.x;//直线ab截距
            //解一元二次方程求圆与直线ab的交点
            ld a0=1.0+k*k;
            ld b0=2.0*k*h;
            ld c0=h*h-r*r;
            c.x=(-b0+sqrt(b0*b0-4.0*a0*c0))/(2.0*a0);
            c.y=k*c.x+h;
            d.x=(-b0-sqrt(b0*b0-4.0*a0*c0))/(2.0*a0);
            d.y=k*d.x+h;
            //交点应为c,d两点中横坐标介于(a.x,b.x)两点之间的那一个
            if(a.x<=c.x&&c.x<=b.x||a.x>=c.x&&c.x>=b.x) e=c;
            else e=d;
        }
        else//直线ab斜率不存在
        {
            c.x=d.x=a.x;
            c.y=-sqrt(r*r-a.x*a.x);
            d.y=sqrt(r*r-a.x*a.x);
            //交点应为c,d两点中纵坐标介于(a.y,b.y)两点之间的那一个
            if(a.y<=c.y&&c.y<=b.y||a.y>=c.y&&c.y>=b.y) e=c;
            else e=d;
        }
        node t1=get_point(a,r);//oa与圆的交点
        ld dist=dis(t1,e);//线段t1e长度
        ld ang=acos(get_angle(r,r,dist));//角eot1
        ld s1=fabs(ang*r*r/2.0);//扇形oet1面积
        ld s2=fabs(multi(b,e))/2.0;//三角形obe面积
        ans=s1+s2;//相交面积即三角形obe面积加上扇形oet1面积
        return sign*ans;
    }
    //第五种情况,b点在圆外,a点在圆内
    else if(ob>=r&&oa<=r)
    {
        node c,d,e;//c,d两点为直线ab与圆的两个交点,e点为介于ab两点之间的那个交点
        if(a.x!=b.x)//直线ab斜率存在
        {
            ld k=(a.y-b.y)/(a.x-b.x);//直线ab斜率
            ld h=a.y-k*a.x;//直线ab截距
            //解一元二次方程求圆与直线ab的交点
            ld a0=1.0+k*k;
            ld b0=2.0*k*h;
            ld c0=h*h-r*r;
            c.x=(-b0+sqrt(b0*b0-4.0*a0*c0))/(2.0*a0);
            c.y=k*c.x+h;
            d.x=(-b0-sqrt(b0*b0-4.0*a0*c0))/(2.0*a0);
            d.y=k*d.x+h;
            //交点应为c,d两点中横坐标介于(a.x,b.x)两点之间的那一个
            if(a.x<=c.x&&c.x<=b.x||a.x>=c.x&&c.x>=b.x) e=c;
            else e=d;
        }
        else//直线ab斜率不存在
        {
            c.x=d.x=a.x;
            c.y=-sqrt(r*r-a.x*c.x);
            d.y=sqrt(r*r-a.x*a.x);
            //交点应为c,d两点中纵坐标介于(a.y,b.y)两点之间的那一个
            if(a.y<=c.y&&c.y<=b.y||a.y>=c.y&&c.y>=b.y) e=c;
            else e=d;
        }
        node t1=get_point(b,r);//oa与圆的交点
        ld dist=dis(t1,e);//线段t1e长度
        ld ang=acos(get_angle(r,r,dist));//角eot1
        ld s1=fabs(ang*r*r/2.0);//扇形oet1面积
        ld s2=fabs(multi(a,e))/2.0;//三角形oae面积
        ans=s1+s2;//相交面积即三角形oae面积加上扇形oet1面积
        return sign*ans;
    }
    return 0;
}
ld solve(node p[],int n,node o,ld r)//p数组为多边形个点坐标,n为多边形顶点数,o为圆心坐标,r为圆的半径
{
    ld ans=0.0;
    for(int i=0;i<n;i++)//移动坐标系使得圆心在原点
    {
        p[i].x-=o.x;
        p[i].y-=o.y;
    }
    o.x=o.y=0;
    for(int i=0;i<n;i++)
    {
        int j=(i+1)%n;
        ans+=get_area(p[i],p[j],o,r);
    }
    return fabs(ans);
}
ld get_S(int n)
{
    ld ans=0;
    for(int i=0;i<n;i++)ans+=multi(p[i],p[(i+1)%n]);
    return fabs(ans)*0.5;
}
int main()
{
    int n;
    scanf("%d",&n);
    for(int i=0;i<n;i++)cin>>p[i].x>>p[i].y;
    ld S=get_S(n);
    int m;
    scanf("%d",&m);
    while(m--)
    {
        ld P,Q;
        cin>>q.x>>q.y>>P>>Q;
        ld area=S*(Q-P)/Q;
        ld l=0,r=1e7,mid,ans=1e7;
        int T=200;
        while(T--)
        {
            mid=0.5*(l+r);
            for(int i=0;i<n;i++)pp[i]=p[i];
            if(sgn(solve(pp,n,q,mid)-area)>=0)r=mid;
            else l=mid;
        }
        cout<<setprecision(10)<<mid<<endl;
    }
    return 0;
}