poj2986A Triangle and a Circle (圆与三角形求交)

题目链接：http://poj.org/problem?id=2986


这题卡了好久TT...之前写错了 斯特瓦尔特(stewart)公式......
公式如下 ：设已知△ABC及其底边上B、C两点间的一点D，则有　　AB^2·DC+AC^2·BD-AD^2·BC=BC·DC·BD。
然后根据圆与三角形相交的4情况来求有向面积.

代码如下：

//============================================================================
// Name        : poj2986.cpp
// Author      : ssslpk
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//============================================================================

#include<iostream>
#include<cstdio>
#include<math.h>
//using namespace std;

#define eps 1e-8
#define PI acos(-1.0)

struct point {double x,y;};

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 tirarea(double a,double b,double c)
{
	double l=(a+b+c)/2;
	return sqrt(l*(l-a)*(l-b)*(l-c));
}

double distance(point p1,point p2)
{
	return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
}
double root(double a,double b,double c)//
{
	if(a==0)return -c/b;
	return (-b+sqrt(b*b-4*a*c))/(2*a);
}
double slove(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()
{
	point p[4],cen;
	double r;
	while(scanf("%lf%lf %lf%lf %lf%lf %lf%lf%lf"
		,&p[0].x,&p[0].y,&p[1].x,&p[1].y,&p[2].x,&p[2].y,&cen.x,&cen.y,&r)!=EOF)
	{
		
		int i;
		
		if(xmult(p[0],p[1],p[2])==0)
		{
			printf("0.00\n");continue;
		}
		p[3]=p[0];
		double res=0;
		for(i=0;i<3;i++)
			res+=slove(p[i],p[i+1],cen,r);
		printf("%.2lf\n",fabs(res));
	}
	return 0;
}