hdu3982 Harry Potter and J.K.Rowling（半平面交 + 圆与多边形求交）

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=3982

题目意思：有一块半径为r的圆形蛋糕，其中心在原点，有一个人在某一个点（x，y），

现把蛋糕按两点所在直线切蛋糕，求切了n次后，这个人所在蛋糕的面积占总面积的百分比。

很容易想到先进行半平面交求出这个人所在位置的区域，再根据这个区域（多边形）求与圆的交，就是就个人得到的蛋糕。

代码如下：

//140MS
#include<iostream>
#include<cstdio>
#include<math.h>
#define eps 1e-8
#define maxn 10000
#define PI acos(-1.0)
struct point{double x,y;};
point p[maxn],save[maxn],temp[maxn];
point points[maxn][2];
point cen,people;
int n,ns,m;
double r;
double xmult(point p1,point p2,point p0)
{
	return (p1.x-p0.x)*(p2.y-p0.y)-(p1.y-p0.y)*(p2.x-p0.x);
}
double dmult(point p1,point p2,point p0)
{
	return (p1.x-p0.x)*(p2.x-p0.x)+(p1.y-p0.y)*(p2.y-p0.y);
}
double distance(point p1,point p2)
{
	return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
}
void mycopy(point *s,int &ns,point *temp,int n)
{
	int i;
	ns=n;
	for(i=0;i<n;i++)s[i]=temp[i];
}
void getline(point p1,point p2,double &a,double &b,double &c)
{
	a=p2.y-p1.y;
	b=p1.x-p2.x;
	c=p2.x*p1.y-p1.x*p2.y;
}
point intersection(point p1,point p2,point p3,point p4)
{
	point ret=p1;
	double t=((p1.x-p3.x)*(p3.y-p4.y)-(p1.y-p3.y)*(p3.x-p4.x))
		/((p1.x-p2.x)*(p3.y-p4.y)-(p1.y-p2.y)*(p3.x-p4.x));
	ret.x+=(p2.x-p1.x)*t;
	ret.y+=(p2.y-p1.y)*t;
	return ret;
}
void cut(point side)
{
	int i,j;
	temp[0].x=-10000;temp[0].y=-10000;
	temp[1].x=-10000;temp[1].y=10000;
	temp[2].x=10000;temp[2].y=10000;
	temp[3].x=10000;temp[3].y=-10000;
	ns=m=4;
	for(i=0;i<n;i++)
	{
		double a,b,c;
		if(xmult(side,points[i][1],points[i][0])>=0)
			getline(points[i][0],points[i][1],a,b,c);
		else  getline(points[i][1],points[i][0],a,b,c);
		int cnt=0;
		for(j=0;j<ns;j++)
		{
			if(a* temp[j].x+b* temp[j].y+c>=0)
			{
				save[cnt++]=temp[j];
			}
			else
			{
				point p1=temp[(j-1+ns)%ns],p2=temp[(j+1)%ns];
				if(a*p1.x+b*p1.y+c>0)
					save[cnt++]=intersection(points[i][0],points[i][1],p1,temp[j]);
				if(a*p2.x+b*p2.y+c>0)
					save[cnt++]=intersection(points[i][0],points[i][1],p2,temp[j]);
			}
		}
		mycopy(temp,ns,save,cnt);
	}
}
double cirtri(point pa,point pb,point po,double r)
{
	double a,b,c,x,y;
	double area=xmult(pa,pb,po)/2;
	a=distance(po,pb);
	b=distance(po,pa);
	c=distance(pa,pb);
	if(a<=r&&b<=r)//1
	{
		return area;
	}
	else if(a<r&&b>=r)//2
	{
		x=(dmult(pa,po,pb)+sqrt(c*c*r*r-xmult(pa,po,pb)*xmult(pa,po,pb)))/c;
		return asin(area*(c-x)*2/c/b/r)*r*r/2+area*x/c;
	}
	else if(a>=r&&b<r)//3
	{
		y=(dmult(pb,po,pa)+sqrt(c*c*r*r-xmult(pb,po,pa)*xmult(pb,po,pa)))/c;
		return asin(area*(c-y)*2/c/a/r)*r*r/2+area*y/c;
	}
	else if(fabs(2*area)>=r*c||dmult(pb,po,pa)<=0
		||dmult(pa,po,pb)<=0)//4
	{
		if(dmult(pa,pb,po)<0)
		{
			if(xmult(pa,pb,po)<0)
			{
				return (-PI-asin(area*2/a/b))*r*r/2;
			}
			else return (PI-asin(area*2/a/b))*r*r/2;
		}
		else return asin(area*2/a/b)*r*r/2;
	}
	else //5
	{
		x=(dmult(pa,po,pb)+sqrt(c*c*r*r-xmult(pa,po,pb)*xmult(pa,po,pb)))/c;
		y=(dmult(pb,po,pa)+sqrt(c*c*r*r-xmult(pb,po,pa)*xmult(pb,po,pa)))/c;
		return (asin(area*(1-x/c)*2/r/b)
			+asin(area*(1-y/c)*2/r/a))*r*r/2
			+area*((y+x)/c-1);
	}
}
int main()
{
	int cas;
	scanf("%d",&cas);
	int i,j,k;
	for(k=1;k<=cas;k++)
	{
		scanf("%lf%d",&r,&n);
		for(i=0;i<n;i++)
			scanf("%lf%lf%lf%lf",&points[i][0].x,&points[i][0].y,
			&points[i][1].x,&points[i][1].y);
		scanf("%lf%lf",&people.x,&people.y);
		cut(people);
		double res=0;
		cen.x=0;cen.y=0;
		for(i=0;i<ns;i++)
			res+=cirtri(save[i],save[(i+1)%ns],cen,r);
		double ans=fabs(100.0*res)/(PI*r*r);
		printf("Case %d: %.5lf%%\n",k,ans);
	}
	return 0;
}