昨晚到现在，终于A掉，上午在J2EE课上想清楚了，不过下午还是拍着数据才调出来BUG的 = =。。。这水平。。。到区域赛遇到计算几何神题肿么办吧。
这题一看就有思路啊，很裸的思路，先求半平面交（交得的面积是樱桃所在的那个大块块），然后求得的多边形和圆形求交。
1、半平面交的话，需要把线段方向改下，都变成有效区域为樱桃所在区域。
2、有一个坑就是如果切痕都切不到蛋糕，是需要输出100%的，直接特判下。
3、半平面交我用的N*LOGN的算法，这个就不多说了。但是记得需要初始化区域，我的区域是(-2*r, -2*r) --- (2*r, 2*r)，这个矩形区域。
关键是求多边形和圆形的交。经过半平面交，交的多边形一定是凸多边形，那么就计算凸多边形和圆形的交。
先求出凸多边形和圆的交点（包括在凸多边形在圆内的顶点），根据樱桃位置极角排序。
然后找到第一个顶点为起始顶点，找相邻顶点和第一个顶点组成的三角形，这个三角形一定是内部的面积。
关键是找弓形面积，如下图，蓝色线表示的是多边形和圆形围城的区域，最终和圆的交点应该是B D E。
当计算到DE这条线段时，做DE的中垂线，交圆与BC（绿色点），判断C或者B是否在多边形内部，如果在的话，说明弓形DE是需要计算的面积。
这里有个计算下面面积和上面面积之分。因为我的计算默认的是逆时针，所以我的计算的有效面积一定是左侧面积，如果需要计算弓形面积，那么计算的一定是DE右侧的面积（这个需要好好想下），然后就需要判断DE是优弧还是劣弧，如果是圆心在DE的右侧，则是优弧，反之为劣弧。这个用叉积判断下就好。
注意计算角度不能用asin - -。。。asin的范围是[-π/2, π/2]。。。（哎，我是挫人。。。）
还有就是，如果最后只是两个交点，注意不要计算两次。。。T T。。可以把线段变换下（因为只有两个点的话，需要保证樱桃在线段的右侧）
#include <set>
#include <map>
#include <queue>
#include <stack>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <limits.h>
#include <string.h>
#include <string>
#include <algorithm>
#define MID(x,y) ( ( x + y ) >> 1 )
#define L(x) ( x << 1 )
#define R(x) ( x << 1 | 1 )
#define FOR(i,s,t) for(int i=s; i<t; i++)
#define BUG puts("here!!!")
using namespace std;
const int MAX = 2010;
struct point {
double x,y;
void P(double xx, double yy){ x = xx; y = yy;}
void get(){ scanf("%lf%lf",&x, &y);}
};
struct line { point a,b; double ang;};
point s[MAX*2], ch, p[MAX*2];
line l[MAX*2],deq[MAX*2];
const double eps = 1e-6;
const double pi = acos(-1.0);
bool dy(double x,double y) { return x > y + eps;} // x > y
bool xy(double x,double y) { return x < y - eps;} // x < y
bool dyd(double x,double y) { return x > y - eps;} // x >= y
bool xyd(double x,double y) { return x < y + eps;} // x <= y
bool dd(double x,double y) { return fabs( x - y ) < eps;} // x == y
double crossProduct(point a,point b,point c)//向量 ac 在 ab 的方向
{
return (c.x - a.x)*(b.y - a.y) - (b.x - a.x)*(c.y - a.y);
}
double disp2p(point a,point b) // a b 两点之间的距离
{
return sqrt( ( a.x - b.x ) * ( a.x - b.x ) + ( a.y - b.y ) * ( a.y - b.y ) );
}
bool parallel(line u,line v)
{
return dd( (u.a.x - u.b.x)*(v.a.y - v.b.y) - (v.a.x - v.b.x)*(u.a.y - u.b.y) , 0.0 );
}
point l2l_inst_p(line l1,line l2)
{
point ans = l1.a;
double t = ((l1.a.x - l2.a.x)*(l2.a.y - l2.b.y) - (l1.a.y - l2.a.y)*(l2.a.x - l2.b.x))/
((l1.a.x - l1.b.x)*(l2.a.y - l2.b.y) - (l1.a.y - l1.b.y)*(l2.a.x - l2.b.x));
ans.x += (l1.b.x - l1.a.x)*t;
ans.y += (l1.b.y - l1.a.y)*t;
return ans;
}
bool equal_ang(line a,line b) // 第一次unique的比较函数
{
return dd(a.ang,b.ang);
}
bool cmphp(line a,line b) // 排序的比较函数
{
if( dd(a.ang,b.ang) ) return xy(crossProduct(b.a,b.b,a.a),0.0);
return xy(a.ang,b.ang);
}
bool equal_p(point a,point b)//第二次unique的比较函数
{
return dd(a.x,b.x) && dd(a.y,b.y);
}
void makeline_hp(double x1,double y1,double x2,double y2,line &l)
{
l.a.x = x1; l.a.y = y1; l.b.x = x2; l.b.y = y2;
l.ang = atan2(y2 - y1,x2 - x1);
}
void inst_hp_nlogn(line *ln,int n,point *s,int &len)
{
len = 0;
sort(ln,ln+n,cmphp);
n = unique(ln,ln+n,equal_ang) - ln;
int bot = 0,top = 1;
deq[0] = ln[0]; deq[1] = ln[1];
for(int i=2; i<n; i++)
{
if( parallel(deq[top],deq[top-1]) || parallel(deq[bot],deq[bot+1]) )
return ;
while( bot < top && dy(crossProduct(ln[i].a,ln[i].b,
l2l_inst_p(deq[top],deq[top-1])),0.0) )
top--;
while( bot < top && dy(crossProduct(ln[i].a,ln[i].b,
l2l_inst_p(deq[bot],deq[bot+1])),0.0) )
bot++;
deq[++top] = ln[i];
}
while( bot < top && dy(crossProduct(deq[bot].a,deq[bot].b,
l2l_inst_p(deq[top],deq[top-1])),0.0) ) top--;
while( bot < top && dy(crossProduct(deq[top].a,deq[top].b,
l2l_inst_p(deq[bot],deq[bot+1])),0.0) ) bot++;
if( top <= bot + 1 ) return ;
for(int i=bot; i<top; i++)
s[len++] = l2l_inst_p(deq[i],deq[i+1]);
if( bot < top + 1 ) s[len++] = l2l_inst_p(deq[bot],deq[top]);
len = unique(s,s+len,equal_p) - s;
}
point l2l_inst_p(point u1,point u2,point v1,point v2)
{
point ans = 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));
ans.x += (u2.x - u1.x)*t;
ans.y += (u2.y - u1.y)*t;
return ans;
}
void l2c_inst_p(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 = l2l_inst_p(p,c,l1,l2);
t = sqrt(r*r - disp2p(p,c)*disp2p(p,c))/disp2p(l1,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 disp2l(point a,point l1,point l2)
{
return fabs( crossProduct(a,l1,l2) )/disp2p(l1,l2);
}
bool onSegment(point a, point b, point c)
{
if( dd(crossProduct(a,b,c),0.0) && dyd(c.x,min(a.x,b.x)) &&
xyd(c.x,max(a.x,b.x)) && dyd(c.y,min(a.y,b.y)) && xyd(c.y,max(a.y,b.y)) )
return true;
return false;
}
bool cmp(point a,point b)
{
double t1 = atan2(a.y - ch.y, a.x - ch.x);
double t2 = atan2(b.y - ch.y, b.x - ch.x);
return xy(t1, t2);
}
double area_triangle(point a,point b,point c)
{
return fabs( crossProduct(a,b,c) )/2.0;
}
double gong_area(point c,point a,point b,double r)
{
if( dd(crossProduct(c,a,b), 0.0) )
return r*r*pi/2;
double area = area_triangle(a, b, c);
double ang = acos(1-disp2p(a,b)*disp2p(a,b)/(r*r*2));
return ang*r*r/2 - area;
}
bool incircle(point a,point c, double r)
{
return xy(disp2p(a, c), r);
}
bool pin_convexh(point *p,int n,point a)
{
p[n] = p[0]; p[n+1] = p[1];
for(int i=0; i<n; i++)
if( xy(crossProduct(p[i],p[i+1],a)*
crossProduct(p[i+1],p[i+2],a),0.0) )
return false;
return true;
}
point foot_line(point a,point l1,point l2)
{
point c;
c.x = a.x - l2.y + l1.y;
c.y = a.y + l2.x - l1.x;
return c;
}
double solve(point *s, int len, double r)
{
point c; c.x = c.y = 0;
point a, b;
int cnt = 0;
s[len] = s[0];
FOR(i, 0, len) // 求得圆形和多边形的交的多边形的点
{
if( xy(disp2l(c, s[i], s[i+1]), r) )
{
l2c_inst_p(c, r, s[i], s[i+1], a, b);
if( onSegment(s[i], s[i+1], a) )
p[cnt++] = a;
if( onSegment(s[i], s[i+1], b) )
p[cnt++] = b;
}
if( incircle(s[i], c, r) )
p[cnt++] = s[i];
if( incircle(s[i+1], c, r) )
p[cnt++] = s[i+1];
}
sort(p, p+cnt, cmp);
cnt = unique(p, p+cnt, equal_p) - p;
p[cnt] = p[0];
double area = 0.0;
if( cnt == 2 ) // 如果有两个点的话 是不能计算两次的
{
if( xy(crossProduct(p[0],p[1],ch), 0.0) )
swap(p[0], p[1]);
cnt--;
}
FOR(i, 0, cnt)
{
area += area_triangle(p[0], p[i], p[i+1]);
if( !dd(disp2p(p[i], c), r) || !dd(disp2p(p[i+1], c), r) )
continue;
point mid;
mid.P((p[i].x + p[i+1].x)/2, (p[i].y + p[i+1].y)/2);
point a, b, k;
k = foot_line(c, p[i], p[i+1]);
l2c_inst_p(c, r, k, c, a, b);
bool flag = false;
// 判断两个交点是否在线段右边而且在多边形内
if( dyd(crossProduct(p[i], p[i+1], a), 0.0) && pin_convexh(s, len, a) )
flag = true;
if( dyd(crossProduct(p[i], p[i+1], b), 0.0) && pin_convexh(s, len, b) )
flag = true;
if( !flag ) continue;
double area_gong = gong_area(c, p[i], p[i+1], r);
if( dy(crossProduct(p[i], p[i+1], c), 0.0) )
area_gong = pi*r*r - area_gong;
area += area_gong;
}
return area/(pi*r*r);
}
int main()
{
int ncases, ind = 1, n;
double r;
scanf("%d", &ncases);
while( ncases-- )
{
scanf("%lf%d",&r, &n);
FOR(i, 0, n)
{
l[i].a.get();
l[i].b.get();
}
ch.get();
bool f = false;
point c; c.P(0,0);
FOR(i, 0, n)
{
if( dy(crossProduct(l[i].a, l[i].b, ch), 0.0) )
swap(l[i].a, l[i].b);
l[i].ang = atan2(l[i].b.y - l[i].a.y, l[i].b.x - l[i].a.x);
if( xy(disp2l(c,l[i].a, l[i].b), r) )
f = true;
}
if( !f ) // 如果输入的切痕都切不到蛋糕，要输出100%
{
printf("Case %d: %.5lf%%\n",ind++, 100.0);
continue;
}
makeline_hp(-2*r, -2*r, 2*r, -2*r, l[n++]);
makeline_hp(2*r, -2*r, 2*r, 2*r, l[n++]);
makeline_hp(2*r, 2*r, -2*r, 2*r, l[n++]);
makeline_hp(-2*r, 2*r, -2*r, -2*r, l[n++]);
int len;
inst_hp_nlogn(l, n, s, len);
double ans = solve(s, len, r);
printf("Case %d: %.5lf%%\n",ind++, ans*100);
}
return 0;
}