計算幾何：模板

叉積：兩個向量的叉積是一個標量，a×b的幾何意義為他們所形成的平行四邊形的有向面積

凸包：把給定點包圍在內部，面積最小的凸多邊形

半平面交：一個半平面指的是由滿足ax+by+c>或ax+by+c>=0的點集組成的二維區域。

一般來說在寫程式碼的時候，我們可以把一個半平面想象成一個向量所在的直線右面的一片區域 

半平面的交可能是一個封閉圖形也可能是沒有邊界的區域或者為空

對踵點對：如果兩個點 p 和 q （屬於 P） 在兩條平行切線上， 那麼他們就形成了一個對踵點對

旋轉卡殼： Shamos提出了一個簡單的 O(n) 時間的演算法來確定一個凸 n 角形的對踵點對

Shamos的演算法就像繞著多邊形旋轉一對卡殼。 因此就有了術語“旋轉卡殼”

pick公式：頂點座標均是整點的簡單多邊形：面積=內部格點數目+邊上格點數目/2-1

/*
計算幾何
目錄 
㈠ 點的基本運算 
1. 平面上兩點之間距離 1 
2. 判斷兩點是否重合 1 
3. 向量叉乘 1 
4. 向量點乘 2 
5. 判斷點是否線上段上 2 
6. 求一點饒某點旋轉後的座標 2 
7. 求向量夾角 2 
㈡ 線段及直線的基本運算 
1. 點與線段的關係 3 
2. 求點到線段所在直線垂線的垂足 4 
3. 點到線段的最近點 4 
4. 點到線段所在直線的距離 4 
5. 點到折線集的最近距離 4 
6. 判斷圓是否在多邊形內 5 
7. 求向量夾角餘弦 5 
8. 求線段之間的夾角 5 
9. 判斷線段是否相交 6 
10.判斷線段是否相交但不交在端點處 6 
11.求線段所在直線的方程 6 
12.求直線的斜率 7 
13.求直線的傾斜角 7 
14.求點關於某直線的對稱點 7 
15.判斷兩條直線是否相交及求直線交點 7 
16.判斷線段是否相交，如果相交返回交點 7 
㈢ 多邊形常用演算法模組 
1. 判斷多邊形是否簡單多邊形 8 
2. 檢查多邊形頂點的凸凹性 9 
3. 判斷多邊形是否凸多邊形 9 
4. 求多邊形面積 9 
5. 判斷多邊形頂點的排列方向，方法一 10 
6. 判斷多邊形頂點的排列方向，方法二 10 
7. 射線法判斷點是否在多邊形內 10 
8. 判斷點是否在凸多邊形內 11 
9. 尋找點集的graham演算法 12 
10.尋找點集凸包的捲包裹法 13 
11.判斷線段是否在多邊形內 14 
12.求簡單多邊形的重心 15 
13.求凸多邊形的重心 17 
14.求肯定在給定多邊形內的一個點 17 
15.求從多邊形外一點出發到該多邊形的切線 18 
16.判斷多邊形的核是否存在 19 
㈣ 圓的基本運算 
1 .點是否在圓內 20 
2 .求不共線的三點所確定的圓 21 
㈤ 矩形的基本運算 
1.已知矩形三點座標，求第4點座標 22 
㈥ 常用演算法的描述 22 
㈦ 補充 
1．兩圓關係： 24 
2．判斷圓是否在矩形內： 24 
3．點到平面的距離： 25 
4．點是否在直線同側： 25 
5．鏡面反射線： 25 
6．矩形包含： 26 
7．兩圓交點： 27 
8．兩圓公共面積： 28 
9. 圓和直線關係： 29 
10. 內切圓： 30 
11. 求切點： 31 
12. 線段的左右旋： 31 
13．公式： 32 
*/
/* 需要包含的標頭檔案 */ 
#include <cmath > 
 
/* 常用的常量定義 */ 
const double    INF        = 1E200    
const double    EP        = 1E-10 
const int        MAXV    = 300 
const double    PI        = 3.14159265 
 
/* 基本幾何結構 */ 
struct POINT 
{ 
    double x; 
    double y; 
    POINT(double a=0, double b=0) { x=a; y=b;} //constructor 
}; 
struct LINESEG 
{ 
    POINT s; 
    POINT e; 
    LINESEG(POINT a, POINT b) { s=a; e=b;} 
    LINESEG() { } 
}; 
struct LINE           // 直線的解析方程 a*x+b*y+c=0  為統一表示，約定 a >= 0 
{ 
   double a; 
   double b; 
   double c; 
   LINE(double d1=1, double d2=-1, double d3=0) {a=d1; b=d2; c=d3;} 
}; 
 
/**********************
 *                    * 
 *   點的基本運算     * 
 *                    * 
 **********************/ 
 
double dist(POINT p1,POINT p2)                // 返回兩點之間歐氏距離 
{ 
    return( sqrt( (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) ) ); 
} 
bool equal_point(POINT p1,POINT p2)           // 判斷兩個點是否重合  
{ 
    return ( (abs(p1.x-p2.x)<EP)&&(abs(p1.y-p2.y)<EP) ); 
} 
/****************************************************************************** 
r=multiply(sp,ep,op),得到(sp-op)和(ep-op)的叉積 
r>0：ep在向量opsp的逆時針方向； 
r=0：opspep三點共線； 
r<0：ep在向量opsp的順時針方向 
*******************************************************************************/ 
double multiply(POINT sp,POINT ep,POINT op) 
{ 
    return((sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y)); 
} 
/* 
r=dotmultiply(p1,p2,op),得到向量(p1-op)和(p2-op)的點積，如果兩個向量都非零向量 
r<0：兩向量夾角為銳角；
r=0：兩向量夾角為直角；
r>0：兩向量夾角為鈍角 
*******************************************************************************/ 
double dotmultiply(POINT p1,POINT p2,POINT p0) 
{ 
    return ((p1.x-p0.x)*(p2.x-p0.x)+(p1.y-p0.y)*(p2.y-p0.y)); 
} 
/****************************************************************************** 
判斷點p是否線上段l上
條件：(p線上段l所在的直線上) && (點p在以線段l為對角線的矩形內)
*******************************************************************************/ 
bool online(LINESEG l,POINT p) 
{ 
    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 ) ) ); 
} 
// 返回點p以點o為圓心逆時針旋轉alpha(單位：弧度)後所在的位置 
POINT rotate(POINT o,double alpha,POINT p) 
{ 
    POINT tp; 
    p.x-=o.x; 
    p.y-=o.y; 
    tp.x=p.x*cos(alpha)-p.y*sin(alpha)+o.x; 
    tp.y=p.y*cos(alpha)+p.x*sin(alpha)+o.y; 
    return tp; 
} 
/* 返回頂角在o點，起始邊為os，終止邊為oe的夾角(單位：弧度) 
    角度小於pi，返回正值 
    角度大於pi，返回負值 
    可以用於求線段之間的夾角 
原理：
    r = dotmultiply(s,e,o) / (dist(o,s)*dist(o,e))
    r'= multiply(s,e,o)
    r >= 1    angle = 0;
    r <= -1    angle = -PI
    -1<r<1 && r'>0    angle = arccos(r)
    -1<r<1 && r'<=0    angle = -arccos(r)
*/ 
double angle(POINT o,POINT s,POINT e) 
{ 
    double cosfi,fi,norm; 
    double dsx = s.x - o.x; 
    double dsy = s.y - o.y; 
    double dex = e.x - o.x; 
    double dey = e.y - o.y; 
 
    cosfi=dsx*dex+dsy*dey; 
    norm=(dsx*dsx+dsy*dsy)*(dex*dex+dey*dey); 
    cosfi /= sqrt( norm ); 
 
    if (cosfi >=  1.0 ) return 0; 
    if (cosfi <= -1.0 ) return -3.1415926; 
 
    fi=acos(cosfi); 
    if (dsx*dey-dsy*dex>0) return fi;      // 說明向量os 在向量 oe的順時針方向 
    return -fi; 
} 
  /*****************************\ 
  *                             * 
  *      線段及直線的基本運算   * 
  *                             * 
  \*****************************/ 
 
/* 判斷點與線段的關係,用途很廣泛 
本函式是根據下面的公式寫的，P是點C到線段AB所在直線的垂足 
                AC dot AB 
        r =     --------- 
                 ||AB||^2 
             (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay) 
          = ------------------------------- 
                          L^2 
    r has the following meaning: 
        r=0      P = A 
        r=1      P = B 
        r<0         P is on the backward extension of AB 
        r>1      P is on the forward extension of AB 
        0<r<1     P is interior to AB 
*/ 
double relation(POINT p,LINESEG l) 
{ 
    LINESEG tl; 
    tl.s=l.s; 
    tl.e=p; 
    return dotmultiply(tl.e,l.e,l.s)/(dist(l.s,l.e)*dist(l.s,l.e)); 
} 
// 求點C到線段AB所在直線的垂足 P 
POINT perpendicular(POINT p,LINESEG l) 
{ 
    double r=relation(p,l); 
    POINT tp; 
    tp.x=l.s.x+r*(l.e.x-l.s.x); 
    tp.y=l.s.y+r*(l.e.y-l.s.y); 
    return tp; 
} 
/* 求點p到線段l的最短距離,並返回線段上距該點最近的點np 
注意：np是線段l上到點p最近的點，不一定是垂足 */ 
double ptolinesegdist(POINT p,LINESEG l,POINT &np) 
{ 
    double r=relation(p,l); 
    if(r<0) 
    { 
        np=l.s; 
        return dist(p,l.s); 
    } 
    if(r>1) 
    { 
        np=l.e; 
        return dist(p,l.e); 
    } 
    np=perpendicular(p,l); 
    return dist(p,np); 
} 
// 求點p到線段l所在直線的距離,請注意本函式與上個函式的區別  
double ptoldist(POINT p,LINESEG l) 
{ 
    return abs(multiply(p,l.e,l.s))/dist(l.s,l.e); 
} 
/* 計算點到折線集的最近距離,並返回最近點. 
注意：呼叫的是ptolineseg()函式 */ 
double ptopointset(int vcount,POINT pointset[],POINT p,POINT &q) 
{ 
    int i; 
    double cd=double(INF),td; 
    LINESEG l; 
    POINT tq,cq; 
 
    for(i=0;i<vcount-1;i++) 
    { 
        l.s=pointset[i]; 
 
        l.e=pointset[i+1]; 
        td=ptolinesegdist(p,l,tq); 
        if(td<cd) 
        { 
            cd=td; 
            cq=tq; 
        } 
    } 
    q=cq; 
    return cd; 
} 
/* 判斷圓是否在多邊形內.ptolineseg()函式的應用2 */ 
bool CircleInsidePolygon(int vcount,POINT center,double radius,POINT polygon[]) 
{ 
    POINT q; 
    double d; 
    q.x=0; 
    q.y=0; 
    d=ptopointset(vcount,polygon,center,q); 
    if(d<radius||fabs(d-radius)<EP) 
        return true; 
    else 
        return false; 
} 
/* 返回兩個向量l1和l2的夾角的餘弦(-1 --- 1)注意：如果想從餘弦求夾角的話，注意反餘弦函式的定義域是從 0到pi */ 
double cosine(LINESEG l1,LINESEG l2) 
{ 
    return (((l1.e.x-l1.s.x)*(l2.e.x-l2.s.x) + 
    (l1.e.y-l1.s.y)*(l2.e.y-l2.s.y))/(dist(l1.e,l1.s)*dist(l2.e,l2.s))) ); 
} 
// 返回線段l1與l2之間的夾角 單位：弧度 範圍(-pi，pi) 
double lsangle(LINESEG l1,LINESEG l2) 
{ 
    POINT o,s,e; 
    o.x=o.y=0; 
    s.x=l1.e.x-l1.s.x; 
    s.y=l1.e.y-l1.s.y; 
    e.x=l2.e.x-l2.s.x; 
    e.y=l2.e.y-l2.s.y; 
    return angle(o,s,e); 
} 
// 如果線段u和v相交(包括相交在端點處)時，返回true 
//
//判斷P1P2跨立Q1Q2的依據是：( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) >= 0。
//判斷Q1Q2跨立P1P2的依據是：( Q1 - P1 ) × ( P2 - P1 ) * ( P2 - P1 ) × ( Q2 - P1 ) >= 0。
bool intersect(LINESEG u,LINESEG v) 
{ 
    return( (max(u.s.x,u.e.x)>=min(v.s.x,v.e.x))&&                     //排斥實驗 
            (max(v.s.x,v.e.x)>=min(u.s.x,u.e.x))&& 
            (max(u.s.y,u.e.y)>=min(v.s.y,v.e.y))&& 
            (max(v.s.y,v.e.y)>=min(u.s.y,u.e.y))&& 
            (multiply(v.s,u.e,u.s)*multiply(u.e,v.e,u.s)>=0)&&         //跨立實驗 
            (multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0)); 
} 
//  (線段u和v相交)&&(交點不是雙方的端點) 時返回true    
bool intersect_A(LINESEG u,LINESEG v) 
{ 
    return    ((intersect(u,v))&& 
            (!online(u,v.s))&& 
            (!online(u,v.e))&& 
            (!online(v,u.e))&& 
            (!online(v,u.s))); 
} 
// 線段v所在直線與線段u相交時返回true；方法：判斷線段u是否跨立線段v  
bool intersect_l(LINESEG u,LINESEG v) 
{ 
    return multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0; 
} 
// 根據已知兩點座標，求過這兩點的直線解析方程： a*x+b*y+c = 0  (a >= 0)  
LINE makeline(POINT p1,POINT p2) 
{ 
    LINE tl; 
    int sign = 1; 
    tl.a=p2.y-p1.y; 
    if(tl.a<0) 
    { 
        sign = -1; 
        tl.a=sign*tl.a; 
    } 
    tl.b=sign*(p1.x-p2.x); 
    tl.c=sign*(p1.y*p2.x-p1.x*p2.y); 
    return tl; 
} 
// 根據直線解析方程返回直線的斜率k,水平線返回 0,豎直線返回 1e200 
double slope(LINE l) 
{ 
    if(abs(l.a) < 1e-20)
        return 0; 
    if(abs(l.b) < 1e-20)
        return INF; 
    return -(l.a/l.b); 
} 
// 返回直線的傾斜角alpha ( 0 - pi) 
double alpha(LINE l) 
{ 
    if(abs(l.a)< EP)
        return 0; 
    if(abs(l.b)< EP)
        return PI/2; 
    double k=slope(l); 
    if(k>0) 
        return atan(k); 
    else 
        return PI+atan(k); 
} 
// 求點p關於直線l的對稱點  
POINT symmetry(LINE l,POINT p) 
{ 
   POINT tp; 
   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); 
   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); 
   return tp; 
} 
// 如果兩條直線 l1(a1*x+b1*y+c1 = 0), l2(a2*x+b2*y+c2 = 0)相交，返回true，且返回交點p  
bool lineintersect(LINE l1,LINE l2,POINT &p) // 是 L1，L2 
{ 
    double d=l1.a*l2.b-l2.a*l1.b; 
    if(abs(d)<EP) // 不相交 
        return false; 
    p.x = (l2.c*l1.b-l1.c*l2.b)/d; 
    p.y = (l2.a*l1.c-l1.a*l2.c)/d; 
    return true; 
} 
// 如果線段l1和l2相交，返回true且交點由(inter)返回，否則返回false 
bool intersection(LINESEG l1,LINESEG l2,POINT &inter) 
{ 
    LINE ll1,ll2; 
    ll1=makeline(l1.s,l1.e); 
    ll2=makeline(l2.s,l2.e); 
    if(lineintersect(ll1,ll2,inter)) 
        return online(l1,inter); 
    else 
        return false; 
} 
 
/******************************\ 
*                              * 
* 多邊形常用演算法模組          * 
*                              * 
\******************************/ 
 
// 如果無特別說明，輸入多邊形頂點要求按逆時針排列 
 
/* 
返回值：輸入的多邊形是簡單多邊形，返回true 
要 求：輸入頂點序列按逆時針排序 
說 明：簡單多邊形定義： 
1：迴圈排序中相鄰線段對的交是他們之間共有的單個點 
2：不相鄰的線段不相交 
本程式預設第一個條件已經滿足 
*/ 
bool issimple(int vcount,POINT polygon[]) 
{ 
    int i,cn; 
    LINESEG l1,l2; 
    for(i=0;i<vcount;i++) 
    { 
        l1.s=polygon[i]; 
        l1.e=polygon[(i+1)%vcount]; 
        cn=vcount-3; 
        while(cn) 
        { 
            l2.s=polygon[(i+2)%vcount]; 
            l2.e=polygon[(i+3)%vcount]; 
            if(intersect(l1,l2)) 
                break; 
            cn--; 
        } 
        if(cn) 
            return false; 
    } 
    return true; 
} 
// 返回值：按輸入順序返回多邊形頂點的凸凹性判斷，bc[i]=1,iff:第i個頂點是凸頂點 
void checkconvex(int vcount,POINT polygon[],bool bc[]) 
{ 
    int i,index=0; 
    POINT tp=polygon[0]; 
    for(i=1;i<vcount;i++) // 尋找第一個凸頂點 
    { 
        if(polygon[i].y<tp.y||(polygon[i].y == tp.y&&polygon[i].x<tp.x)) 
        { 
            tp=polygon[i]; 
            index=i; 
        } 
    } 
    int count=vcount-1; 
    bc[index]=1; 
    while(count) // 判斷凸凹性 
    { 
        if(multiply(polygon[(index+1)%vcount],polygon[(index+2)%vcount],polygon[index])>=0 ) 
            bc[(index+1)%vcount]=1; 
        else 
            bc[(index+1)%vcount]=0; 
        index++; 
        count--; 
    } 
}
// 返回值：多邊形polygon是凸多邊形時，返回true  
bool isconvex(int vcount,POINT polygon[]) 
{ 
    bool bc[MAXV]; 
    checkconvex(vcount,polygon,bc); 
    for(int i=0;i<vcount;i++) // 逐一檢查頂點，是否全部是凸頂點 
        if(!bc[i]) 
            return false; 
    return true; 
} 
// 返回多邊形面積(signed)；輸入頂點按逆時針排列時，返回正值；否則返回負值 
double area_of_polygon(int vcount,POINT polygon[]) 
{ 
    int i; 
    double s; 
    if (vcount<3) 
        return 0; 
    s=polygon[0].y*(polygon[vcount-1].x-polygon[1].x); 
    for (i=1;i<vcount;i++) 
        s+=polygon[i].y*(polygon[(i-1)].x-polygon[(i+1)%vcount].x); 
    return s/2; 
} 
// 如果輸入頂點按逆時針排列，返回true 
bool isconterclock(int vcount,POINT polygon[]) 
{ 
    return area_of_polygon(vcount,polygon)>0; 
} 
// 另一種判斷多邊形頂點排列方向的方法  
bool isccwize(int vcount,POINT polygon[]) 
{ 
    int i,index; 
    POINT a,b,v; 
    v=polygon[0]; 
    index=0; 
    for(i=1;i<vcount;i++) // 找到最低且最左頂點，肯定是凸頂點 
    { 
        if(polygon[i].y<v.y||polygon[i].y == v.y && polygon[i].x<v.x) 
        { 
            index=i; 
        } 
    } 
    a=polygon[(index-1+vcount)%vcount]; // 頂點v的前一頂點 
    b=polygon[(index+1)%vcount]; // 頂點v的後一頂點 
    return multiply(v,b,a)>0; 
} 
/********************************************************************************************   
射線法判斷點q與多邊形polygon的位置關係，要求polygon為簡單多邊形，頂點逆時針排列 
   如果點在多邊形內：   返回0 
   如果點在多邊形邊上： 返回1 
   如果點在多邊形外：    返回2 
*********************************************************************************************/ 
int insidepolygon(int vcount,POINT Polygon[],POINT q) 
{ 
    int c=0,i,n; 
    LINESEG l1,l2; 
    bool bintersect_a,bonline1,bonline2,bonline3; 
    double r1,r2; 
 
    l1.s=q; 
    l1.e=q; 
    l1.e.x=double(INF); 
    n=vcount; 
    for (i=0;i<vcount;i++) 
    { 
        l2.s=Polygon[i]; 
        l2.e=Polygon[(i+1)%n]; 
        if(online(l2,q))
            return 1; // 如果點在邊上，返回1 
        if ( (bintersect_a=intersect_A(l1,l2))|| // 相交且不在端點 
        ( (bonline1=online(l1,Polygon[(i+1)%n]))&& // 第二個端點在射線上 
        ( (!(bonline2=online(l1,Polygon[(i+2)%n])))&& /* 前一個端點和後一個端點在射線兩側 */ 
        ((r1=multiply(Polygon[i],Polygon[(i+1)%n],l1.s)*multiply(Polygon[(i+1)%n],Polygon[(i+2)%n],l1.s))>0) ||    
        (bonline3=online(l1,Polygon[(i+2)%n]))&&     /* 下一條邊是水平線，前一個端點和後一個端點在射線兩側  */ 
            ((r2=multiply(Polygon[i],Polygon[(i+2)%n],l1.s)*multiply(Polygon[(i+2)%n], 
        Polygon[(i+3)%n],l1.s))>0) 
                ) 
            ) 
        ) c++; 
    } 
    if(c%2 == 1) 
        return 0; 
    else 
        return 2; 
} 
//點q是凸多邊形polygon內時，返回true；注意：多邊形polygon一定要是凸多邊形  
bool InsideConvexPolygon(int vcount,POINT polygon[],POINT q) // 可用於三角形！ 
{ 
    POINT p; 
    LINESEG l; 
    int i; 
    p.x=0;p.y=0; 
    for(i=0;i<vcount;i++) // 尋找一個肯定在多邊形polygon內的點p：多邊形頂點平均值 
    { 
        p.x+=polygon[i].x; 
        p.y+=polygon[i].y; 
    } 
    p.x /= vcount; 
    p.y /= vcount; 
 
    for(i=0;i<vcount;i++) 
    { 
        l.s=polygon[i];l.e=polygon[(i+1)%vcount]; 
        if(multiply(p,l.e,l.s)*multiply(q,l.e,l.s)<0) /* 點p和點q在邊l的兩側，說明點q肯定在多邊形外 */ 
        break; 
    } 
    return (i==vcount); 
} 
/********************************************** 
尋找凸包的graham 掃描法 
PointSet為輸入的點集； 
ch為輸出的凸包上的點集，按照逆時針方向排列; 
n為PointSet中的點的數目 
len為輸出的凸包上的點的個數 
**********************************************/ 
void Graham_scan(POINT PointSet[],POINT ch[],int n,int &len) 
{ 
    int i,j,k=0,top=2; 
    POINT tmp; 
    // 選取PointSet中y座標最小的點PointSet[k]，如果這樣的點有多個，則取最左邊的一個 
    for(i=1;i<n;i++) 
        if ( PointSet[i].y<PointSet[k].y || (PointSet[i].y==PointSet[k].y) && (PointSet[i].x<PointSet[k].x) ) 
            k=i; 
    tmp=PointSet[0]; 
    PointSet[0]=PointSet[k]; 
    PointSet[k]=tmp; // 現在PointSet中y座標最小的點在PointSet[0] 
    for (i=1;i<n-1;i++) /* 對頂點按照相對PointSet[0]的極角從小到大進行排序，極角相同的按照距離PointSet[0]從近到遠進行排序 */ 
    { 
        k=i; 
        for (j=i+1;j<n;j++) 
            if ( multiply(PointSet[j],PointSet[k],PointSet[0])>0 ||  // 極角更小    
                (multiply(PointSet[j],PointSet[k],PointSet[0])==0) && /* 極角相等，距離更短 */        
                dist(PointSet[0],PointSet[j])<dist(PointSet[0],PointSet[k])
               ) 
                k=j; 
        tmp=PointSet[i]; 
        PointSet[i]=PointSet[k]; 
        PointSet[k]=tmp; 
    } 
    ch[0]=PointSet[0]; 
    ch[1]=PointSet[1]; 
    ch[2]=PointSet[2]; 
    for (i=3;i<n;i++) 
    { 
        while (multiply(PointSet[i],ch[top],ch[top-1])>=0) 
            top--; 
        ch[++top]=PointSet[i]; 
    } 
    len=top+1; 
} 
// 捲包裹法求點集凸殼，引數說明同graham演算法    
void ConvexClosure(POINT PointSet[],POINT ch[],int n,int &len) 
{ 
    int top=0,i,index,first; 
    double curmax,curcos,curdis; 
    POINT tmp; 
    LINESEG l1,l2; 
    bool use[MAXV]; 
    tmp=PointSet[0]; 
    index=0; 
    // 選取y最小點，如果多於一個，則選取最左點 
    for(i=1;i<n;i++) 
    { 
        if(PointSet[i].y<tmp.y||PointSet[i].y == tmp.y&&PointSet[i].x<tmp.x) 
        { 
            index=i; 
        } 
        use[i]=false; 
    } 
    tmp=PointSet[index]; 
    first=index; 
    use[index]=true; 
 
    index=-1; 
    ch[top++]=tmp; 
    tmp.x-=100; 
    l1.s=tmp; 
    l1.e=ch[0]; 
    l2.s=ch[0]; 
 
    while(index!=first) 
    { 
        curmax=-100; 
        curdis=0; 
        // 選取與最後一條確定邊夾角最小的點，即餘弦值最大者 
        for(i=0;i<n;i++) 
        { 
            if(use[i])continue; 
            l2.e=PointSet[i]; 
            curcos=cosine(l1,l2); // 根據cos值求夾角餘弦，範圍在 （-1 -- 1 ） 
            if(curcos>curmax || fabs(curcos-curmax)<1e-6 && dist(l2.s,l2.e)>curdis) 
            { 
                curmax=curcos; 
                index=i; 
                curdis=dist(l2.s,l2.e); 
            } 
        } 
        use[first]=false;            //清空第first個頂點標誌，使最後能形成封閉的hull 
        use[index]=true; 
        ch[top++]=PointSet[index]; 
        l1.s=ch[top-2]; 
        l1.e=ch[top-1]; 
        l2.s=ch[top-1]; 
    } 
    len=top-1; 
} 
/*********************************************************************************************  
    判斷線段是否在簡單多邊形內(注意：如果多邊形是凸多邊形，下面的演算法可以化簡) 
       必要條件一：線段的兩個端點都在多邊形內； 
    必要條件二：線段和多邊形的所有邊都不內交； 
    用途：    1. 判斷折線是否在簡單多邊形內 
            2. 判斷簡單多邊形是否在另一個簡單多邊形內 
**********************************************************************************************/ 
bool LinesegInsidePolygon(int vcount,POINT polygon[],LINESEG l) 
{ 
    // 判斷線端l的端點是否不都在多邊形內 
    if(!insidepolygon(vcount,polygon,l.s)||!insidepolygon(vcount,polygon,l.e)) 
        return false; 
    int top=0,i,j; 
    POINT PointSet[MAXV],tmp; 
    LINESEG s; 
 
    for(i=0;i<vcount;i++) 
    { 
        s.s=polygon[i]; 
        s.e=polygon[(i+1)%vcount]; 
        if(online(s,l.s)) //線段l的起始端點線上段s上 
            PointSet[top++]=l.s; 
        else if(online(s,l.e)) //線段l的終止端點線上段s上 
            PointSet[top++]=l.e; 
        else 
        { 
            if(online(l,s.s)) //線段s的起始端點線上段l上 
                PointSet[top++]=s.s; 
            else if(online(l,s.e)) // 線段s的終止端點線上段l上 
                PointSet[top++]=s.e; 
            else 
            { 
                if(intersect(l,s)) // 這個時候如果相交，肯定是內交，返回false 
                return false; 
            } 
        } 
    } 
 
    for(i=0;i<top-1;i++) /* 氣泡排序，x座標小的排在前面；x座標相同者，y座標小的排在前面 */ 
    { 
        for(j=i+1;j<top;j++) 
        { 
            if( PointSet[i].x>PointSet[j].x || fabs(PointSet[i].x-PointSet[j].x)<EP && PointSet[i].y>PointSet[j].y ) 
            { 
                tmp=PointSet[i]; 
                PointSet[i]=PointSet[j]; 
                PointSet[j]=tmp; 
            } 
        } 
    } 
 
    for(i=0;i<top-1;i++) 
    { 
        tmp.x=(PointSet[i].x+PointSet[i+1].x)/2; //得到兩個相鄰交點的中點 
        tmp.y=(PointSet[i].y+PointSet[i+1].y)/2; 
        if(!insidepolygon(vcount,polygon,tmp)) 
            return false; 
    } 
    return true; 
} 
/*********************************************************************************************  
求任意簡單多邊形polygon的重心 
需要呼叫下面幾個函式： 
    void AddPosPart(); 增加右邊區域的面積 
    void AddNegPart(); 增加左邊區域的面積 
    void AddRegion(); 增加區域面積 
在使用該程式時，如果把xtr,ytr,wtr,xtl,ytl,wtl設成全域性變數就可以使這些函式的形式得到化簡,
但要注意函式的宣告和呼叫要做相應變化 
**********************************************************************************************/ 
void AddPosPart(double x, double y, double w, double &xtr, double &ytr, double &wtr) 
{ 
    if (abs(wtr + w)<1e-10 ) return; // detect zero regions 
    xtr = ( wtr*xtr + w*x ) / ( wtr + w ); 
    ytr = ( wtr*ytr + w*y ) / ( wtr + w ); 
    wtr = w + wtr; 
    return; 
} 
void AddNegPart(double x, ouble y, double w, double &xtl, double &ytl, double &wtl) 
{ 
    if ( abs(wtl + w)<1e-10 ) 
        return; // detect zero regions 
 
    xtl = ( wtl*xtl + w*x ) / ( wtl + w ); 
    ytl = ( wtl*ytl + w*y ) / ( wtl + w ); 
    wtl = w + wtl; 
    return; 
} 
void AddRegion ( double x1, double y1, double x2, double y2, double &xtr, double &ytr, 
        double &wtr, double &xtl, double &ytl, double &wtl ) 
{ 
    if ( abs (x1 - x2)< 1e-10 ) 
        return; 
 
    if ( x2 > x1 ) 
    { 
        AddPosPart ((x2+x1)/2, y1/2, (x2-x1) * y1,xtr,ytr,wtr); /* rectangle 全域性變數變化處 */ 
        AddPosPart ((x1+x2+x2)/3, (y1+y1+y2)/3, (x2-x1)*(y2-y1)/2,xtr,ytr,wtr);    
        // triangle 全域性變數變化處 
    } 
    else 
    { 
        AddNegPart ((x2+x1)/2, y1/2, (x2-x1) * y1,xtl,ytl,wtl);  
        // rectangle 全域性變數變化處 
        AddNegPart ((x1+x2+x2)/3, (y1+y1+y2)/3, (x2-x1)*(y2-y1)/2,xtl,ytl,wtl);  
        // triangle  全域性變數變化處 
    } 
} 
POINT cg_simple(int vcount,POINT polygon[]) 
{ 
    double xtr,ytr,wtr,xtl,ytl,wtl;        
    //注意： 如果把xtr,ytr,wtr,xtl,ytl,wtl改成全域性變數後這裡要刪去 
    POINT p1,p2,tp; 
    xtr = ytr = wtr = 0.0; 
    xtl = ytl = wtl = 0.0; 
    for(int i=0;i<vcount;i++) 
    { 
        p1=polygon[i]; 
        p2=polygon[(i+1)%vcount]; 
        AddRegion(p1.x,p1.y,p2.x,p2.y,xtr,ytr,wtr,xtl,ytl,wtl); //全域性變數變化處 
    } 
    tp.x = (wtr*xtr + wtl*xtl) / (wtr + wtl); 
    tp.y = (wtr*ytr + wtl*ytl) / (wtr + wtl); 
    return tp; 
} 
// 求凸多邊形的重心,要求輸入多邊形按逆時針排序 
POINT gravitycenter(int vcount,POINT polygon[]) 
{ 
    POINT tp; 
    double x,y,s,x0,y0,cs,k; 
    x=0;y=0;s=0; 
    for(int i=1;i<vcount-1;i++) 
    { 
        x0=(polygon[0].x+polygon[i].x+polygon[i+1].x)/3; 
        y0=(polygon[0].y+polygon[i].y+polygon[i+1].y)/3; //求當前三角形的重心 
        cs=multiply(polygon[i],polygon[i+1],polygon[0])/2; 
        //三角形面積可以直接利用該公式求解 
        if(abs(s)<1e-20) 
        { 
            x=x0;y=y0;s+=cs;continue; 
        } 
        k=cs/s; //求面積比例 
        x=(x+k*x0)/(1+k); 
        y=(y+k*y0)/(1+k); 
        s += cs; 
    } 
    tp.x=x; 
    tp.y=y; 
    return tp; 
} 
 
/************************************************
給定一簡單多邊形，找出一個肯定在該多邊形內的點 
定理1    ：每個多邊形至少有一個凸頂點 
定理2    ：頂點數>=4的簡單多邊形至少有一條對角線 
結論    ： x座標最大，最小的點肯定是凸頂點 
    y座標最大，最小的點肯定是凸頂點            
************************************************/ 
POINT a_point_insidepoly(int vcount,POINT polygon[]) 
{ 
    POINT v,a,b,r; 
    int i,index; 
    v=polygon[0]; 
    index=0; 
    for(i=1;i<vcount;i++) //尋找一個凸頂點 
    { 
        if(polygon[i].y<v.y) 
        { 
            v=polygon[i]; 
            index=i; 
        } 
    } 
    a=polygon[(index-1+vcount)%vcount]; //得到v的前一個頂點 
    b=polygon[(index+1)%vcount]; //得到v的後一個頂點 
    POINT tri[3],q; 
    tri[0]=a;tri[1]=v;tri[2]=b; 
    double md=INF; 
    int in1=index; 
    bool bin=false; 
    for(i=0;i<vcount;i++) //尋找在三角形avb內且離頂點v最近的頂點q 
    { 
        if(i == index)continue; 
        if(i == (index-1+vcount)%vcount)continue; 
        if(i == (index+1)%vcount)continue; 
        if(!InsideConvexPolygon(3,tri,polygon[i]))continue; 
        bin=true; 
        if(dist(v,polygon[i])<md) 
        { 
            q=polygon[i]; 
            md=dist(v,q); 
        } 
    } 
    if(!bin) //沒有頂點在三角形avb內，返回線段ab中點 
    { 
        r.x=(a.x+b.x)/2; 
        r.y=(a.y+b.y)/2; 
        return r; 
    } 
    r.x=(v.x+q.x)/2; //返回線段vq的中點 
    r.y=(v.y+q.y)/2; 
    return r; 
} 
/***********************************************************************************************
求從多邊形外一點p出發到一個簡單多邊形的切線,如果存在返回切點,其中rp點是右切點,lp是左切點 
注意：p點一定要在多邊形外 ,輸入頂點序列是逆時針排列 
原 理：    如果點在多邊形內肯定無切線;凸多邊形有唯一的兩個切點,凹多邊形就可能有多於兩個的切點; 
        如果polygon是凸多邊形，切點只有兩個只要找到就可以,可以化簡此演算法 
        如果是凹多邊形還有一種演算法可以求解:先求凹多邊形的凸包,然後求凸包的切線 
/***********************************************************************************************/ 
void pointtangentpoly(int vcount,POINT polygon[],POINT p,POINT &rp,POINT &lp) 
{ 
    LINESEG ep,en; 
    bool blp,bln; 
    rp=polygon[0]; 
    lp=polygon[0]; 
    for(int i=1;i<vcount;i++) 
    { 
        ep.s=polygon[(i+vcount-1)%vcount]; 
        ep.e=polygon[i]; 
        en.s=polygon[i]; 
        en.e=polygon[(i+1)%vcount]; 
        blp=multiply(ep.e,p,ep.s)>=0;                // p is to the left of pre edge 
        bln=multiply(en.e,p,en.s)>=0;                // p is to the left of next edge 
        if(!blp&&bln) 
        { 
            if(multiply(polygon[i],rp,p)>0)           // polygon[i] is above rp 
            rp=polygon[i]; 
        } 
        if(blp&&!bln) 
        { 
            if(multiply(lp,polygon[i],p)>0)           // polygon[i] is below lp 
            lp=polygon[i]; 
        } 
    } 
    return ; 
} 
// 如果多邊形polygon的核存在，返回true，返回核上的一點p.頂點按逆時針方向輸入  
bool core_exist(int vcount,POINT polygon[],POINT &p) 
{ 
    int i,j,k; 
    LINESEG l; 
    LINE lineset[MAXV]; 
    for(i=0;i<vcount;i++) 
    { 
        lineset[i]=makeline(polygon[i],polygon[(i+1)%vcount]); 
    } 
    for(i=0;i<vcount;i++) 
    { 
        for(j=0;j<vcount;j++) 
        { 
            if(i == j)continue; 
            if(lineintersect(lineset[i],lineset[j],p)) 
            { 
                for(k=0;k<vcount;k++) 
                { 
                    l.s=polygon[k]; 
                    l.e=polygon[(k+1)%vcount]; 
                    if(multiply(p,l.e,l.s)>0)      
                    //多邊形頂點按逆時針方向排列，核肯定在每條邊的左側或邊上 
                    break; 
                } 
                if(k == vcount)             //找到了一個核上的點 
                break; 
            } 
        } 
        if(j<vcount) break; 
    } 
    if(i<vcount) 
        return true; 
    else 
        return false; 
} 
/*************************\ 
*                         * 
* 圓的基本運算           * 
*                         * 
\*************************/ 
/******************************************************************************
返回值    ： 點p在圓內(包括邊界)時，返回true 
用途    ： 因為圓為凸集，所以判斷點集，折線，多邊形是否在圓內時，
    只需要逐一判斷點是否在圓內即可。 
*******************************************************************************/ 
bool point_in_circle(POINT o,double r,POINT p) 
{ 
    double d2=(p.x-o.x)*(p.x-o.x)+(p.y-o.y)*(p.y-o.y); 
    double r2=r*r; 
    return d2<r2||abs(d2-r2)<EP; 
} 
/******************************************************************************
用 途    ：求不共線的三點確定一個圓 
輸 入    ：三個點p1,p2,p3 
返回值    ：如果三點共線，返回false；反之，返回true。圓心由q返回，半徑由r返回 
*******************************************************************************/ 
bool cocircle(POINT p1,POINT p2,POINT p3,POINT &q,double &r) 
{ 
    double x12=p2.x-p1.x; 
    double y12=p2.y-p1.y; 
    double x13=p3.x-p1.x; 
    double y13=p3.y-p1.y; 
    double z2=x12*(p1.x+p2.x)+y12*(p1.y+p2.y); 
    double z3=x13*(p1.x+p3.x)+y13*(p1.y+p3.y); 
    double d=2.0*(x12*(p3.y-p2.y)-y12*(p3.x-p2.x)); 
    if(abs(d)<EP) //共線，圓不存在 
        return false; 
    q.x=(y13*z2-y12*z3)/d; 
    q.y=(x12*z3-x13*z2)/d; 
    r=dist(p1,q); 
    return true; 
} 
int line_circle(LINE l,POINT o,double r,POINT &p1,POINT &p2) 
{ 
    return true; 
} 
 
/**************************\ 
*                          * 
* 矩形的基本運算          * 
*                         * 
\**************************/ 
/* 
說明：因為矩形的特殊性，常用演算法可以化簡： 
1.判斷矩形是否包含點 
只要判斷該點的橫座標和縱座標是否夾在矩形的左右邊和上下邊之間。 
2.判斷線段、折線、多邊形是否在矩形中 
因為矩形是個凸集，所以只要判斷所有端點是否都在矩形中就可以了。 
3.判斷圓是否在矩形中 
圓在矩形中的充要條件是：圓心在矩形中且圓的半徑小於等於圓心到矩形四邊的距離的最小值。 
*/ 
// 已知矩形的三個頂點(a,b,c)，計算第四個頂點d的座標. 注意：已知的三個頂點可以是無序的 
POINT rect4th(POINT a,POINT b,POINT c) 
{ 
    POINT d; 
    if(abs(dotmultiply(a,b,c))<EP) // 說明c點是直角拐角處 
    { 
        d.x=a.x+b.x-c.x; 
        d.y=a.y+b.y-c.y; 
    } 
    if(abs(dotmultiply(a,c,b))<EP) // 說明b點是直角拐角處 
    { 
        d.x=a.x+c.x-b.x; 
        d.y=a.y+c.y-b.x; 
    } 
    if(abs(dotmultiply(c,b,a))<EP) // 說明a點是直角拐角處 
    { 
        d.x=c.x+b.x-a.x; 
        d.y=c.y+b.y-a.y; 
    } 
    return d; 
} 
 
/*************************\ 
*                        * 
* 常用演算法的描述        * 
*                        * 
\*************************/ 
/* 
尚未實現的演算法： 
1. 求包含點集的最小圓 
2. 求多邊形的交 
3. 簡單多邊形的三角剖分 
4. 尋找包含點集的最小矩形 
5. 折線的化簡 
6. 判斷矩形是否在矩形中 
7. 判斷矩形能否放在矩形中 
8. 矩形並的面積與周長 
9. 矩形並的輪廓 
10.矩形並的閉包 
11.矩形的交 
12.點集中的最近點對 
13.多邊形的並 
14.圓的交與並 
15.直線與圓的關係 
16.線段與圓的關係 
17.求多邊形的核監視攝象機 
18.求點集中不相交點對 railwai 
*//* 
尋找包含點集的最小矩形 
原理：該矩形至少一條邊與點集的凸殼的某條邊共線 
First take the convex hull of the points. Let the resulting convex 
polygon be P. It has been known for some time that the minimum 
area rectangle enclosing P must have one rectangle side flush with 
(i.e., collinear with and overlapping) one edge of P. This geometric 
fact was used by Godfried Toussaint to develop the "rotating calipers" 
algorithm in a hard-to-find 1983 paper, "Solving Geometric Problems 
with the Rotating Calipers" (Proc. IEEE MELECON). The algorithm 
rotates a surrounding rectangle from one flush edge to the next, 
keeping track of the minimum area for each edge. It achieves O(n) 
time (after hull computation). See the "Rotating Calipers Homepage" 
http://www.cs.mcgill.ca/~orm/rotcal.frame.html for a description 
and applet. 
*//* 
折線的化簡 偽碼如下： 
Input: tol = the approximation tolerance 
L = {V0,V1,...,Vn-1} is any n-vertex polyline 
Set start = 0; 
Set k = 0; 
Set W0 = V0; 
for each vertex Vi (i=1,n-1) 
{ 
if Vi is within tol from Vstart 
then ignore it, and continue with the next vertex 
Vi is further than tol away from Vstart 
so add it as a new vertex of the reduced polyline 
Increment k++; 
Set Wk = Vi; 
Set start = i; as the new initial vertex 
} 
Output: W = {W0,W1,...,Wk-1} = the k-vertex simplified polyline 
*/ 
/********************\ 
*                    * 
* 補充                * 
*                    * 
\********************/ 
 
//兩圓關係： 
/* 兩圓： 
相離： return 1； 
外切： return 2； 
相交： return 3； 
內切： return 4； 
內含： return 5； 
*/ 
int CircleRelation(POINT p1, double r1, POINT p2, double r2) 
{ 
    double d = sqrt( (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) ); 
 
    if( fabs(d-r1-r2) < EP ) // 必須保證前兩個if先被判定！ 
        return 2; 
    if( fabs(d-fabs(r1-r2)) < EP ) 
        return 4; 
    if( d > r1+r2 ) 
        return 1; 
    if( d < fabs(r1-r2) ) 
        return 5; 
    if( fabs(r1-r2) < d && d < r1+r2 ) 
        return 3; 
    return 0; // indicate an error! 
} 
//判斷圓是否在矩形內：
// 判定圓是否在矩形內，是就返回true（設矩形水平，且其四個頂點由左上開始按順時針排列） 
// 呼叫ptoldist函式，在第4頁 
bool CircleRecRelation(POINT pc, double r, POINT pr1, POINT pr2, POINT pr3, POINT pr4) 
{ 
    if( pr1.x < pc.x && pc.x < pr2.x && pr3.y < pc.y && pc.y < pr2.y ) 
    { 
        LINESEG line1(pr1, pr2); 
        LINESEG line2(pr2, pr3); 
        LINESEG line3(pr3, pr4); 
        LINESEG line4(pr4, pr1); 
        if( r<ptoldist(pc,line1) && r<ptoldist(pc,line2) &&    r<ptoldist(pc,line3) && r<ptoldist(pc,line4) ) 
            return true; 
    } 
    return false; 
} 
//點到平面的距離： 
//點到平面的距離,平面用一般式表示ax+by+cz+d=0 
double P2planeDist(double x, double y, double z, double a, double b, double c, double d) 
{ 
    return fabs(a*x+b*y+c*z+d) / sqrt(a*a+b*b+c*c); 
} 
//點是否在直線同側：
//兩個點是否在直線同側，是則返回true 
bool SameSide(POINT p1, POINT p2, LINE line) 
{ 
    return (line.a * p1.x + line.b * p1.y + line.c) * 
    (line.a * p2.x + line.b * p2.y + line.c) > 0; 
}
//鏡面反射線：
// 已知入射線、鏡面，求反射線。 
// a1,b1,c1為鏡面直線方程(a1 x + b1 y + c1 = 0 ,下同)係數;  
//a2,b2,c2為入射光直線方程係數;  
//a,b,c為反射光直線方程係數. 
// 光是有方向的，使用時注意：入射光向量:<-b2,a2>；反射光向量:<b,-a>. 
// 不要忘記結果中可能會有"negative zeros" 
void reflect(double a1,double b1,double c1,double a2,double b2,double c2,double &a,double &b,double &c) 
{ 
    double n,m; 
    double tpb,tpa; 
    tpb=b1*b2+a1*a2; 
    tpa=a2*b1-a1*b2; 
    m=(tpb*b1+tpa*a1)/(b1*b1+a1*a1); 
    n=(tpa*b1-tpb*a1)/(b1*b1+a1*a1); 
    if(fabs(a1*b2-a2*b1)<1e-20) 
    { 
        a=a2;b=b2;c=c2; 
        return; 
    } 
    double xx,yy; //(xx,yy)是入射線與鏡面的交點。 
    xx=(b1*c2-b2*c1)/(a1*b2-a2*b1); 
    yy=(a2*c1-a1*c2)/(a1*b2-a2*b1); 
    a=n; 
    b=-m; 
    c=m*yy-xx*n; 
} 
//矩形包含： 
// 矩形2（C，D）是否在1（A，B）內
bool r2inr1(double A,double B,double C,double D)  
{ 
    double X,Y,L,K,DMax; 
    if (A < B) 
    { 
        double tmp = A; 
        A = B; 
        B = tmp; 
    } 
    if (C < D) 
    { 
        double tmp = C; 
        C = D; 
        D = tmp; 
    } 
    if (A > C && B > D)                 // trivial case  
        return true; 
    else 
        if (D >= B) 
            return false; 
        else 
        { 
            X = sqrt(A * A + B * B);         // outer rectangle's diagonal  
            Y = sqrt(C * C + D * D);         // inner rectangle's diagonal  
            if (Y < B) // check for marginal conditions 
                return true; // the inner rectangle can freely rotate inside 
            else 
                if (Y > X) 
                    return false; 
                else 
                { 
                    L = (B - sqrt(Y * Y - A * A)) / 2; 
                    K = (A - sqrt(Y * Y - B * B)) / 2; 
                    DMax = sqrt(L * L + K * K); 
                    if (D >= DMax) 
                    return false; 
                    else 
                    return true; 
                } 
        } 
} 
//兩圓交點： 
// 兩圓已經相交（相切） 
void  c2point(POINT p1,double r1,POINT p2,double r2,POINT &rp1,POINT &rp2) 
{ 
    double a,b,r; 
    a=p2.x-p1.x; 
    b=p2.y-p1.y; 
    r=(a*a+b*b+r1*r1-r2*r2)/2; 
    if(a==0&&b!=0) 
    { 
        rp1.y=rp2.y=r/b; 
        rp1.x=sqrt(r1*r1-rp1.y*rp1.y); 
        rp2.x=-rp1.x; 
    } 
    else if(a!=0&&b==0) 
    { 
        rp1.x=rp2.x=r/a; 
        rp1.y=sqrt(r1*r1-rp1.x*rp2.x); 
        rp2.y=-rp1.y; 
    } 
    else if(a!=0&&b!=0) 
    { 
        double delta; 
        delta=b*b*r*r-(a*a+b*b)*(r*r-r1*r1*a*a); 
        rp1.y=(b*r+sqrt(delta))/(a*a+b*b); 
        rp2.y=(b*r-sqrt(delta))/(a*a+b*b); 
        rp1.x=(r-b*rp1.y)/a; 
        rp2.x=(r-b*rp2.y)/a; 
    } 
 
    rp1.x+=p1.x; 
    rp1.y+=p1.y; 
    rp2.x+=p1.x; 
    rp2.y+=p1.y; 
} 
//兩圓公共面積：
// 必須保證相交 
double c2area(POINT p1,double r1,POINT p2,double r2) 
{ 
    POINT rp1,rp2; 
    c2point(p1,r1,p2,r2,rp1,rp2); 
 
    if(r1>r2) //保證r2>r1 
    { 
        swap(p1,p2); 
        swap(r1,r2); 
    } 
    double a,b,rr; 
    a=p1.x-p2.x; 
    b=p1.y-p2.y; 
    rr=sqrt(a*a+b*b); 
 
    double dx1,dy1,dx2,dy2; 
    double sita1,sita2; 
    dx1=rp1.x-p1.x; 
    dy1=rp1.y-p1.y; 
    dx2=rp2.x-p1.x; 
    dy2=rp2.y-p1.y; 
    sita1=acos((dx1*dx2+dy1*dy2)/r1/r1); 
 
    dx1=rp1.x-p2.x; 
    dy1=rp1.y-p2.y; 
    dx2=rp2.x-p2.x; 
    dy2=rp2.y-p2.y; 
    sita2=acos((dx1*dx2+dy1*dy2)/r2/r2); 
    double s=0; 
    if(rr<r2)//相交弧為優弧 
        s=r1*r1*(PI-sita1/2+sin(sita1)/2)+r2*r2*(sita2-sin(sita2))/2; 
    else//相交弧為劣弧 
        s=(r1*r1*(sita1-sin(sita1))+r2*r2*(sita2-sin(sita2)))/2; 
 
    return s; 
} 
//圓和直線關係： 
//0----相離 1----相切 2----相交 
int clpoint(POINT p,double r,double a,double b,double c,POINT &rp1,POINT &rp2) 
{ 
    int res=0; 
 
    c=c+a*p.x+b*p.y; 
    double tmp; 
    if(a==0&&b!=0) 
    { 
        tmp=-c/b; 
        if(r*r<tmp*tmp) 
            res=0; 
        else if(r*r==tmp*tmp) 
        { 
            res=1; 
            rp1.y=tmp; 
            rp1.x=0; 
        } 
        else 
        { 
            res=2; 
            rp1.y=rp2.y=tmp; 
            rp1.x=sqrt(r*r-tmp*tmp); 
            rp2.x=-rp1.x; 
        } 
    } 
    else if(a!=0&&b==0) 
    { 
        tmp=-c/a; 
        if(r*r<tmp*tmp) 
            res=0; 
        else if(r*r==tmp*tmp) 
        { 
            res=1; 
            rp1.x=tmp; 
            rp1.y=0; 
        } 
        else 
        { 
            res=2; 
            rp1.x=rp2.x=tmp; 
            rp1.y=sqrt(r*r-tmp*tmp); 
            rp2.y=-rp1.y; 
        } 
    } 
    else if(a!=0&&b!=0) 
    { 
        double delta; 
        delta=b*b*c*c-(a*a+b*b)*(c*c-a*a*r*r); 
        if(delta<0) 
            res=0; 
        else if(delta==0) 
        { 
            res=1; 
            rp1.y=-b*c/(a*a+b*b); 
            rp1.x=(-c-b*rp1.y)/a; 
        } 
        else 
        { 
            res=2; 
            rp1.y=(-b*c+sqrt(delta))/(a*a+b*b); 
            rp2.y=(-b*c-sqrt(delta))/(a*a+b*b); 
            rp1.x=(-c-b*rp1.y)/a; 
            rp2.x=(-c-b*rp2.y)/a; 
        } 
    } 
    rp1.x+=p.x; 
    rp1.y+=p.y; 
    rp2.x+=p.x; 
    rp2.y+=p.y; 
    return res; 
} 
//內切圓： 
void incircle(POINT p1,POINT p2,POINT p3,POINT &rp,double &r) 
{ 
    double dx31,dy31,dx21,dy21,d31,d21,a1,b1,c1; 
    dx31=p3.x-p1.x; 
    dy31=p3.y-p1.y; 
    dx21=p2.x-p1.x; 
    dy21=p2.y-p1.y; 
 
    d31=sqrt(dx31*dx31+dy31*dy31); 
    d21=sqrt(dx21*dx21+dy21*dy21); 
    a1=dx31*d21-dx21*d31; 
    b1=dy31*d21-dy21*d31; 
    c1=a1*p1.x+b1*p1.y; 
 
    double dx32,dy32,dx12,dy12,d32,d12,a2,b2,c2; 
    dx32=p3.x-p2.x; 
    dy32=p3.y-p2.y; 
    dx12=-dx21; 
    dy12=-dy21; 
 
    d32=sqrt(dx32*dx32+dy32*dy32); 
    d12=d21; 
    a2=dx12*d32-dx32*d12; 
    b2=dy12*d32-dy32*d12; 
    c2=a2*p2.x+b2*p2.y; 
 
    rp.x=(c1*b2-c2*b1)/(a1*b2-a2*b1); 
    rp.y=(c2*a1-c1*a2)/(a1*b2-a2*b1); 
    r=fabs(dy21*rp.x-dx21*rp.y+dx21*p1.y-dy21*p1.x)/d21; 
} 
//求切點： 
// p---圓心座標， r---圓半徑， sp---圓外一點， rp1,rp2---切點座標 
void cutpoint(POINT p,double r,POINT sp,POINT &rp1,POINT &rp2) 
{ 
    POINT p2; 
    p2.x=(p.x+sp.x)/2; 
    p2.y=(p.y+sp.y)/2; 
 
    double dx2,dy2,r2; 
    dx2=p2.x-p.x; 
    dy2=p2.y-p.y; 
    r2=sqrt(dx2*dx2+dy2*dy2); 
    c2point(p,r,p2,r2,rp1,rp2); 
} 
//線段的左右旋： 
/* l2在l1的左/右方向（l1為基準線） 
返回    0    ： 重合； 
返回    1    ： 右旋； 
返回    –1 ： 左旋； 
*/ 
int rotat(LINESEG l1,LINESEG l2) 
{ 
    double dx1,dx2,dy1,dy2; 
    dx1=l1.s.x-l1.e.x; 
    dy1=l1.s.y-l1.e.y; 
    dx2=l2.s.x-l2.e.x; 
    dy2=l2.s.y-l2.e.y; 
 
    double d; 
    d=dx1*dy2-dx2*dy1; 
    if(d==0) 
        return 0; 
    else if(d>0) 
        return -1; 
    else 
        return 1; 
} 
/*
公式： 
球座標公式： 
直角座標為 P(x, y, z) 時，對應的球座標是(rsinφcosθ, rsinφsinθ, rcosφ),其中φ是向量OP與Z軸的夾角，範圍[0，π]；是OP在XOY面上的投影到X軸的旋角，範圍[0，2π]  
直線的一般方程轉化成向量方程： 
ax+by+c=0 
x-x0     y-y0 
   ------ = ------- // (x0,y0)為直線上一點，m,n為向量 
m        n 
轉換關係： 
a=n；b=-m；c=m·y0-n·x0； 
m=-b; n=a; 
三點平面方程： 
三點為P1，P2，P3 
設向量  M1=P2-P1; M2=P3-P1; 
平面法向量：  M=M1 x M2 （） 
平面方程：    M.i(x-P1.x)+M.j(y-P1.y)+M.k(z-P1.z)=0 
*/