POJ 4720 Naive and Silly Muggles -

题目地址：http://acm.hdu.edu.cn/showproblem.php?pid=4720

求围住一个三角形的最小圆

如果是钝角，那么圆的直径就是最长的那条边
其他情况圆是三角形的外接圆

可用以下方法求的外心：

给定三角形三个顶点的坐标，如何求三角形的外心的坐标呢？

例如 ：给定a（x1,y1） b（x2,y2） c（x3,y3）求外接圆心坐标O（x，y）

1. 首先，外接圆的圆心是三角形三条边的垂直平分线的交点，我们根据圆心到顶点的距离相等，可以列出以下方程：

(x1-x)*(x1-x)-(y1-y)*(y1-y)=(x2-x)*(x2-x)+(y2-y)*(y2-y);

(x2-x)*(x2-x)+(y2-y)*(y2-y)=(x3-x)*(x3-x)+(y3-y)*(y3-y);

2.化简得到：

2*(x2-x1)*x+2*(y2-y1)y=x2^2+y2^2-x1^2-y1^2;

2*(x3-x2)*x+2*(y3-y2)y=x3^2+y3^2-x2^2-y2^2;

令A1=2*(x2-x1)；

B1=2*(y2-y1)；

C1=x2^2+y2^2-x1^2-y1^2;

A2=2*(x3-x2)；

B2=2*(y3-y2)；

C2=x3^2+y3^2-x2^2-y2^2;

即

A1*x+B1y=C1;

A2*x+B2y=C2;

3.最后根据克拉默法则：

x=((C1*B2)-(C2*B1))/((A1*B2)-(A2*B1))；

y=((A1*C2)-(A2*C1))/((A1*B2)-(A2*B1))；

因此，x，y为最终结果；

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
#define REP(i,a,b) for(int i=a;i<=(int)(b);++i)
#define REPD(i,a,b) for(int i=a;i>=(int)(b);--i)
/*Point的模板*/
struct Point{
double x,y;
Point(double x=0,double y=0):x(x),y(y){}
};
typedef Point Vector;
//向量+向量=向量，点+向量=点
Vector operator + (Vector A, Vector B) { return Vector(A.x+B.x, A.y+B.y); }
//点-点=向量
Vector operator - (Vector A, Vector B) { return Vector(A.x-B.x, A.y-B.y); }
//向量*数=向量
Vector operator * (Vector A, double p) { return Vector(A.x*p, A.y*p); }
//向量/数=向量
Vector operator / (Vector A, double p) { return Vector(A.x/p, A.y/p); }

bool operator < (const Point& a, const Point& b) {
return a.x < b.x || (a.x == b.x && a.y < b.y);
}

//比较
const double eps = 1e-10;
int dcmp(double x){
if(fabs(x) < eps) return 0;
else return x < 0 ? -1: 1;
}

bool operator == (const Point& a,const Point& b) {
return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}

//基本计算
double Dot(Vector A, Vector B) { return A.x*B.x + A.y*B.y; }
double Length(Vector A) { return sqrt(Dot(A,A)); }
double Angle(Vector A, Vector B) { return acos(Dot(A, B)/Length(A)/Length(B)); } //两向量的夹角
double Cross(Vector A, Vector B) { return A.x*B.y - A.y*B.x; } //叉积
double Area2(Point A, Point B, Point C) { return Cross(B-A,C-A); } //有向面积

//A向量逆时针旋转α rad
//x'=xcosα-ysinα;
//y'=xsinα+ycosα;
Vector Rotate(Vector A, double rad){
return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}

//A的单位法线,也就是逆时针90°，长度变为1，注意A要非零向量
Vector Normal(Vector A){
double L=Length(A);
return Vector(-A.y/L,A.x/L);
}
/*Point的模板*/

/*Circle的模板*/
const double PI=acos(-1);
struct Circle
{
Point c;
double r;
Circle(Point c, double r):c(c), r(r) {}
Point point(double a) { //根据极角给出圆上的点
return Point(c.x + cos(a)*r, c.y + sin(a)*r);
}
};
struct Line
{
Vector v; //参数方程L=p+v*t
Point p;
Line(Point p,Vector v):v(v), p(p) {}
Point point(double t){ //根据参数方程L=v+(p-v)*t中t求出线上的点
return p+v*t;
}
};
//求直线L与圆C的交点，两个交点参数方程值t1,t2，sol保存交点本身，函数返回交点个数
int getLCI(Line L, Circle C, double& t1, double& t2, vector<Point>& sol){
double a = L.v.x, b = L.p.x-C.c.x, c = L.v.y, d = L.p.y-C.c.y;
double e = a*a+c*c, f = 2*(a*b+c*d), g = b*b+d*d-C.r*C.r;
double delta = f*f-4*e*g;
if(dcmp(delta) < 0) return 0;
if(dcmp(delta) == 0) {
t1 = t2 = -f / (2*e) ; sol.push_back(L.point(t1));
return 1;
}

t1 = (-f-sqrt(delta)) / (2*e); sol.push_back(L.point(t1));
t2 = (-f+sqrt(delta)) / (2*e); sol.push_back(L.point(t2));
return 2;
}
double Angle(Vector v) { return atan2(v.y,v.x); } //求v的极角
//求圆与圆的交点，函数返回交点数，sol保存交点
int getCCI(Circle C1, Circle C2, vector<Point>& sol){
double d = Length(C1.c-C2.c);
if(dcmp(d)==0){
if(dcmp(C1.r-C2.r)==0) return -1; //两圆重合
return 0;
}
if(dcmp(C1.r+C2.r-d)<0) return 0; //相离
if(dcmp(fabs(C1.r-C2.r)-d)>0) return 0; //内含

double a=Angle(C2.c-C1.c); //求极角
double da=acos((C1.r*C1.r+d*d-C2.r*C2.r) / (2*C1.r*d));
Point p1=C1.point(a-da), p2=C1.point(a+da);

sol.push_back(p1);
if(p1==p2) return 1;
sol.push_back(p2);
return 2;
}
//过点p到圆C做切线，函数返回切线条数，v[i]保存第i条切线的向量
int getT(Point p, Circle C, Vector* v){
Vector u=C.c-p;
double dist=Length(u);
if(dist<C.r) return 0;
else if(dcmp(dist-C.r)==0) {
v[0]=Rotate(u,PI/2);
return 1;
}
else {
double ang=asin(C.r/dist);
v[0]=Rotate(u, -ang);
v[1]=Rotate(u, +ang);
return 2;
}
}
//求两圆的公切线，返回切线条数，-1表示无穷条，a[i]和b[i]表示第i切线在圆A和圆B上的切点
int getT(Circle A, Circle B, Point* a, Point* b){
int cnt=0;
if(A.r<B.r) { swap(A, B); swap(a, b); } //保证圆A比圆B大
int d2=(A.c.x-B.c.x)*(A.c.x-B.c.x)+(A.c.y-B.c.y)*(A.c.y-B.c.y);
int rdiff=A.r-B.r;
int rsum=A.r+B.r;
if(d2<rdiff*rdiff) return 0; //内含

double base=Angle(B.c-A.c); // atan2(B.c.y-A.c.y,B.c.x-A.c.x);

if(d2==0&&A.r==B.r) return -1; //两圆重合，有无数条
if(d2==rdiff*rdiff){ //内切，一条切线
a[cnt]=A.point(base); b[cnt]=B.point(base); cnt++;
return 1;
}
//有外共切线
double ang=acos((A.r-B.r)/sqrt(d2));
a[cnt]=A.point(base+ang); b[cnt]=B.point(base+ang); cnt++;
a[cnt]=A.point(base-ang); b[cnt]=B.point(base-ang); cnt++;
if(d2==rsum*rsum){ //存在一条
a[cnt]=A.point(base); b[cnt]=B.point(PI+base); cnt++;
}
else if(d2>rsum*rsum){ //存在两条
double ang=acos((A.r+B.r)/sqrt(d2));
a[cnt]=A.point(base+ang); b[cnt]=B.point(PI+base+ang); cnt++;
a[cnt]=A.point(base-ang); b[cnt]=B.point(PI+base-ang); cnt++;
}
return cnt;
}
/*Circle的模板*/
inline Point ReadPoint(){
double a,b;
scanf("%lf%lf",&a,&b);
return Point(a,b);
}

bool isDun(Point p1, Point p2, Point p3){
double ang1=Angle(p2-p1,p3-p1),ang2=Angle(p3-p2,p1-p2),ang3=Angle(p1-p3,p2-p3);
// printf("%lf %lf %lf %lf\n", ang1,ang2,ang3,PI/2);
if(dcmp(ang1-PI/2)>0||dcmp(ang2-PI/2)>0||dcmp(ang3-PI/2)>0) return true;
else return false;

}
#define max3(a,b,c) max(a,max(b,c))
Circle getCircle(Point p1, Point p2, Point p3){
Point c0; double r;
double x1=p1.x,x2=p2.x,x3=p3.x,y1=p1.y,y2=p2.y,y3=p3.y;
if(isDun(p1,p2,p3)){
double a=Length(p2-p3),b=Length(p3-p1),c=Length(p2-p1);
if(dcmp(a-max3(a,b,c))==0) c0=Point((x2+x3)/2,(y2+y3)/2),r=a/2;
else if(dcmp(b-max3(a,b,c))==0) c0=Point((x3+x1)/2.0,(y3+y1)/2),r=b/2;
else c0=Point((x2+x1)/2,(y2+y1)/2),r=c/2;
}
else {
double A1=2*(x2-x1),B1=2*(y2-y1),C1=x2*x2+y2*y2-x1*x1-y1*y1;
double A2=2*(x3-x2),B2=2*(y3-y2),C2=x3*x3+y3*y3-x2*x2-y2*y2;
double x=((C1*B2)-(C2*B1)) / ((A1*B2)-(A2*B1));
double y=((A1*C2)-(A2*C1)) / ((A1*B2)-(A2*B1));
c0=Point(x,y); r=Length(c0-p1);
}

return Circle(c0,r);
}

int main(int argc, char const *argv[])
{
// freopen("input.in","r",stdin);
int T,kase=0; scanf("%d",&T);
while(T--){
Point p1,p2,p3,p;
p1=ReadPoint();
p2=ReadPoint();
p3=ReadPoint();
p=ReadPoint();
Circle c=getCircle(p1,p2,p3);
// cout<<c.c.x<<' '<<c.c.y<<' '<<c.r<<' '<<Length(c.c-p)<<endl;
printf("Case #%d: ", ++kase);
if(dcmp(Length(c.c-p)-c.r)<=0) printf("Danger\n");
else printf("Safe\n");
}
return 0;
}