计算几何模板--更新中

weeping的私人模板，需要的话就拿去用吧（虽然可能会有小错，如果有请提醒我）

这是现在日常用的模板

2017.10.11 update:

　　把白书上的半平面交模板换成别的了，原来的那个好像有问题，一直wa。

　　把求三角形外心的模板改了精度更高的模板。

2017.10.08 update:

　　原来的模板太啰嗦，今天全部重新简化了，但同时保留以前模板。

 2017.09.04 update:

　　1.发现射线法判断点在多边形内部的代码有误，已修正

2017.08.06 update:

　　1.卷包裹法求凸包 

　　2.三角形和圆的相关模板

　　3.O(logn)判断点在凸包内

现在模板：

  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 class Polygon
273 {
274 public:
275     const static int maxpn = 5e4+7;
276     Point pt[maxpn];//点（顺时针或逆时针）
277     Line dq[maxpn]; //求半平面交打开注释
278     int n;//点的个数
279 
280 
281     //求多边形面积，多边形内点必须顺时针或逆时针
282     double area()
283     {
284         double ans = 0.0;
285         for(int i = 0; i < n; i ++)
286         {
287             int nt = (i + 1) % n;
288             ans += pt[i].x * pt[nt].y - pt[nt].x * pt[i].y;
289         }
290         return fabs( ans / 2.0);
291     }
292     //求多边形重心，多边形内点必须顺时针或逆时针
293     Point gravity()
294     {
295         Point ans;
296         ans.x = ans.y = 0.0;
297         double area = 0.0;
298         for(int i = 0; i < n; i ++)
299         {
300             int nt = (i + 1) % n;
301             double tp = pt[i].x * pt[nt].y - pt[nt].x * pt[i].y;
302             area += tp;
303             ans.x += tp * (pt[i].x + pt[nt].x);
304             ans.y += tp * (pt[i].y + pt[nt].y);
305         }
306         ans.x /= 3 * area;
307         ans.y /= 3 * area;
308         return ans;
309     }
310     //判断点是否在任意多边形内[射线法]，O(n)
311     bool ahas( Point &_Off)
312     {
313         int ret = 0;
314         double infv = 1e20;//坐标系最大范围
315         Line l = Line( _Off, Point( -infv ,_Off.y));
316         for(int i = 0; i < n; i ++)
317         {
318             Line ln = Line( pt[i], pt[(i + 1) % n]);
319             if(fabs(ln.s.y - ln.e.y) > eps)
320             {
321                 Point tp = (ln.s.y > ln.e.y)? ln.s: ln.e;
322                 if( ( fabs( tp.y - _Off.y) < eps && tp.x < _Off.x + eps) || Line::isCrossSS( ln, l))
323                     ret++;
324             }
325             else if( Line::isCrossSS( ln, l))
326                 ret++;
327         }
328         return ret&1;
329     }
330 
331     //判断任意点是否在凸包内，O(logn)
332     bool bhas( Point & p)
333     {
334         if( n < 3)
335             return false;
336         if( xmult( pt[0], p, pt[1]) > eps)
337             return false;
338         if( xmult( pt[0], p, pt[n-1]) < -eps)
339             return false;
340         int l = 2,r = n-1;
341         int line = -1;
342         while( l <= r)
343         {
344             int mid = ( l + r) >> 1;
345             if( xmult( pt[0], p, pt[mid]) >= 0)
346                 line = mid,r = mid - 1;
347             else l = mid + 1;
348         }
349         return xmult( pt[line-1], p, pt[line]) <= eps;
350     }
351 
352 
353 
354     //凸多边形被直线分割
355     Polygon split( Line &_Off)
356     {
357         //注意确保多边形能被分割
358         Polygon ret;
359         Point spt[2];
360         double tp = 0.0, np;
361         bool flag = true;
362         int i, pn = 0, spn = 0;
363         for(i = 0; i < n; i ++)
364         {
365             if(flag)
366                 pt[pn ++] = pt[i];
367             else
368                 ret.pt[ret.n ++] = pt[i];
369             np = xmult( _Off.s, _Off.e, pt[(i + 1) % n]);
370             if(tp * np < -eps)
371             {
372                 flag = !flag;
373                 Line::crossLPt( _Off, Line(pt[i], pt[(i + 1) % n]), spt[spn++]);
374             }
375             tp = (fabs(np) > eps)?np: tp;
376         }
377         ret.pt[ret.n ++] = spt[0];
378         ret.pt[ret.n ++] = spt[1];
379         n = pn;
380         return ret;
381     }
382 
383 
384     /** 卷包裹法求点集凸包，_p为输入点集，_n为点的数量 **/
385     void ConvexClosure( Point _p[], int _n)
386     {
387         sort( _p, _p + _n);
388         n = 0;
389         for(int i = 0; i < _n; i++)
390         {
391             while( n > 1 && sgn( xmult( pt[n-2], pt[n-1], _p[i]), 0) <= 0)
392                 n--;
393             pt[n++] = _p[i];
394         }
395         int _key = n;
396         for(int i = _n - 2; i >= 0; i--)
397         {
398             while( n > _key && sgn( xmult( pt[n-2], pt[n-1], _p[i]), 0) <= 0)
399                 n--;
400             pt[n++] = _p[i];
401         }
402         if(n>1)   n--;//除去重复的点，该点已是凸包凸包起点
403     }
404     /****** 寻找凸包的graham 扫描法********************/
405     /****** _p为输入的点集,_n为点的数量****************/
406 
407     void graham( Point _p[], int _n)
408     {
409         int cur=0;
410         for(int i = 1; i < _n; i++)
411             if( sgn( _p[cur].y, _p[i].y) > 0 || ( sgn( _p[cur].y, _p[i].y) == 0 && sgn( _p[cur].x, _p[i].x) > 0) )
412                 cur = i;
413         swap( _p[cur], _p[0]);
414         n = 0, gsort = pt[n++] = _p[0];
415         if( _n <= 1)   return;
416         sort( _p + 1, _p+_n ,gcmp);
417         pt[n++] = _p[1];
418         for(int i = 2; i < _n; i++)
419         {
420             while(n>1 && sgn( xmult( pt[n-2], pt[n-1], _p[i]), 0) <= 0)// 当凸包退化成直线时需特别注意n
421                 n--;
422             pt[n++] = _p[i];
423         }
424     }
425     //凸包旋转卡壳(注意点必须顺时针或逆时针排列)
426     //返回值凸包直径的平方（最远两点距离的平方）
427     pair<Point,Point> rotating_calipers()
428     {
429         int i = 1 % n;
430         double ret = 0.0;
431         pt[n] = pt[0];
432         pair<Point,Point>ans=make_pair(pt[0],pt[0]);
433         for(int j = 0; j < n; j ++)
434         {
435             while( fabs( xmult( pt[i+1], pt[j], pt[j + 1])) > fabs( xmult( pt[i], pt[j], pt[j + 1])) + eps)
436                 i = (i + 1) % n;
437             //pt[i]和pt[j],pt[i + 1]和pt[j + 1]可能是对踵点
438             if(ret < getdis2(pt[i],pt[j]))  ret = getdis2(pt[i],pt[j]), ans = make_pair(pt[i],pt[j]);
439             if(ret < getdis2(pt[i+1],pt[j+1]))  ret = getdis(pt[i+1],pt[j+1]), ans = make_pair(pt[i+1],pt[j+1]);
440         }
441         return ans;
442     }
443 
444     //凸包旋转卡壳(注意点必须逆时针排列)
445     //返回值两凸包的最短距离
446     double rotating_calipers( Polygon &_Off)
447     {
448         int i = 0;
449         double ret = 1e10;//inf
450         pt[n] = pt[0];
451         _Off.pt[_Off.n] = _Off.pt[0];
452         //注意凸包必须逆时针排列且pt[0]是左下角点的位置
453         while( _Off.pt[i + 1].y > _Off.pt[i].y)
454             i = (i + 1) % _Off.n;
455         for(int j = 0; j < n; j ++)
456         {
457             double tp;
458             //逆时针时为 >,顺时针则相反
459             while((tp = xmult(_Off.pt[i + 1],pt[j], pt[j + 1]) - xmult(_Off.pt[i], pt[j], pt[j + 1])) > eps)
460                 i = (i + 1) % _Off.n;
461             //(pt[i],pt[i+1])和(_Off.pt[j],_Off.pt[j + 1])可能是最近线段
462             ret = min(ret, Line(pt[j], pt[j + 1]).dis(_Off.pt[i], true));
463             ret = min(ret, Line(_Off.pt[i], _Off.pt[i + 1]).dis(pt[j + 1], true));
464             if(tp > -eps)//如果不考虑TLE问题最好不要加这个判断
465             {
466                 ret = min(ret, Line(pt[j], pt[j + 1]).dis(_Off.pt[i + 1], true));
467                 ret = min(ret, Line(_Off.pt[i], _Off.pt[i + 1]).dis(pt[j], true));
468             }
469         }
470         return ret;
471     }
472 
473     //-----------半平面交-------------
474     //复杂度:O(nlog2(n))
475     //获取半平面交的多边形（多边形的核）
476     //参数：向量集合[l]，向量数量[ln];(半平面方向在向量左边）
477     //函数运行后如果n[即返回多边形的点数量]为0则不存在半平面交的多边形（不存在区域或区域面积无穷大）
478     int judege( Line &_lx, Line &_ly, Line &_lz)
479     {
480         Point tmp;
481         Line::crossLPt(_lx,_ly,tmp);
482         return sgn(xmult(_lz.s,tmp,_lz.e),0);
483     }
484     int halfPanelCross(Line L[], int ln)
485     {
486         int i, tn, bot, top;
487         for(int i = 0; i < ln; i++)
488             L[i].get_angle();
489         sort(L, L + ln);
490         //平面在向量左边的筛选
491         for(i = tn = 1; i < ln; i ++)
492             if(fabs(L[i].angle - L[i - 1].angle) > eps)
493                 L[tn ++] = L[i];
494         ln = tn, n = 0, bot = 0, top = 1;
495         dq[0] = L[0], dq[1] = L[1];
496         for(i = 2; i < ln; i ++)
497         {
498             while(bot < top &&  judege(dq[top],dq[top-1],L[i]) > 0)
499                 top --;
500             while(bot < top &&  judege(dq[bot],dq[bot+1],L[i]) > 0)
501                 bot ++;
502             dq[++ top] = L[i];
503         }
504         while(bot < top && judege(dq[top],dq[top-1],dq[bot]) > 0)
505             top --;
506         while(bot < top && judege(dq[bot],dq[bot+1],dq[top]) > 0)
507             bot ++;
508         //若半平面交退化为点或线
509         //        if(top <= bot + 1)
510         //            return 0;
511         dq[++top] = dq[bot];
512         for(i = bot; i < top; i ++)
513             Line::crossLPt(dq[i],dq[i + 1],pt[n++]);
514         return n;
515     }
516 };
517 
518 
519 class Circle
520 {
521 public:
522     Point c;//圆心
523     double r;//半径
524     double db, de;//圆弧度数起点， 圆弧度数终点(逆时针0-360)
525 
526     //-------圆---------
527 
528     //判断圆在多边形内
529     bool inside( Polygon &_Off)
530     {
531         if(_Off.ahas(c) == false)
532             return false;
533         for(int i = 0; i < _Off.n; i ++)
534         {
535             Line l = Line(_Off.pt[i], _Off.pt[(i + 1) % _Off.n]);
536             if(l.dis(c, true) < r - eps)
537                 return false;
538         }
539         return true;
540     }
541 
542     //判断多边形在圆内（线段和折线类似）
543     bool has( Polygon &_Off)
544     {
545         for(int i = 0; i < _Off.n; i ++)
546             if( getdis2(_Off.pt[i],c) > r * r - eps)
547                 return false;
548         return true;
549     }
550 
551     //-------圆弧-------
552     //圆被其他圆截得的圆弧，参数：圆[_Off]
553     Circle operator-(Circle &_Off) const
554     {
555         //注意圆必须相交，圆心不能重合
556         double d2 = getdis2(c,_Off.c);
557         double d = getdis(c,_Off.c);
558         double ans = acos((d2 + r * r - _Off.r * _Off.r) / (2 * d * r));
559         Point py = _Off.c - c;
560         double oans = atan2(py.y, py.x);
561         Circle res;
562         res.c = c;
563         res.r = r;
564         res.db = oans + ans;
565         res.de = oans - ans + 2 * PI;
566         return res;
567     }
568     //圆被其他圆截得的圆弧，参数：圆[_Off]
569     Circle operator+(Circle &_Off) const
570     {
571         //注意圆必须相交，圆心不能重合
572         double d2 = getdis2(c,_Off.c);
573         double d = getdis(c,_Off.c);
574         double ans = acos((d2 + r * r - _Off.r * _Off.r) / (2 * d * r));
575         Point py = _Off.c - c;
576         double oans = atan2(py.y, py.x);
577         Circle res;
578         res.c = c;
579         res.r = r;
580         res.db = oans - ans;
581         res.de = oans + ans;
582         return res;
583     }
584 
585     //过圆外一点的两条切线
586     //参数：点[_Off](必须在圆外),返回：两条切线(切线的s点为_Off,e点为切点)
587     pair<Line, Line>  tangent( Point &_Off)
588     {
589         double d = getdis(c,_Off);
590         //计算角度偏移的方式
591         double angp = acos(r / d), ango = atan2(_Off.y - c.y, _Off.x - c.x);
592         Point pl = Point(c.x + r * cos(ango + angp), c.y + r * sin(ango + angp)),
593             pr = Point(c.x + r * cos(ango - angp), c.y + r * sin(ango - angp));
594         return make_pair(Line(_Off, pl), Line(_Off, pr));
595     }
596 
597     //计算直线和圆的两个交点
598     //参数：直线[_Off](两点式)，返回两个交点，注意直线必须和圆有两个交点
599     pair<Point, Point> cross(Line _Off)
600     {
601         _Off.pton();
602         //到直线垂足的距离
603         double td = fabs(_Off.a * c.x + _Off.b * c.y + _Off.c) / sqrt(_Off.a * _Off.a + _Off.b * _Off.b);
604 
605         //计算垂足坐标
606         double xp = (_Off.b * _Off.b * c.x - _Off.a * _Off.b * c.y - _Off.a * _Off.c) / ( _Off.a * _Off.a + _Off.b * _Off.b);
607         double yp = (- _Off.a * _Off.b * c.x + _Off.a * _Off.a * c.y - _Off.b * _Off.c) / (_Off.a * _Off.a + _Off.b * _Off.b);
608 
609         double ango = atan2(yp - c.y, xp - c.x);
610         double angp = acos(td / r);
611 
612         return make_pair(Point(c.x + r * cos(ango + angp), c.y + r * sin(ango + angp)),
613             Point(c.x + r * cos(ango - angp), c.y + r * sin(ango - angp)));
614     }
615 };
616 
617 class triangle
618 {
619 public:
620     Point a, b, c;//顶点
621     triangle(){}
622     triangle(Point a, Point b, Point c): a(a), b(b), c(c){}
623 
624     //计算三角形面积
625     double area()
626     {
627         return fabs( xmult(a, b, c)) / 2.0;
628     }
629 
630     //计算三角形外心
631     //返回：外接圆圆心
632     Point circumcenter()
633     {
634         double pa = a.x * a.x + a.y * a.y;
635         double pb = b.x * b.x + b.y * b.y;
636         double pc = c.x * c.x + c.y * c.y;
637         double ta = pa * ( b.y - c.y) - pb * ( a.y - c.y) + pc * ( a.y - b.y);
638         double tb = -pa * ( b.x - c.x) + pb * ( a.x - c.x) - pc * ( a.x - b.x);
639         double tc = a.x * ( b.y - c.y) - b.x * ( a.y - c.y) + c.x * ( a.y - b.y);
640         return Point( ta / 2.0 / tc, tb / 2.0 / tc);
641     }
642 
643     //计算三角形内心
644     //返回：内接圆圆心
645     Point incenter()
646     {
647         Line u, v;
648         double m, n;
649         u.s = a;
650         m = atan2(b.y - a.y, b.x - a.x);
651         n = atan2(c.y - a.y, c.x - a.x);
652         u.e.x = u.s.x + cos((m + n) / 2);
653         u.e.y = u.s.y + sin((m + n) / 2);
654         v.s = b;
655         m = atan2(a.y - b.y, a.x - b.x);
656         n = atan2(c.y - b.y, c.x - b.x);
657         v.e.x = v.s.x + cos((m + n) / 2);
658         v.e.y = v.s.y + sin((m + n) / 2);
659         Point ret;
660         Line::crossLPt(u,v,ret);
661         return ret;
662     }
663 
664     //计算三角形垂心
665     //返回：高的交点
666     Point perpencenter()
667     {
668         Line u,v;
669         u.s = c;
670         u.e.x = u.s.x - a.y + b.y;
671         u.e.y = u.s.y + a.x - b.x;
672         v.s = b;
673         v.e.x = v.s.x - a.y + c.y;
674         v.e.y = v.s.y + a.x - c.x;
675         Point ret;
676         Line::crossLPt(u,v,ret);
677         return ret;
678     }
679 
680     //计算三角形重心
681     //返回：重心
682     //到三角形三顶点距离的平方和最小的点
683     //三角形内到三边距离之积最大的点
684     Point barycenter()
685     {
686         Line u,v;
687         u.s.x = (a.x + b.x) / 2;
688         u.s.y = (a.y + b.y) / 2;
689         u.e = c;
690         v.s.x = (a.x + c.x) / 2;
691         v.s.y = (a.y + c.y) / 2;
692         v.e = b;
693         Point ret;
694         Line::crossLPt(u,v,ret);
695         return ret;
696     }
697 
698     //计算三角形费马点
699     //返回：到三角形三顶点距离之和最小的点
700     Point fermentPoint()
701     {
702         Point u, v;
703         double step = fabs(a.x) + fabs(a.y) + fabs(b.x) + fabs(b.y) + fabs(c.x) + fabs(c.y);
704         int i, j, k;
705         u.x = (a.x + b.x + c.x) / 3;
706         u.y = (a.y + b.y + c.y) / 3;
707         while (step > eps)
708         {
709             for (k = 0; k < 10; step /= 2, k ++)
710             {
711                 for (i = -1; i <= 1; i ++)
712                 {
713                     for (j =- 1; j <= 1; j ++)
714                     {
715                         v.x = u.x + step * i;
716                         v.y = u.y + step * j;
717                         if (getdis(u,a) + getdis(u,b) + getdis(u,c) > getdis(v,a) + getdis(v,b) + getdis(v,c))
718                             u = v;
719                     }
720                 }
721             }
722         }
723         return u;
724     }
725 };
726 
727 int main(void)
728 {
729 
730     return 0;
731 }

以前模板：

  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 class Point
 12 {
 13 public:
 14     double x, y;
 15 
 16     Point(){}
 17     Point(double x, double y):x(x),y(y){}
 18 
 19     bool operator < (const Point &_se) const
 20     {
 21         return x<_se.x || (x==_se.x && y<_se.y);
 22     }
 23     /*******判断ta与tb的大小关系*******/
 24     static int sgn(double ta,double tb)
 25     {
 26         if(fabs(ta-tb)<eps)return 0;
 27         if(ta<tb)   return -1;
 28         return 1;
 29     }
 30     static double xmult(const Point &po, const Point &ps, const Point &pe)
 31     {
 32         return (ps.x - po.x) * (pe.y - po.y) - (pe.x - po.x) * (ps.y - po.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  Point::sgn(x, _off.x) == 0 && Point::sgn(y, _off.y) == 0;
 46     }
 47     //点位置不同(double类型)
 48     bool operator != (const Point &_Off) const
 49     {
 50         return ((*this) == _Off) == false;
 51     }
 52     //两点间距离的平方
 53     static double dis2(const Point &_st,const Point &_se)
 54     {
 55         return (_st.x - _se.x) * (_st.x - _se.x) + (_st.y - _se.y) * (_st.y - _se.y);
 56     }
 57     //两点间距离
 58     static double dis(const Point &_st, const Point &_se)
 59     {
 60         return sqrt((_st.x - _se.x) * (_st.x - _se.x) + (_st.y - _se.y) * (_st.y - _se.y));
 61     }
 62 };
 63 //两点表示的向量
 64 class Line
 65 {
 66 public:
 67     Point s, e;//两点表示，起点[s]，终点[e]
 68     double a, b, c;//一般式,ax+by+c=0
 69     double angle;//向量的角度，[-pi,pi]
 70     Line(){}
 71     Line(const Point &s, const Point &e):s(s),e(e){get_angle(1);}
 72     Line(double _a,double _b,double _c):a(_a),b(_b),c(_c){}
 73 
 74     //向量与点的叉乘,参数：点[_Off]
 75     //[点相对向量位置判断]
 76     double operator /(const Point &_Off) const
 77     {
 78         return (_Off.y - s.y) * (e.x - s.x) - (_Off.x - s.x) * (e.y - s.y);
 79     }
 80     //向量与向量的叉乘,参数：向量[_Off]
 81     friend double operator /(const Line &_st,const Line &_se)
 82     {
 83         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);
 84     }
 85     friend double operator *(const Line &_st,const Line &_se)
 86     {
 87         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);
 88     }
 89     //从两点表示转换为一般表示
 90     //a=y2-y1,b=x1-x2,c=x2*y1-x1*y2
 91     bool pton()
 92     {
 93         a = e.y - s.y;
 94         b = s.x - e.x;
 95         c = e.x * s.y - e.y * s.x;
 96         return true;
 97     }
 98     //求直线或向量的角度
 99     double get_angle(bool isVector)
100     {
101         angle=atan2(e.y-s.y,e.x-s.x);
102         if(!isVector && angle<0)
103             angle+=PI;
104         return angle;
105     }
106 
107     //
108     bool operator < (const Line &ta)const
109     {
110         return angle<ta.angle;
111     }
112     //-----------点和直线（向量）-----------
113     //点在向量左边（右边的小于号改成大于号即可,在对应直线上则加上=号）
114     //参数：点[_Off],向量[_Ori]
115     friend bool operator<(const Point &_Off, const Line &_Ori)
116     {
117         return (_Ori.e.y - _Ori.s.y) * (_Off.x - _Ori.s.x)
118             < (_Off.y - _Ori.s.y) * (_Ori.e.x - _Ori.s.x);
119     }
120 
121     //点在直线上,参数：点[_Off]
122     bool lhas(const Point &_Off) const
123     {
124         return Point::sgn((*this) / _Off, 0) == 0;
125     }
126     //点在线段上,参数：点[_Off]
127     bool shas(const Point &_Off) const
128     {
129         return lhas(_Off)
130             && Point::sgn(_Off.x - min(s.x, e.x), 0) > 0 && Point::sgn(_Off.x - max(s.x, e.x), 0) < 0
131             && Point::sgn(_Off.y - min(s.y, e.y), 0) > 0 && Point::sgn(_Off.y - max(s.y, e.y), 0) < 0;
132     }
133 
134     //点到直线/线段的距离
135     //参数： 点[_Off], 是否是线段[isSegment](默认为直线)
136     double dis(const Point &_Off, bool isSegment = false)
137     {
138         ///化为一般式
139         pton();
140 
141         //到直线垂足的距离
142         double td = (a * _Off.x + b * _Off.y + c) / sqrt(a * a + b * b);
143 
144         //如果是线段判断垂足
145         if(isSegment)
146         {
147             double xp = (b * b * _Off.x - a * b * _Off.y - a * c) / ( a * a + b * b);
148             double yp = (-a * b * _Off.x + a * a * _Off.y - b * c) / (a * a + b * b);
149             double xb = max(s.x, e.x);
150             double yb = max(s.y, e.y);
151             double xs = s.x + e.x - xb;
152             double ys = s.y + e.y - yb;
153             if(xp > xb + eps || xp < xs - eps || yp > yb + eps || yp < ys - eps)
154                 td = min(Point::dis(_Off,s), Point::dis(_Off,e));
155         }
156 
157         return fabs(td);
158     }
159 
160     //关于直线对称的点
161     Point mirror(const Point &_Off) const
162     {
163         ///注意先转为一般式
164         Point ret;
165         double d = a * a + b * b;
166         ret.x = (b * b * _Off.x - a * a * _Off.x - 2 * a * b * _Off.y - 2 * a * c) / d;
167         ret.y = (a * a * _Off.y - b * b * _Off.y - 2 * a * b * _Off.x - 2 * b * c) / d;
168         return ret;
169     }
170     //计算两点的中垂线
171     static Line ppline(const Point &_a, const Point &_b)
172     {
173         Line ret;
174         ret.s.x = (_a.x + _b.x) / 2;
175         ret.s.y = (_a.y + _b.y) / 2;
176         //一般式
177         ret.a = _b.x - _a.x;
178         ret.b = _b.y - _a.y;
179         ret.c = (_a.y - _b.y) * ret.s.y + (_a.x - _b.x) * ret.s.x;
180         //两点式
181         if(std::fabs(ret.a) > eps)
182         {
183             ret.e.y = 0.0;
184             ret.e.x = - ret.c / ret.a;
185             if(ret.e == ret. s)
186             {
187                 ret.e.y = 1e10;
188                 ret.e.x = - (ret.c - ret.b * ret.e.y) / ret.a;
189             }
190         }
191         else
192         {
193             ret.e.x = 0.0;
194             ret.e.y = - ret.c / ret.b;
195             if(ret.e == ret. s)
196             {
197                 ret.e.x = 1e10;
198                 ret.e.y = - (ret.c - ret.a * ret.e.x) / ret.b;
199             }
200         }
201         return ret;
202     }
203 
204     //------------直线和直线（向量）-------------
205     //直线向左边平移t的距离
206     Line& moveLine(double t)
207     {
208         Point of;
209         of=Point(-(e.y-s.y),e.x-s.x);
210         double dis=sqrt(of.x*of.x+of.y*of.y);
211         of.x=of.x*t/dis,of.y=of.y*t/dis;
212         s=s+of,e=e+of;
213         return *this;
214     }
215     //直线重合,参数：直线向量[_st],[_se]
216     static bool equal(const Line &_st, const Line &_se)
217     {
218         return _st.lhas(_se.e) && _se.lhas(_se.s);
219     }
220     //直线平行，参数：直线向量[_st],[_se]
221     static bool parallel(const Line &_st,const Line &_se)
222     {
223         return Point::sgn(_st / _se, 0) == 0;
224     }
225     //两直线（线段）交点，参数：直线向量[_st],[_se],交点
226     //返回-1代表平行，0代表重合，1代表相交
227     static bool crossLPt(const Line &_st,const Line &_se,Point &ret)
228     {
229         if(Line::parallel(_st,_se))
230         {
231             if(Line::equal(_st,_se)) return 0;
232             return -1;
233         }
234         ret = _st.s;
235         double t = (Line(_st.s,_se.s)/_se)/(_st/_se);
236         ret.x += (_st.e.x - _st.s.x) * t;
237         ret.y += (_st.e.y - _st.s.y) * t;
238         return 1;
239     }
240     //------------线段和直线（向量）----------
241     //线段和直线交
242     //参数：直线[_st],线段[_se]
243     friend bool crossSL(const Line &_st,const Line &_se)
244     {
245         return Point::sgn((_st / _se.s) * (_st / _se.e) ,0) <= 0;
246     }
247 
248     //------------线段和线段（向量）----------
249     //判断线段是否相交(注意添加eps)，参数：线段[_st],线段[_se]
250     static bool isCrossSS(const Line &_st,const Line &_se)
251     {
252         //1.快速排斥试验判断以两条线段为对角线的两个矩形是否相交
253         //2.跨立试验（等于0时端点重合）
254         return
255             max(_st.s.x, _st.e.x) >= min(_se.s.x, _se.e.x) &&
256             max(_se.s.x, _se.e.x) >= min(_st.s.x, _st.e.x) &&
257             max(_st.s.y, _st.e.y) >= min(_se.s.y, _se.e.y) &&
258             max(_se.s.y, _se.e.y) >= min(_st.s.y, _st.e.y) &&
259             Point::sgn((_st / Line(_st.s, _se.s)) * (_st / Line(_st.s, _se.e)), 0) <= 0 &&
260             Point::sgn((_se / Line(_se.s, _st.s)) * (_se / Line(_se.s, _st.e)), 0) <= 0;
261     }
262 };
263 
264 class Polygon
265 {
266 public:
267     const static int maxpn = 5e4+7;
268     Point pt[maxpn];//点（顺时针或逆时针）
269     int n;//点的个数
270 
271 
272     //求多边形面积，多边形内点必须顺时针或逆时针
273     double area() const
274     {
275         double ans = 0.0;
276         for(int i = 0; i < n; i ++)
277         {
278             int nt = (i + 1) % n;
279             ans += pt[i].x * pt[nt].y - pt[nt].x * pt[i].y;
280         }
281         return fabs(ans / 2.0);
282     }
283     //求多边形重心，多边形内点必须顺时针或逆时针
284     Point gravity() const
285     {
286         Point ans;
287         ans.x = ans.y = 0.0;
288         double area = 0.0;
289         for(int i = 0; i < n; i ++)
290         {
291             int nt = (i + 1) % n;
292             double tp = pt[i].x * pt[nt].y - pt[nt].x * pt[i].y;
293             area += tp;
294             ans.x += tp * (pt[i].x + pt[nt].x);
295             ans.y += tp * (pt[i].y + pt[nt].y);
296         }
297         ans.x /= 3 * area;
298         ans.y /= 3 * area;
299         return ans;
300     }
301     //判断点在凸多边形内，参数：点[_Off]
302     bool chas(const Point &_Off) const
303     {
304         double tp = 0, np;
305         for(int i = 0; i < n; i ++)
306         {
307             np = Line(pt[i], pt[(i + 1) % n]) / _Off;
308             if(tp * np < -eps)
309                 return false;
310             tp = (fabs(np) > eps)?np: tp;
311         }
312         return true;
313     }
314     //判断点是否在任意多边形内[射线法]，O(n)
315     bool ahas(const Point &_Off) const
316     {
317         int ret = 0;
318         double infv = 1e20;//坐标系最大范围
319         Line l = Line(_Off, Point( -infv ,_Off.y));
320         for(int i = 0; i < n; i ++)
321         {
322             Line ln = Line(pt[i], pt[(i + 1) % n]);
323             if(fabs(ln.s.y - ln.e.y) > eps)
324             {
325                 Point tp = (ln.s.y > ln.e.y)? ln.s: ln.e;
326                 if((fabs(tp.y - _Off.y) < eps && tp.x < _Off.x + eps)||Line::isCrossSS(ln,l))
327                     ret ++;
328             }
329             else if(Line::isCrossSS(ln,l))
330                 ret ++;
331         }
332         return ret&1;
333     }
334 
335     //判断任意点是否在凸包内，O(lgn)
336     bool bhas(const Point & p)
337     {
338         if(n<3)
339             return false;
340         if(Point::xmult(pt[0],p,pt[1])>eps)
341             return false;
342         if(Point::xmult(pt[0],p,pt[n-1])<-eps)
343             return false;
344         int l=2,r=n-1;
345         int line=-1;
346         while(l<=r)
347         {
348             int mid=(l+r)>>1;
349             if(Point::xmult(pt[0],p,pt[mid])>=0)
350                 line=mid,r=mid-1;
351             else l=mid+1;
352         }
353         return Point::xmult(pt[line-1],p,pt[line])<=eps;
354     }
355 
356 
357 
358     //凸多边形被直线分割,参数：直线[_Off]
359     Polygon split(Line _Off)
360     {
361         //注意确保多边形能被分割
362         Polygon ret;
363         Point spt[2];
364         double tp = 0.0, np;
365         bool flag = true;
366         int i, pn = 0, spn = 0;
367         for(i = 0; i < n; i ++)
368         {
369             if(flag)
370                 pt[pn ++] = pt[i];
371             else
372                 ret.pt[ret.n ++] = pt[i];
373             np = _Off / pt[(i + 1) % n];
374             if(tp * np < -eps)
375             {
376                 flag = !flag;
377                 Line::crossLPt(_Off,Line(pt[i], pt[(i + 1) % n]),spt[spn++]);
378             }
379             tp = (fabs(np) > eps)?np: tp;
380         }
381         ret.pt[ret.n ++] = spt[0];
382         ret.pt[ret.n ++] = spt[1];
383         n = pn;
384         return ret;
385     }
386 
387 
388     /** 卷包裹法求点集凸包，_p为输入点集，_n为点的数量 **/
389     void ConvexClosure(Point _p[],int _n)
390     {
391         sort(_p,_p+_n);
392         n=0;
393         for(int i=0;i<_n;i++)
394         {
395             while(n>1&&Point::sgn(Line(pt[n-2],pt[n-1])/Line(pt[n-2],_p[i]),0)<=0)
396                 n--;
397             pt[n++]=_p[i];
398         }
399         int _key=n;
400         for(int i=_n-2;i>=0;i--)
401         {
402             while(n>_key&&Point::sgn(Line(pt[n-2],pt[n-1])/Line(pt[n-2],_p[i]),0)<=0)
403                 n--;
404             pt[n++]=_p[i];
405         }
406         if(n>1)   n--;//除去重复的点，该点已是凸包凸包起点
407     }
408 //    /****** 寻找凸包的graham 扫描法********************/
409 //    /****** _p为输入的点集,_n为点的数量****************/
410 //    /**使用时需把gmp函数放在Polygon类上面L，ine类下面，并且看情况修改pt[0]**/
411 //    bool gcmp(const Point &ta,const Point &tb)/// 选取与最后一条确定边夹角最小的点，即余弦值最大者
412 //    {
413 //        double tmp=Line(pt[0],ta)/Line(pt[0],tb);
414 //        if(Point::sgn(tmp,0)==0)
415 //            return Point::dis(pt[0],ta)<Point::dis(pt[0],tb);
416 //        else if(tmp>0)
417 //            return 1;
418 //        return 0;
419 //    }
420 //    void graham(Point _p[],int _n)
421 //    {
422 //        int cur=0;
423 //        for(int i=1;i<_n;i++)
424 //            if(Point::sgn(_p[cur].y,_p[i].y)>0 || (Point::sgn(_p[cur].y,_p[i].y)==0 && Point::sgn(_p[cur].x,_p[i].x)>0))
425 //                cur=i;
426 //        swap(_p[cur],_p[0]);
427 //        n=0,pt[n++]=_p[0];
428 //        if(_n==1)   return;
429 //        sort(_p+1,_p+_n,Polygon::gcmp);
430 //        pt[n++]=_p[1],pt[n++]=_p[2];
431 //        for(int i=3;i<_n;i++)
432 //        {
433 //            while(n>1 && Point::sgn(Line(pt[n-2],pt[n-1])/Line(pt[n-2],_p[i]),0)<0)// 当凸包退化成直线时需特别注意n
434 //                n--;
435 //            pt[n++]=_p[i];
436 //        }
437 //    }
438     //凸包旋转卡壳(注意点必须顺时针或逆时针排列)
439     //返回值凸包直径的平方（最远两点距离的平方）
440     double rotating_calipers()
441     {
442         int i = 1;
443         double ret = 0.0;
444         pt[n] = pt[0];
445         for(int j = 0; j < n; j ++)
446         {
447             while(fabs(Point::xmult(pt[i+1],pt[j], pt[j + 1])) > fabs(Point::xmult(pt[i],pt[j], pt[j + 1])) + eps)
448                 i = (i + 1) % n;
449             //pt[i]和pt[j],pt[i + 1]和pt[j + 1]可能是对踵点
450             ret = max(ret, max(Point::dis(pt[i],pt[j]), Point::dis(pt[i + 1],pt[j + 1])));
451         }
452         return ret;
453     }
454 
455     //凸包旋转卡壳(注意点必须逆时针排列)
456     //返回值两凸包的最短距离
457     double rotating_calipers(Polygon &_Off)
458     {
459         int i = 0;
460         double ret = 1e10;//inf
461         pt[n] = pt[0];
462         _Off.pt[_Off.n] = _Off.pt[0];
463         //注意凸包必须逆时针排列且pt[0]是左下角点的位置
464         while(_Off.pt[i + 1].y > _Off.pt[i].y)
465             i = (i + 1) % _Off.n;
466         for(int j = 0; j < n; j ++)
467         {
468             double tp;
469             //逆时针时为 >,顺时针则相反
470             while((tp = Point::xmult(_Off.pt[i + 1],pt[j], pt[j + 1]) - Point::xmult(_Off.pt[i], pt[j], pt[j + 1])) > eps)
471                 i = (i + 1) % _Off.n;
472             //(pt[i],pt[i+1])和(_Off.pt[j],_Off.pt[j + 1])可能是最近线段
473             ret = min(ret, Line(pt[j], pt[j + 1]).dis(_Off.pt[i], true));
474             ret = min(ret, Line(_Off.pt[i], _Off.pt[i + 1]).dis(pt[j + 1], true));
475             if(tp > -eps)//如果不考虑TLE问题最好不要加这个判断
476             {
477                 ret = min(ret, Line(pt[j], pt[j + 1]).dis(_Off.pt[i + 1], true));
478                 ret = min(ret, Line(_Off.pt[i], _Off.pt[i + 1]).dis(pt[j], true));
479             }
480         }
481         return ret;
482     }
483 
484     //-----------半平面交-------------
485     //复杂度:O(nlog2(n))
486     //#include <algorithm>
487 
488     //获取半平面交的多边形（多边形的核）
489     //参数：向量集合[l]，向量数量[ln];(半平面方向在向量左边）
490     //函数运行后如果n[即返回多边形的点数量]为0则不存在半平面交的多边形（不存在区域或区域面积无穷大）
491     int halfPanelCross(Line *L, int ln)
492     {
493         sort(L,L + ln);     //极角排序
494         int first,last;     //双端队列的队首和队尾的下标
495         Point *p = new Point[ln];   //p[i]为q[i]和q[i+1]的交点
496         Line *q = new Line[ln];     //双端队列
497         q[first=last=0] = L[0];     //双端队列初始只有一个半平面L[0]
498         for(int i = 1; i < ln; i++)
499         {
500             while(first < last && !(p[last-1] < L[i])) last--;
501             while(first < last && !(p[first] < L[i])) first++;
502             q[++last] = L[i];
503             if(fabs(q[last]/q[last-1]) < eps)
504             {
505                 if(L[i].s < q[--last]) q[last] = L[i];  //两向量平行，取内侧的一个
506             }
507             if(first < last) Line::crossLPt(q[last-1], q[last], p[last-1]);
508         }
509         while(first < last && !(p[last-1] < q[first])) last--; //删除无用平面
510         if(last - first <= 1) return 0; //空集
511         Line::crossLPt(q[last], q[first], p[last]);   //计算首尾两个半平面的交点
512         int m=0;
513         for(int i = 1; i <= last; i++) pt[m++] = p[i];  //复制到pt中
514         return m;
515     }
516 };
517 
518 class Circle
519 {
520 public:
521     Point c;//圆心
522     double r;//半径
523     double db, de;//圆弧度数起点， 圆弧度数终点(逆时针0-360)
524 
525     //-------圆---------
526 
527     //判断圆在多边形内
528     bool inside(const Polygon &_Off) const
529     {
530         if(_Off.ahas(c) == false)
531             return false;
532         for(int i = 0; i < _Off.n; i ++)
533         {
534             Line l = Line(_Off.pt[i], _Off.pt[(i + 1) % _Off.n]);
535             if(l.dis(c, true) < r - eps)
536                 return false;
537         }
538         return true;
539     }
540 
541     //判断多边形在圆内（线段和折线类似）
542     bool has(const Polygon &_Off) const
543     {
544         for(int i = 0; i < _Off.n; i ++)
545             if(Point::dis2(_Off.pt[i],c) > r * r - eps)
546                 return false;
547         return true;
548     }
549 
550     //-------圆弧-------
551     //圆被其他圆截得的圆弧，参数：圆[_Off]
552     Circle operator-(Circle &_Off) const
553     {
554         //注意圆必须相交，圆心不能重合
555         double d2 = Point::dis2(c,_Off.c);
556         double d = Point::dis(c,_Off.c);
557         double ans = acos((d2 + r * r - _Off.r * _Off.r) / (2 * d * r));
558         Point py = _Off.c - c;
559         double oans = atan2(py.y, py.x);
560         Circle res;
561         res.c = c;
562         res.r = r;
563         res.db = oans + ans;
564         res.de = oans - ans + 2 * PI;
565         return res;
566     }
567     //圆被其他圆截得的圆弧，参数：圆[_Off]
568     Circle operator+(Circle &_Off) const
569     {
570         //注意圆必须相交，圆心不能重合
571         double d2 = Point::dis2(c,_Off.c);
572         double d = Point::dis(c,_Off.c);
573         double ans = acos((d2 + r * r - _Off.r * _Off.r) / (2 * d * r));
574         Point py = _Off.c - c;
575         double oans = atan2(py.y, py.x);
576         Circle res;
577         res.c = c;
578         res.r = r;
579         res.db = oans - ans;
580         res.de = oans + ans;
581         return res;
582     }
583 
584     //过圆外一点的两条切线
585     //参数：点[_Off](必须在圆外),返回：两条切线(切线的s点为_Off,e点为切点)
586     std::pair<Line, Line>  tangent(const Point &_Off) const
587     {
588         double d = Point::dis(c,_Off);
589         //计算角度偏移的方式
590         double angp = std::acos(r / d), ango = std::atan2(_Off.y - c.y, _Off.x - c.x);
591         Point pl = Point(c.x + r * std::cos(ango + angp), c.y + r * std::sin(ango + angp)),
592             pr = Point(c.x + r * std::cos(ango - angp), c.y + r * std::sin(ango - angp));
593         return std::make_pair(Line(_Off, pl), Line(_Off, pr));
594     }
595 
596     //计算直线和圆的两个交点
597     //参数：直线[_Off](两点式)，返回两个交点，注意直线必须和圆有两个交点
598     std::pair<Point, Point> cross(Line _Off) const
599     {
600         _Off.pton();
601         //到直线垂足的距离
602         double td = fabs(_Off.a * c.x + _Off.b * c.y + _Off.c) / sqrt(_Off.a * _Off.a + _Off.b * _Off.b);
603 
604         //计算垂足坐标
605         double xp = (_Off.b * _Off.b * c.x - _Off.a * _Off.b * c.y - _Off.a * _Off.c) / ( _Off.a * _Off.a + _Off.b * _Off.b);
606         double yp = (- _Off.a * _Off.b * c.x + _Off.a * _Off.a * c.y - _Off.b * _Off.c) / (_Off.a * _Off.a + _Off.b * _Off.b);
607 
608         double ango = std::atan2(yp - c.y, xp - c.x);
609         double angp = std::acos(td / r);
610 
611         return std::make_pair(Point(c.x + r * std::cos(ango + angp), c.y + r * std::sin(ango + angp)),
612             Point(c.x + r * std::cos(ango - angp), c.y + r * std::sin(ango - angp)));
613     }
614 };
615 
616 class triangle
617 {
618 public:
619     Point a, b, c;//顶点
620     triangle(){}
621     triangle(Point a, Point b, Point c): a(a), b(b), c(c){}
622 
623     //计算三角形面积
624     double area()
625     {
626         return fabs(Point::xmult(a, b, c)) / 2.0;
627     }
628 
629     //计算三角形外心
630     //返回：外接圆圆心
631     Point circumcenter()
632     {
633         Line u,v;
634         u.s.x = (a.x + b.x) / 2;
635         u.s.y = (a.y + b.y) / 2;
636         u.e.x = u.s.x - a.y + b.y;
637         u.e.y = u.s.y + a.x - b.x;
638         v.s.x = (a.x + c.x) / 2;
639         v.s.y = (a.y + c.y) / 2;
640         v.e.x = v.s.x - a.y + c.y;
641         v.e.y = v.s.y + a.x - c.x;
642         Point ret;
643         Line::crossLPt(u,v,ret);
644         return ret;
645     }
646 
647     //计算三角形内心
648     //返回：内接圆圆心
649     Point incenter()
650     {
651         Line u, v;
652         double m, n;
653         u.s = a;
654         m = atan2(b.y - a.y, b.x - a.x);
655         n = atan2(c.y - a.y, c.x - a.x);
656         u.e.x = u.s.x + cos((m + n) / 2);
657         u.e.y = u.s.y + sin((m + n) / 2);
658         v.s = b;
659         m = atan2(a.y - b.y, a.x - b.x);
660         n = atan2(c.y - b.y, c.x - b.x);
661         v.e.x = v.s.x + cos((m + n) / 2);
662         v.e.y = v.s.y + sin((m + n) / 2);
663         Point ret;
664         Line::crossLPt(u,v,ret);
665         return ret;
666     }
667 
668     //计算三角形垂心
669     //返回：高的交点
670     Point perpencenter()
671     {
672         Line u,v;
673         u.s = c;
674         u.e.x = u.s.x - a.y + b.y;
675         u.e.y = u.s.y + a.x - b.x;
676         v.s = b;
677         v.e.x = v.s.x - a.y + c.y;
678         v.e.y = v.s.y + a.x - c.x;
679         Point ret;
680         Line::crossLPt(u,v,ret);
681         return ret;
682     }
683 
684     //计算三角形重心
685     //返回：重心
686     //到三角形三顶点距离的平方和最小的点
687     //三角形内到三边距离之积最大的点
688     Point barycenter()
689     {
690         Line u,v;
691         u.s.x = (a.x + b.x) / 2;
692         u.s.y = (a.y + b.y) / 2;
693         u.e = c;
694         v.s.x = (a.x + c.x) / 2;
695         v.s.y = (a.y + c.y) / 2;
696         v.e = b;
697         Point ret;
698         Line::crossLPt(u,v,ret);
699         return ret;
700     }
701 
702     //计算三角形费马点
703     //返回：到三角形三顶点距离之和最小的点
704     Point fermentPoint()
705     {
706         Point u, v;
707         double step = fabs(a.x) + fabs(a.y) + fabs(b.x) + fabs(b.y) + fabs(c.x) + fabs(c.y);
708         int i, j, k;
709         u.x = (a.x + b.x + c.x) / 3;
710         u.y = (a.y + b.y + c.y) / 3;
711         while (step > eps)
712         {
713             for (k = 0; k < 10; step /= 2, k ++)
714             {
715                 for (i = -1; i <= 1; i ++)
716                 {
717                     for (j =- 1; j <= 1; j ++)
718                     {
719                         v.x = u.x + step * i;
720                         v.y = u.y + step * j;
721                         if (Point::dis(u,a) + Point::dis(u,b) + Point::dis(u,c) > Point::dis(v,a) + Point::dis(v,b) + Point::dis(v,c))
722                             u = v;
723                     }
724                 }
725             }
726         }
727         return u;
728     }
729 };
730 
731 int main(void)
732 {
733 
734     return 0;
735 }
【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/weeping/p/6363742.html