地址:http://poj.org/problem?id=1329
题目:
Circle Through Three Points
Time Limit: 1000MS | Memory Limit: 10000K | |
Total Submissions: 3970 | Accepted: 1667 |
Description
Your team is to write a program that, given the Cartesian coordinates of three points on a plane, will find the equation of the circle through them all. The three points will not be on a straight line.
The solution is to be printed as an equation of the form
and an equation of the form
The solution is to be printed as an equation of the form
(x - h)^2 + (y - k)^2 = r^2 (1)
and an equation of the form
x^2 + y^2 + cx + dy - e = 0 (2)
Input
Each line of input to your program will contain the x and y coordinates of three points, in the order Ax, Ay, Bx, By, Cx, Cy. These coordinates will be real numbers separated from each other by one or more spaces.
Output
Your program must print the required equations on two lines using the format given in the sample below. Your computed values for h, k, r, c, d, and e in Equations 1 and 2 above are to be printed with three digits after the decimal point. Plus and minus signs in the equations should be changed as needed to avoid multiple signs before a number. Plus, minus, and equal signs must be separated from the adjacent characters by a single space on each side. No other spaces are to appear in the equations. Print a single blank line after each equation pair.
Sample Input
7.0 -5.0 -1.0 1.0 0.0 -6.0 1.0 7.0 8.0 6.0 7.0 -2.0
Sample Output
(x - 3.000)^2 + (y + 2.000)^2 = 5.000^2 x^2 + y^2 - 6.000x + 4.000y - 12.000 = 0 (x - 3.921)^2 + (y - 2.447)^2 = 5.409^2 x^2 + y^2 - 7.842x - 4.895y - 7.895 = 0
Source
Southern California 1989,UVA 190
思路:
维基有公式。
1 #include <iostream> 2 #include <cstdio> 3 #include <cmath> 4 #include <algorithm> 5 6 7 using namespace std; 8 const double PI = acos(-1.0); 9 const double eps = 1e-10; 10 11 /****************常用函数***************/ 12 //判断ta与tb的大小关系 13 int sgn( double ta, double tb) 14 { 15 if(fabs(ta-tb)<eps)return 0; 16 if(ta<tb) return -1; 17 return 1; 18 } 19 20 //点 21 class Point 22 { 23 public: 24 25 double x, y; 26 27 Point(){} 28 Point( double tx, double ty){ x = tx, y = ty;} 29 30 bool operator < (const Point &_se) const 31 { 32 return x<_se.x || (x==_se.x && y<_se.y); 33 } 34 friend Point operator + (const Point &_st,const Point &_se) 35 { 36 return Point(_st.x + _se.x, _st.y + _se.y); 37 } 38 friend Point operator - (const Point &_st,const Point &_se) 39 { 40 return Point(_st.x - _se.x, _st.y - _se.y); 41 } 42 //点位置相同(double类型) 43 bool operator == (const Point &_off)const 44 { 45 return sgn(x, _off.x) == 0 && sgn(y, _off.y) == 0; 46 } 47 48 }; 49 50 /****************常用函数***************/ 51 //点乘 52 double dot(const Point &po,const Point &ps,const Point &pe) 53 { 54 return (ps.x - po.x) * (pe.x - po.x) + (ps.y - po.y) * (pe.y - po.y); 55 } 56 //叉乘 57 double xmult(const Point &po,const Point &ps,const Point &pe) 58 { 59 return (ps.x - po.x) * (pe.y - po.y) - (pe.x - po.x) * (ps.y - po.y); 60 } 61 //两点间距离的平方 62 double getdis2(const Point &st,const Point &se) 63 { 64 return (st.x - se.x) * (st.x - se.x) + (st.y - se.y) * (st.y - se.y); 65 } 66 //两点间距离 67 double getdis(const Point &st,const Point &se) 68 { 69 return sqrt((st.x - se.x) * (st.x - se.x) + (st.y - se.y) * (st.y - se.y)); 70 } 71 72 //两点表示的向量 73 class Line 74 { 75 public: 76 77 Point s, e;//两点表示,起点[s],终点[e] 78 double a, b, c;//一般式,ax+by+c=0 79 double angle;//向量的角度,[-pi,pi] 80 81 Line(){} 82 Line( Point ts, Point te):s(ts),e(te){}//get_angle();} 83 Line(double _a,double _b,double _c):a(_a),b(_b),c(_c){} 84 85 //排序用 86 bool operator < (const Line &ta)const 87 { 88 return angle<ta.angle; 89 } 90 //向量与向量的叉乘 91 friend double operator / ( const Line &_st, const Line &_se) 92 { 93 return (_st.e.x - _st.s.x) * (_se.e.y - _se.s.y) - (_st.e.y - _st.s.y) * (_se.e.x - _se.s.x); 94 } 95 //向量间的点乘 96 friend double operator *( const Line &_st, const Line &_se) 97 { 98 return (_st.e.x - _st.s.x) * (_se.e.x - _se.s.x) - (_st.e.y - _st.s.y) * (_se.e.y - _se.s.y); 99 } 100 //从两点表示转换为一般表示 101 //a=y2-y1,b=x1-x2,c=x2*y1-x1*y2 102 bool pton() 103 { 104 a = e.y - s.y; 105 b = s.x - e.x; 106 c = e.x * s.y - e.y * s.x; 107 return true; 108 } 109 //半平面交用 110 //点在向量左边(右边的小于号改成大于号即可,在对应直线上则加上=号) 111 friend bool operator < (const Point &_Off, const Line &_Ori) 112 { 113 return (_Ori.e.y - _Ori.s.y) * (_Off.x - _Ori.s.x) 114 < (_Off.y - _Ori.s.y) * (_Ori.e.x - _Ori.s.x); 115 } 116 //求直线或向量的角度 117 double get_angle( bool isVector = true) 118 { 119 angle = atan2( e.y - s.y, e.x - s.x); 120 if(!isVector && angle < 0) 121 angle += PI; 122 return angle; 123 } 124 125 //点在线段或直线上 1:点在直线上 2点在s,e所在矩形内 126 bool has(const Point &_Off, bool isSegment = false) const 127 { 128 bool ff = sgn( xmult( s, e, _Off), 0) == 0; 129 if( !isSegment) return ff; 130 return ff 131 && sgn(_Off.x - min(s.x, e.x), 0) >= 0 && sgn(_Off.x - max(s.x, e.x), 0) <= 0 132 && sgn(_Off.y - min(s.y, e.y), 0) >= 0 && sgn(_Off.y - max(s.y, e.y), 0) <= 0; 133 } 134 135 //点到直线/线段的距离 136 double dis(const Point &_Off, bool isSegment = false) 137 { 138 ///化为一般式 139 pton(); 140 //到直线垂足的距离 141 double td = (a * _Off.x + b * _Off.y + c) / sqrt(a * a + b * b); 142 //如果是线段判断垂足 143 if(isSegment) 144 { 145 double xp = (b * b * _Off.x - a * b * _Off.y - a * c) / ( a * a + b * b); 146 double yp = (-a * b * _Off.x + a * a * _Off.y - b * c) / (a * a + b * b); 147 double xb = max(s.x, e.x); 148 double yb = max(s.y, e.y); 149 double xs = s.x + e.x - xb; 150 double ys = s.y + e.y - yb; 151 if(xp > xb + eps || xp < xs - eps || yp > yb + eps || yp < ys - eps) 152 td = min( getdis(_Off,s), getdis(_Off,e)); 153 } 154 return fabs(td); 155 } 156 157 //关于直线对称的点 158 Point mirror(const Point &_Off) 159 { 160 ///注意先转为一般式 161 Point ret; 162 double d = a * a + b * b; 163 ret.x = (b * b * _Off.x - a * a * _Off.x - 2 * a * b * _Off.y - 2 * a * c) / d; 164 ret.y = (a * a * _Off.y - b * b * _Off.y - 2 * a * b * _Off.x - 2 * b * c) / d; 165 return ret; 166 } 167 //计算两点的中垂线 168 static Line ppline(const Point &_a,const Point &_b) 169 { 170 Line ret; 171 ret.s.x = (_a.x + _b.x) / 2; 172 ret.s.y = (_a.y + _b.y) / 2; 173 //一般式 174 ret.a = _b.x - _a.x; 175 ret.b = _b.y - _a.y; 176 ret.c = (_a.y - _b.y) * ret.s.y + (_a.x - _b.x) * ret.s.x; 177 //两点式 178 if(fabs(ret.a) > eps) 179 { 180 ret.e.y = 0.0; 181 ret.e.x = - ret.c / ret.a; 182 if(ret.e == ret. s) 183 { 184 ret.e.y = 1e10; 185 ret.e.x = - (ret.c - ret.b * ret.e.y) / ret.a; 186 } 187 } 188 else 189 { 190 ret.e.x = 0.0; 191 ret.e.y = - ret.c / ret.b; 192 if(ret.e == ret. s) 193 { 194 ret.e.x = 1e10; 195 ret.e.y = - (ret.c - ret.a * ret.e.x) / ret.b; 196 } 197 } 198 return ret; 199 } 200 201 //------------直线和直线(向量)------------- 202 //向量向左边平移t的距离 203 Line& moveLine( double t) 204 { 205 Point of; 206 of = Point( -( e.y - s.y), e.x - s.x); 207 double dis = sqrt( of.x * of.x + of.y * of.y); 208 of.x= of.x * t / dis, of.y = of.y * t / dis; 209 s = s + of, e = e + of; 210 return *this; 211 } 212 //直线重合 213 static bool equal(const Line &_st,const Line &_se) 214 { 215 return _st.has( _se.e) && _se.has( _st.s); 216 } 217 //直线平行 218 static bool parallel(const Line &_st,const Line &_se) 219 { 220 return sgn( _st / _se, 0) == 0; 221 } 222 //两直线(线段)交点 223 //返回-1代表平行,0代表重合,1代表相交 224 static bool crossLPt(const Line &_st,const Line &_se, Point &ret) 225 { 226 if(parallel(_st,_se)) 227 { 228 if(Line::equal(_st,_se)) return 0; 229 return -1; 230 } 231 ret = _st.s; 232 double t = ( Line(_st.s,_se.s) / _se) / ( _st / _se); 233 ret.x += (_st.e.x - _st.s.x) * t; 234 ret.y += (_st.e.y - _st.s.y) * t; 235 return 1; 236 } 237 //------------线段和直线(向量)---------- 238 //直线和线段相交 239 //参数:直线[_st],线段[_se] 240 friend bool crossSL( Line &_st, Line &_se) 241 { 242 return sgn( xmult( _st.s, _se.s, _st.e) * xmult( _st.s, _st.e, _se.e), 0) >= 0; 243 } 244 245 //判断线段是否相交(注意添加eps) 246 static bool isCrossSS( const Line &_st, const Line &_se) 247 { 248 //1.快速排斥试验判断以两条线段为对角线的两个矩形是否相交 249 //2.跨立试验(等于0时端点重合) 250 return 251 max(_st.s.x, _st.e.x) >= min(_se.s.x, _se.e.x) && 252 max(_se.s.x, _se.e.x) >= min(_st.s.x, _st.e.x) && 253 max(_st.s.y, _st.e.y) >= min(_se.s.y, _se.e.y) && 254 max(_se.s.y, _se.e.y) >= min(_st.s.y, _st.e.y) && 255 sgn( xmult( _se.s, _st.s, _se.e) * xmult( _se.s, _se.e, _st.s), 0) >= 0 && 256 sgn( xmult( _st.s, _se.s, _st.e) * xmult( _st.s, _st.e, _se.s), 0) >= 0; 257 } 258 }; 259 260 //寻找凸包的graham 扫描法所需的排序函数 261 Point gsort; 262 bool gcmp( const Point &ta, const Point &tb)/// 选取与最后一条确定边夹角最小的点,即余弦值最大者 263 { 264 double tmp = xmult( gsort, ta, tb); 265 if( fabs( tmp) < eps) 266 return getdis( gsort, ta) < getdis( gsort, tb); 267 else if( tmp > 0) 268 return 1; 269 return 0; 270 } 271 272 273 class triangle 274 { 275 public: 276 Point a, b, c;//顶点 277 triangle(){} 278 triangle(Point a, Point b, Point c): a(a), b(b), c(c){} 279 280 //计算三角形面积 281 double area() 282 { 283 return fabs( xmult(a, b, c)) / 2.0; 284 } 285 286 //计算三角形外心 287 //返回:外接圆圆心 288 Point circumcenter() 289 { 290 double pa = a.x * a.x + a.y * a.y; 291 double pb = b.x * b.x + b.y * b.y; 292 double pc = c.x * c.x + c.y * c.y; 293 double ta = pa * ( b.y - c.y) - pb * ( a.y - c.y) + pc * ( a.y - b.y); 294 double tb = -pa * ( b.x - c.x) + pb * ( a.x - c.x) - pc * ( a.x - b.x); 295 double tc = a.x * ( b.y - c.y) - b.x * ( a.y - c.y) + c.x * ( a.y - b.y); 296 return Point( ta / 2.0 / tc, tb / 2.0 / tc); 297 } 298 299 //计算三角形内心 300 //返回:内接圆圆心 301 Point incenter() 302 { 303 Line u, v; 304 double m, n; 305 u.s = a; 306 m = atan2(b.y - a.y, b.x - a.x); 307 n = atan2(c.y - a.y, c.x - a.x); 308 u.e.x = u.s.x + cos((m + n) / 2); 309 u.e.y = u.s.y + sin((m + n) / 2); 310 v.s = b; 311 m = atan2(a.y - b.y, a.x - b.x); 312 n = atan2(c.y - b.y, c.x - b.x); 313 v.e.x = v.s.x + cos((m + n) / 2); 314 v.e.y = v.s.y + sin((m + n) / 2); 315 Point ret; 316 Line::crossLPt(u,v,ret); 317 return ret; 318 } 319 320 //计算三角形垂心 321 //返回:高的交点 322 Point perpencenter() 323 { 324 Line u,v; 325 u.s = c; 326 u.e.x = u.s.x - a.y + b.y; 327 u.e.y = u.s.y + a.x - b.x; 328 v.s = b; 329 v.e.x = v.s.x - a.y + c.y; 330 v.e.y = v.s.y + a.x - c.x; 331 Point ret; 332 Line::crossLPt(u,v,ret); 333 return ret; 334 } 335 336 //计算三角形重心 337 //返回:重心 338 //到三角形三顶点距离的平方和最小的点 339 //三角形内到三边距离之积最大的点 340 Point barycenter() 341 { 342 Line u,v; 343 u.s.x = (a.x + b.x) / 2; 344 u.s.y = (a.y + b.y) / 2; 345 u.e = c; 346 v.s.x = (a.x + c.x) / 2; 347 v.s.y = (a.y + c.y) / 2; 348 v.e = b; 349 Point ret; 350 Line::crossLPt(u,v,ret); 351 return ret; 352 } 353 354 //计算三角形费马点 355 //返回:到三角形三顶点距离之和最小的点 356 Point fermentPoint() 357 { 358 Point u, v; 359 double step = fabs(a.x) + fabs(a.y) + fabs(b.x) + fabs(b.y) + fabs(c.x) + fabs(c.y); 360 int i, j, k; 361 u.x = (a.x + b.x + c.x) / 3; 362 u.y = (a.y + b.y + c.y) / 3; 363 while (step > eps) 364 { 365 for (k = 0; k < 10; step /= 2, k ++) 366 { 367 for (i = -1; i <= 1; i ++) 368 { 369 for (j =- 1; j <= 1; j ++) 370 { 371 v.x = u.x + step * i; 372 v.y = u.y + step * j; 373 if (getdis(u,a) + getdis(u,b) + getdis(u,c) > getdis(v,a) + getdis(v,b) + getdis(v,c)) 374 u = v; 375 } 376 } 377 } 378 } 379 return u; 380 } 381 }; 382 383 triangle tr; 384 385 char cgn(double x) 386 { 387 return x>0?'+':'-'; 388 } 389 int main(void) 390 { 391 392 while(~scanf("%lf%lf%lf%lf%lf%lf",&tr.a.x,&tr.a.y,&tr.b.x,&tr.b.y,&tr.c.x,&tr.c.y)) 393 { 394 Point pp = tr.circumcenter(); 395 double r = getdis( tr.a, pp),d = pp.x*pp.x+pp.y*pp.y-r*r; 396 pp.x = -pp.x, pp.y = -pp.y; 397 printf("(x %c %.3f)^2 + (y %c %.3f)^2 = %.3f^2 ",cgn(pp.x),fabs(pp.x),cgn(pp.y), fabs(pp.y),r); 398 printf("x^2 + y^2 %c %.3fx %c %.3fy %c %.3f = 0 ", cgn(pp.x),fabs(pp.x*2),cgn(pp.y),fabs(pp.y*2),cgn(d), fabs(d)); 399 400 } 401 return 0; 402 }