poj1329 Circle Through Three Points

地址: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 
	(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

 
思路:
  维基有公式。
  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\n",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\n\n", cgn(pp.x),fabs(pp.x*2),cgn(pp.y),fabs(pp.y*2),cgn(d), fabs(d));
399 
400     }
401     return 0;
402 }

 

posted @ 2017-10-11 20:37  weeping  阅读(216)  评论(0编辑  收藏  举报