HDOJ 5130 Signal Interference（圆与多边形面积交）

通过化简，发现p点轨迹是圆，于是问题转化为圆与多边形面积交。

设k1 = k / (1.0 + k)，k2 = k / (1.0 - k)，则半径r = l * (k1 + k2) * 0.5（l为AB两点间距离）。

ox = 0.5 * ((ax + (1 - k1) * (bx - ax)) + (bx + k2 * (bx - ax)))，

oy = 0.5 * ((ay + (1 - k1) * (by - ay)) + (by + k2 * (by - ay)))。

模板大部分参考了这篇博文：计算几何--简单多边形与圆面积交 - lxglbk - 博客园

#include <algorithm>
#include <iostream>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
using namespace std;
#define LL long long
#define MAXN 510
#define eps 1e-8

int sig(double x)
{
    return (x > eps) - (x < -eps);
}

typedef struct Point
{
        double x, y;
        Point() {}
        Point(double _x, double _y):
                x(_x), y(_y) {}
        bool operator <(const Point &argu) const
        {
                if(sig(x - argu.x) == 0) return y < argu.y;
                return x < argu.x;
        }
        double dis(const Point &argu) const
        {
            return sqrt((x - argu.x) * (x - argu.x) + (y - argu.y) * (y - argu.y));
        }
        double dis2(const Point &argu) const
        {
            return (x - argu.x) * (x - argu.x) + (y - argu.y) * (y - argu.y);
        }
        double operator ^(const Point &argu) const
        {
            return x * argu.y - y * argu.x;
        }
        double operator *(const Point &argu) const
        {
            return x * argu.x + y * argu.y;
        }
        Point operator -(const Point &argu) const
        {
            return Point(x - argu.x, y - argu.y);
        }
        double len2() const
        {
            return x * x + y * y;
        }
        double len() const
        {
            return sqrt(x * x + y * y);
        }
        void in()
        {
            scanf("%lf%lf", &x, &y);
        }
        void out()
        {
            printf("%.10lf %.10lf\n", x, y);
        }
}Vector;

Point Get_Intersection(Point u1, Point u2, Point v1, Point v2)
{
    Point ret = u1;
    double t = ((u1.x - v1.x) * (v1.y - v2.y) - (u1.y-v1.y) * (v1.x - v2.x))
            / ((u1.x - u2.x) * (v1.y - v2.y) - (u1.y - u2.y) * (v1.x - v2.x));
    ret.x += (u2.x - u1.x) * t;
    ret.y += (u2.y - u1.y) * t;
    return ret;
}

void Intersection_Line_Circle(Point c, double r, Point l1, Point l2, Point& p1, Point& p2)
{
    Point p = c;
    double t;
    p.x += l1.y - l2.y;
    p.y += l2.x - l1.x;
    p = Get_Intersection(p, c, l1, l2);
    t = sqrt(r * r - p.dis2(c)) / l1.dis(l2);
    p1.x = p.x + (l2.x - l1.x) * t;
    p1.y = p.y + (l2.y - l1.y) * t;
    p2.x = p.x - (l2.x - l1.x) * t;
    p2.y = p.y - (l2.y - l1.y) * t;
}

double Dis_Line(Point a, Point b, Point p)
{
    double area = fabs((a - p) ^ (b - p));
    if(sig(area) == 0)
        return min(a.dis(p), b.dis(p));
    return area / a.dis(b);
}

Point FootPoint(Point p, Point l1, Point l2)
{
    Point t = p;
    t.x += l1.y - l2.y;
    t.y += l2.x - l1.x;
    if (((l1 - p) ^ (t - p)) * ((l2 - p) ^ (t - p)) > eps)
        return p.dis2(l1) < p.dis2(l2) ? l1 : l2;
    return Get_Intersection(p, t, l1, l2);
}

double Direct_Triangle_Circle_Area(Point a,Point b,Point o,double r)
{
    double sign = 1.0;
    a = a - o, b = b - o;
    o = Point(0.0, 0.0);

    if(sig(a ^ b) == 0)
        return 0.0;

    if(sig(a.len2() - b.len2()) > 0)
        swap(a, b), sign = -1.0;

    Point pa, pb;
    Intersection_Line_Circle(o, r, a, b, pa, pb);
    if(pb.dis2(b) > pa.dis2(b))
        swap(pa, pb);

    if(a.len2() < r * r + eps)
    {
        if(b.len2() < r * r + eps)
            return (a ^ b) * 0.5 * sign;
        double ang = acos(pb * b / pb.len() / b.len());
        double ret = (ang * r * r + fabs(a ^ pb)) * 0.5;
        if(sig((a ^ b) * sign) < 0)
            ret = -ret;
        return ret;
    }

    Point f = FootPoint(o, a, b);
    if(f.len2() > r * r)
    {
        double ang = acos(a * b / (a.len() * b.len()));
        double ret = ang * r * r * 0.5;
        if(sig((a ^ b) * sign) < 0)
            ret = -ret;
        return ret;

    }

    double cm = r / (a.len() - r);
    Point m = Point((o.x + cm * a.x) / (1 + cm), (o.y + cm * a.y) / (1 + cm));
    double cn = r / (b.len() - r);
    Point n = Point((o.x + cn * b.x) / (1 + cn), (o.y + cn * b.y) / (1 + cn));
    double ang1 = acos( m * n / m.len()/ n.len());
    double ang2 = acos(pa * pb / pa.len()/ pb.len());
    double ret = ((ang1 - ang2) * r * r + fabs(pa ^ pb)) *0.5;
    if(sig((a ^ b) * sign) < 0)
        ret = -ret;
    return ret;
}

int n;
double k;
Point al[MAXN], a, b, o;

double work()
{
    double ret = 0;

    for(int i = 0; i < n; i++)
        al[i].in();
    a.in(), b.in();
    double l = a.dis(b);

    double k1 = k / (1.0 + k);
    double k2 = k / (1.0 - k);
    double r = l * (k1 + k2) * 0.5;

    o.x = 0.5 * ((a.x + (1 - k1) * (b.x - a.x)) + (b.x + k2 * (b.x - a.x)));
    o.y = 0.5 * ((a.y + (1 - k1) * (b.y - a.y)) + (b.y + k2 * (b.y - a.y)));

    double ans = 0.0;
    for(int i = 0; i < n - 1; i++)
        ans += Direct_Triangle_Circle_Area(al[i], al[i + 1], o, r);
    ans += Direct_Triangle_Circle_Area(al[n - 1], al[0], o, r);
    return fabs(ans);
}

int main()
{
//    freopen("5130.in", "r", stdin);
//    freopen("5130.out", "w", stdout);

    int cas = 1;
    while(~scanf("%d%lf\n", &n, &k))
    {
        printf("Case %d: %.10f\n", cas++, work());
    }
    return 0;
}