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*tPoint 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;}