HDU5130 Signal Interference

/*
 HDU5130 Signal Interference
 http://acm.hdu.edu.cn/showproblem.php?pid=5130
 计算几何 圆与多边形面积交
 *
 */


#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
const double Pi=acos(-1.0);
const double eps = 1e-8;
const int Nmax=1005;
double k;
int sgn(double x)
{ 
    if(x<-eps)
        return -1;
    if(x>eps)
        return 1;
    return 0;
}
double Abs(double x)
{
    if(sgn(x)<0)
        x=-x;
    else if(sgn(x)==0)
        x=0.0;
    return x;
}
double sqr(double x)
{
    return x*x;
}
double Sqrt(double x) 
{ 
    return max(0.0,sqrt(x)); 
}
struct Pt 
{
    double x,y;
    Pt() { }
    Pt(double x, double y) : x(x), y(y) { }
    
    Pt operator - (const Pt &b)
    {
        return Pt(x-b.x,y-b.y);
    }
    Pt operator + (const Pt &b)
    {
        return Pt(x+b.x,y+b.y);
    }
    friend Pt operator * (double a,Pt p)
    {
        return Pt(p.x*a,p.y*a);
    }
    friend Pt operator * (Pt p,double a)
    {
        return Pt(p.x*a,p.y*a);
    }
    friend Pt operator / (Pt p,double a)
    {
        return Pt(p.x/a,p.y/a);
    }
    double norm()
    { 
        return Sqrt(x*x + y*y); 
    }
    double len()
    {
        return norm();
    }
    void print() 
    { 
        printf("(%f, %f)\n", x, y); 
    }
};
int n;
Pt pl[Nmax];
Pt pr,A,B,p0,p1;
double dist (Pt a, Pt b) 
{    
    return (a-b).norm(); 
}
double dot(Pt a,Pt b)//点乘 
{
    return a.x*b.x+a.y*b.y; 
}
double det(Pt a,Pt b)
{
    return a.x*b.y-a.y*b.x;
}
Pt rotate(Pt p,double a)//P点绕原点逆时针旋转a弧度
{
    return Pt(p.x*cos(a)-p.y*sin(a),p.x*sin(a)+p.y*cos(a));
}
struct Sg//线段 
{
    Pt s, t;
    Sg() { }
    Sg(Pt s, Pt t) : s(s), t(t) { }
    Sg(double a, double b, double c, double d) : s(a, b), t(c, d) { }
};
bool PtOnSegment(Pt p, Pt a, Pt b) //p是否在线段ab上，把<=改成<就能实现不含线段端点的点在线段上的判断。
{
    return !sgn(det(p-a, b-a)) && sgn(dot(p-a, p-b)) <= 0;
}
bool PtOnLine(Pt p, Pt a, Pt b) //p是否在直线ab上
{
    return !sgn(det(p-a, b-a));
}
Pt PtLineProj(Pt s, Pt t, Pt p) //p到直线st的投影
{
    double r = dot(p-s, t-s) / (t - s).norm();
    return s + (t - s) * r;
}
bool parallel(Pt a, Pt b, Pt s, Pt t) 
{
    return !sgn(det(a-b, s-t));
}
Pt triangleMassCenter(Pt a, Pt b, Pt c) 
{
    return (a+b+c) / 3.0;
}
double polygon_area(Pt poly[],int n)
{
    double ans=0.0;
    if(n<3)
        return ans;
    for(int i=1;i<=n;i++)
        ans+=det(poly[i],poly[i%n+1]);
    return ans*0.5;
}

struct Circle
{
    Pt c;
    double r;
    Circle(Pt _c,double _r)
    {
        c=_c;
        r=_r;
    }
};
bool PointInCircle(Circle a,Pt p)
{
    Pt c=a.c;
    return sgn( (p-c).norm()-a.r )<0;
}

bool PointOnCircle(Circle a,Pt p)
{
    Pt c=a.c;
    return sgn( (p-c).norm()-a.r )==0;
}

bool PointIOCircle(Circle a,Pt p)
{
    Pt c=a.c;
    return sgn( (p-c).norm()-a.r )<=0;
}

void circle_cross_line(Circle c,Pt a, Pt b, Pt ans[], int &num)
{
    double x1=a.x,y1=a.y,x2=b.x,y2=b.y;
    double dx=x2-x1,dy=y2-y1;
    double tmpx=x1-c.c.x,tmpy=y1-c.c.y;
    double A=dx*dx+dy*dy;
    double B=2.0*( dx*tmpx+dy*tmpy );
    double C=tmpx*tmpx+tmpy*tmpy-c.r*c.r;
    double delta=B*B-4.0*A*C;
    num=0;
    if(sgn(delta)<0)
        return;
    double t1=(-B-Sqrt(delta))/(2.0*A);
    double t2=(-B+Sqrt(delta))/(2.0*A);
    if(sgn(t1-1.0)<=0 && sgn(t1)>=0)
        ans[++num]=Pt(x1+t1*dx,y1+t1*dy);
    if(sgn(delta)==0)
        return;
    if(sgn(t2-1.0)<=0 && sgn(t2)>=0)
        ans[++num]=Pt(x1+t2*dx,y1+t2*dy);
}   

double sector_area(Circle c,Pt a,Pt b)
{
    a=a-c.c,b=b-c.c;
    double theta = atan2(a.y, a.x) - atan2(b.y, b.x);
    while (sgn(theta) <= 0) theta += 2.0*Pi;
    while (sgn(theta-2.0*Pi)>0) theta -= 2.0*Pi;
    theta = min(theta, 2.0*Pi - theta);
    return c.r*c.r*theta*0.5;
}

double CirclePolyArea(Circle c,Pt poly[],int n)
{
    double ans=0.0;
    Pt p[3];
    int num;
    for(int i=1;i<=n;i++)
    {
        Pt a=poly[i],b=poly[i%n+1];
        int ina=PointInCircle(c,a);
        int inb=PointInCircle(c,b);
        Pt da=a-c.c,db=b-c.c;
        int s=sgn(det(da,db));
        double part=0.0;
        if(ina)
        {
            if(inb)
            {
                part=Abs(det(da,db))*0.5;
            }
            else
            {
                circle_cross_line(c,a,b,p,num);
                part=sector_area(c,p[1],b)+Abs(det(da,p[1]-c.c))*0.5;
            }
        }
        else
        {
            if(inb)
            {
                circle_cross_line(c,a,b,p,num);
                part=sector_area(c,p[1],a)+Abs(det(db,p[1]-c.c))*0.5;
            }
            else
            {
                circle_cross_line(c,a,b,p,num);
                if(num==2)
                {
                    part=sector_area(c,a,p[1])+sector_area(c,b,p[2])+Abs(det(p[1]-c.c,p[2]-c.c))*0.5;
                }
                else
                {
                    part=sector_area(c,a,b);
                }
            }
        }
        if(s!=0)
            ans+=1.0*s*part;
    }
    return ans;
}
int main()
{
   int t=0;
    while(scanf("%d%lf",&n,&k)==2)
    {
        t++;
        for(int i=1;i<=n;i++)
            scanf("%lf%lf",&pl[i].x,&pl[i].y);
        scanf("%lf%lf%lf%lf",&A.x,&A.y,&B.x,&B.y);
        double r;
        double E=(2.0*B.x-2.0*sqr(k)*A.x)/(1.0-sqr(k));  
        double F=(2.0*B.y-2.0*sqr(k)*A.y)/(1.0-sqr(k));  
        double G=(sqr(k*A.x)+sqr(k*A.y)-sqr(B.x)-sqr(B.y))/(1.0-sqr(k));  
        pr.x=E/2.0,pr.y=F/2.0;  
        r=Sqrt(G+sqr(E)/4.0+sqr(F)/4.0);  
        Circle C(pr,r);
        double ans=Abs( CirclePolyArea(C,pl,n) );
        printf("Case %d: ",t);
        printf("%.10f\n",ans+eps);
    }
    return 0;
}