计算几何：模板

叉积：两个向量的叉积是一个标量，a×b的几何意义为他们所形成的平行四边形的有向面积
凸包：把给定点包围在内部，面积最小的凸多边形
半平面交：一个半平面指的是由满足 $a x + b y + c > 0$
$a x + b y + c > 0$
$a x + b y + c > 0$ $a x + b y + c > 0$
$a x + b y + c > 0$
   1 /*
   2 计算几何
   3 目录 
   4 ㈠ 点的基本运算 
   5 1. 平面上两点之间距离 1 
   6 2. 判断两点是否重合 1 
   7 3. 矢量叉乘 1 
   8 4. 矢量点乘 2 
   9 5. 判断点是否在线段上 2 
  10 6. 求一点饶某点旋转后的坐标 2 
  11 7. 求矢量夹角 2 
  12 ㈡ 线段及直线的基本运算 
  13 1. 点与线段的关系 3 
  14 2. 求点到线段所在直线垂线的垂足 4 
  15 3. 点到线段的最近点 4 
  16 4. 点到线段所在直线的距离 4 
  17 5. 点到折线集的最近距离 4 
  18 6. 判断圆是否在多边形内 5 
  19 7. 求矢量夹角余弦 5 
  20 8. 求线段之间的夹角 5 
  21 9. 判断线段是否相交 6 
  22 10.判断线段是否相交但不交在端点处 6 
  23 11.求线段所在直线的方程 6 
  24 12.求直线的斜率 7 
  25 13.求直线的倾斜角 7 
  26 14.求点关于某直线的对称点 7 
  27 15.判断两条直线是否相交及求直线交点 7 
  28 16.判断线段是否相交，如果相交返回交点 7 
  29 ㈢ 多边形常用算法模块 
  30 1. 判断多边形是否简单多边形 8 
  31 2. 检查多边形顶点的凸凹性 9 
  32 3. 判断多边形是否凸多边形 9 
  33 4. 求多边形面积 9 
  34 5. 判断多边形顶点的排列方向，方法一 10 
  35 6. 判断多边形顶点的排列方向，方法二 10 
  36 7. 射线法判断点是否在多边形内 10 
  37 8. 判断点是否在凸多边形内 11 
  38 9. 寻找点集的graham算法 12 
  39 10.寻找点集凸包的卷包裹法 13 
  40 11.判断线段是否在多边形内 14 
  41 12.求简单多边形的重心 15 
  42 13.求凸多边形的重心 17 
  43 14.求肯定在给定多边形内的一个点 17 
  44 15.求从多边形外一点出发到该多边形的切线 18 
  45 16.判断多边形的核是否存在 19 
  46 ㈣ 圆的基本运算 
  47 1 .点是否在圆内 20 
  48 2 .求不共线的三点所确定的圆 21 
  49 ㈤ 矩形的基本运算 
  50 1.已知矩形三点坐标，求第4点坐标 22 
  51 ㈥ 常用算法的描述 22 
  52 ㈦ 补充 
  53 1．两圆关系： 24 
  54 2．判断圆是否在矩形内： 24 
  55 3．点到平面的距离： 25 
  56 4．点是否在直线同侧： 25 
  57 5．镜面反射线： 25 
  58 6．矩形包含： 26 
  59 7．两圆交点： 27 
  60 8．两圆公共面积： 28 
  61 9. 圆和直线关系： 29 
  62 10. 内切圆： 30 
  63 11. 求切点： 31 
  64 12. 线段的左右旋： 31 
  65 13．公式： 32 
  66 */
  67 /* 需要包含的头文件 */ 
  68 #include <cmath > 
  69  
  70 /* 常用的常量定义 */ 
  71 const double    INF        = 1E200    
  72 const double    EP        = 1E-10 
  73 const int        MAXV    = 300 
  74 const double    PI        = 3.14159265 
  75  
  76 /* 基本几何结构 */ 
  77 struct POINT 
  78 { 
  79     double x; 
  80     double y; 
  81     POINT(double a=0, double b=0) { x=a; y=b;} //constructor 
  82 }; 
  83 struct LINESEG 
  84 { 
  85     POINT s; 
  86     POINT e; 
  87     LINESEG(POINT a, POINT b) { s=a; e=b;} 
  88     LINESEG() { } 
  89 }; 
  90 struct LINE           // 直线的解析方程 a*x+b*y+c=0  为统一表示，约定 a >= 0 
  91 { 
  92    double a; 
  93    double b; 
  94    double c; 
  95    LINE(double d1=1, double d2=-1, double d3=0) {a=d1; b=d2; c=d3;} 
  96 }; 
  97  
  98 /**********************
  99  *                    * 
 100  *   点的基本运算     * 
 101  *                    * 
 102  **********************/ 
 103  
 104 double dist(POINT p1,POINT p2)                // 返回两点之间欧氏距离 
 105 { 
 106     return( sqrt( (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) ) ); 
 107 } 
 108 bool equal_point(POINT p1,POINT p2)           // 判断两个点是否重合  
 109 { 
 110     return ( (abs(p1.x-p2.x)<EP)&&(abs(p1.y-p2.y)<EP) ); 
 111 } 
 112 /****************************************************************************** 
 113 r=multiply(sp,ep,op),得到(sp-op)和(ep-op)的叉积 
 114 r>0：ep在矢量opsp的逆时针方向； 
 115 r=0：opspep三点共线； 
 116 r<0：ep在矢量opsp的顺时针方向 
 117 *******************************************************************************/ 
 118 double multiply(POINT sp,POINT ep,POINT op) 
 119 { 
 120     return((sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y)); 
 121 } 
 122 /* 
 123 r=dotmultiply(p1,p2,op),得到矢量(p1-op)和(p2-op)的点积，如果两个矢量都非零矢量 
 124 r<0：两矢量夹角为锐角；
 125 r=0：两矢量夹角为直角；
 126 r>0：两矢量夹角为钝角 
 127 *******************************************************************************/ 
 128 double dotmultiply(POINT p1,POINT p2,POINT p0) 
 129 { 
 130     return ((p1.x-p0.x)*(p2.x-p0.x)+(p1.y-p0.y)*(p2.y-p0.y)); 
 131 } 
 132 /****************************************************************************** 
 133 判断点p是否在线段l上
 134 条件：(p在线段l所在的直线上) && (点p在以线段l为对角线的矩形内)
 135 *******************************************************************************/ 
 136 bool online(LINESEG l,POINT p) 
 137 { 
 138     return( (multiply(l.e,p,l.s)==0) &&( ( (p.x-l.s.x)*(p.x-l.e.x)<=0 )&&( (p.y-l.s.y)*(p.y-l.e.y)<=0 ) ) ); 
 139 } 
 140 // 返回点p以点o为圆心逆时针旋转alpha(单位：弧度)后所在的位置 
 141 POINT rotate(POINT o,double alpha,POINT p) 
 142 { 
 143     POINT tp; 
 144     p.x-=o.x; 
 145     p.y-=o.y; 
 146     tp.x=p.x*cos(alpha)-p.y*sin(alpha)+o.x; 
 147     tp.y=p.y*cos(alpha)+p.x*sin(alpha)+o.y; 
 148     return tp; 
 149 } 
 150 /* 返回顶角在o点，起始边为os，终止边为oe的夹角(单位：弧度) 
 151     角度小于pi，返回正值 
 152     角度大于pi，返回负值 
 153     可以用于求线段之间的夹角 
 154 原理：
 155     r = dotmultiply(s,e,o) / (dist(o,s)*dist(o,e))
 156     r'= multiply(s,e,o)
 157     r >= 1    angle = 0;
 158     r <= -1    angle = -PI
 159     -1<r<1 && r'>0    angle = arccos(r)
 160     -1<r<1 && r'<=0    angle = -arccos(r)
 161 */ 
 162 double angle(POINT o,POINT s,POINT e) 
 163 { 
 164     double cosfi,fi,norm; 
 165     double dsx = s.x - o.x; 
 166     double dsy = s.y - o.y; 
 167     double dex = e.x - o.x; 
 168     double dey = e.y - o.y; 
 169  
 170     cosfi=dsx*dex+dsy*dey; 
 171     norm=(dsx*dsx+dsy*dsy)*(dex*dex+dey*dey); 
 172     cosfi /= sqrt( norm ); 
 173  
 174     if (cosfi >=  1.0 ) return 0; 
 175     if (cosfi <= -1.0 ) return -3.1415926; 
 176  
 177     fi=acos(cosfi); 
 178     if (dsx*dey-dsy*dex>0) return fi;      // 说明矢量os 在矢量 oe的顺时针方向 
 179     return -fi; 
 180 } 
 181   /***************************** 
 182   *                             * 
 183   *      线段及直线的基本运算   * 
 184   *                             * 
 185   *****************************/ 
 186  
 187 /* 判断点与线段的关系,用途很广泛 
 188 本函数是根据下面的公式写的，P是点C到线段AB所在直线的垂足 
 189                 AC dot AB 
 190         r =     --------- 
 191                  ||AB||^2 
 192              (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay) 
 193           = ------------------------------- 
 194                           L^2 
 195     r has the following meaning: 
 196         r=0      P = A 
 197         r=1      P = B 
 198         r<0         P is on the backward extension of AB 
 199         r>1      P is on the forward extension of AB 
 200         0<r<1     P is interior to AB 
 201 */ 
 202 double relation(POINT p,LINESEG l) 
 203 { 
 204     LINESEG tl; 
 205     tl.s=l.s; 
 206     tl.e=p; 
 207     return dotmultiply(tl.e,l.e,l.s)/(dist(l.s,l.e)*dist(l.s,l.e)); 
 208 } 
 209 // 求点C到线段AB所在直线的垂足 P 
 210 POINT perpendicular(POINT p,LINESEG l) 
 211 { 
 212     double r=relation(p,l); 
 213     POINT tp; 
 214     tp.x=l.s.x+r*(l.e.x-l.s.x); 
 215     tp.y=l.s.y+r*(l.e.y-l.s.y); 
 216     return tp; 
 217 } 
 218 /* 求点p到线段l的最短距离,并返回线段上距该点最近的点np 
 219 注意：np是线段l上到点p最近的点，不一定是垂足 */ 
 220 double ptolinesegdist(POINT p,LINESEG l,POINT &np) 
 221 { 
 222     double r=relation(p,l); 
 223     if(r<0) 
 224     { 
 225         np=l.s; 
 226         return dist(p,l.s); 
 227     } 
 228     if(r>1) 
 229     { 
 230         np=l.e; 
 231         return dist(p,l.e); 
 232     } 
 233     np=perpendicular(p,l); 
 234     return dist(p,np); 
 235 } 
 236 // 求点p到线段l所在直线的距离,请注意本函数与上个函数的区别  
 237 double ptoldist(POINT p,LINESEG l) 
 238 { 
 239     return abs(multiply(p,l.e,l.s))/dist(l.s,l.e); 
 240 } 
 241 /* 计算点到折线集的最近距离,并返回最近点. 
 242 注意：调用的是ptolineseg()函数 */ 
 243 double ptopointset(int vcount,POINT pointset[],POINT p,POINT &q) 
 244 { 
 245     int i; 
 246     double cd=double(INF),td; 
 247     LINESEG l; 
 248     POINT tq,cq; 
 249  
 250     for(i=0;i<vcount-1;i++) 
 251     { 
 252         l.s=pointset[i]; 
 253  
 254         l.e=pointset[i+1]; 
 255         td=ptolinesegdist(p,l,tq); 
 256         if(td<cd) 
 257         { 
 258             cd=td; 
 259             cq=tq; 
 260         } 
 261     } 
 262     q=cq; 
 263     return cd; 
 264 } 
 265 /* 判断圆是否在多边形内.ptolineseg()函数的应用2 */ 
 266 bool CircleInsidePolygon(int vcount,POINT center,double radius,POINT polygon[]) 
 267 { 
 268     POINT q; 
 269     double d; 
 270     q.x=0; 
 271     q.y=0; 
 272     d=ptopointset(vcount,polygon,center,q); 
 273     if(d<radius||fabs(d-radius)<EP) 
 274         return true; 
 275     else 
 276         return false; 
 277 } 
 278 /* 返回两个矢量l1和l2的夹角的余弦(-1 --- 1)注意：如果想从余弦求夹角的话，注意反余弦函数的定义域是从 0到pi */ 
 279 double cosine(LINESEG l1,LINESEG l2) 
 280 { 
 281     return (((l1.e.x-l1.s.x)*(l2.e.x-l2.s.x) + 
 282     (l1.e.y-l1.s.y)*(l2.e.y-l2.s.y))/(dist(l1.e,l1.s)*dist(l2.e,l2.s))) ); 
 283 } 
 284 // 返回线段l1与l2之间的夹角 单位：弧度 范围(-pi，pi) 
 285 double lsangle(LINESEG l1,LINESEG l2) 
 286 { 
 287     POINT o,s,e; 
 288     o.x=o.y=0; 
 289     s.x=l1.e.x-l1.s.x; 
 290     s.y=l1.e.y-l1.s.y; 
 291     e.x=l2.e.x-l2.s.x; 
 292     e.y=l2.e.y-l2.s.y; 
 293     return angle(o,s,e); 
 294 } 
 295 // 如果线段u和v相交(包括相交在端点处)时，返回true 
 296 //
 297 //判断P1P2跨立Q1Q2的依据是：( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) >= 0。
 298 //判断Q1Q2跨立P1P2的依据是：( Q1 - P1 ) × ( P2 - P1 ) * ( P2 - P1 ) × ( Q2 - P1 ) >= 0。
 299 bool intersect(LINESEG u,LINESEG v) 
 300 { 
 301     return( (max(u.s.x,u.e.x)>=min(v.s.x,v.e.x))&&                     //排斥实验 
 302             (max(v.s.x,v.e.x)>=min(u.s.x,u.e.x))&& 
 303             (max(u.s.y,u.e.y)>=min(v.s.y,v.e.y))&& 
 304             (max(v.s.y,v.e.y)>=min(u.s.y,u.e.y))&& 
 305             (multiply(v.s,u.e,u.s)*multiply(u.e,v.e,u.s)>=0)&&         //跨立实验 
 306             (multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0)); 
 307 } 
 308 //  (线段u和v相交)&&(交点不是双方的端点) 时返回true    
 309 bool intersect_A(LINESEG u,LINESEG v) 
 310 { 
 311     return    ((intersect(u,v))&& 
 312             (!online(u,v.s))&& 
 313             (!online(u,v.e))&& 
 314             (!online(v,u.e))&& 
 315             (!online(v,u.s))); 
 316 } 
 317 // 线段v所在直线与线段u相交时返回true；方法：判断线段u是否跨立线段v  
 318 bool intersect_l(LINESEG u,LINESEG v) 
 319 { 
 320     return multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0; 
 321 } 
 322 // 根据已知两点坐标，求过这两点的直线解析方程： a*x+b*y+c = 0  (a >= 0)  
 323 LINE makeline(POINT p1,POINT p2) 
 324 { 
 325     LINE tl; 
 326     int sign = 1; 
 327     tl.a=p2.y-p1.y; 
 328     if(tl.a<0) 
 329     { 
 330         sign = -1; 
 331         tl.a=sign*tl.a; 
 332     } 
 333     tl.b=sign*(p1.x-p2.x); 
 334     tl.c=sign*(p1.y*p2.x-p1.x*p2.y); 
 335     return tl; 
 336 } 
 337 // 根据直线解析方程返回直线的斜率k,水平线返回 0,竖直线返回 1e200 
 338 double slope(LINE l) 
 339 { 
 340     if(abs(l.a) < 1e-20)
 341         return 0; 
 342     if(abs(l.b) < 1e-20)
 343         return INF; 
 344     return -(l.a/l.b); 
 345 } 
 346 // 返回直线的倾斜角alpha ( 0 - pi) 
 347 double alpha(LINE l) 
 348 { 
 349     if(abs(l.a)< EP)
 350         return 0; 
 351     if(abs(l.b)< EP)
 352         return PI/2; 
 353     double k=slope(l); 
 354     if(k>0) 
 355         return atan(k); 
 356     else 
 357         return PI+atan(k); 
 358 } 
 359 // 求点p关于直线l的对称点  
 360 POINT symmetry(LINE l,POINT p) 
 361 { 
 362    POINT tp; 
 363    tp.x=((l.b*l.b-l.a*l.a)*p.x-2*l.a*l.b*p.y-2*l.a*l.c)/(l.a*l.a+l.b*l.b); 
 364    tp.y=((l.a*l.a-l.b*l.b)*p.y-2*l.a*l.b*p.x-2*l.b*l.c)/(l.a*l.a+l.b*l.b); 
 365    return tp; 
 366 } 
 367 // 如果两条直线 l1(a1*x+b1*y+c1 = 0), l2(a2*x+b2*y+c2 = 0)相交，返回true，且返回交点p  
 368 bool lineintersect(LINE l1,LINE l2,POINT &p) // 是 L1，L2 
 369 { 
 370     double d=l1.a*l2.b-l2.a*l1.b; 
 371     if(abs(d)<EP) // 不相交 
 372         return false; 
 373     p.x = (l2.c*l1.b-l1.c*l2.b)/d; 
 374     p.y = (l2.a*l1.c-l1.a*l2.c)/d; 
 375     return true; 
 376 } 
 377 // 如果线段l1和l2相交，返回true且交点由(inter)返回，否则返回false 
 378 bool intersection(LINESEG l1,LINESEG l2,POINT &inter) 
 379 { 
 380     LINE ll1,ll2; 
 381     ll1=makeline(l1.s,l1.e); 
 382     ll2=makeline(l2.s,l2.e); 
 383     if(lineintersect(ll1,ll2,inter)) 
 384         return online(l1,inter); 
 385     else 
 386         return false; 
 387 } 
 388  
 389 /****************************** 
 390 *                              * 
 391 * 多边形常用算法模块          * 
 392 *                              * 
 393 ******************************/ 
 394  
 395 // 如果无特别说明，输入多边形顶点要求按逆时针排列 
 396  
 397 /* 
 398 返回值：输入的多边形是简单多边形，返回true 
 399 要 求：输入顶点序列按逆时针排序 
 400 说 明：简单多边形定义： 
 401 1：循环排序中相邻线段对的交是他们之间共有的单个点 
 402 2：不相邻的线段不相交 
 403 本程序默认第一个条件已经满足 
 404 */ 
 405 bool issimple(int vcount,POINT polygon[]) 
 406 { 
 407     int i,cn; 
 408     LINESEG l1,l2; 
 409     for(i=0;i<vcount;i++) 
 410     { 
 411         l1.s=polygon[i]; 
 412         l1.e=polygon[(i+1)%vcount]; 
 413         cn=vcount-3; 
 414         while(cn) 
 415         { 
 416             l2.s=polygon[(i+2)%vcount]; 
 417             l2.e=polygon[(i+3)%vcount]; 
 418             if(intersect(l1,l2)) 
 419                 break; 
 420             cn--; 
 421         } 
 422         if(cn) 
 423             return false; 
 424     } 
 425     return true; 
 426 } 
 427 // 返回值：按输入顺序返回多边形顶点的凸凹性判断，bc[i]=1,iff:第i个顶点是凸顶点 
 428 void checkconvex(int vcount,POINT polygon[],bool bc[]) 
 429 { 
 430     int i,index=0; 
 431     POINT tp=polygon[0]; 
 432     for(i=1;i<vcount;i++) // 寻找第一个凸顶点 
 433     { 
 434         if(polygon[i].y<tp.y||(polygon[i].y == tp.y&&polygon[i].x<tp.x)) 
 435         { 
 436             tp=polygon[i]; 
 437             index=i; 
 438         } 
 439     } 
 440     int count=vcount-1; 
 441     bc[index]=1; 
 442     while(count) // 判断凸凹性 
 443     { 
 444         if(multiply(polygon[(index+1)%vcount],polygon[(index+2)%vcount],polygon[index])>=0 ) 
 445             bc[(index+1)%vcount]=1; 
 446         else 
 447             bc[(index+1)%vcount]=0; 
 448         index++; 
 449         count--; 
 450     } 
 451 }
 452 // 返回值：多边形polygon是凸多边形时，返回true  
 453 bool isconvex(int vcount,POINT polygon[]) 
 454 { 
 455     bool bc[MAXV]; 
 456     checkconvex(vcount,polygon,bc); 
 457     for(int i=0;i<vcount;i++) // 逐一检查顶点，是否全部是凸顶点 
 458         if(!bc[i]) 
 459             return false; 
 460     return true; 
 461 } 
 462 // 返回多边形面积(signed)；输入顶点按逆时针排列时，返回正值；否则返回负值 
 463 double area_of_polygon(int vcount,POINT polygon[]) 
 464 { 
 465     int i; 
 466     double s; 
 467     if (vcount<3) 
 468         return 0; 
 469     s=polygon[0].y*(polygon[vcount-1].x-polygon[1].x); 
 470     for (i=1;i<vcount;i++) 
 471         s+=polygon[i].y*(polygon[(i-1)].x-polygon[(i+1)%vcount].x); 
 472     return s/2; 
 473 } 
 474 // 如果输入顶点按逆时针排列，返回true 
 475 bool isconterclock(int vcount,POINT polygon[]) 
 476 { 
 477     return area_of_polygon(vcount,polygon)>0; 
 478 } 
 479 // 另一种判断多边形顶点排列方向的方法  
 480 bool isccwize(int vcount,POINT polygon[]) 
 481 { 
 482     int i,index; 
 483     POINT a,b,v; 
 484     v=polygon[0]; 
 485     index=0; 
 486     for(i=1;i<vcount;i++) // 找到最低且最左顶点，肯定是凸顶点 
 487     { 
 488         if(polygon[i].y<v.y||polygon[i].y == v.y && polygon[i].x<v.x) 
 489         { 
 490             index=i; 
 491         } 
 492     } 
 493     a=polygon[(index-1+vcount)%vcount]; // 顶点v的前一顶点 
 494     b=polygon[(index+1)%vcount]; // 顶点v的后一顶点 
 495     return multiply(v,b,a)>0; 
 496 } 
 497 /********************************************************************************************   
 498 射线法判断点q与多边形polygon的位置关系，要求polygon为简单多边形，顶点逆时针排列 
 499    如果点在多边形内：   返回0 
 500    如果点在多边形边上： 返回1 
 501    如果点在多边形外：    返回2 
 502 *********************************************************************************************/ 
 503 int insidepolygon(int vcount,POINT Polygon[],POINT q) 
 504 { 
 505     int c=0,i,n; 
 506     LINESEG l1,l2; 
 507     bool bintersect_a,bonline1,bonline2,bonline3; 
 508     double r1,r2; 
 509  
 510     l1.s=q; 
 511     l1.e=q; 
 512     l1.e.x=double(INF); 
 513     n=vcount; 
 514     for (i=0;i<vcount;i++) 
 515     { 
 516         l2.s=Polygon[i]; 
 517         l2.e=Polygon[(i+1)%n]; 
 518         if(online(l2,q))
 519             return 1; // 如果点在边上，返回1 
 520         if ( (bintersect_a=intersect_A(l1,l2))|| // 相交且不在端点 
 521         ( (bonline1=online(l1,Polygon[(i+1)%n]))&& // 第二个端点在射线上 
 522         ( (!(bonline2=online(l1,Polygon[(i+2)%n])))&& /* 前一个端点和后一个端点在射线两侧 */ 
 523         ((r1=multiply(Polygon[i],Polygon[(i+1)%n],l1.s)*multiply(Polygon[(i+1)%n],Polygon[(i+2)%n],l1.s))>0) ||    
 524         (bonline3=online(l1,Polygon[(i+2)%n]))&&     /* 下一条边是水平线，前一个端点和后一个端点在射线两侧  */ 
 525             ((r2=multiply(Polygon[i],Polygon[(i+2)%n],l1.s)*multiply(Polygon[(i+2)%n], 
 526         Polygon[(i+3)%n],l1.s))>0) 
 527                 ) 
 528             ) 
 529         ) c++; 
 530     } 
 531     if(c%2 == 1) 
 532         return 0; 
 533     else 
 534         return 2; 
 535 } 
 536 //点q是凸多边形polygon内时，返回true；注意：多边形polygon一定要是凸多边形  
 537 bool InsideConvexPolygon(int vcount,POINT polygon[],POINT q) // 可用于三角形！ 
 538 { 
 539     POINT p; 
 540     LINESEG l; 
 541     int i; 
 542     p.x=0;p.y=0; 
 543     for(i=0;i<vcount;i++) // 寻找一个肯定在多边形polygon内的点p：多边形顶点平均值 
 544     { 
 545         p.x+=polygon[i].x; 
 546         p.y+=polygon[i].y; 
 547     } 
 548     p.x /= vcount; 
 549     p.y /= vcount; 
 550  
 551     for(i=0;i<vcount;i++) 
 552     { 
 553         l.s=polygon[i];l.e=polygon[(i+1)%vcount]; 
 554         if(multiply(p,l.e,l.s)*multiply(q,l.e,l.s)<0) /* 点p和点q在边l的两侧，说明点q肯定在多边形外 */ 
 555         break; 
 556     } 
 557     return (i==vcount); 
 558 } 
 559 /********************************************** 
 560 寻找凸包的graham 扫描法 
 561 PointSet为输入的点集； 
 562 ch为输出的凸包上的点集，按照逆时针方向排列; 
 563 n为PointSet中的点的数目 
 564 len为输出的凸包上的点的个数 
 565 **********************************************/ 
 566 void Graham_scan(POINT PointSet[],POINT ch[],int n,int &len) 
 567 { 
 568     int i,j,k=0,top=2; 
 569     POINT tmp; 
 570     // 选取PointSet中y坐标最小的点PointSet[k]，如果这样的点有多个，则取最左边的一个 
 571     for(i=1;i<n;i++) 
 572         if ( PointSet[i].y<PointSet[k].y || (PointSet[i].y==PointSet[k].y) && (PointSet[i].x<PointSet[k].x) ) 
 573             k=i; 
 574     tmp=PointSet[0]; 
 575     PointSet[0]=PointSet[k]; 
 576     PointSet[k]=tmp; // 现在PointSet中y坐标最小的点在PointSet[0] 
 577     for (i=1;i<n-1;i++) /* 对顶点按照相对PointSet[0]的极角从小到大进行排序，极角相同的按照距离PointSet[0]从近到远进行排序 */ 
 578     { 
 579         k=i; 
 580         for (j=i+1;j<n;j++) 
 581             if ( multiply(PointSet[j],PointSet[k],PointSet[0])>0 ||  // 极角更小    
 582                 (multiply(PointSet[j],PointSet[k],PointSet[0])==0) && /* 极角相等，距离更短 */        
 583                 dist(PointSet[0],PointSet[j])<dist(PointSet[0],PointSet[k])
 584                ) 
 585                 k=j; 
 586         tmp=PointSet[i]; 
 587         PointSet[i]=PointSet[k]; 
 588         PointSet[k]=tmp; 
 589     } 
 590     ch[0]=PointSet[0]; 
 591     ch[1]=PointSet[1]; 
 592     ch[2]=PointSet[2]; 
 593     for (i=3;i<n;i++) 
 594     { 
 595         while (multiply(PointSet[i],ch[top],ch[top-1])>=0) 
 596             top--; 
 597         ch[++top]=PointSet[i]; 
 598     } 
 599     len=top+1; 
 600 } 
 601 // 卷包裹法求点集凸壳，参数说明同graham算法    
 602 void ConvexClosure(POINT PointSet[],POINT ch[],int n,int &len) 
 603 { 
 604     int top=0,i,index,first; 
 605     double curmax,curcos,curdis; 
 606     POINT tmp; 
 607     LINESEG l1,l2; 
 608     bool use[MAXV]; 
 609     tmp=PointSet[0]; 
 610     index=0; 
 611     // 选取y最小点，如果多于一个，则选取最左点 
 612     for(i=1;i<n;i++) 
 613     { 
 614         if(PointSet[i].y<tmp.y||PointSet[i].y == tmp.y&&PointSet[i].x<tmp.x) 
 615         { 
 616             index=i; 
 617         } 
 618         use[i]=false; 
 619     } 
 620     tmp=PointSet[index]; 
 621     first=index; 
 622     use[index]=true; 
 623  
 624     index=-1; 
 625     ch[top++]=tmp; 
 626     tmp.x-=100; 
 627     l1.s=tmp; 
 628     l1.e=ch[0]; 
 629     l2.s=ch[0]; 
 630  
 631     while(index!=first) 
 632     { 
 633         curmax=-100; 
 634         curdis=0; 
 635         // 选取与最后一条确定边夹角最小的点，即余弦值最大者 
 636         for(i=0;i<n;i++) 
 637         { 
 638             if(use[i])continue; 
 639             l2.e=PointSet[i]; 
 640             curcos=cosine(l1,l2); // 根据cos值求夹角余弦，范围在 （-1 -- 1 ） 
 641             if(curcos>curmax || fabs(curcos-curmax)<1e-6 && dist(l2.s,l2.e)>curdis) 
 642             { 
 643                 curmax=curcos; 
 644                 index=i; 
 645                 curdis=dist(l2.s,l2.e); 
 646             } 
 647         } 
 648         use[first]=false;            //清空第first个顶点标志，使最后能形成封闭的hull 
 649         use[index]=true; 
 650         ch[top++]=PointSet[index]; 
 651         l1.s=ch[top-2]; 
 652         l1.e=ch[top-1]; 
 653         l2.s=ch[top-1]; 
 654     } 
 655     len=top-1; 
 656 } 
 657 /*********************************************************************************************  
 658     判断线段是否在简单多边形内(注意：如果多边形是凸多边形，下面的算法可以化简) 
 659        必要条件一：线段的两个端点都在多边形内； 
 660     必要条件二：线段和多边形的所有边都不内交； 
 661     用途：    1. 判断折线是否在简单多边形内 
 662             2. 判断简单多边形是否在另一个简单多边形内 
 663 **********************************************************************************************/ 
 664 bool LinesegInsidePolygon(int vcount,POINT polygon[],LINESEG l) 
 665 { 
 666     // 判断线端l的端点是否不都在多边形内 
 667     if(!insidepolygon(vcount,polygon,l.s)||!insidepolygon(vcount,polygon,l.e)) 
 668         return false; 
 669     int top=0,i,j; 
 670     POINT PointSet[MAXV],tmp; 
 671     LINESEG s; 
 672  
 673     for(i=0;i<vcount;i++) 
 674     { 
 675         s.s=polygon[i]; 
 676         s.e=polygon[(i+1)%vcount]; 
 677         if(online(s,l.s)) //线段l的起始端点在线段s上 
 678             PointSet[top++]=l.s; 
 679         else if(online(s,l.e)) //线段l的终止端点在线段s上 
 680             PointSet[top++]=l.e; 
 681         else 
 682         { 
 683             if(online(l,s.s)) //线段s的起始端点在线段l上 
 684                 PointSet[top++]=s.s; 
 685             else if(online(l,s.e)) // 线段s的终止端点在线段l上 
 686                 PointSet[top++]=s.e; 
 687             else 
 688             { 
 689                 if(intersect(l,s)) // 这个时候如果相交，肯定是内交，返回false 
 690                 return false; 
 691             } 
 692         } 
 693     } 
 694  
 695     for(i=0;i<top-1;i++) /* 冒泡排序，x坐标小的排在前面；x坐标相同者，y坐标小的排在前面 */ 
 696     { 
 697         for(j=i+1;j<top;j++) 
 698         { 
 699             if( PointSet[i].x>PointSet[j].x || fabs(PointSet[i].x-PointSet[j].x)<EP && PointSet[i].y>PointSet[j].y ) 
 700             { 
 701                 tmp=PointSet[i]; 
 702                 PointSet[i]=PointSet[j]; 
 703                 PointSet[j]=tmp; 
 704             } 
 705         } 
 706     } 
 707  
 708     for(i=0;i<top-1;i++) 
 709     { 
 710         tmp.x=(PointSet[i].x+PointSet[i+1].x)/2; //得到两个相邻交点的中点 
 711         tmp.y=(PointSet[i].y+PointSet[i+1].y)/2; 
 712         if(!insidepolygon(vcount,polygon,tmp)) 
 713             return false; 
 714     } 
 715     return true; 
 716 } 
 717 /*********************************************************************************************  
 718 求任意简单多边形polygon的重心 
 719 需要调用下面几个函数： 
 720     void AddPosPart(); 增加右边区域的面积 
 721     void AddNegPart(); 增加左边区域的面积 
 722     void AddRegion(); 增加区域面积 
 723 在使用该程序时，如果把xtr,ytr,wtr,xtl,ytl,wtl设成全局变量就可以使这些函数的形式得到化简,
 724 但要注意函数的声明和调用要做相应变化 
 725 **********************************************************************************************/ 
 726 void AddPosPart(double x, double y, double w, double &xtr, double &ytr, double &wtr) 
 727 { 
 728     if (abs(wtr + w)<1e-10 ) return; // detect zero regions 
 729     xtr = ( wtr*xtr + w*x ) / ( wtr + w ); 
 730     ytr = ( wtr*ytr + w*y ) / ( wtr + w ); 
 731     wtr = w + wtr; 
 732     return; 
 733 } 
 734 void AddNegPart(double x, ouble y, double w, double &xtl, double &ytl, double &wtl) 
 735 { 
 736     if ( abs(wtl + w)<1e-10 ) 
 737         return; // detect zero regions 
 738  
 739     xtl = ( wtl*xtl + w*x ) / ( wtl + w ); 
 740     ytl = ( wtl*ytl + w*y ) / ( wtl + w ); 
 741     wtl = w + wtl; 
 742     return; 
 743 } 
 744 void AddRegion ( double x1, double y1, double x2, double y2, double &xtr, double &ytr, 
 745         double &wtr, double &xtl, double &ytl, double &wtl ) 
 746 { 
 747     if ( abs (x1 - x2)< 1e-10 ) 
 748         return; 
 749  
 750     if ( x2 > x1 ) 
 751     { 
 752         AddPosPart ((x2+x1)/2, y1/2, (x2-x1) * y1,xtr,ytr,wtr); /* rectangle 全局变量变化处 */ 
 753         AddPosPart ((x1+x2+x2)/3, (y1+y1+y2)/3, (x2-x1)*(y2-y1)/2,xtr,ytr,wtr);    
 754         // triangle 全局变量变化处 
 755     } 
 756     else 
 757     { 
 758         AddNegPart ((x2+x1)/2, y1/2, (x2-x1) * y1,xtl,ytl,wtl);  
 759         // rectangle 全局变量变化处 
 760         AddNegPart ((x1+x2+x2)/3, (y1+y1+y2)/3, (x2-x1)*(y2-y1)/2,xtl,ytl,wtl);  
 761         // triangle  全局变量变化处 
 762     } 
 763 } 
 764 POINT cg_simple(int vcount,POINT polygon[]) 
 765 { 
 766     double xtr,ytr,wtr,xtl,ytl,wtl;        
 767     //注意： 如果把xtr,ytr,wtr,xtl,ytl,wtl改成全局变量后这里要删去 
 768     POINT p1,p2,tp; 
 769     xtr = ytr = wtr = 0.0; 
 770     xtl = ytl = wtl = 0.0; 
 771     for(int i=0;i<vcount;i++) 
 772     { 
 773         p1=polygon[i]; 
 774         p2=polygon[(i+1)%vcount]; 
 775         AddRegion(p1.x,p1.y,p2.x,p2.y,xtr,ytr,wtr,xtl,ytl,wtl); //全局变量变化处 
 776     } 
 777     tp.x = (wtr*xtr + wtl*xtl) / (wtr + wtl); 
 778     tp.y = (wtr*ytr + wtl*ytl) / (wtr + wtl); 
 779     return tp; 
 780 } 
 781 // 求凸多边形的重心,要求输入多边形按逆时针排序 
 782 POINT gravitycenter(int vcount,POINT polygon[]) 
 783 { 
 784     POINT tp; 
 785     double x,y,s,x0,y0,cs,k; 
 786     x=0;y=0;s=0; 
 787     for(int i=1;i<vcount-1;i++) 
 788     { 
 789         x0=(polygon[0].x+polygon[i].x+polygon[i+1].x)/3; 
 790         y0=(polygon[0].y+polygon[i].y+polygon[i+1].y)/3; //求当前三角形的重心 
 791         cs=multiply(polygon[i],polygon[i+1],polygon[0])/2; 
 792         //三角形面积可以直接利用该公式求解 
 793         if(abs(s)<1e-20) 
 794         { 
 795             x=x0;y=y0;s+=cs;continue; 
 796         } 
 797         k=cs/s; //求面积比例 
 798         x=(x+k*x0)/(1+k); 
 799         y=(y+k*y0)/(1+k); 
 800         s += cs; 
 801     } 
 802     tp.x=x; 
 803     tp.y=y; 
 804     return tp; 
 805 } 
 806  
 807 /************************************************
 808 给定一简单多边形，找出一个肯定在该多边形内的点 
 809 定理1    ：每个多边形至少有一个凸顶点 
 810 定理2    ：顶点数>=4的简单多边形至少有一条对角线 
 811 结论    ： x坐标最大，最小的点肯定是凸顶点 
 812     y坐标最大，最小的点肯定是凸顶点            
 813 ************************************************/ 
 814 POINT a_point_insidepoly(int vcount,POINT polygon[]) 
 815 { 
 816     POINT v,a,b,r; 
 817     int i,index; 
 818     v=polygon[0]; 
 819     index=0; 
 820     for(i=1;i<vcount;i++) //寻找一个凸顶点 
 821     { 
 822         if(polygon[i].y<v.y) 
 823         { 
 824             v=polygon[i]; 
 825             index=i; 
 826         } 
 827     } 
 828     a=polygon[(index-1+vcount)%vcount]; //得到v的前一个顶点 
 829     b=polygon[(index+1)%vcount]; //得到v的后一个顶点 
 830     POINT tri[3],q; 
 831     tri[0]=a;tri[1]=v;tri[2]=b; 
 832     double md=INF; 
 833     int in1=index; 
 834     bool bin=false; 
 835     for(i=0;i<vcount;i++) //寻找在三角形avb内且离顶点v最近的顶点q 
 836     { 
 837         if(i == index)continue; 
 838         if(i == (index-1+vcount)%vcount)continue; 
 839         if(i == (index+1)%vcount)continue; 
 840         if(!InsideConvexPolygon(3,tri,polygon[i]))continue; 
 841         bin=true; 
 842         if(dist(v,polygon[i])<md) 
 843         { 
 844             q=polygon[i]; 
 845             md=dist(v,q); 
 846         } 
 847     } 
 848     if(!bin) //没有顶点在三角形avb内，返回线段ab中点 
 849     { 
 850         r.x=(a.x+b.x)/2; 
 851         r.y=(a.y+b.y)/2; 
 852         return r; 
 853     } 
 854     r.x=(v.x+q.x)/2; //返回线段vq的中点 
 855     r.y=(v.y+q.y)/2; 
 856     return r; 
 857 } 
 858 /***********************************************************************************************
 859 求从多边形外一点p出发到一个简单多边形的切线,如果存在返回切点,其中rp点是右切点,lp是左切点 
 860 注意：p点一定要在多边形外 ,输入顶点序列是逆时针排列 
 861 原 理：    如果点在多边形内肯定无切线;凸多边形有唯一的两个切点,凹多边形就可能有多于两个的切点; 
 862         如果polygon是凸多边形，切点只有两个只要找到就可以,可以化简此算法 
 863         如果是凹多边形还有一种算法可以求解:先求凹多边形的凸包,然后求凸包的切线 
 864 /***********************************************************************************************/ 
 865 void pointtangentpoly(int vcount,POINT polygon[],POINT p,POINT &rp,POINT &lp) 
 866 { 
 867     LINESEG ep,en; 
 868     bool blp,bln; 
 869     rp=polygon[0]; 
 870     lp=polygon[0]; 
 871     for(int i=1;i<vcount;i++) 
 872     { 
 873         ep.s=polygon[(i+vcount-1)%vcount]; 
 874         ep.e=polygon[i]; 
 875         en.s=polygon[i]; 
 876         en.e=polygon[(i+1)%vcount]; 
 877         blp=multiply(ep.e,p,ep.s)>=0;                // p is to the left of pre edge 
 878         bln=multiply(en.e,p,en.s)>=0;                // p is to the left of next edge 
 879         if(!blp&&bln) 
 880         { 
 881             if(multiply(polygon[i],rp,p)>0)           // polygon[i] is above rp 
 882             rp=polygon[i]; 
 883         } 
 884         if(blp&&!bln) 
 885         { 
 886             if(multiply(lp,polygon[i],p)>0)           // polygon[i] is below lp 
 887             lp=polygon[i]; 
 888         } 
 889     } 
 890     return ; 
 891 } 
 892 // 如果多边形polygon的核存在，返回true，返回核上的一点p.顶点按逆时针方向输入  
 893 bool core_exist(int vcount,POINT polygon[],POINT &p) 
 894 { 
 895     int i,j,k; 
 896     LINESEG l; 
 897     LINE lineset[MAXV]; 
 898     for(i=0;i<vcount;i++) 
 899     { 
 900         lineset[i]=makeline(polygon[i],polygon[(i+1)%vcount]); 
 901     } 
 902     for(i=0;i<vcount;i++) 
 903     { 
 904         for(j=0;j<vcount;j++) 
 905         { 
 906             if(i == j)continue; 
 907             if(lineintersect(lineset[i],lineset[j],p)) 
 908             { 
 909                 for(k=0;k<vcount;k++) 
 910                 { 
 911                     l.s=polygon[k]; 
 912                     l.e=polygon[(k+1)%vcount]; 
 913                     if(multiply(p,l.e,l.s)>0)      
 914                     //多边形顶点按逆时针方向排列，核肯定在每条边的左侧或边上 
 915                     break; 
 916                 } 
 917                 if(k == vcount)             //找到了一个核上的点 
 918                 break; 
 919             } 
 920         } 
 921         if(j<vcount) break; 
 922     } 
 923     if(i<vcount) 
 924         return true; 
 925     else 
 926         return false; 
 927 } 
 928 /************************* 
 929 *                         * 
 930 * 圆的基本运算           * 
 931 *                         * 
 932 *************************/ 
 933 /******************************************************************************
 934 返回值    ： 点p在圆内(包括边界)时，返回true 
 935 用途    ： 因为圆为凸集，所以判断点集，折线，多边形是否在圆内时，
 936     只需要逐一判断点是否在圆内即可。 
 937 *******************************************************************************/ 
 938 bool point_in_circle(POINT o,double r,POINT p) 
 939 { 
 940     double d2=(p.x-o.x)*(p.x-o.x)+(p.y-o.y)*(p.y-o.y); 
 941     double r2=r*r; 
 942     return d2<r2||abs(d2-r2)<EP; 
 943 } 
 944 /******************************************************************************
 945 用 途    ：求不共线的三点确定一个圆 
 946 输 入    ：三个点p1,p2,p3 
 947 返回值    ：如果三点共线，返回false；反之，返回true。圆心由q返回，半径由r返回 
 948 *******************************************************************************/ 
 949 bool cocircle(POINT p1,POINT p2,POINT p3,POINT &q,double &r) 
 950 { 
 951     double x12=p2.x-p1.x; 
 952     double y12=p2.y-p1.y; 
 953     double x13=p3.x-p1.x; 
 954     double y13=p3.y-p1.y; 
 955     double z2=x12*(p1.x+p2.x)+y12*(p1.y+p2.y); 
 956     double z3=x13*(p1.x+p3.x)+y13*(p1.y+p3.y); 
 957     double d=2.0*(x12*(p3.y-p2.y)-y12*(p3.x-p2.x)); 
 958     if(abs(d)<EP) //共线，圆不存在 
 959         return false; 
 960     q.x=(y13*z2-y12*z3)/d; 
 961     q.y=(x12*z3-x13*z2)/d; 
 962     r=dist(p1,q); 
 963     return true; 
 964 } 
 965 int line_circle(LINE l,POINT o,double r,POINT &p1,POINT &p2) 
 966 { 
 967     return true; 
 968 } 
 969  
 970 /************************** 
 971 *                          * 
 972 * 矩形的基本运算          * 
 973 *                         * 
 974 **************************/ 
 975 /* 
 976 说明：因为矩形的特殊性，常用算法可以化简： 
 977 1.判断矩形是否包含点 
 978 只要判断该点的横坐标和纵坐标是否夹在矩形的左右边和上下边之间。 
 979 2.判断线段、折线、多边形是否在矩形中 
 980 因为矩形是个凸集，所以只要判断所有端点是否都在矩形中就可以了。 
 981 3.判断圆是否在矩形中 
 982 圆在矩形中的充要条件是：圆心在矩形中且圆的半径小于等于圆心到矩形四边的距离的最小值。 
 983 */ 
 984 // 已知矩形的三个顶点(a,b,c)，计算第四个顶点d的坐标. 注意：已知的三个顶点可以是无序的 
 985 POINT rect4th(POINT a,POINT b,POINT c) 
 986 { 
 987     POINT d; 
 988     if(abs(dotmultiply(a,b,c))<EP) // 说明c点是直角拐角处 
 989     { 
 990         d.x=a.x+b.x-c.x; 
 991         d.y=a.y+b.y-c.y; 
 992     } 
 993     if(abs(dotmultiply(a,c,b))<EP) // 说明b点是直角拐角处 
 994     { 
 995         d.x=a.x+c.x-b.x; 
 996         d.y=a.y+c.y-b.x; 
 997     } 
 998     if(abs(dotmultiply(c,b,a))<EP) // 说明a点是直角拐角处 
 999     { 
1000         d.x=c.x+b.x-a.x; 
1001         d.y=c.y+b.y-a.y; 
1002     } 
1003     return d; 
1004 } 
1005  
1006 /************************* 
1007 *                        * 
1008 * 常用算法的描述        * 
1009 *                        * 
1010 *************************/ 
1011 /* 
1012 尚未实现的算法： 
1013 1. 求包含点集的最小圆 
1014 2. 求多边形的交 
1015 3. 简单多边形的三角剖分 
1016 4. 寻找包含点集的最小矩形 
1017 5. 折线的化简 
1018 6. 判断矩形是否在矩形中 
1019 7. 判断矩形能否放在矩形中 
1020 8. 矩形并的面积与周长 
1021 9. 矩形并的轮廓 
1022 10.矩形并的闭包 
1023 11.矩形的交 
1024 12.点集中的最近点对 
1025 13.多边形的并 
1026 14.圆的交与并 
1027 15.直线与圆的关系 
1028 16.线段与圆的关系 
1029 17.求多边形的核监视摄象机 
1030 18.求点集中不相交点对 railwai 
1031 *//* 
1032 寻找包含点集的最小矩形 
1033 原理：该矩形至少一条边与点集的凸壳的某条边共线 
1034 First take the convex hull of the points. Let the resulting convex 
1035 polygon be P. It has been known for some time that the minimum 
1036 area rectangle enclosing P must have one rectangle side flush with 
1037 (i.e., collinear with and overlapping) one edge of P. This geometric 
1038 fact was used by Godfried Toussaint to develop the "rotating calipers" 
1039 algorithm in a hard-to-find 1983 paper, "Solving Geometric Problems 
1040 with the Rotating Calipers" (Proc. IEEE MELECON). The algorithm 
1041 rotates a surrounding rectangle from one flush edge to the next, 
1042 keeping track of the minimum area for each edge. It achieves O(n) 
1043 time (after hull computation). See the "Rotating Calipers Homepage" 
1044 http://www.cs.mcgill.ca/~orm/rotcal.frame.html for a description 
1045 and applet. 
1046 *//* 
1047 折线的化简 伪码如下： 
1048 Input: tol = the approximation tolerance 
1049 L = {V0,V1,...,Vn-1} is any n-vertex polyline 
1050 Set start = 0; 
1051 Set k = 0; 
1052 Set W0 = V0; 
1053 for each vertex Vi (i=1,n-1) 
1054 { 
1055 if Vi is within tol from Vstart 
1056 then ignore it, and continue with the next vertex 
1057 Vi is further than tol away from Vstart 
1058 so add it as a new vertex of the reduced polyline 
1059 Increment k++; 
1060 Set Wk = Vi; 
1061 Set start = i; as the new initial vertex 
1062 } 
1063 Output: W = {W0,W1,...,Wk-1} = the k-vertex simplified polyline 
1064 */ 
1065 /******************** 
1066 *                    * 
1067 * 补充                * 
1068 *                    * 
1069 ********************/ 
1070  
1071 //两圆关系： 
1072 /* 两圆： 
1073 相离： return 1； 
1074 外切： return 2； 
1075 相交： return 3； 
1076 内切： return 4； 
1077 内含： return 5； 
1078 */ 
1079 int CircleRelation(POINT p1, double r1, POINT p2, double r2) 
1080 { 
1081     double d = sqrt( (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) ); 
1082  
1083     if( fabs(d-r1-r2) < EP ) // 必须保证前两个if先被判定！ 
1084         return 2; 
1085     if( fabs(d-fabs(r1-r2)) < EP ) 
1086         return 4; 
1087     if( d > r1+r2 ) 
1088         return 1; 
1089     if( d < fabs(r1-r2) ) 
1090         return 5; 
1091     if( fabs(r1-r2) < d && d < r1+r2 ) 
1092         return 3; 
1093     return 0; // indicate an error! 
1094 } 
1095 //判断圆是否在矩形内：
1096 // 判定圆是否在矩形内，是就返回true（设矩形水平，且其四个顶点由左上开始按顺时针排列） 
1097 // 调用ptoldist函数，在第4页 
1098 bool CircleRecRelation(POINT pc, double r, POINT pr1, POINT pr2, POINT pr3, POINT pr4) 
1099 { 
1100     if( pr1.x < pc.x && pc.x < pr2.x && pr3.y < pc.y && pc.y < pr2.y ) 
1101     { 
1102         LINESEG line1(pr1, pr2); 
1103         LINESEG line2(pr2, pr3); 
1104         LINESEG line3(pr3, pr4); 
1105         LINESEG line4(pr4, pr1); 
1106         if( r<ptoldist(pc,line1) && r<ptoldist(pc,line2) &&    r<ptoldist(pc,line3) && r<ptoldist(pc,line4) ) 
1107             return true; 
1108     } 
1109     return false; 
1110 } 
1111 //点到平面的距离： 
1112 //点到平面的距离,平面用一般式表示ax+by+cz+d=0 
1113 double P2planeDist(double x, double y, double z, double a, double b, double c, double d) 
1114 { 
1115     return fabs(a*x+b*y+c*z+d) / sqrt(a*a+b*b+c*c); 
1116 } 
1117 //点是否在直线同侧：
1118 //两个点是否在直线同侧，是则返回true 
1119 bool SameSide(POINT p1, POINT p2, LINE line) 
1120 { 
1121     return (line.a * p1.x + line.b * p1.y + line.c) * 
1122     (line.a * p2.x + line.b * p2.y + line.c) > 0; 
1123 }
1124 //镜面反射线：
1125 // 已知入射线、镜面，求反射线。 
1126 // a1,b1,c1为镜面直线方程(a1 x + b1 y + c1 = 0 ,下同)系数;  
1127 //a2,b2,c2为入射光直线方程系数;  
1128 //a,b,c为反射光直线方程系数. 
1129 // 光是有方向的，使用时注意：入射光向量:<-b2,a2>；反射光向量:<b,-a>. 
1130 // 不要忘记结果中可能会有"negative zeros" 
1131 void reflect(double a1,double b1,double c1,double a2,double b2,double c2,double &a,double &b,double &c) 
1132 { 
1133     double n,m; 
1134     double tpb,tpa; 
1135     tpb=b1*b2+a1*a2; 
1136     tpa=a2*b1-a1*b2; 
1137     m=(tpb*b1+tpa*a1)/(b1*b1+a1*a1); 
1138     n=(tpa*b1-tpb*a1)/(b1*b1+a1*a1); 
1139     if(fabs(a1*b2-a2*b1)<1e-20) 
1140     { 
1141         a=a2;b=b2;c=c2; 
1142         return; 
1143     } 
1144     double xx,yy; //(xx,yy)是入射线与镜面的交点。 
1145     xx=(b1*c2-b2*c1)/(a1*b2-a2*b1); 
1146     yy=(a2*c1-a1*c2)/(a1*b2-a2*b1); 
1147     a=n; 
1148     b=-m; 
1149     c=m*yy-xx*n; 
1150 } 
1151 //矩形包含： 
1152 // 矩形2（C，D）是否在1（A，B）内
1153 bool r2inr1(double A,double B,double C,double D)  
1154 { 
1155     double X,Y,L,K,DMax; 
1156     if (A < B) 
1157     { 
1158         double tmp = A; 
1159         A = B; 
1160         B = tmp; 
1161     } 
1162     if (C < D) 
1163     { 
1164         double tmp = C; 
1165         C = D; 
1166         D = tmp; 
1167     } 
1168     if (A > C && B > D)                 // trivial case  
1169         return true; 
1170     else 
1171         if (D >= B) 
1172             return false; 
1173         else 
1174         { 
1175             X = sqrt(A * A + B * B);         // outer rectangle's diagonal  
1176             Y = sqrt(C * C + D * D);         // inner rectangle's diagonal  
1177             if (Y < B) // check for marginal conditions 
1178                 return true; // the inner rectangle can freely rotate inside 
1179             else 
1180                 if (Y > X) 
1181                     return false; 
1182                 else 
1183                 { 
1184                     L = (B - sqrt(Y * Y - A * A)) / 2; 
1185                     K = (A - sqrt(Y * Y - B * B)) / 2; 
1186                     DMax = sqrt(L * L + K * K); 
1187                     if (D >= DMax) 
1188                     return false; 
1189                     else 
1190                     return true; 
1191                 } 
1192         } 
1193 } 
1194 //两圆交点： 
1195 // 两圆已经相交（相切） 
1196 void  c2point(POINT p1,double r1,POINT p2,double r2,POINT &rp1,POINT &rp2) 
1197 { 
1198     double a,b,r; 
1199     a=p2.x-p1.x; 
1200     b=p2.y-p1.y; 
1201     r=(a*a+b*b+r1*r1-r2*r2)/2; 
1202     if(a==0&&b!=0) 
1203     { 
1204         rp1.y=rp2.y=r/b; 
1205         rp1.x=sqrt(r1*r1-rp1.y*rp1.y); 
1206         rp2.x=-rp1.x; 
1207     } 
1208     else if(a!=0&&b==0) 
1209     { 
1210         rp1.x=rp2.x=r/a; 
1211         rp1.y=sqrt(r1*r1-rp1.x*rp2.x); 
1212         rp2.y=-rp1.y; 
1213     } 
1214     else if(a!=0&&b!=0) 
1215     { 
1216         double delta; 
1217         delta=b*b*r*r-(a*a+b*b)*(r*r-r1*r1*a*a); 
1218         rp1.y=(b*r+sqrt(delta))/(a*a+b*b); 
1219         rp2.y=(b*r-sqrt(delta))/(a*a+b*b); 
1220         rp1.x=(r-b*rp1.y)/a; 
1221         rp2.x=(r-b*rp2.y)/a; 
1222     } 
1223  
1224     rp1.x+=p1.x; 
1225     rp1.y+=p1.y; 
1226     rp2.x+=p1.x; 
1227     rp2.y+=p1.y; 
1228 } 
1229 //两圆公共面积：
1230 // 必须保证相交 
1231 double c2area(POINT p1,double r1,POINT p2,double r2) 
1232 { 
1233     POINT rp1,rp2; 
1234     c2point(p1,r1,p2,r2,rp1,rp2); 
1235  
1236     if(r1>r2) //保证r2>r1 
1237     { 
1238         swap(p1,p2); 
1239         swap(r1,r2); 
1240     } 
1241     double a,b,rr; 
1242     a=p1.x-p2.x; 
1243     b=p1.y-p2.y; 
1244     rr=sqrt(a*a+b*b); 
1245  
1246     double dx1,dy1,dx2,dy2; 
1247     double sita1,sita2; 
1248     dx1=rp1.x-p1.x; 
1249     dy1=rp1.y-p1.y; 
1250     dx2=rp2.x-p1.x; 
1251     dy2=rp2.y-p1.y; 
1252     sita1=acos((dx1*dx2+dy1*dy2)/r1/r1); 
1253  
1254     dx1=rp1.x-p2.x; 
1255     dy1=rp1.y-p2.y; 
1256     dx2=rp2.x-p2.x; 
1257     dy2=rp2.y-p2.y; 
1258     sita2=acos((dx1*dx2+dy1*dy2)/r2/r2); 
1259     double s=0; 
1260     if(rr<r2)//相交弧为优弧 
1261         s=r1*r1*(PI-sita1/2+sin(sita1)/2)+r2*r2*(sita2-sin(sita2))/2; 
1262     else//相交弧为劣弧 
1263         s=(r1*r1*(sita1-sin(sita1))+r2*r2*(sita2-sin(sita2)))/2; 
1264  
1265     return s; 
1266 } 
1267 //圆和直线关系： 
1268 //0----相离 1----相切 2----相交 
1269 int clpoint(POINT p,double r,double a,double b,double c,POINT &rp1,POINT &rp2) 
1270 { 
1271     int res=0; 
1272  
1273     c=c+a*p.x+b*p.y; 
1274     double tmp; 
1275     if(a==0&&b!=0) 
1276     { 
1277         tmp=-c/b; 
1278         if(r*r<tmp*tmp) 
1279             res=0; 
1280         else if(r*r==tmp*tmp) 
1281         { 
1282             res=1; 
1283             rp1.y=tmp; 
1284             rp1.x=0; 
1285         } 
1286         else 
1287         { 
1288             res=2; 
1289             rp1.y=rp2.y=tmp; 
1290             rp1.x=sqrt(r*r-tmp*tmp); 
1291             rp2.x=-rp1.x; 
1292         } 
1293     } 
1294     else if(a!=0&&b==0) 
1295     { 
1296         tmp=-c/a; 
1297         if(r*r<tmp*tmp) 
1298             res=0; 
1299         else if(r*r==tmp*tmp) 
1300         { 
1301             res=1; 
1302             rp1.x=tmp; 
1303             rp1.y=0; 
1304         } 
1305         else 
1306         { 
1307             res=2; 
1308             rp1.x=rp2.x=tmp; 
1309             rp1.y=sqrt(r*r-tmp*tmp); 
1310             rp2.y=-rp1.y; 
1311         } 
1312     } 
1313     else if(a!=0&&b!=0) 
1314     { 
1315         double delta; 
1316         delta=b*b*c*c-(a*a+b*b)*(c*c-a*a*r*r); 
1317         if(delta<0) 
1318             res=0; 
1319         else if(delta==0) 
1320         { 
1321             res=1; 
1322             rp1.y=-b*c/(a*a+b*b); 
1323             rp1.x=(-c-b*rp1.y)/a; 
1324         } 
1325         else 
1326         { 
1327             res=2; 
1328             rp1.y=(-b*c+sqrt(delta))/(a*a+b*b); 
1329             rp2.y=(-b*c-sqrt(delta))/(a*a+b*b); 
1330             rp1.x=(-c-b*rp1.y)/a; 
1331             rp2.x=(-c-b*rp2.y)/a; 
1332         } 
1333     } 
1334     rp1.x+=p.x; 
1335     rp1.y+=p.y; 
1336     rp2.x+=p.x; 
1337     rp2.y+=p.y; 
1338     return res; 
1339 } 
1340 //内切圆： 
1341 void incircle(POINT p1,POINT p2,POINT p3,POINT &rp,double &r) 
1342 { 
1343     double dx31,dy31,dx21,dy21,d31,d21,a1,b1,c1; 
1344     dx31=p3.x-p1.x; 
1345     dy31=p3.y-p1.y; 
1346     dx21=p2.x-p1.x; 
1347     dy21=p2.y-p1.y; 
1348  
1349     d31=sqrt(dx31*dx31+dy31*dy31); 
1350     d21=sqrt(dx21*dx21+dy21*dy21); 
1351     a1=dx31*d21-dx21*d31; 
1352     b1=dy31*d21-dy21*d31; 
1353     c1=a1*p1.x+b1*p1.y; 
1354  
1355     double dx32,dy32,dx12,dy12,d32,d12,a2,b2,c2; 
1356     dx32=p3.x-p2.x; 
1357     dy32=p3.y-p2.y; 
1358     dx12=-dx21; 
1359     dy12=-dy21; 
1360  
1361     d32=sqrt(dx32*dx32+dy32*dy32); 
1362     d12=d21; 
1363     a2=dx12*d32-dx32*d12; 
1364     b2=dy12*d32-dy32*d12; 
1365     c2=a2*p2.x+b2*p2.y; 
1366  
1367     rp.x=(c1*b2-c2*b1)/(a1*b2-a2*b1); 
1368     rp.y=(c2*a1-c1*a2)/(a1*b2-a2*b1); 
1369     r=fabs(dy21*rp.x-dx21*rp.y+dx21*p1.y-dy21*p1.x)/d21; 
1370 } 
1371 //求切点： 
1372 // p---圆心坐标， r---圆半径， sp---圆外一点， rp1,rp2---切点坐标 
1373 void cutpoint(POINT p,double r,POINT sp,POINT &rp1,POINT &rp2) 
1374 { 
1375     POINT p2; 
1376     p2.x=(p.x+sp.x)/2; 
1377     p2.y=(p.y+sp.y)/2; 
1378  
1379     double dx2,dy2,r2; 
1380     dx2=p2.x-p.x; 
1381     dy2=p2.y-p.y; 
1382     r2=sqrt(dx2*dx2+dy2*dy2); 
1383     c2point(p,r,p2,r2,rp1,rp2); 
1384 } 
1385 //线段的左右旋： 
1386 /* l2在l1的左/右方向（l1为基准线） 
1387 返回    0    ： 重合； 
1388 返回    1    ： 右旋； 
1389 返回    –1 ： 左旋； 
1390 */ 
1391 int rotat(LINESEG l1,LINESEG l2) 
1392 { 
1393     double dx1,dx2,dy1,dy2; 
1394     dx1=l1.s.x-l1.e.x; 
1395     dy1=l1.s.y-l1.e.y; 
1396     dx2=l2.s.x-l2.e.x; 
1397     dy2=l2.s.y-l2.e.y; 
1398  
1399     double d; 
1400     d=dx1*dy2-dx2*dy1; 
1401     if(d==0) 
1402         return 0; 
1403     else if(d>0) 
1404         return -1; 
1405     else 
1406         return 1; 
1407 } 
1408 /*
1409 公式： 
1410 球坐标公式： 
1411 直角坐标为 P(x, y, z) 时，对应的球坐标是(rsinφcosθ, rsinφsinθ, rcosφ),其中φ是向量OP与Z轴的夹角，范围[0，π]；是OP在XOY面上的投影到X轴的旋角，范围[0，2π]  
1412 直线的一般方程转化成向量方程： 
1413 ax+by+c=0 
1414 x-x0     y-y0 
1415    ------ = ------- // (x0,y0)为直线上一点，m,n为向量 
1416 m        n 
1417 转换关系： 
1418 a=n；b=-m；c=m·y0-n·x0； 
1419 m=-b; n=a; 
1420 三点平面方程： 
1421 三点为P1，P2，P3 
1422 设向量  M1=P2-P1; M2=P3-P1; 
1423 平面法向量：  M=M1 x M2 （） 
1424 平面方程：    M.i(x-P1.x)+M.j(y-P1.y)+M.k(z-P1.z)=0 
1425 */