计算几何基础（模板）

1     double S(Point p[],int n)  
2     {  
3         double ans = 0;  
4         p[n] = p[0];  
5         for(int i=1;i<n;i++)  
6            ans += cross(p[0],p[i],p[i+1]);  
7         if(ans < 0) ans = -ans;  
8         return ans / 2.0;  
9     }  

 1     bool cmp(Point A,Point B)  
 2     {  
 3         double k = cross(MinA,A,B);  
 4         if(k<0) return 0;  
 5         if(k>0) return 1;  
 6         return dist(MinA,A)<dist(MinA,B);  
 7     }  
 8       
 9     void Graham(Point p[],int n)  
10     {  
11         for(int i=1;i<n;i++)  
12            if(p[i].y<p[0].y || (p[i].y == p[0].y && p[i].x < p[0].x))  
13                 swap(p[i],p[0]);  
14         MinA = p[0];  
15         p[n] = p[0];  
16         sort(p+1,p+n,cmp);  
17         stack[0] = p[0];  
18         stack[1] = p[1];  
19         top = 2;  
20         for(int i=2;i<n;i++)  
21         {  
22             while(top >= 2 && cross(stack[top-2],stack[top-1],p[i])<=0) top--;  
23             stack[top++] = p[i];  
24         }  
25     }  

 1     Point Gravity(Point p[],int n)  
 2     {  
 3         Point O,t;  
 4         O.x = O.y = 0;  
 5         t.x = t.y = 0;  
 6         p[n] = p[0];  
 7         double A = 0;  
 8         for(int i=0; i<n; i++)  
 9             A += cross(O,p[i],p[i+1]);  
10         A /= 2.0;  
11         for(int i=0; i<n; i++)  
12         {  
13             t.x += (p[i].x + p[i+1].x) * cross(O,p[i],p[i+1]);  
14             t.y += (p[i].y + p[i+1].y) * cross(O,p[i],p[i+1]);  
15         }  
16         t.x /= 6*A;  
17         t.y /= 6*A;  
18         return t;  
19     }  

 1     bool Segment_crossing(Segment u,Segment v)   /** 判断线段是否相交 */  
 2     {  
 3              return((max(u.a.x,u.b.x)>=min(v.a.x,v.b.x))&&  
 4                (max(v.a.x,v.b.x)>=min(u.a.x,u.b.x))&&  
 5                (max(u.a.y,u.b.y)>=min(v.a.y,v.b.y))&&  
 6                (max(v.a.y,v.b.y)>=min(u.a.y,u.b.y))&&  
 7                (cross(v.a,u.b,u.a)*cross(u.b,v.b,u.a)>=0)&&  
 8                (cross(u.a,v.b,v.a)*cross(v.b,u.b,v.a)>=0));  
 9     }  
10       
11     /**求直线交点的坐标，如果没有交点返回NULL，否则返回交点p的地址*/  
12     Point* CrossPoint(Segment u,Segment v)  
13     {  
14         Point p;  
15         if(Segment_crossing(u,v))  
16         {  
17             p.x=(cross(v.b,u.b,u.a)*v.a.x-cross(v.a,u.b,u.a)*v.b.x)/(cross(v.b,u.b,u.a)-cross(v.a,u.b,u.a));  
18             p.y=(cross(v.b,u.b,u.a)*v.a.y-cross(v.a,u.b,u.a)*v.b.y)/(cross(v.b,u.b,u.a)-cross(v.a,u.b,u.a));  
19             return &p;  
20         }  
21         return NULL;  
22     }  

 1     Point Circle_Point(Point A,Point B,Point C)  
 2     {  
 3         double a=dist(B,C);  
 4         double b=dist(A,C);  
 5         double c=dist(A,B);  
 6         double p=(a+b+c)/2.0;  
 7         double S=sqrt(p*(p-a)*(p-b)*(p-c));  
 8         R=(a*b*c)/(4*S);    //三角形外接圆的半径为R  
 9       
10         double t1=(A.x*A.x+A.y*A.y-B.x*B.x-B.y*B.y)/2;  
11         double t2=(A.x*A.x+A.y*A.y-C.x*C.x-C.y*C.y)/2;  
12         Point center;  
13         center.x=(t1*(A.y-C.y)-t2*(A.y-B.y))/((A.x-B.x)*(A.y-C.y)-(A.x-C.x)*(A.y-B.y));  
14         center.y=(t1*(A.x-C.x)-t2*(A.x-B.x))/((A.y-B.y)*(A.x-C.x)-(A.y-C.y)*(A.x-B.x));  
15         return center;  
16     }  

 1     double rotating_calipers(Point p[],int n)  
 2     {  
 3         int k = 1;  
 4         double ans = 0;  
 5         p[n] = p[0];  
 6         for(int i=0;i<n;i++)  
 7         {  
 8             while(fabs(cross(p[i],p[i+1],p[k])) < fabs(cross(p[i],p[i+1],p[k+1])))  
 9                  k = (k+1) % n;  
10             ans = max(ans, max(dist(p[i],p[k]),dist(p[i+1],p[k])));  
11         }  
12         return ans;  
13     }  

 1     double rotating_calipers(Point p[],int n)  
 2     {  
 3         int k = 1;  
 4         double ans = 0x7FFFFFFF;  
 5         p[n] = p[0];  
 6         for(int i=0;i<n;i++)  
 7         {  
 8             while(fabs(cross(p[i],p[i+1],p[k])) < fabs(cross(p[i],p[i+1],p[k+1])))  
 9                  k = (k+1) % n;  
10             double tmp = fabs(cross(p[i],p[i+1],p[k]));  
11             double d   = dist(p[i],p[i+1]);  
12             ans = min(ans,tmp/d);  
13         }  
14         return ans;  
15     }