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 }
扫一扫
847
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
